Log: a zero value gets the error of x, as in WMA

mathics/builtin/numbers/exp.py

@@ -15,14 +15,21 @@
 
 import mpmath
 
-from mathics.core.atoms import Real
+from mathics.core.atoms import Integer0, Integer1, Real
 from mathics.core.attributes import A_LISTABLE, A_NUMERIC_FUNCTION, A_PROTECTED
 from mathics.core.builtin import Builtin, MPMathFunction
 from mathics.core.convert.python import from_python
 from mathics.core.expression import Expression
 from mathics.core.number import min_prec
 from mathics.core.symbols import SymbolPower
-from mathics.core.systemsymbols import SymbolCos, SymbolE, SymbolSin
+from mathics.core.systemsymbols import (
+    SymbolComplexInfinity,
+    SymbolCos,
+    SymbolE,
+    SymbolIndeterminate,
+    SymbolSin,
+)
+from mathics.eval.numeric import covers_point
 
 
 class Fold:
@@ -228,6 +235,25 @@ def get_mpmath_function(self, args):
             return mpmath.log
         return lambda base, x: mpmath.log(x, base)
 
+    def special_point_value(self, args):
+        # As in WMA, an inexact 0 gives Indeterminate, as Log[2.`20, 0``20],
+        # and so does a base that is 1, or whose error covers 1, when x is
+        # or covers 1 too: Log[1.`20, 1.`20]. With any other x, Log[x] / 0 is
+        # ComplexInfinity, as Log[1.`20, 2]. Any other base gives Log[b, 1]
+        # an exact 0.
+        if any(isinstance(arg, Real) and arg.is_zero for arg in args):
+            return SymbolIndeterminate
+        if len(args) == 1:
+            return None
+        base, x = args
+        if _is_or_covers_one(base):
+            return (
+                SymbolIndeterminate if _is_or_covers_one(x) else SymbolComplexInfinity
+            )
+        if x.sameQ(Integer1):
+            return Integer0
+        return None
+
     @staticmethod
     def condition_numbers(*args: mpmath.mpf) -> tuple[mpmath.mpf, ...]:
         # The arguments, then the value. Log[b, x] is Log[x] / Log[b]. As in
@@ -239,6 +265,10 @@ def condition_numbers(*args: mpmath.mpf) -> tuple[mpmath.mpf, ...]:
         return (1 / abs(mpmath.log(base)), 1 / abs(mpmath.log(x)))
 
 
