arithmetic refactor

mathics/builtin/arithfns/basic.py

@@ -6,6 +6,8 @@
 
 """
 
+import sympy
+
 from mathics.builtin.arithmetic import _MPMathFunction, create_infix
 from mathics.builtin.base import BinaryOperator, Builtin, PrefixOperator, SympyFunction
 from mathics.core.atoms import (
@@ -38,7 +40,6 @@
     Symbol,
     SymbolDivide,
     SymbolHoldForm,
-    SymbolNull,
     SymbolPower,
     SymbolTimes,
 )
@@ -49,10 +50,17 @@
     SymbolInfix,
     SymbolLeft,
     SymbolMinus,
+    SymbolOverflow,
     SymbolPattern,
-    SymbolSequence,
 )
-from mathics.eval.arithmetic import eval_Plus, eval_Times
+from mathics.eval.arithmetic import (
+    associate_powers,
+    eval_Exponential,
+    eval_Plus,
+    eval_Power_inexact,
+    eval_Power_number,
+    eval_Times,
+)
 from mathics.eval.nevaluator import eval_N
 from mathics.eval.numerify import numerify
 
@@ -520,6 +528,8 @@ class Power(BinaryOperator, _MPMathFunction):
     rules = {
         "Power[]": "1",
         "Power[x_]": "x",
+        "Power[I,-1]": "-I",
+        "Power[-1, 1/2]": "I",
     }
 
     summary_text = "exponentiate"
@@ -528,15 +538,15 @@ class Power(BinaryOperator, _MPMathFunction):
     # Remember to up sympy doc link when this is corrected
     sympy_name = "Pow"
 
+    def eval_exp(self, x, evaluation):
+        "Power[E, x]"
+        return eval_Exponential(x)
+
     def eval_check(self, x, y, evaluation):
         "Power[x_, y_]"
-
-        # Power uses _MPMathFunction but does some error checking first
-        if isinstance(x, Number) and x.is_zero:
-            if isinstance(y, Number):
-                y_err = y
-            else:
-                y_err = eval_N(y, evaluation)
+        # if x is zero
+        if x.is_zero:
+            y_err = y if isinstance(y, Number) else eval_N(y, evaluation)
             if isinstance(y_err, Number):
                 py_y = y_err.round_to_float(permit_complex=True).real
                 if py_y > 0:
@@ -550,17 +560,47 @@ def eval_check(self, x, y, evaluation):
                     evaluation.message(
                         "Power", "infy", Expression(SymbolPower, x, y_err)
                     )
-                    return SymbolComplexInfinity
-        if isinstance(x, Complex) and x.real.is_zero:
-            yhalf = Expression(SymbolTimes, y, RationalOneHalf)
-            factor = self.eval(Expression(SymbolSequence, x.imag, y), evaluation)
-            return Expression(
-                SymbolTimes, factor, Expression(SymbolPower, IntegerM1, yhalf)
-            )
-
-        result = self.eval(Expression(SymbolSequence, x, y), evaluation)
-        if result is None or result != SymbolNull:
-            return result
+                return SymbolComplexInfinity
+
+        # If x and y are inexact numbers, use the numerical function
+
+        if x.is_inexact() and y.is_inexact():
+            try:
+                return eval_Power_inexact(x, y)
+            except OverflowError:
+                evaluation.message("General", "ovfl")
+                return Expression(SymbolOverflow)
+
+        # Tries to associate powers a^b^c-> a^(b*c)
+        assoc = associate_powers(x, y)
+        if not assoc.has_form("Power", 2):
+            return assoc
+
+        assoc = numerify(assoc, evaluation)
+        x, y = assoc.elements
+        # If x and y are numbers
+        if isinstance(x, Number) and isinstance(y, Number):
+            try:
+                return eval_Power_number(x, y)
+            except OverflowError:
+                evaluation.message("General", "ovfl")
+                return Expression(SymbolOverflow)
+
+        # if x or y are inexact, leave the expression
+        # as it is:
+        if x.is_inexact() or y.is_inexact():
+            return assoc
+
+        # Finally, try to convert to sympy
+        base_sp, exp_sp = x.to_sympy(), y.to_sympy()
+        if base_sp is None or exp_sp is None:
+            # If base or exp can not be converted to sympy,
+            # returns the result of applying the associative
+            # rule.
+            return assoc
+
+        result = from_sympy(sympy.Pow(base_sp, exp_sp))
+        return result.evaluate_elements(evaluation)
 
 
 class Sqrt(SympyFunction):

mathics/builtin/arithmetic.py

@@ -629,10 +629,16 @@ class DirectedInfinity(SympyFunction):
     rules = {
         "DirectedInfinity[args___] ^ -1": "0",
         # Special arguments:
-        "DirectedInfinity[DirectedInfinity[args___]]": "DirectedInfinity[args]",
         "DirectedInfinity[Indeterminate]": "Indeterminate",
         "DirectedInfinity[Alternatives[0, 0.]]": "DirectedInfinity[]",
         # Plus
+        "DirectedInfinity[DirectedInfinity[args___]]": "DirectedInfinity[args]",
+        # "DirectedInfinity[a_?NumericQ] /; N[Abs[a]] != 1": "DirectedInfinity[a / Abs[a]]",
+        # "DirectedInfinity[a_] * DirectedInfinity[b_]": "DirectedInfinity[a*b]",
+        # "DirectedInfinity[] * DirectedInfinity[args___]": "DirectedInfinity[]",
+        # Rules already implemented in Times.eval
+        #        "z_?NumberQ * DirectedInfinity[]": "DirectedInfinity[]",
+        #        "z_?NumberQ * DirectedInfinity[a_]": "DirectedInfinity[z * a]",
         "DirectedInfinity[a_] + DirectedInfinity[b_] /; b == -a": (
             "Message[Infinity::indet,"
             "  Unevaluated[DirectedInfinity[a] + DirectedInfinity[b]]];"
@@ -644,23 +650,27 @@ class DirectedInfinity(SympyFunction):
             "Indeterminate"
         ),
         "DirectedInfinity[args___] + _?NumberQ": "DirectedInfinity[args]",
-        # Times. See if can be reinstalled in eval_Times
+        "DirectedInfinity[Alternatives[0, 0.]]": "DirectedInfinity[]",
+        "DirectedInfinity[0.]": (
+            "Message[Infinity::indet,"
+            "  Unevaluated[DirectedInfinity[0.]]];"
+            "Indeterminate"
+        ),
         "Alternatives[0, 0.] DirectedInfinity[z___]": (
             "Message[Infinity::indet,"
             "  Unevaluated[0 DirectedInfinity[z]]];"
             "Indeterminate"
         ),
-        "a_?NumericQ * DirectedInfinity[b_]": "DirectedInfinity[a * b]",
         "a_ DirectedInfinity[]": "DirectedInfinity[]",
         "DirectedInfinity[a_] * DirectedInfinity[b_]": "DirectedInfinity[a * b]",
+        "a_?NumericQ * DirectedInfinity[b_]": "DirectedInfinity[a * b]",
     }
 
     formats = {
         "DirectedInfinity[1]": "HoldForm[Infinity]",
-        "DirectedInfinity[-1]": "HoldForm[-Infinity]",
+        "DirectedInfinity[-1]": "PrecedenceForm[-HoldForm[Infinity], 390]",
         "DirectedInfinity[]": "HoldForm[ComplexInfinity]",
-        "DirectedInfinity[DirectedInfinity[z_]]": "DirectedInfinity[z]",
-        "DirectedInfinity[z_?NumericQ]": "HoldForm[z Infinity]",
+        "DirectedInfinity[z_]": "PrecedenceForm[z HoldForm[Infinity], 390]",
     }
 
     def eval_complex_infinity(self, evaluation: Evaluation):

mathics/builtin/base.py

@@ -799,10 +799,14 @@ def get_operator_display(self) -> Optional[str]:
 
 
 class Predefined(Builtin):
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.symbol = Symbol(self.get_name())
+
     def get_functions(self, prefix="eval", is_pymodule=False) -> List[Callable]:
         functions = list(super().get_functions(prefix))
         if prefix == "eval" and hasattr(self, "evaluate"):
-            functions.append((Symbol(self.get_name()), self.evaluate))
+            functions.append((self.symbol, self.evaluate))
         return functions
 
 

mathics/builtin/numbers/constants.py

@@ -9,25 +9,19 @@
 # This tells documentation how to sort this module
 sort_order = "mathics.builtin.mathematical-constants"
 
-
 import math
+from typing import Optional
 
 import mpmath
 import numpy
 import sympy
 
 from mathics.builtin.base import Builtin, Predefined, SympyObject
-from mathics.core.atoms import MachineReal, PrecisionReal
+from mathics.core.atoms import NUMERICAL_CONSTANTS, MachineReal, PrecisionReal
 from mathics.core.attributes import A_CONSTANT, A_PROTECTED, A_READ_PROTECTED
+from mathics.core.element import BaseElement
 from mathics.core.evaluation import Evaluation
-from mathics.core.number import (
-    MACHINE_DIGITS,
-    MAX_MACHINE_NUMBER,
-    MIN_MACHINE_NUMBER,
-    PrecisionValueError,
-    get_precision,
-    prec,
-)
+from mathics.core.number import MACHINE_DIGITS, PrecisionValueError, get_precision, prec
 from mathics.core.symbols import Atom, Symbol, strip_context
 from mathics.core.systemsymbols import SymbolIndeterminate
 
@@ -89,28 +83,33 @@ def eval_N(self, precision, evaluation):
     def is_constant(self) -> bool:
         return True
 
-    def get_constant(self, precision, evaluation):
+    def get_constant(
+        self,
+        precision: Optional[BaseElement] = None,
+        evaluation: Optional[Evaluation] = None,
+    ):
         # first, determine the precision
         d = None
-        if precision:
-            try:
-                d = get_precision(precision, evaluation)
-            except PrecisionValueError:
-                pass
+        preference = None
+        if evaluation:
+            if precision:
+                try:
+                    d = get_precision(precision, evaluation)
+                except PrecisionValueError:
+                    pass
+
+            preflist = evaluation._preferred_n_method.copy()
+            while preflist:
+                pref_method = preflist.pop()
+                if pref_method in ("numpy", "mpmath", "sympy"):
+                    preference = pref_method
+                    break
 
         if d is None:
             d = MACHINE_DIGITS
 
         # If preference not especified, determine it
         # from the precision.
-        preference = None
-        preflist = evaluation._preferred_n_method.copy()
-        while preflist:
-            pref_method = preflist.pop()
-            if pref_method in ("numpy", "mpmath", "sympy"):
-                preference = pref_method
-                break
-
         if preference is None:
             if d <= MACHINE_DIGITS:
                 preference = "numpy"
@@ -131,10 +130,16 @@ def get_constant(self, precision, evaluation):
                 preference = "mpmath"
             else:
                 preference = ""
+
         if preference == "numpy":
-            value = numpy_constant(self.numpy_name)
             if d == MACHINE_DIGITS:
-                return MachineReal(value)
+                try:
+                    return NUMERICAL_CONSTANTS[self.symbol]
+                except KeyError:
+                    value = MachineReal(numpy_constant(self.numpy_name))
+                    NUMERICAL_CONSTANTS[self.symbol] = value
+                    return value
+            value = numpy_constant(self.numpy_name)
         if preference == "sympy":
             value = sympy_constant(self.sympy_name, d + 2)
         if preference == "mpmath":
@@ -177,13 +182,16 @@ class _NumpyConstant(_Constant_Common):
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
         if self.numpy_name is None:
-            self.numpy_name = strip_context(self.get_name()).lower()
+            self.numpy_name = strip_context(self.symbol.name).lower()
         self.mathics_to_numpy[self.__class__.__name__] = self.numpy_name
+        try:
+            value_float = numpy_constant(self.numpy_name)
+        except AttributeError:
+            value_float = self.to_numpy(self.symbol)
+        NUMERICAL_CONSTANTS[self.symbol] = MachineReal(value_float)
 
     def to_numpy(self, args):
-        if self.numpy_name is None or len(args) != 0:
-            return None
-        return self.get_constant()
+        return NUMERICAL_CONSTANTS[self.symbol]
 
 
 class _SympyConstant(_Constant_Common, SympyObject):
@@ -608,7 +616,7 @@ class MaxMachineNumber(Predefined):
     summary_text = "largest normalized positive machine number"
 
     def evaluate(self, evaluation: Evaluation) -> MachineReal:
-        return MachineReal(MAX_MACHINE_NUMBER)
+        return NUMERICAL_CONSTANTS[self.symbol]
 
 
 class MinMachineNumber(Predefined):
@@ -635,10 +643,10 @@ class MinMachineNumber(Predefined):
     summary_text = "smallest normalized positive machine number"
 
     def evaluate(self, evaluation: Evaluation) -> MachineReal:
-        return MachineReal(MIN_MACHINE_NUMBER)
+        return NUMERICAL_CONSTANTS[self.symbol]
 
 
-class Pi(_MPMathConstant, _SympyConstant):
+class Pi(_MPMathConstant, _NumpyConstant, _SympyConstant):
     """
     <url>
     :Pi, \u03c0: https://en.wikipedia.org/wiki/Pi</url> (<url>
