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

@@ -7,7 +7,7 @@
 Basic arithmetic functions, including complex number arithmetic.
 """
 
-from mathics.eval.numerify import numerify
+import sympy
 
 # This tells documentation how to sort this module
 sort_order = "mathics.builtin.mathematical-functions"
@@ -16,7 +16,6 @@
 from functools import lru_cache
 
 import mpmath
-import sympy
 
 from mathics.builtin.base import (
     Builtin,
@@ -55,11 +54,9 @@
 from mathics.core.symbols import (
     Atom,
     Symbol,
-    SymbolAbs,
     SymbolFalse,
     SymbolList,
     SymbolPlus,
-    SymbolPower,
     SymbolTimes,
     SymbolTrue,
 )
@@ -72,8 +69,10 @@
     SymbolTable,
     SymbolUndefined,
 )
-from mathics.eval.arithmetic import eval_mpmath_function
-from mathics.eval.nevaluator import eval_N
+from mathics.eval.arithmetic import eval_abs, eval_mpmath_function, eval_sign
+from mathics.eval.numerify import numerify
+
+ExpressionComplexInfinity = Expression(SymbolDirectedInfinity)
 
 
 class _MPMathFunction(SympyFunction):
@@ -208,6 +207,13 @@ class Abs(_MPMathFunction):
     summary_text = "absolute value of a number"
     sympy_name = "Abs"
 
+    def eval(self, x, evaluation):
+        "%(name)s[x_]"
+        result = eval_abs(x)
+        if result is not None:
+            return result
+        return super(Abs, self).eval(x, evaluation)
+
 
 class Arg(_MPMathFunction):
     """
@@ -590,8 +596,7 @@ class DirectedInfinity(SympyFunction):
      = Indeterminate
 
     >> DirectedInfinity[0]
-     : Indeterminate expression 0 Infinity encountered.
-     = Indeterminate
+     = ComplexInfinity
 
     #> DirectedInfinity[1+I]+DirectedInfinity[2+I]
      = (2 / 5 + I / 5) Sqrt[5] Infinity + (1 / 2 + I / 2) Sqrt[2] Infinity
@@ -604,7 +609,7 @@ class DirectedInfinity(SympyFunction):
     rules = {
         "DirectedInfinity[Indeterminate]": "Indeterminate",
         "DirectedInfinity[args___] ^ -1": "0",
-        "0 * DirectedInfinity[args___]": "Message[Infinity::indet, Unevaluated[0 DirectedInfinity[args]]]; Indeterminate",
+        "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[]",
@@ -622,39 +627,44 @@ class DirectedInfinity(SympyFunction):
             "Indeterminate"
         ),
         "DirectedInfinity[args___] + _?NumberQ": "DirectedInfinity[args]",
