Merge pull request #775 from Mathics3/from_mpmath_precision

from_mpmath precision twaks

mathics/builtin/arithfns/basic.py

@@ -446,24 +446,31 @@ def append_last():
                     last_item = rest
                     last_count = count
         append_last()
-
         if numbers:
+            # TODO: reorganize de conditions to avoid compute unnecesary
+            # quantities. In particular, is we check mathine_precision,
+            # we do not need to evaluate prec.
             if prec is not None:
                 if is_machine_precision:
                     numbers = [item.to_mpmath() for item in numbers]
                     number = mpmath.fsum(numbers)
                     number = from_mpmath(number)
                 else:
-                    # For a sum, what is relevant is the minimum accuracy of the terms
-                    acc = (
-                        Expression(SymbolAccuracy, ListExpression(items))
-                        .evaluate(evaluation)
-                        .to_python()
-                    )
+                    # TODO: If there are Complex numbers in `numbers`,
+                    # and we are not working in machine precision, compute the sum of the real and imaginary
+                    # parts separately, to preserve precision. For example,
+                    # 1.`2 + 1.`3 I should produce
+                    # Complex[1.`2, 1.`3]
+                    # but with this implementation returns
+                    # Complex[1.`2, 1.`2]
+                    #
+                    # TODO: if the precision are not equal for each number,
+                    # we should estimate the result precision by computing the sum of individual errors
+                    # prec =  sum(abs(n.value) * 2**(-n.value._prec)    for n in number if n.value._prec is not None)/sum(abs(n))
                     with mpmath.workprec(prec):
                         numbers = [item.to_mpmath() for item in numbers]
                         number = mpmath.fsum(numbers)
-                        number = from_mpmath(number, acc=acc)
+                        number = from_mpmath(number, precision=prec)
             else:
                 number = from_sympy(sum(item.to_sympy() for item in numbers))
         else:
@@ -936,7 +943,8 @@ def eval(self, items, evaluation):
         elements = []
         numbers = []
         infinity_factor = False
-
+        # These quantities only have sense if there are numeric terms.
+        # Also, prec is only needed if is_machine_precision is not True.
         prec = min_prec(*items)
         is_machine_precision = any(item.is_machine_precision() for item in items)
 
@@ -1000,7 +1008,7 @@ def eval(self, items, evaluation):
                     with mpmath.workprec(prec):
                         numbers = [item.to_mpmath() for item in numbers]
                         number = mpmath.fprod(numbers)
-                        number = from_mpmath(number, dps(prec))
+                        number = from_mpmath(number, precision=prec)
             else:
                 number = sympy.Mul(*[item.to_sympy() for item in numbers])
                 number = from_sympy(number)

mathics/builtin/arithmetic.py

@@ -178,7 +178,7 @@ def eval(self, z, evaluation):
                     return
                 result = call_mpmath(mpmath_function, tuple(mpmath_args))
                 if isinstance(result, (mpmath.mpc, mpmath.mpf)):
-                    result = from_mpmath(result, d)
+                    result = from_mpmath(result, precision=prec)
         return result
 
 

mathics/builtin/atomic/numbers.py

@@ -177,10 +177,10 @@ 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 the smaller of the accuracies of its \
+    # real and imaginary parts:
+    # >> Accuracy[1.00`2 + 2.00`2 I]
+    # = 1.
 
     Accuracy of expressions is given by the minimum accuracy of its elements:
     >> Accuracy[F[1, Pi, A]]

mathics/builtin/binary/io.py

@@ -193,7 +193,7 @@ def _Real128_reader(s):
             else:
                 result = mpmath.fdiv(core, 2**-exp)
 
-            return from_mpmath(result, dps(112))
+            return from_mpmath(result, precision=112)
 
     @staticmethod
     def _TerminatedString_reader(s):

mathics/builtin/specialfns/bessel.py

@@ -157,7 +157,7 @@ def eval_N(self, k, precision, evaluation: Evaluation):
 
         with mpmath.workprec(p):
             result = mpmath.airyaizero(k_int)
-            return from_mpmath(result, d)
+            return from_mpmath(result, precision=p)
 
 
 class AiryBi(_MPMathFunction):
@@ -283,7 +283,7 @@ def eval_N(self, k: Integer, precision, evaluation: Evaluation):
 
         with mpmath.workprec(p):
             result = mpmath.airybizero(k_int)
-            return from_mpmath(result, d)
+            return from_mpmath(result, precision=p)
 
 
 class AngerJ(_Bessel):

mathics/builtin/specialfns/gamma.py

@@ -133,7 +133,7 @@ def eval_with_z(self, z, a, b, evaluation):
                     return
                 result = call_mpmath(mpmath_function, tuple(mpmath_args))
                 if isinstance(result, (mpmath.mpc, mpmath.mpf)):
-                    result = from_mpmath(result, d)
+                    result = from_mpmath(result, precision=prec)
         return result
 
 

mathics/core/convert/mpmath.py

@@ -13,38 +13,34 @@
 @lru_cache(maxsize=1024)
 def from_mpmath(
     value: Union[mpmath.mpf, mpmath.mpc],
-    prec: Optional[float] = None,
-    acc: Optional[float] = None,
+    precision: Optional[int] = None,
 ) -> Atom:
-    "Converts mpf or mpc to Number."
+    """
+    Converts mpf or mpc to Number.
+    The optional parameter `precision` represents
+    the binary precision.
+    """
     if mpmath.isnan(value):
         return SymbolIndeterminate
     if isinstance(value, mpmath.mpf):
         if mpmath.isinf(value):
             direction = Integer1 if value > 0 else IntegerM1
             return Expression(SymbolDirectedInfinity, direction)
-        # if accuracy is given, override
-        # prec:
-        if acc is not None:
-            prec = acc
-            if value != 0.0:
-                offset = mpmath.log(-value if value < 0.0 else value, 10)
-                prec += offset
-        if prec is None:
+        if precision is None:
             return MachineReal(float(value))
         # If the error if of the order of the number, the number
         # is compatible with 0.
-        if prec < 1.0:
+        if precision < 1:
             return MachineReal0
         # HACK: use str here to prevent loss of precision
-        return PrecisionReal(sympy.Float(str(value), prec))
+        return PrecisionReal(sympy.Float(str(value), precision=precision - 1))
     elif isinstance(value, mpmath.mpc):
         if mpmath.isinf(value):
             return SymbolComplexInfinity
         if value.imag == 0.0:
-            return from_mpmath(value.real, prec, acc)
-        real = from_mpmath(value.real, prec, acc)
-        imag = from_mpmath(value.imag, prec, acc)
+            return from_mpmath(value.real, precision=precision)
+        real = from_mpmath(value.real, precision=precision)
+        imag = from_mpmath(value.imag, precision=precision)
         return Complex(real, imag)
     else:
         raise TypeError(type(value))