From 5fd3a7465c3a45e23e3624a33857dca80bd7fdfd Mon Sep 17 00:00:00 2001 From: mmatera Date: Sat, 4 Feb 2023 11:15:51 -0300 Subject: [PATCH 1/3] fix accuracy and precision for complex numbers --- mathics/builtin/atomic/numbers.py | 9 +++++---- mathics/core/atoms.py | 6 ++++-- mathics/core/convert/mpmath.py | 8 ++++++-- mathics/eval/numbers.py | 12 +++++++++--- test/builtin/atomic/test_numbers.py | 2 +- 5 files changed, 25 insertions(+), 12 deletions(-) diff --git a/mathics/builtin/atomic/numbers.py b/mathics/builtin/atomic/numbers.py index 2737c5401..d1d9cbe00 100644 --- a/mathics/builtin/atomic/numbers.py +++ b/mathics/builtin/atomic/numbers.py @@ -177,10 +177,11 @@ class Accuracy(Builtin): >> Accuracy[A] = Infinity - For Complex numbers, the accuracy is the smaller of the accuracies of its \ - real and imaginary parts: - >> Accuracy[1.00`2 + 2.00`2 I] - = 1. + For Complex numbers, the accuracy is estimated as (minus) the base-10 log + of the square root of the squares of the errors on the real and complex parts: + >> z=Complex[3.00``2, 4..00``2]; + >> Accuracy[z] == -Log[10, Sqrt[10^(-2 Accuracy[Re[z]]) + 10^(-2 Accuracy[Im[z]])]] + = True Accuracy of expressions is given by the minimum accuracy of its elements: >> Accuracy[F[1, Pi, A]] diff --git a/mathics/core/atoms.py b/mathics/core/atoms.py index 99cf37ab7..68f0e8542 100644 --- a/mathics/core/atoms.py +++ b/mathics/core/atoms.py @@ -297,10 +297,11 @@ class Real(Number): # __new__ rather than __init__ is used here because the kind of # object created differs based on contents of "value". - def __new__(cls, value, p=None) -> "Real": + def __new__(cls, value, p: int = None) -> "Real": """ Return either a MachineReal or a PrecisionReal object. - Or raise a TypeError + Or raise a TypeError. + p is the number of binary digits of precision. """ if isinstance(value, str): value = str(value) @@ -325,6 +326,7 @@ def __new__(cls, value, p=None) -> "Real": if p is None or p == FP_MANTISA_BINARY_DIGITS: return MachineReal.__new__(MachineReal, value) else: + # TODO: check where p is set in value: return PrecisionReal.__new__(PrecisionReal, value) def __eq__(self, other) -> bool: diff --git a/mathics/core/convert/mpmath.py b/mathics/core/convert/mpmath.py index 87739f714..28d4ed4e2 100644 --- a/mathics/core/convert/mpmath.py +++ b/mathics/core/convert/mpmath.py @@ -7,6 +7,7 @@ import sympy from mathics.core.atoms import Complex, MachineReal, MachineReal0, PrecisionReal +from mathics.core.number import LOG2_10 from mathics.core.symbols import Atom @@ -36,8 +37,11 @@ def from_mpmath( # is compatible with 0. if prec < 1.0: return MachineReal0 - # HACK: use str here to prevent loss of precision - return PrecisionReal(sympy.Float(str(value), prec)) + # HACK: use str here to prevent loss of precision. + # + # sympy.Float internally work with binary precision instead of + # decimal precision. TODO: avoid the irrelevant conversion + return PrecisionReal(sympy.Float(str(value), precision=(LOG2_10 * prec + 2))) elif isinstance(value, mpmath.mpc): if mpmath.isinf(value): return SymbolComplexInfinity diff --git a/mathics/eval/numbers.py b/mathics/eval/numbers.py index 69a8371c4..7389ac8d3 100644 --- a/mathics/eval/numbers.py +++ b/mathics/eval/numbers.py @@ -53,7 +53,8 @@ def eval_Accuracy(z: BaseElement) -> Optional[float]: return acc_imag if acc_imag is None: return acc_real - return min(acc_real, acc_imag) + + return -mpmath.log(10 ** (-2 * acc_real) + 10 ** (-2 * acc_imag), 10.0) * 0.5 if isinstance(z, Expression): elem_accuracies = (eval_Accuracy(z_elem) for z_elem in z.elements) @@ -91,11 +92,16 @@ def eval_Precision(z: BaseElement) -> Optional[float]: if isinstance(z, Complex): prec_real = eval_Precision(z.real) prec_imag = eval_Precision(z.imag) - if prec_real is None: + if prec_real is None or prec_imag == prec_real: return prec_imag if prec_imag is None: return prec_real - return min(prec_real, prec_imag) + # both numbers have different precision. + # Evaluate the accuracy and add the log of + # the module. + 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 isinstance(z, Expression): elem_prec = (eval_Precision(z_elem) for z_elem in z.elements) diff --git a/test/builtin/atomic/test_numbers.py b/test/builtin/atomic/test_numbers.py index da539617c..ebc15da91 100644 --- a/test/builtin/atomic/test_numbers.py +++ b/test/builtin/atomic/test_numbers.py @@ -213,7 +213,7 @@ def test_n(): # For some reason, the following test # would fail in WMA ("1. I", "Accuracy[1.]"), - (" 0.4 + 2.4 I", "Accuracy[2.4]"), + (" 0.4 + 2.4 I", "$MachinePrecision-Log[10, Abs[.4+2.4 I]]"), ("2 + 3 I", "Infinity"), ('"abc"', "Infinity"), # Returns the accuracy of ``` 3.2`3 ``` From aedc2657b6781a2a61389c56341cff9f07818493 Mon Sep 17 00:00:00 2001 From: mmatera Date: Sat, 4 Feb 2023 11:22:55 -0300 Subject: [PATCH 2/3] improving comment --- mathics/core/convert/mpmath.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/mathics/core/convert/mpmath.py b/mathics/core/convert/mpmath.py index 28d4ed4e2..03bce491c 100644 --- a/mathics/core/convert/mpmath.py +++ b/mathics/core/convert/mpmath.py @@ -39,8 +39,9 @@ def from_mpmath( return MachineReal0 # HACK: use str here to prevent loss of precision. # - # sympy.Float internally work with binary precision instead of - # decimal precision. TODO: avoid the irrelevant conversion + # sympy.Float internally works with (rounded) binary precision + # instead of decimal precision. Let's use it to avoid rounding errors. + # TODO: avoid the irrelevant conversion return PrecisionReal(sympy.Float(str(value), precision=(LOG2_10 * prec + 2))) elif isinstance(value, mpmath.mpc): if mpmath.isinf(value): From 42699d168b8f30558be62e037c45eb1ba81597f9 Mon Sep 17 00:00:00 2001 From: Juan Mauricio Matera Date: Sun, 5 Feb 2023 00:17:08 -0300 Subject: [PATCH 3/3] Update mpmath.py --- mathics/core/convert/mpmath.py | 1 - 1 file changed, 1 deletion(-) diff --git a/mathics/core/convert/mpmath.py b/mathics/core/convert/mpmath.py index 6f6a3331e..f6e2a885e 100644 --- a/mathics/core/convert/mpmath.py +++ b/mathics/core/convert/mpmath.py @@ -7,7 +7,6 @@ import sympy from mathics.core.atoms import Complex, MachineReal, MachineReal0, PrecisionReal -from mathics.core.number import LOG2_10 from mathics.core.symbols import Atom