@@ -740,3 +748,10 @@ class Underflow(Builtin):
         "Underflow[] * x_Real": "0.",
     }
     summary_text = "underflow in numeric evaluation"
+
+
+# Constants that are not numpy constants,
+for cls in (Catalan, Degree, Glaisher, GoldenRatio, Khinchin):
+    instance = cls(expression=False)
+    val = instance.get_constant()
+    NUMERICAL_CONSTANTS[instance.symbol] = MachineReal(val.value)

mathics/builtin/testing_expressions/numerical_properties.py

@@ -13,6 +13,7 @@
 from mathics.core.expression import Expression
 from mathics.core.symbols import BooleanType, SymbolFalse, SymbolTrue
 from mathics.core.systemsymbols import SymbolExpandAll, SymbolSimplify
+from mathics.eval.arithmetic import test_zero_arithmetic_expr
 from mathics.eval.nevaluator import eval_N
 
 
@@ -460,6 +461,10 @@ def eval(self, expr, evaluation):
         "%(name)s[expr_]"
         from sympy.matrices.utilities import _iszero
 
+        # This handles most of the arithmetic cases
+        if test_zero_arithmetic_expr(expr):
+            return SymbolTrue
+
         sympy_expr = expr.to_sympy()
         result = _iszero(sympy_expr)
         if result is None:

