From c1c378a016201943b4fadddfafd008ed79d624fc Mon Sep 17 00:00:00 2001 From: Stanislav Date: Sun, 25 Oct 2026 11:36:12 +0400 Subject: [PATCH] Log: a zero value gets the error of x, as in WMA MATH-258. Log[1.`20] is 0``20. and Log[b, 1.`20] is 0 with the accuracy of x plus Log10[Abs[Log[b]]], as in WMA: the change of f within the error of each argument, the others held fixed, added up. The base does not change a zero, so Log[2.`20, 1.`25] is 0``24.84. - A base that is 1, or whose error covers 1, gives ComplexInfinity, and Indeterminate when x is or covers 1 too; these raised ZeroDivisionError. An inexact 0 gives Indeterminate, and Log[b, 1] with an exact 1 is an exact 0. - The change within an error is the leading term c delta^p of f, as in WMA, so Log[1.`3] is 0``3. and Sin[0``1.5] is 0``1.5 rather than showing the next terms. - MPMathFunction.special_point_value takes the Indeterminate check of indeterminate_at, and Log adds its own points there. Co-Authored-By: Claude Opus 5.5 --- mathics/builtin/numbers/exp.py | 34 ++++++++++- mathics/core/builtin.py | 15 ++++- mathics/eval/numeric.py | 95 +++++++++++++++++++++++------ test/builtin/atomic/test_numbers.py | 38 +++++++++++- 4 files changed, 160 insertions(+), 22 deletions(-) diff --git a/mathics/builtin/numbers/exp.py b/mathics/builtin/numbers/exp.py index 87242c440..5777e7bab 100644 --- a/mathics/builtin/numbers/exp.py +++ b/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): """ :WMA link:https://reference.wolfram.com/language/ref/Log2.html diff --git a/mathics/core/builtin.py b/mathics/core/builtin.py index 2fd98afa9..12807ef07 100644 --- a/mathics/core/builtin.py +++ b/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( diff --git a/mathics/eval/numeric.py b/mathics/eval/numeric.py index 290043b7e..b815369d4 100644 --- a/mathics/eval/numeric.py +++ b/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( diff --git a/test/builtin/atomic/test_numbers.py b/test/builtin/atomic/test_numbers.py index 033229e0c..6b4a859d5 100644 --- a/test/builtin/atomic/test_numbers.py +++ b/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),