From 564c05f91498d05d7c0c33f964191fc16e8a96fd Mon Sep 17 00:00:00 2001 From: Stanislav Date: Wed, 28 Oct 2026 19:56:00 +0400 Subject: [PATCH] Add a page on how precision and accuracy are implemented How the kernel tracks the error of arbitrary-precision numbers, as WMA does: PrecisionReal's decimal precision and 64-bit value, the error rules for sums and products, condition numbers for functions, zero values, special points, complex parts, the second pass of N, and SetPrecision. Ends with how to add precision to a new function. The user-level description is in the Mathics3 Manual (Mathics3/mathics-core#2069). Co-Authored-By: Claude Opus 5.5 --- docs/code-overview.rst | 1 + docs/code-overview/precision-and-accuracy.rst | 276 ++++++++++++++++++ 2 files changed, 277 insertions(+) create mode 100644 docs/code-overview/precision-and-accuracy.rst diff --git a/docs/code-overview.rst b/docs/code-overview.rst index 30214c4..f98a03f 100644 --- a/docs/code-overview.rst +++ b/docs/code-overview.rst @@ -36,6 +36,7 @@ Each of the above steps can be involved, so we break these down below. code-overview/evaluation code-overview/pattern-matching code-overview/rules + code-overview/precision-and-accuracy code-overview/boxes-forms-and-formats code-overview/python-modules code-overview/object-classes diff --git a/docs/code-overview/precision-and-accuracy.rst b/docs/code-overview/precision-and-accuracy.rst new file mode 100644 index 0000000..ea2df7d --- /dev/null +++ b/docs/code-overview/precision-and-accuracy.rst @@ -0,0 +1,276 @@ +.. index:: precision, accuracy, PrecisionReal + +====================== +Precision and Accuracy +====================== + +As in the Wolfram Language, every arbitrary-precision number in +Mathics3 carries its error, and every computation works out the error +of its result from the errors of its inputs. How this looks to a user +is described in the section "Precision and Accuracy" of the Mathics3 +Manual (``mathics/doc/documentation/1-Manual.mdoc``). This page +describes how it is implemented. + +The aim throughout is to give the same precision and accuracy as WMA. +The tests that check this are in ``test/builtin/atomic/test_numbers.py`` +(``test_real_valued_precision``, ``test_power_real_valued_precision``, +``test_trig_real_valued_precision``, ``test_complex_precision``, +``test_log_zero_precision`` and others), and their expected values +come from WMA 15. When you change this code, take new expected values +from WMA, not from what Mathics3 gives. + +Precision is the relative error in decimal digits, ``-Log10`` of the +relative error, and accuracy is the absolute error in decimal digits. +So for a nonzero ``x``, ``Precision[x] == Accuracy[x] + Log10[Abs[x]]``. + +Numbers +======= + +``MachineReal`` is a Python float. It has no error to track: an +expression with a machine number anywhere in it has ``MachinePrecision`` +(``has_machine_number`` in ``mathics/eval/numbers/numbers.py``), and +arithmetic with a machine number gives a machine number. + +``PrecisionReal`` (``mathics/core/atoms/numerics.py``) holds a +``sympy.Float``, whose precision is a whole number of binary digits. +Two slots add what WMA keeps and the ``sympy.Float`` does not: + +``_precision10`` + The precision in decimal digits, which can be any real number, as + 20.3 in ``1.5`20.3``. It is ``None`` when it is the precision that the + binary digits of the ``sympy.Float`` already give, so that most + numbers carry nothing extra. For a zero, it is the accuracy. + +``_extended`` + The value with ``MIN_STORED_BITS`` (64) bits + (``mathics/core/number.py``), kept when the precision needs fewer. + WMA keeps about this many bits of a low-precision number, and the + difference shows when numbers cancel: ``1.4`2`` would otherwise be + stored as 1.400390625. + +They are read through these properties and methods: + +* ``precision`` is the binary precision, as it was before these slots + existed. Most builtins use it, and it is right for them. +* ``precision10`` is the decimal precision (or a zero's accuracy). + ``Precision`` and ``Accuracy`` return it, and the error-tracking code + works with it. +* ``value`` is the value rounded to the precision. ``extended_value`` + is the value with all the bits kept, and it is what the error-tracking + code computes with. +* ``PrecisionReal.from_value(value, precision10)`` makes a number from + an mpmath or sympy value with enough bits: it rounds the value to the + precision and keeps the 64-bit value if that has more bits. Use it + rather than the constructor when making a number from a computed + value. + +The object cache is keyed by the ``sympy.Float``, ``_precision10`` and +``_extended``. Two numbers that differ in any of these are different +objects, though they can still be ``SameQ`` and ``Equal``. The cache +key of ``Complex`` includes the parts' ``_precision10`` for the same +reason. + +``OutputForm`` shows as many digits as the rounded precision. +``InputForm`` and ``FullForm`` show the shortest decimal string that +reads back as the kept value, followed by the stored precision, as +``3.1415926535897932385`3.``. + +Where errors are worked out +=========================== + +Sums and products +----------------- + +``eval_add_numbers`` and ``eval_multiply_numbers`` are in +``mathics/eval/numeric.py``. ``Plus`` and ``Times`` use them, and so +do ``Abs``, ``Sign`` and other builtins. When every number is exact or +a ``PrecisionReal`` (``_exact_or_precision_reals``): + +* In a sum, the absolute errors add up: the accuracy of the result is + ``-Log10[Total[10^-acc_i]]``. In a product, the relative errors add up + in the same way. ``_combine_digits`` does this sum. +* The value is computed from ``extended_value`` with ``_stored_bits`` + bits: the most bits any input keeps, and at least 64. +* ``_precision_real`` makes the result. When no digit is left, it gives + a zero whose accuracy is the error. +* An inexact zero counts in a sum with its absolute error. In a + product it gives a zero (``_multiply_by_zeros``), and the size of the + other factors goes into the zero's accuracy. An exact ``0`` times + anything gives an exact ``0``. + +Complex numbers with exact or ``PrecisionReal`` parts +(``_precision_complexes``) are added and multiplied part by part with +the same rules, so each part keeps its own error. ``Times`` multiplies +as a balanced tree, ``z^n`` squares starting from the top bit of ``n``, +and ``1/z`` is ``(a - b I)/(a^2 + b^2)``. This is the order WMA uses, +and it shows in the last digits of the precision. + +Any other combination falls back to the old rule: the result gets the +lowest precision of the inputs (``min_prec``). + +Exact constants +--------------- + +``Plus``, ``Times`` and ``Power`` first call ``numerify`` +(``mathics/eval/numerify.py``), which turns exact numeric elements such +as ``E`` or ``Sqrt[2]`` into numbers. It gives them ``GUARD_DIGITS`` +(16) digits more than the least precise inexact element. Their error +then vanishes next to the error of the inexact elements, so +``N[Pi, 30] E`` keeps precision 30, as in WMA. With a machine number +present, no digits are added. + +Functions: the condition number +------------------------------- + +The relative error of ``f[x]`` is the relative error of ``x`` times the +condition number ``Abs[x f'[x] / f[x]]``. For several arguments the +errors add up. A function opts in by giving ``MPMathFunction`` a static +method ``condition_numbers`` (``mathics/core/builtin.py``). The method +is given the arguments and then the value, all as mpmath numbers, and +returns one factor per argument. For ``Sin``: + +.. code-block:: python + + @staticmethod + def condition_numbers(x: mpmath.mpf, value: mpmath.mpf) -> tuple[mpmath.mpf]: + return (abs(x * mpmath.cos(x) / value),) + +and for ``Power``, ``(Abs[e], Abs[e Log[b]])``. The factors must stay +mpmath numbers, because a tiny factor would become ``0.`` as a float. + +``eval_mpmath_function_by_condition`` (``mathics/eval/numeric.py``) +does the rest. It computes the value with the stored bits plus a few +guard bits. If the result has more digits than the inputs, as for +``Sqrt``, it computes the value again with enough bits. It takes these +class attributes into account: + +``has_poles`` + When no digit is left, the result is ``ComplexInfinity``, as for + ``Tan[1.5707963`5]``. + +``value_bound`` + The error is no larger than this bound, as for ``Sin[1.`20*^30]``. + +``complex_values`` + When ``True`` (the default), a complex value at real arguments is + handled by ``_complex_value_by_parts``, described below. + +Arguments that are not exact reals or ``PrecisionReal`` numbers, for +example machine numbers or complex arguments, take the older path in +``MPMathFunction.eval``: the arguments are rounded to the lowest +precision, and the function is computed at that precision. + +Zero values and zero arguments +------------------------------ + +Where the condition number does not exist, the error is how much ``f`` +changes within the error of its argument (``_change_within``). That is +the case when the value is 0, when the argument is an inexact zero +(which has no relative error), or at the end of the domain, as for +``ArcCos[1.`20]``, where ``f'[x]`` is infinite. For example, +``Cos[0``20]`` is 1 with precision 40.301, since ``Cos`` changes only by +``(10^-20)^2 / 2``. + +As in WMA, the change is the leading Taylor term ``c δ^p`` evaluated at +the size of the error ``δ``, even where ``δ`` is not small. The change +is measured at ``δ 2^-40`` and at half of that, and ``p`` is read off +the two and rounded to a multiple of 1/2. ``c δ^p`` is then +extrapolated from them. Measuring at ``δ`` itself is wrong at low +precision. On each side the smaller change is taken: where ``f`` is +continuous the two sides agree, and a jump, as ``ArcCot`` has at 0, is +on one side only. + +With several arguments (``_change_within_args``), each inexact argument +is varied within its own error while the others are held fixed, and +the changes add up. So ``Log[2.`20, 1.`20]`` gets its error only from +``x``. + +Special points +-------------- + +``MPMathFunction.special_point_value(args)`` is called before anything +is computed. It returns a value where ``f`` has none of its own, or +``None`` otherwise. By default it returns ``Indeterminate`` when the error +of a single argument covers one of ``indeterminate_at``. For example, +``ArcTanh`` has ``indeterminate_at = (1, -1)``, so ``ArcTanh[1.0000001`5]`` +is ``Indeterminate``. ``covers_point`` makes this test. A machine number +has no error, so it covers a point only if it is that point. + +A function can override ``special_point_value``. ``Log`` +(``mathics/builtin/numbers/exp.py``) does: an inexact 0 in either +argument gives ``Indeterminate``, a base that is or covers 1 gives +``ComplexInfinity`` (or ``Indeterminate`` when ``x`` is or covers 1 too), +and ``Log[b, 1]`` with an exact 1 is an exact 0. + +Complex values at real arguments +-------------------------------- + +For ``ArcSin[2.`20]``, ``Log[-1.`20]`` and similar values, +``_complex_value_by_parts`` gives each part of the value its own error. +For each argument, the part's error is the argument's relative error +times that part of ``x f'[x]``, which is taken from a central +difference. A part that does not change with the arguments, as the +``Pi`` of ``Log[-1.`20]``, still gets ``CONSTANT_PART_ERROR`` (``2^-54``) +times ``Abs[x f'[x]]``, as in WMA. A part whose value is 0 and does not +change is an exact 0. + +N +- + +``N[expr, p]`` (``eval_NValues_to_precision`` in +``mathics/eval/nevaluator.py``) first evaluates at ``p``. If a nonzero +number in the result has fewer than ``p`` digits, it evaluates again +with the lost digits added. It stops at ``$MaxExtraPrecision`` (50 by +default) digits more than ``p``, and then gives ``N::meprec``. It goes +on while some number that is short of ``p`` gains digits. A number that +gains nothing comes from an input with fewer than ``p`` digits, and +keeps its precision. Finally, every number with more than ``p`` digits +is set to ``p``. + +N rules see the working precision. So a retry is used only if +``_same_form`` holds: the retry must have the same form as the first +result, up to the values of its numbers. Otherwise, a rule given for +``p`` alone could be lost. + +``N[expr, {p, a}]`` gives each number the precision +``min(p, a + Log10[Abs[x]])``. + +SetPrecision and SetAccuracy +---------------------------- + +They are in ``mathics/eval/numbers/precision.py``. An inexact number +keeps its binary value, so ``SetPrecision[0.1, 20]`` shows the digits of +the double. With precision ``Infinity`` a number becomes the exact +rational it stores. The numbers are set first and the expression is +evaluated afterwards, so exact exponents stay exact. + +Adding precision to a new function +================================== + +For an ``MPMathFunction`` with a real value: + +#. Give it ``condition_numbers``: the factors ``Abs[x D[f, x] / f]``, as + mpmath numbers. Where ``f`` is small, use the value you are given + rather than computing ``f`` again, so that the digits lost near a zero + are measured at full precision. +#. Set ``has_poles`` or ``value_bound`` if they apply. +#. If an inexact argument should give ``Indeterminate`` or another value + at some points, set ``indeterminate_at`` or override + ``special_point_value``. +#. Take the expected values for the tests from WMA, at several points, + near zeros and poles of the function, and for an inexact zero + argument. + +Not done yet +============ + +These cases still take the older path, and their precision can differ +from WMA: + +* functions of a complex argument, such as ``Sqrt[1.`20 + I]`` and + ``Exp[1.`20 I]``; +* a negative base with an inexact exponent; +* ``N`` of an exact complex number whose two parts are nonzero (WMA splits + the precision between the parts); +* ``Power`` of an inexact zero, such as ``Sqrt[0``3]``; +* a precision below 1, as in ``Log[1.`0.5]``.