Arithmetic refactor

CHANGES.rst

@@ -9,7 +9,7 @@ New Builtins
 
 
 * `Elements`
-
+* `RealAbs` and `RealSign`
 
 Compatibility
 -------------
@@ -24,7 +24,7 @@ Internals
 * ``eval_abs`` and ``eval_sign`` extracted from ``Abs`` and ``Sign`` and added to ``mathics.eval.arithmetic``.
 * Maximum number of digits allowed in a string set to 7000 and can be adjusted using environment variable
   ``MATHICS_MAX_STR_DIGITS`` on Python versions that don't adjust automatically (like pyston).
-
+* Real number comparisons implemented is based now in the internal implementation of `RealSign`.
 
 Bugs
 ----

SYMBOLS_MANIFEST.txt

@@ -824,8 +824,10 @@ System`Read
 System`ReadList
 System`ReadProtected
 System`Real
+System`RealAbs
 System`RealDigits
 System`RealNumberQ
+System`RealSign
 System`Reals
 System`Reap
 System`Record

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

@@ -77,7 +77,13 @@
     SymbolTable,
     SymbolUndefined,
 )
-from mathics.eval.arithmetic import eval_Abs, eval_mpmath_function, eval_Sign
+from mathics.eval.arithmetic import (
+    eval_Abs,
+    eval_mpmath_function,
+    eval_negate_number,
+    eval_RealSign,
+    eval_Sign,
+)
 from mathics.eval.nevaluator import eval_N
 from mathics.eval.numerify import numerify
 
@@ -823,7 +829,9 @@ def evaluate(self, evaluation: Evaluation):
 
 class Im(SympyFunction):
     """
-    <url>:WMA link:https://reference.wolfram.com/language/ref/Im.html</url>
+    <url>
+    :WMA link:
+    https://reference.wolfram.com/language/ref/Im.html</url>
 
     <dl>
       <dt>'Im[$z$]'
@@ -1207,6 +1215,49 @@ class Real_(Builtin):
     name = "Real"
 
 
+class RealAbs(Builtin):
+    """
+    <url>
+        :Abs (Real):
+        https://en.wikipedia.org/wiki/Absolute_value</url> (<url>
+        :WMA link:
+        https://reference.wolfram.com/language/ref/RealAbs.html
+        </url>)
+
+    <dl>
+      <dt>'RealAbs[$x$]'
+      <dd>returns the absolute value of a real number $x$.
+    </dl>
+    'RealAbs' is also known as modulus. It is evaluated if $x$ can be compared \
+    with $0$.
+
+    >> RealAbs[-3.]
+     = 3.
+    'RealAbs[$z$]' is left unevaluated for complex $z$:
+    >> RealAbs[2. + 3. I]
+     = RealAbs[2. + 3. I]
+    >> D[RealAbs[x ^ 2], x]
+     = 2 x ^ 3 / RealAbs[x ^ 2]
+    """
+
+    attributes = A_LISTABLE | A_NUMERIC_FUNCTION | A_PROTECTED
+    rules = {
+        "D[RealAbs[x_],x_]": "x/RealAbs[x]",
+        "Integrate[RealAbs[x_],x_]": "1/2 x RealAbs[x]",
+        "Integrate[RealAbs[u_],{u_,a_,b_}]": "1/2 b RealAbs[b]-1/2 a RealAbs[a]",
+    }
+    summary_text = "real absolute value"
+
+    def eval(self, x: BaseElement, evaluation: Evaluation):
+        """RealAbs[x_]"""
+        real_sign = eval_RealSign(x)
+        if real_sign is IntegerM1:
+            return eval_negate_number(x)
+        if real_sign is None:
+            return
+        return x
+
+
 class RealNumberQ(Test):
     """
     ## Not found in WMA
@@ -1237,9 +1288,54 @@ def test(self, expr) -> bool:
         return isinstance(expr, (Integer, Rational, Real))
 
 
+class RealSign(Builtin):
+    """
+    <url>
+        :Signum:
+        https://en.wikipedia.org/wiki/Sign_function</url> (<url>
+        :WMA link:
+        https://reference.wolfram.com/language/ref/RealSign.html
+        </url>)
+    <dl>
+      <dt>'RealSign[$x$]'
+      <dd>returns $-1$, $0$ or $1$ depending on whether $x$ is negative,
+    zero or positive.
+    </dl>
+    'RealSign' is also known as $sgn$ or $signum$ function.
+
+    >> RealSign[-3.]
+     = -1
+    'RealSign[$z$]' is left unevaluated for complex $z$:
+    >> RealSign[2. + 3. I]
+     = RealSign[2. + 3. I]
+
+    >> D[RealSign[x^2],x]
+     = 2 x Piecewise[{{0, x ^ 2 != 0}}, Indeterminate]
+    >> Integrate[RealSign[u],{u,0,x}]
+     = RealAbs[x]
+    """
+
+    attributes = A_LISTABLE | A_NUMERIC_FUNCTION | A_PROTECTED
+    rules = {
+        "D[RealSign[x_],x_]": "Piecewise[{{0, x!=0}}, Indeterminate]",
+        "Integrate[RealSign[x_],x_]": "RealAbs[x]",
+        "Integrate[RealSign[u_],{u_, a_, b_}]": "RealAbs[b]-RealSign[a]",
+    }
+    summary_text = "real sign"
+
+    def eval(self, x: Number, evaluation: Evaluation) -> Optional[Integer]:
+        """RealSign[x_]"""
+        return eval_RealSign(x)
+
+
 class Sign(SympyFunction):
     """
-    <url>:WMA link:https://reference.wolfram.com/language/ref/Sign.html</url>
+    <url>
+        :Sign:
+        https://en.wikipedia.org/wiki/Sign_function</url> (<url>
+        :WMA link:
+        https://reference.wolfram.com/language/ref/Sign.html
+        </url>)
 
     <dl>
       <dt>'Sign[$x$]'

mathics/builtin/testing_expressions/equality_inequality.py

@@ -69,8 +69,7 @@ def numerify_args(items, evaluation) -> list:
             n_items = []
             for item in items:
                 if not isinstance(item, Number):
-                    # TODO: use $MaxExtraPrecision insterad of hard-coded 50
-                    item = eval_N(item, evaluation, SymbolMaxPrecision)
+                    item = eval_N(item, evaluation, SymbolMaxExtraPrecision)
                 n_items.append(item)
             items = n_items
         else:

mathics/core/systemsymbols.py

@@ -200,7 +200,9 @@
 SymbolRational = Symbol("System`Rational")
 SymbolRe = Symbol("System`Re")
 SymbolReal = Symbol("System`Real")
+SymbolRealAbs = Symbol("System`RealAbs")
 SymbolRealDigits = Symbol("System`RealDigits")
+SymbolRealSign = Symbol("System`RealSign")
 SymbolRepeated = Symbol("System`Repeated")
 SymbolRepeatedNull = Symbol("System`RepeatedNull")
 SymbolReturn = Symbol("System`Return")
@@ -223,7 +225,7 @@
 SymbolSlot = Symbol("System`Slot")
 SymbolSparseArray = Symbol("System`SparseArray")
 SymbolSplit = Symbol("System`Split")
-SymbolSqrt = Symbol("System'Sqrt")
+SymbolSqrt = Symbol("System`Sqrt")
 SymbolSqrtBox = Symbol("System`SqrtBox")
 SymbolStandardDeviation = Symbol("System`StandardDeviation")
 SymbolStandardForm = Symbol("System`StandardForm")