Complex numbers: the parts keep their own errors, as in WMA

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(

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)))
 
 

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)&",
     }
 

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:

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

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(

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)