+def _is_or_covers_one(number) -> bool:
+    return number.sameQ(Integer1) or covers_point(number, (1,))
+
+
 class Log2(Builtin):
     """
     <url>:WMA link:https://reference.wolfram.com/language/ref/Log2.html</url>

mathics/core/builtin.py

@@ -839,6 +839,16 @@ def get_mpmath_function(self, args):
             return None
         return getattr(mpmath, self.mpmath_name)
 
+    def special_point_value(self, args: Sequence[Number]) -> Optional[BaseElement]:
+        """
+        The value at numbers, not all exact, where f has no value of its
+        own, or None: Indeterminate where the error of the argument covers
+        one of ``indeterminate_at``.
+        """
+        if len(args) == 1 and covers_point(args[0], self.indeterminate_at):
+            return SymbolIndeterminate
+        return None
+
     def eval(self, z, evaluation: Evaluation):
         "%(name)s[z__]"
 
@@ -864,8 +874,9 @@ def eval(self, z, evaluation: Evaluation):
         if mpmath_function is None:
             return
 
-        if len(args) == 1 and covers_point(args[0], self.indeterminate_at):
-            return SymbolIndeterminate
+        special = self.special_point_value(args)
+        if special is not None:
+            return special
 
         if self.condition_numbers is not None:
             result = eval_mpmath_function_by_condition(

mathics/eval/numeric.py

@@ -728,8 +728,16 @@ def eval_mpmath_function_by_condition(
         if factors is not None and not all(mpmath.isfinite(f) for f in factors):
             factors = None
 
+    if factors is None and len(args) > 1:
+        if isinstance(value, mpmath.mpc) or value != 0:
+            return None
+        error = _change_within_args(mpmath_function, args, values)
+        if error is None:
+            return None
+        return _value_with_error(value, error, has_poles, bound)
+
     if factors is None:
-        if len(args) != 1 or args[0].precision10 < 1:
+        if args[0].precision10 < 1:
             return None
         x = values[0]
         digits = args[0].precision10 - _log10_abs(x)
@@ -889,16 +897,24 @@ def _change_within(
     (10^-20)^2 / 2, and ArcCos[1.`20] is 0``9.849, since ArcCos changes by
     Sqrt[2 10^-20].
 
-    The smaller change of the two sides is taken: where f is continuous they
-    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.
+    As in WMA, the change is the leading term of f at x, c delta^p, even
+    where delta is not small: Log[1.`3] is 0``3., not the 0``3.0002 of
+    Log[1.001]. The power p is read off the change within two much smaller
+    steps.
+
+    On each side, the smaller change is taken: where f is continuous the
+    sides 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)):
+    # Enough digits for a change of second order within the small steps.
+    shrink_bits = 40
+    with mpmath.workprec(
+        number_prec(math.ceil(2 * (max(digits, 1) + shrink_bits * math.log10(2))) + 20)
+    ):
         x = mpmath.mpf(x)
         try:
             value = mpmath_function(x)
@@ -911,17 +927,62 @@ def _change_within(
         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, kinds) and mpmath.isfinite(other):
-                changes.append(abs(other - value))
-        if not changes or min(changes) == 0:
+
+        def change(step: mpmath.mpf) -> Optional[mpmath.mpf]:
+            changes = []
+            for side in (x + step, x - step):
+                try:
+                    other = mpmath_function(side)
+                except (ZeroDivisionError, ValueError):
+                    continue
+                if isinstance(other, kinds) and mpmath.isfinite(other):
+                    changes.append(abs(other - value))
+            return min(changes) if changes else None
+
+        whole = change(delta)
+        if not whole:
             return None
-        return value, min(changes)
+        step = mpmath.ldexp(delta, -shrink_bits)
+        small, smaller = change(step), change(step / 2)
+        if not small or not smaller:
+            return value, whole
+        # A power like 1, 2 or 1/2 is taken as exact.
+        power = mpmath.log(small / smaller, 2)
+        if abs(2 * power - mpmath.nint(2 * power)) < 0.01:
+            power = mpmath.nint(2 * power) / 2
+        return value, small * mpmath.power(2, shrink_bits * power)
+
+
+def _change_within_args(
+    mpmath_function: Callable, args: tuple[Number, ...], values: list[mpmath.mpf]
+) -> Optional[mpmath.mpf]:
+    """
+    How much f changes within the errors of its arguments: the change
+    within the error of each PrecisionReal argument, the others held fixed,
+    added up. So Log[2.`20, 1.`20] is 0``19.84: Log[2, x] changes by
+    10^-20 / Log[2] at x = 1, and Log[b, 1] not at all. Returns None if
+    no argument changes f.
+    """
+    total = mpmath.mpf(0)
+    for k, arg in enumerate(args):
+        if not isinstance(arg, PrecisionReal):
+            continue
+        if arg.precision10 < 1:
+            return None
+        x = values[k]
+
+        def partial(t, k=k):
+            return mpmath_function(*values[:k], t, *values[k + 1 :])
+
+        change = _change_within(
+            partial,
+            x,
+            abs(x) * mpmath.mpf(10) ** -arg.precision10,
+            arg.precision10 - _log10_abs(x),
+        )
+        if change is not None:
+            total += change[1]
+    return total or None
 
 
 def _value_with_error(

test/builtin/atomic/test_numbers.py

@@ -391,6 +391,42 @@ def test_real_valued_precision(str_expr, str_expected, msg):
     check_evaluation(str_expr, str_expected, failure_message=msg)
 
 
+@pytest.mark.parametrize(
+    ("str_expr", "str_expected", "msg"),
+    [
+        ("Accuracy[Log[1.`20]]", "20.", "the error of x, divided by x"),
+        ("Accuracy[Log[1.`3]]", "3.", "the leading term, also at low precision"),
+        ("Accuracy[Log[2, 1.`20]]", "19.8408", "and by Log[b]"),
+        ("Accuracy[Log[10, 1.`20]]", "20.3622", None),
+        ("Accuracy[Log10[1.`20]]", "20.3622", None),
+        ("Accuracy[Log2[1.`20]]", "19.8408", None),
+        ("Accuracy[Log[2.`20, 1.`25]]", "24.8408", "the base does not change 0"),
+        ("Accuracy[Log[2.`25, 1.`20]]", "19.8408", None),
+        ("Accuracy[Log[1.0001`20, 1.`20]]", "16.", None),
+        ("Log[2.`20, 1]", "0", "an exact 1 gives an exact 0"),
+        ("Log[1.`20, 1.`20]", "Indeterminate", "the base covers 1, and x too"),
+        ("Log[1.0001`3, 1.`20]", "Indeterminate", None),
+        ("Log[1, 1.`20]", "Indeterminate", None),
+        ("Log[1.`20, 1]", "Indeterminate", None),
+        ("Log[1., 1.]", "Indeterminate", None),
+        ("Log[1.`20, 2]", "ComplexInfinity", "the base covers 1"),
+        ("Log[1.0001`3, 2]", "ComplexInfinity", None),
+        ("Log[1, 2.`20]", "ComplexInfinity", None),
+        ("Log[1.`20, 2 + I]", "ComplexInfinity", None),
+        ("Precision[Log[1.0001`5, 2]]", "0.999978", "the base does not cover 1"),
+        ("Log[2.`20, 0``20]", "Indeterminate", "an inexact 0"),
+        ("Log[0``20]", "Indeterminate", None),
+        ("Log[0., 2.`20]", "Indeterminate", None),
+        ("Accuracy[Sin[0``1.5]]", "1.5", "the leading term of f"),
+        ("Accuracy[Cos[0``2]]", "4.30103", None),
+        ("Accuracy[ArcCos[1.`3]]", "1.34949", None),
+    ],
+)
+def test_log_zero_precision(str_expr, str_expected, msg):
+    """Compared with WMA 15"""
+    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}"
@@ -823,7 +859,7 @@ def test_trig_zero(str_expr, str_expected, msg):
         ("Precision[ArcSin[0.0`20]]", "MachinePrecision", "0.0`20 is 0."),
         ("Accuracy[ArcCos[1.`20]]", "9.84949", "ArcCos changes by Sqrt[2 10^-20]"),
         ("Accuracy[ArcSec[1.`20]]", "9.84949", None),
-        ("Accuracy[ArcCos[1.`5]]", "2.34948", None),
+        ("Accuracy[ArcCos[1.`5]]", "2.34949", None),
         ("Precision[ArcSec[-1.`20]]", "10.3466", None),
         ("Precision[ArcCos[0.99999999`20]]", "12.301", "next to the end"),
         ("Precision[ArcSin[0.99999999`20]]", "16.3466", None),