-        "DirectedInfinity[0]": (
+        "DirectedInfinity[Alternatives[0, 0.]]": "DirectedInfinity[]",
+        "DirectedInfinity[0.]": (
             "Message[Infinity::indet,"
-            "  Unevaluated[DirectedInfinity[0]]];"
+            "  Unevaluated[DirectedInfinity[0.]]];"
             "Indeterminate"
         ),
-        "DirectedInfinity[0.]": (
+        "Alternatives[0, 0.] DirectedInfinity[z___]": (
             "Message[Infinity::indet,"
-            "  Unevaluated[DirectedInfinity[0.]]];"
+            "  Unevaluated[0 DirectedInfinity[z]]];"
             "Indeterminate"
         ),
-        "DirectedInfinity[DirectedInfinity[x___]]": "DirectedInfinity[x]",
+        "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(self, z, evaluation):
-        """DirectedInfinity[z_]"""
-        if z in (Integer1, IntegerM1):
+    def eval_directed_infinity(self, direction, evaluation):
+        """DirectedInfinity[direction_]"""
+        if direction in (Integer1, IntegerM1):
             return None
-        if isinstance(z, Number) or isinstance(eval_N(z, evaluation), Number):
-            direction = (z / Expression(SymbolAbs, z)).evaluate(evaluation)
-            return Expression(
-                SymbolDirectedInfinity,
-                direction,
-                elements_properties=ElementsProperties(True, True, True),
-            )
-        return None
+        if direction.is_zero:
+            return ExpressionComplexInfinity
+
+        normalized_direction = eval_sign(direction)
+        if normalized_direction is None:
+            return None
+        return Expression(
+            SymbolDirectedInfinity,
+            normalized_direction.evaluate(evaluation),
+            elements_properties=ElementsProperties(True, False, False),
+        )
 
     def to_sympy(self, expr, **kwargs):
         if len(expr.elements) == 1:
@@ -719,8 +729,8 @@ class Im(SympyFunction):
 
     def eval_complex(self, number, evaluation):
         "Im[number_Complex]"
-
-        return number.imag
+        if isinstance(number, Complex):
+            return number.imag
 
     def eval_number(self, number, evaluation):
         "Im[number_?NumberQ]"
@@ -990,8 +1000,8 @@ class Re(SympyFunction):
 
     def eval_complex(self, number, evaluation):
         "Re[number_Complex]"
-
-        return number.real
+        if isinstance(number, Complex):
+            return number.real
 
     def eval_number(self, number, evaluation):
         "Re[number_?NumberQ]"
@@ -1000,7 +1010,6 @@ def eval_number(self, number, evaluation):
 
     def eval(self, number, evaluation):
         "Re[number_]"
-
         return from_sympy(sympy.re(number.to_sympy().expand(complex=True)))
 
 
@@ -1150,18 +1159,16 @@ class Sign(SympyFunction):
 
     def eval(self, x, evaluation):
         "%(name)s[x_]"
-        # Sympy and mpmath do not give the desired form of complex number
-        if isinstance(x, Complex):
-            return Expression(
-                SymbolTimes,
-                x,
-                Expression(SymbolPower, Expression(SymbolAbs, x), IntegerM1),
-            )
+        result = eval_sign(x)
+        if result is not None:
+            return result
+        return None
 
         sympy_x = x.to_sympy()
         if sympy_x is None:
             return None
-        return super().eval(x, evaluation)
+        # Unhandled cases. Use sympy
+        return super(Sign, self).eval(x, evaluation)
 
     def eval_error(self, x, seqs, evaluation):
         "Sign[x_, seqs__]"

mathics/builtin/atomic/numbers.py

@@ -179,7 +179,7 @@ class Accuracy(Builtin):
 
     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];
+    >> z=Complex[3.00``2, 4.00``2];
     >> Accuracy[z] == -Log[10, Sqrt[10^(-2 Accuracy[Re[z]]) + 10^(-2 Accuracy[Im[z]])]]
      = True
 

mathics/builtin/drawing/plot.py

@@ -2431,6 +2431,7 @@ class ParametricPlot(_Plot):
      = -Graphics-
 
     >> ParametricPlot[{Cos[u] / u, Sin[u] / u}, {u, 0, 50}, PlotRange->0.5]
+     :  Indeterminate expression 0 ComplexInfinity encountered.
      = -Graphics-
 
     >> ParametricPlot[{{Sin[u], Cos[u]},{0.6 Sin[u], 0.6 Cos[u]}, {0.2 Sin[u], 0.2 Cos[u]}}, {u, 0, 2 Pi}, PlotRange->1, AspectRatio->1]

mathics/builtin/layout.py

@@ -288,6 +288,19 @@ def eval(self, expr, evaluation) -> Real:
         return Real(precedence)
 
 
+class PrecedenceForm(Builtin):
+    """
+    <url>:WMA link:https://reference.wolfram.com/language/ref/PrecedenceForm.html</url>
+
+    <dl>
+      <dt>'PrecedenceForm'[$expr$, $prec$]
+      <dd> format $expr$ parenthesized as it would be if it contained an operator of precedence $prec$.
+    </dl>
+    """
+
+    summary_text = "parenthesize with a precedence"
+
+
 class Prefix(BinaryOperator):
     """
     <url>:WMA link:https://reference.wolfram.com/language/ref/Prefix.html</url>

mathics/builtin/makeboxes.py

@@ -343,7 +343,7 @@ class MakeBoxes(Builtin):
         "MathMLForm)[expr_], StandardForm|TraditionalForm]": ("MakeBoxes[expr, form]"),
         "MakeBoxes[(form:StandardForm|OutputForm|MathMLForm|TeXForm)[expr_], OutputForm]": "MakeBoxes[expr, form]",
         "MakeBoxes[(form:FullForm|InputForm)[expr_], StandardForm|TraditionalForm|OutputForm]": "StyleBox[MakeBoxes[expr, form], ShowStringCharacters->True]",
-        "MakeBoxes[PrecedenceForm[expr_, prec_], f_]": "MakeBoxes[expr, f]",
+        "MakeBoxes[PrecedenceForm[expr_, prec_], f_]": 'Print["PrecedenceForm", f];MakeBoxes[expr, f]',
         "MakeBoxes[Style[expr_, OptionsPattern[Style]], f_]": (
             "StyleBox[MakeBoxes[expr, f], "
             "ImageSizeMultipliers -> OptionValue[ImageSizeMultipliers]]"
@@ -400,12 +400,12 @@ def eval_general(self, expr, f, evaluation):
             result.append(to_boxes(String(right), evaluation))
             return RowBox(*result)
 
-    def eval_outerprecedenceform(self, expr, precedence, evaluation):
+    def eval_outerprecedenceform(self, expr, precedence, form, evaluation):
         """MakeBoxes[PrecedenceForm[expr_, precedence_],
-        StandardForm|TraditionalForm|OutputForm|InputForm]"""
+        form:StandardForm|TraditionalForm|OutputForm|InputForm]"""
 
         py_precedence = precedence.get_int_value()
-        boxes = MakeBoxes(expr)
+        boxes = MakeBoxes(expr, form)
         return parenthesize(py_precedence, expr, boxes, True)
 
     def eval_postprefix(self, p, expr, h, precedence, form, evaluation):