diff --git a/mathics/builtin/arithfns/basic.py b/mathics/builtin/arithfns/basic.py index 756628ce4..9c9b3a7c9 100644 --- a/mathics/builtin/arithfns/basic.py +++ b/mathics/builtin/arithfns/basic.py @@ -357,6 +357,7 @@ class Power(InfixOperator, MPMathFunction): grouping = "Right" mpmath_name = "power" + complex_values = False @staticmethod def condition_numbers( diff --git a/mathics/builtin/numbers/exp.py b/mathics/builtin/numbers/exp.py index 19aa98803..87242c440 100644 --- a/mathics/builtin/numbers/exp.py +++ b/mathics/builtin/numbers/exp.py @@ -203,7 +203,7 @@ class Log(MPMathFunction): """ summary_text = "logarithm function" - nargs = {2} + nargs = {1, 2} mpmath_name = "log" sympy_name = "log" @@ -216,7 +216,6 @@ class Log(MPMathFunction): "Log[E^x_Integer]": "x", "Log[Overflow[]]": "Overflow[]", "Log[Undefined]": "Undefined", - "Log[x_?InexactNumberQ]": "Log[E, x]", } def prepare_sympy(self, elements): @@ -225,13 +224,18 @@ def prepare_sympy(self, elements): return elements def get_mpmath_function(self, args): + if len(args) == 1: + return mpmath.log return lambda base, x: mpmath.log(x, base) @staticmethod - def condition_numbers( - base: mpmath.mpf, x: mpmath.mpf, value: mpmath.mpf - ) -> tuple[mpmath.mpf, mpmath.mpf]: - # Log[b, x] is Log[x] / Log[b]. + def condition_numbers(*args: mpmath.mpf) -> tuple[mpmath.mpf, ...]: + # The arguments, then the value. Log[b, x] is Log[x] / Log[b]. As in + # WMA, Log[x] has no base E that would add an error of its own. + *arguments, value = args + if len(arguments) == 1: + return (1 / abs(value),) + base, x = arguments return (1 / abs(mpmath.log(base)), 1 / abs(mpmath.log(x))) diff --git a/mathics/builtin/numbers/hyperbolic.py b/mathics/builtin/numbers/hyperbolic.py index efc26236b..62d26c4ff 100644 --- a/mathics/builtin/numbers/hyperbolic.py +++ b/mathics/builtin/numbers/hyperbolic.py @@ -111,7 +111,6 @@ def condition_numbers(x: mpmath.mpf, value: mpmath.mpf) -> tuple[mpmath.mpf]: rules = { "ArcCoth[Undefined]": "Undefined", - "ArcCoth[z:0.0]": "N[I / 2 Pi, Precision[1+z]]", "Derivative[1][ArcCoth]": "1/(1-#^2)&", } diff --git a/mathics/core/atoms/numerics.py b/mathics/core/atoms/numerics.py index 8b23cb1fe..c0c04543f 100644 --- a/mathics/core/atoms/numerics.py +++ b/mathics/core/atoms/numerics.py @@ -1157,8 +1157,12 @@ def precision(self) -> Optional[int]: return self._precision def round(self, d=None) -> "Complex": - real = self._real.round(d) - imag = self._imag.round(d) + # As in WMA, N[I, 30] keeps the exact zero real part, and N[I] is + # 0. + 1. I. + real, imag = ( + part if d is not None and part.sameQ(Integer0) else part.round(d) + for part in (self._real, self._imag) + ) return Complex(real, imag) def sameQ(self, rhs) -> bool: diff --git a/mathics/core/builtin.py b/mathics/core/builtin.py index 9d3562bb4..2fd98afa9 100644 --- a/mathics/core/builtin.py +++ b/mathics/core/builtin.py @@ -829,6 +829,10 @@ class MPMathFunction(SympyFunction): # Real points where, as in WMA, an inexact argument whose error covers # one of them gives Indeterminate, like ArcTanh at 1. indeterminate_at: tuple[int, ...] = () + # As in WMA, a complex value at real arguments, like ArcSin[2.`20], gets + # an error for each part. Power does not follow that rule in WMA: see + # eval_Power_by_parts for a negative base to a rational power. + complex_values: bool = True def get_mpmath_function(self, args): if self.mpmath_name is None or len(args) not in self.nargs: @@ -870,6 +874,7 @@ def eval(self, z, evaluation: Evaluation): *args, has_poles=self.has_poles, bound=self.value_bound, + complex_values=self.complex_values, ) if result is not None: return result diff --git a/mathics/eval/arithfns/basic.py b/mathics/eval/arithfns/basic.py index 09871884b..6f70971ad 100644 --- a/mathics/eval/arithfns/basic.py +++ b/mathics/eval/arithfns/basic.py @@ -44,7 +44,11 @@ segregate_numbers_from_sorted_list, ) from mathics.eval.nevaluator import eval_N -from mathics.eval.numeric import eval_add_numbers, eval_multiply_numbers +from mathics.eval.numeric import ( + eval_add_numbers, + eval_multiply_numbers, + eval_Power_by_parts, +) RealM0p5 = Real(-0.5) RealOne = Real(1.0) @@ -169,6 +173,11 @@ def eval_Power(base, exponent, power_fn_eval: Callable, evaluation: Evaluation): pass return result + if isinstance(base, Number) and isinstance(exponent, Number): + result = eval_Power_by_parts(base, exponent) + if result is not None: + return result + if isinstance(base, Complex) and base.real.is_zero: yhalf = Expression(SymbolTimes, exponent, RationalOneHalf) factor = power_fn_eval( diff --git a/mathics/eval/numbers/numbers.py b/mathics/eval/numbers/numbers.py index 34489d00b..0784738a4 100644 --- a/mathics/eval/numbers/numbers.py +++ b/mathics/eval/numbers/numbers.py @@ -3,6 +3,7 @@ Implementation of number-handling functions. """ +import math from functools import cache from typing import Optional @@ -160,7 +161,11 @@ def eval_Accuracy(z: BaseElement) -> Optional[float]: if acc_imag is None: return acc_real - return -mpmath.log(10 ** (-2 * acc_real) + 10 ** (-2 * acc_imag), 10.0) * 0.5 + # As in WMA, the errors of the parts add in quadrature. + least = min(acc_real, acc_imag) + return least - 0.5 * math.log10( + 10.0 ** (2 * (least - acc_real)) + 10.0 ** (2 * (least - acc_imag)) + ) if isinstance(z, Expression): elem_accuracies = (eval_Accuracy(z_elem) for z_elem in z.elements) @@ -209,18 +214,18 @@ def eval_Precision(z: BaseElement) -> Optional[float]: return z.precision10 if isinstance(z, Complex): - prec_real = eval_Precision(z.real) - prec_imag = eval_Precision(z.imag) - if prec_real is None or prec_imag == prec_real: - return prec_imag - if prec_imag is None: - return prec_real - # both numbers have different precision. - # Evaluate the accuracy and add the log of - # the module. + if z.is_machine_precision: + return MACHINE_PRECISION_VALUE acc = eval_Accuracy(z) - abs_sq = z.real.value**2 + z.imag.value**2 - return acc + mpmath.log(abs_sq, 10.0) * 0.5 + if acc is None: + return None + # As in WMA, the accuracy plus the log of the modulus, to which an + # exact part also counts: Precision[1.`20 + 3 I] is 20.5. + with mpmath.workprec(53): + modulus = abs(z.to_mpmath()) + if modulus == 0: + return 0.0 + return acc + float(mpmath.log10(modulus)) if isinstance(z, Expression): elem_prec = (eval_Precision(z_elem) for z_elem in z.elements) diff --git a/mathics/eval/numeric.py b/mathics/eval/numeric.py index 69dee5ff3..290043b7e 100644 --- a/mathics/eval/numeric.py +++ b/mathics/eval/numeric.py @@ -9,7 +9,7 @@ import math from fractions import Fraction -from typing import Callable, Optional +from typing import Callable, Final, Optional import mpmath import sympy @@ -678,13 +678,13 @@ def eval_mpmath_function_by_condition( *args: Number, has_poles: bool = False, bound: Optional[float] = None, + complex_values: bool = True, ) -> Optional[BaseElement]: """ Applies ``mpmath_function`` to ``args``, when they are exact real - numbers or PrecisionReals and the value is a real number. As in WMA, - the relative error of each argument is multiplied by its factor in - ``condition_numbers``, and these errors add up. Returns None otherwise, - or if no argument adds an error. + numbers or PrecisionReals. As in WMA, the relative error of each + argument is multiplied by its factor in ``condition_numbers``, and these + errors add up. Returns None otherwise, or if no argument adds an error. Where the factor does not exist, the error is how much f changes within the error of its single argument (see ``_change_within``): for a zero @@ -694,6 +694,9 @@ def eval_mpmath_function_by_condition( As in WMA, a function with poles (``has_poles``) gives ComplexInfinity when no digit is left, as Tan[1.5707963`5]; and a function bounded by ``bound`` has an error no larger than that, as Sin[1.`20*^30]. + + A complex value, as ArcSin[2.`20], gets the errors of its parts from + ``_complex_value_by_parts``, unless ``complex_values`` is False. """ if len(args) == 1 and isinstance(args[0], PrecisionReal) and args[0].is_zero: accuracy10 = args[0].precision10 @@ -712,7 +715,11 @@ def eval_mpmath_function_by_condition( with mpmath.workprec(input_bits + 10): values = [_full_mpmath(arg) for arg in args] value = mpmath_function(*values) - if not isinstance(value, mpmath.mpf) or not mpmath.isfinite(value): + if not isinstance(value, (mpmath.mpf, mpmath.mpc)) or not mpmath.isfinite( + value + ): + return None + if isinstance(value, mpmath.mpc) and not complex_values: return None try: factors = condition_numbers(*values, value) if value != 0 else None @@ -733,6 +740,9 @@ def eval_mpmath_function_by_condition( return None return _value_with_error(*change, has_poles, bound) + if isinstance(value, mpmath.mpc): + return _complex_value_by_parts(mpmath_function, args, input_bits) + digits = [ arg.precision10 - _log10_abs(factor) for arg, factor in zip(args, factors) @@ -754,6 +764,101 @@ def eval_mpmath_function_by_condition( return _precision_real(value, precision10, accuracy10) +# As in WMA, a part of a complex value that does not change with the +# arguments, like the Pi of Log[-2.`20], still gets this fraction of the +# error of the whole: Log[-1.`20] is 0``20 + Pi with precision 36.75. +CONSTANT_PART_ERROR: Final = mpmath.mpf(2) ** -54 + + +def _complex_value_by_parts( + mpmath_function: Callable, args: tuple[Number, ...], input_bits: int +) -> Optional[BaseElement]: + """ + The complex value of ``mpmath_function`` at the real ``args``. As in + WMA, each part has its own error: the sum, over the arguments, of the + relative error of the argument times that part of ``x D[f, x]``. A part + that does not change gets CONSTANT_PART_ERROR times ``Abs[x D[f, x]]``, + and is an exact 0 if its value is 0, as the real part of ArcSec[0.5`20]. + + ``x D[f, x]`` is taken from a central difference. Returns None if f + has no complex value near ``args``. + """ + errors = [mpmath.mpf(0), mpmath.mpf(0)] + changes = [False, False] + bits = input_bits + 64 + with mpmath.workprec(bits): + values = [_full_mpmath(arg) for arg in args] + step = mpmath.mpf(2) ** -(bits // 3) + for k, arg in enumerate(args): + if not isinstance(arg, PrecisionReal): + continue + up, down = list(values), list(values) + up[k] *= 1 + step + down[k] *= 1 - step + try: + difference = mpmath_function(*up) - mpmath_function(*down) + except (ZeroDivisionError, ValueError): + return None + if not mpmath.isfinite(difference): + return None + sensitivity = mpmath.mpc(difference) / (2 * step) + error = mpmath.mpf(10) ** -arg.precision10 + floor = abs(sensitivity) * CONSTANT_PART_ERROR + for i, part in enumerate((sensitivity.real, sensitivity.imag)): + errors[i] += error * max(abs(part), floor) + changes[i] = changes[i] or abs(part) > floor + rough = mpmath.mpc(mpmath_function(*values)) + + if errors[0] == 0: + # f does not change at all with its arguments. + return None + accuracies = [-_log10_abs(error) for error in errors] + precisions = [ + accuracy + _log10_abs(part) if part != 0 else 0 + for accuracy, part in zip(accuracies, (rough.real, rough.imag)) + ] + needed_bits = number_prec(max(1, math.ceil(max(precisions)))) + 10 + value = rough + if needed_bits > bits: + # The constant part has more digits than the arguments. + with mpmath.workprec(needed_bits): + value = mpmath.mpc(mpmath_function(*(_full_mpmath(arg) for arg in args))) + return Complex( + *( + ( + Integer0 + if part == 0 and not change + else _precision_real(part, precision, accuracy) + ) + for part, change, precision, accuracy in zip( + (value.real, value.imag), changes, precisions, accuracies + ) + ) + ) + + +def _complex_with_error(value: mpmath.mpc, error: mpmath.mpf) -> BaseElement: + """ + ``value`` with the absolute error ``error`` of the whole number. As in + WMA, a part that is 0 is an exact 0, as the real part of ArcCosh[0``20], + and if no part is 0, each gets the error divided by Sqrt[2]. + """ + parts = (value.real, value.imag) + if all(part != 0 for part in parts): + error = error / mpmath.sqrt(2) + accuracy10 = -_log10_abs(error) + return Complex( + *( + ( + Integer0 + if part == 0 + else _precision_real(part, accuracy10 + _log10_abs(part), accuracy10) + ) + for part in parts + ) + ) + + def covers_point(number: Number, points: tuple[int, ...]) -> bool: """ True if ``number`` is a nonzero inexact real number and one of @@ -788,6 +893,9 @@ def _change_within( agree, and a branch cut, as for ArcCot at 0, jumps on one side only. A side where f has no real value, as above 1 for ArcCos, is left out. Returns None if f has no real value at x or on either side. + + If f has a complex value at x, as ArcCosh[0], the sides with a complex + value count as well. """ # Enough digits for a change of second order. with mpmath.workprec(number_prec(math.ceil(2 * max(digits, 1)) + 20)): @@ -796,15 +904,20 @@ def _change_within( value = mpmath_function(x) except (ZeroDivisionError, ValueError): return None - if not isinstance(value, mpmath.mpf) or not mpmath.isfinite(value): + if not isinstance(value, (mpmath.mpf, mpmath.mpc)) or not mpmath.isfinite( + value + ): return None + kinds = ( + (mpmath.mpf, mpmath.mpc) if isinstance(value, mpmath.mpc) else mpmath.mpf + ) changes = [] for side in (x + delta, x - delta): try: other = mpmath_function(side) except (ZeroDivisionError, ValueError): continue - if isinstance(other, mpmath.mpf) and mpmath.isfinite(other): + if isinstance(other, kinds) and mpmath.isfinite(other): changes.append(abs(other - value)) if not changes or min(changes) == 0: return None @@ -815,6 +928,8 @@ def _value_with_error( value: mpmath.mpf, error: mpmath.mpf, has_poles: bool, bound: Optional[float] ) -> BaseElement: """``value`` with the absolute error ``error``.""" + if isinstance(value, mpmath.mpc): + return _complex_with_error(value, error) accuracy10 = -_log10_abs(error) if bound is not None: accuracy10 = max(accuracy10, -math.log10(bound)) @@ -842,6 +957,9 @@ def eval_add_numbers( number = mpmath.fsum(terms) return from_mpmath(number) + if _precision_complexes(numbers): + return _add_complex_numbers(numbers) + if _exact_or_precision_reals(numbers, allow_zeros=True): # As in WMA, the absolute errors of the terms add up. For a zero, # precision10 is its accuracy. @@ -909,6 +1027,9 @@ def eval_multiply_numbers(*numbers: Number) -> Number: number = mpmath.fprod(factors) return from_mpmath(number) + if _precision_complexes(numbers): + return _multiply_complex_numbers(numbers) + if _exact_or_precision_reals(numbers, allow_zeros=True): if any(number is Integer0 for number in numbers): # As in WMA, an exact zero gives an exact zero. @@ -954,6 +1075,196 @@ def _multiply_by_zeros(zeros: list[Number], numbers: tuple[Number, ...]) -> Numb return _precision_real(mpmath.mpf(0), 0, accuracy10) +def _precision_complexes(numbers: tuple[Number, ...]) -> bool: + """ + True if ``numbers`` has a Complex, every element is an exact real + number, a PrecisionReal, or a Complex with such parts, and one of them + is a PrecisionReal. As in WMA, the parts of these complex numbers keep + their own errors, and sums and products of them are worked out from + sums and products of the parts. + """ + has_complex = has_precision_real = False + for number in numbers: + if isinstance(number, Complex): + has_complex = True + parts: tuple[Number, ...] = (number.real, number.imag) + else: + parts = (number,) + for part in parts: + if isinstance(part, PrecisionReal): + has_precision_real = True + elif not isinstance(part, (Integer, Rational)): + return False + return has_complex and has_precision_real + + +def _real_and_imag(number: Number) -> tuple[Number, Number]: + if isinstance(number, Complex): + return number.real, number.imag + return number, Integer0 + + +def _add_complex_numbers(numbers: tuple[Number, ...]) -> Number: + """The sum of ``numbers``: the real parts and the imaginary parts add up.""" + parts = [_real_and_imag(number) for number in numbers] + real = eval_add_numbers(*(part[0] for part in parts)) + imag = eval_add_numbers(*(part[1] for part in parts)) + return Complex(real, imag) + + +def _times_complex(x: Number, y: Number) -> Number: + """``x y``, with x or y a Complex: (a + b I)(c + d I) is ac - bd + (ad + bc) I""" + if not isinstance(x, Complex): + return Complex( + eval_multiply_numbers(x, y.real), eval_multiply_numbers(x, y.imag) + ) + if not isinstance(y, Complex): + return Complex( + eval_multiply_numbers(x.real, y), eval_multiply_numbers(x.imag, y) + ) + a, b = x.real, x.imag + c, d = y.real, y.imag + real = eval_add_numbers( + eval_multiply_numbers(a, c), eval_negate_number(eval_multiply_numbers(b, d)) + ) + imag = eval_add_numbers(eval_multiply_numbers(a, d), eval_multiply_numbers(b, c)) + return Complex(real, imag) + + +def _square_complex(z: Complex) -> Number: + """``z^2``: as in WMA, (a + b I)^2 is (a + b)(a - b) + 2 a b I.""" + a, b = z.real, z.imag + real = eval_multiply_numbers( + eval_add_numbers(a, b), eval_add_numbers(a, eval_negate_number(b)) + ) + imag = eval_multiply_numbers(Integer(2), eval_multiply_numbers(a, b)) + return Complex(real, imag) + + +def _multiply_complex_numbers(numbers: tuple[Number, ...]) -> Number: + """ + The product of ``numbers``. As in WMA, the factors are multiplied as a + balanced tree, the first half times the second half, and a pair of the + same complex number is squared. The errors depend on the order: + ``(z w) u`` and ``z (w u)`` differ in their last digits. + """ + if len(numbers) == 1: + return numbers[0] + if ( + len(numbers) == 2 + and isinstance(numbers[0], Complex) + and numbers[0].sameQ(numbers[1]) + ): + return _square_complex(numbers[0]) + half = (len(numbers) + 1) // 2 + x = _multiply_complex_numbers(numbers[:half]) + y = _multiply_complex_numbers(numbers[half:]) + if isinstance(x, Complex) or isinstance(y, Complex): + return _times_complex(x, y) + return eval_multiply_numbers(x, y) + + +def _inverse_real(x: Number) -> Optional[Number]: + """``1/x`` for an exact real number or a nonzero PrecisionReal.""" + if not isinstance(x, PrecisionReal): + return eval_inverse_number(x) + if x.is_zero: + return None + with mpmath.workprec(_stored_bits((x,))): + value = 1 / _full_mpmath(x) + return _precision_real(value, x.precision10, 0) + + +def _inverse_complex(z: Complex) -> Optional[Number]: + """``1/z``: as in WMA, (a - b I) / (a^2 + b^2).""" + a, b = z.real, z.imag + norm = eval_add_numbers(eval_multiply_numbers(a, a), eval_multiply_numbers(b, b)) + inverse_norm = _inverse_real(norm) + if inverse_norm is None: + return None + return Complex( + eval_multiply_numbers(a, inverse_norm), + eval_negate_number(eval_multiply_numbers(b, inverse_norm)), + ) + + +def eval_Power_by_parts(base: Number, exponent: Number) -> Optional[Number]: + """ + ``base^exponent`` where, as in WMA, the parts of the complex result get + their errors one by one: a Complex with exact or PrecisionReal parts to + an Integer power, and a negative PrecisionReal to a Rational power. + Returns None otherwise. + """ + if isinstance(exponent, Integer): + if isinstance(base, Complex) and _precision_complexes((base,)): + return _complex_integer_power(base, exponent.value) + return None + if ( + isinstance(exponent, Rational) + and isinstance(base, PrecisionReal) + and not base.is_zero + and base.value < 0 + ): + # As in WMA, x^(7/3) is x^2 x^(1/3), and x^(2/3) is (x^(1/3))^2. + p, q = (int(n) for n in exponent.value.as_numer_denom()) + sign = 1 if p > 0 else -1 + whole, rest = divmod(abs(p), q) + root = _negative_real_root(base, sign, q) + if rest > 1: + root = _complex_integer_power(root, rest) + if whole == 0 or root is None: + return root + with mpmath.workprec(_stored_bits((base,))): + power = _full_mpmath(base) ** (sign * whole) + factor = _precision_real(power, base.precision10 - math.log10(whole), 0) + return _times_complex(factor, root) + return None + + +def _negative_real_root(x: PrecisionReal, sign: int, q: int) -> BaseElement: + """ + ``x^(sign/q)`` for a negative ``x``. As in WMA, its error is that of the + whole number, Abs[x^(1/q)] 10^-precision / q, shared by the parts as in + ``_complex_with_error``. + """ + precision10 = x.precision10 + math.log10(q) + bits = max(_stored_bits((x,)), number_prec(math.ceil(precision10) + 1)) + 10 + with mpmath.workprec(bits): + magnitude = mpmath.root(-_full_mpmath(x), q) ** sign + if q == 2: + value = mpmath.mpc(0, sign * magnitude) + else: + angle = sign * mpmath.pi / q + value = mpmath.mpc( + magnitude * mpmath.cos(angle), magnitude * mpmath.sin(angle) + ) + return _complex_with_error(value, magnitude * mpmath.mpf(10) ** -precision10) + + +def _complex_integer_power(z: Complex, n: int) -> Optional[Number]: + """ + ``z^n`` for a Complex ``z`` whose parts are exact or PrecisionReals. + As in WMA, z^n is worked out by squaring from the highest bit of n + down, multiplying by z at each bit that is set, and z^-n is 1/z^n. + Returns None if z^n is 0 within its error, so that 1/z^n has no value. + """ + if n == 0: + return Integer1 + result: Number = z + for bit in bin(abs(n))[3:]: + if isinstance(result, Complex): + result = _square_complex(result) + else: + result = eval_multiply_numbers(result, result) + if bit == "1": + result = _times_complex(result, z) + if n > 0: + return result + if isinstance(result, Complex): + return _inverse_complex(result) + return _inverse_real(result) + + def eval_negate_number(n: Number) -> Number: """ Changes the sign of n diff --git a/mathics/format/box/formatvalues.py b/mathics/format/box/formatvalues.py index 262679cc4..5b73eed34 100644 --- a/mathics/format/box/formatvalues.py +++ b/mathics/format/box/formatvalues.py @@ -205,7 +205,9 @@ def do_format_complex( return None parts: List[Any] = [] - if element.is_machine_precision or not element.real.is_zero: + # As in WMA, only an exact zero real part is left out: 0``19.4 + 2.`19.7 I + # shows its zero. + if not element.real.sameQ(Integer(0)): parts.append(element.real) if element.imag.sameQ(Integer(1)): parts.append(SymbolI) diff --git a/test/builtin/atomic/test_numbers.py b/test/builtin/atomic/test_numbers.py index f711f1d23..033229e0c 100644 --- a/test/builtin/atomic/test_numbers.py +++ b/test/builtin/atomic/test_numbers.py @@ -391,6 +391,177 @@ def test_real_valued_precision(str_expr, str_expected, msg): check_evaluation(str_expr, str_expected, failure_message=msg) +# Arguments of Through: the precision and the accuracy of each part. +PARTS = "{Precision@*Re, Accuracy@*Re, Precision@*Im, Accuracy@*Im}" +ACCURACIES = "{Accuracy@*Re, Accuracy@*Im}" +Z = "Complex[1.`20, 2.`25]" +W = "Complex[3.`22, -1.`18]" + + +@pytest.mark.parametrize( + ("str_expr", "str_expected", "msg"), + [ + ( + "Precision[1.5`20.3 + 2.5`20.3 I]", + "20.3", + "the errors of the parts add in quadrature", + ), + ("Accuracy[1.5`20.3 + 2.5`20.3 I]", "19.8353", None), + ("Precision[Complex[1.`20, 1.`10]]", "10.1505", None), + ("Accuracy[Complex[1000.`20, 1.`20]]", "17.", None), + ("Precision[Complex[1.`20, 3]]", "20.5", "an exact part counts to Abs"), + ("Precision[2.`20 + I]", "20.0485", None), + ("{Re[2.`20 I], Im[2.`20 + I]}", "{0, 1}", "exact parts stay exact"), + ( + f"Through[{PARTS}[Complex[1.`20, 1.`20] + 1]]", + "{20.301, 20., 20., 20.}", + "Plus: each part adds its own errors", + ), + ( + f"Through[{PARTS}[Complex[1.`20, 1.`20] + 1.`10]]", + "{10.301, 10., 20., 20.}", + None, + ), + ( + f"Through[{PARTS}[{Z} + {W}]]", + "{20.5892, 19.9872, 18., 18.}", + None, + ), + ( + f"Through[{PARTS}[Complex[1.`20, 1.`20] Complex[1.`20, -1.`20]]]", + "{19.699, 19.3979, 0., 19.3979}", + "Times: (a + b I)(c + d I) is ac - bd + (ad + bc) I", + ), + ( + f"Through[{PARTS}[Complex[1.`20, 1.`20] 1.`10]]", + "{10., 10., 10., 10.}", + None, + ), + (f"Through[{ACCURACIES}[{Z} {W}]]", "{17.6924, 17.9954}", None), + ( + f"Through[{ACCURACIES}[{Z} {Z} {Z} {Z}]]", + "{18.2518, 18.2518}", + "the factors are multiplied as a balanced tree, (z z)(z z)", + ), + ( + f"Through[{ACCURACIES}[{Z} / {W}]]", + "{18.5184, 18.6176}", + "z / w is z (1 / w)", + ), + ( + f"Through[{PARTS}[Complex[1.`20, 1.`20]^2]]", + "{0., 19.3979, 19.699, 19.3979}", + "z^2 is (a + b)(a - b) + 2 a b I", + ), + ( + f"Through[{ACCURACIES}[Complex[1.`20, 2]^2]]", + "{19.3979, 19.3979}", + None, + ), + ( + f"Through[{ACCURACIES}[{Z}^3]]", + "{18.8239, 18.7959}", + "z^n by squaring, from the highest bit of n", + ), + (f"Through[{ACCURACIES}[{Z}^4]]", "{18.1938, 18.2518}", None), + (f"Through[{ACCURACIES}[{Z}^5]]", "{17.7375, 17.6819}", None), + (f"Through[{ACCURACIES}[{Z}^6]]", "{17.1662, 17.3851}", None), + ( + f"Through[{PARTS}[1 / Complex[1.`20, 1.`20]]]", + "{19.5229, 19.8239, 19.5229, 19.8239}", + "1 / z is (a - b I) / (a^2 + b^2)", + ), + (f"Through[{ACCURACIES}[{Z}^-2]]", "{20.3677, 20.2853}", None), + ( + f"Through[{ACCURACIES}[1 / Complex[0``10, 1.`20]]]", + "{10., 19.3979}", + None, + ), + ], +) +def test_complex_precision(str_expr, str_expected, msg): + """Complex numbers whose parts keep their own errors: compared with WMA 15""" + check_evaluation(str_expr, str_expected, failure_message=msg) + + +@pytest.mark.parametrize( + ("str_expr", "str_expected", "msg"), + [ + ( + "{Re[#], Precision[#]}&[ArcSec[0.5`20]]", + "{0, 20.0571}", + "a complex value at a real argument", + ), + ("{Re[#], Precision[#]}&[Sqrt[-2.`20]]", "{0, 20.301}", None), + ( + f"Through[{PARTS}[Log[-1.`20]]]", + "{0., 20., 36.7528, 36.2556}", + "a part that does not change gets 2^-54 of the error", + ), + (f"Through[{PARTS}[Log[-1.`2]]]", "{0., 2., 18.7528, 18.2556}", None), + ("Precision[Re[Log[-2.`20]]]", "19.8408", None), + ( + f"Through[{PARTS}[ArcSin[2.`20]]]", + "{36.3893, 36.1932, 20.0571, 19.9375}", + None, + ), + ("Precision[ArcSin[2.`20]]", "20.2492", None), + ("Precision[ArcCos[-2.`20]]", "20.4698", None), + ( + "Through[{Precision@*Re, Precision@*Im}[ArcTanh[2.`20]]]", + "{19.9159, 36.6278}", + None, + ), + ( + "Through[{Precision@*Re, Precision@*Im}[ArcCoth[0.5`20]]]", + "{19.9159, 36.6278}", + None, + ), + ("Precision[Im[ArcCosh[-2.`20]]]", "36.6903", None), + ("{Re[#], Precision[#]}&[ArcCosh[0.5`20]]", "{0, 20.2586}", None), + ("Precision[Im[ArcSin[1.`20*^10]]]", "21.3751", None), + ("Precision[ArcCos[1.5`5]]", "4.85573", None), + ( + "{Re[#], Precision[#]}&[ArcCosh[0``20]]", + "{0, 20.1961}", + "at a zero, the change of f within the error", + ), + ("{Re[#], Precision[#]}&[ArcCoth[0``20]]", "{0, 20.1961}", None), + ( + "{Re[#], Precision[#]}&[ArcCosh[-1.`20]]", + "{0, 10.3466}", + "at a branch point", + ), + ( + f"Through[{PARTS}[(-8.`20)^(1/3)]]", + "{20.3266, 20.3266, 20.5652, 20.3266}", + "(-x)^(1/q): the error of the whole, shared by the parts", + ), + ("Accuracy[Re[(-8.`20)^(-1/3)]]", "20.9287", None), + ( + "Through[{Precision@*Re, Precision@*Im}[(-2.`20)^(1/6)]]", + "{20.8662, 20.6276}", + None, + ), + ("{Re[#], Precision[#]}&[(-2.`20)^(-1/2)]", "{0, 20.301}", None), + ( + f"Through[{ACCURACIES}[(-8.`20)^(2/3)]]", + "{19.486, 19.5891}", + "(-x)^(2/3) is ((-x)^(1/3))^2", + ), + ( + "Through[{Precision@*Re, Precision@*Im}[(-8.`20)^(4/3)]]", + "{19.8323, 19.8955}", + "(-x)^(4/3) is (-x) (-x)^(1/3)", + ), + ("{Re[#], Precision[#]}&[(-8.`20)^(5/2)]", "{0, 19.6021}", None), + ], +) +def test_complex_valued_precision(str_expr, str_expected, msg): + """Functions with a complex value at a real argument: compared with WMA 15""" + check_evaluation(str_expr, str_expected, failure_message=msg) + + @pytest.mark.parametrize( ("str_expr", "str_expected", "msg"), [