mathics/core/atoms.py

@@ -13,6 +13,8 @@
 from mathics.core.number import (
     FP_MANTISA_BINARY_DIGITS,
     MACHINE_PRECISION_VALUE,
+    MAX_MACHINE_NUMBER,
+    MIN_MACHINE_NUMBER,
     dps,
     min_prec,
     prec,
@@ -946,6 +948,14 @@ def is_zero(self) -> bool:
 MATHICS3_COMPLEX_I = Complex(Integer0, Integer1)
 MATHICS3_COMPLEX_I_NEG = Complex(Integer0, IntegerM1)
 
+# Numerical constants
+# These constants are populated by the `Predefined`
+# classes. See `mathics.builtin.numbers.constants`
+NUMERICAL_CONSTANTS = {
+    Symbol("System`$MaxMachineNumber"): MachineReal(MAX_MACHINE_NUMBER),
+    Symbol("System`$MinMachineNumber"): MachineReal(MIN_MACHINE_NUMBER),
+}
+
 
 class String(Atom, BoxElementMixin):
     value: str

mathics/core/number.py

@@ -68,7 +68,9 @@ def _get_float_inf(value, evaluation) -> Optional[float]:
     return value.round_to_float(evaluation)
 
 
-def get_precision(value, evaluation, show_messages=True) -> Optional[float]:
+def get_precision(
+    value: BaseElement, evaluation, show_messages: bool = True
+) -> Optional[float]:
     """
     Returns the ``float`` in the interval     [``$MinPrecision``, ``$MaxPrecision``] closest to ``value``.
     If ``value`` does not belongs to that interval, and ``show_messages`` is True, a Message warning is shown.