refactor arithmetic power

Author: mmatera

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):
@@ -788,7 +828,6 @@ def inverse(item):
                 and isinstance(item.elements[1], (Integer, Rational, Real))
                 and item.elements[1].to_sympy() < 0
             ):  # nopep8
-
                 negative.append(inverse(item))
             elif isinstance(item, Rational):
                 numerator = item.numerator()

mathics/eval/arithmetic.py

@@ -1,7 +1,9 @@
 # -*- coding: utf-8 -*-
 
 """
-arithmetic-related evaluation functions.
+helper functions for arithmetic evaluation, which do not
+depends on the evaluation context. Conversions to Sympy are
+used just as a last resource.
 
 Many of these do do depend on the evaluation context. Conversions to Sympy are
 used just as a last resource.
@@ -319,48 +321,27 @@ def eval_complex_sign(n: BaseElement) -> Optional[BaseElement]:
     sign = eval_RealSign(expr)
     return sign or eval_complex_sign(expr)
 
-    if expr.has_form("Power", 2):
-        base, exp = expr.elements
-        if exp.is_zero:
-            return Integer1
-        if isinstance(exp, (Integer, Real, Rational)):
-            sign = eval_Sign(base) or Expression(SymbolSign, base)
-            return Expression(SymbolPower, sign, exp)
-        if isinstance(exp, Complex):
-            sign = eval_Sign(base) or Expression(SymbolSign, base)
-            return Expression(SymbolPower, sign, exp.real)
-        if test_arithmetic_expr(exp):
-            sign = eval_Sign(base) or Expression(SymbolSign, base)
-            return Expression(SymbolPower, sign, exp)
-        return None
-    if expr.get_head() is SymbolTimes:
-        abs_value = eval_Abs(eval_multiply_numbers(*expr.elements))
-        if abs_value is Integer1:
-            return expr
-        if abs_value is None:
-            return None
-        criteria = eval_add_numbers(abs_value, IntegerM1)
-        if test_zero_arithmetic_expr(criteria, numeric=True):
-            return expr
-        return None
-    if expr.get_head() is SymbolPlus:
-        abs_value = eval_Abs(eval_add_numbers(*expr.elements))
-        if abs_value is Integer1:
-            return expr
-        if abs_value is None:
-            return None
-        criteria = eval_add_numbers(abs_value, IntegerM1)
-        if test_zero_arithmetic_expr(criteria, numeric=True):
-            return expr
-        return None
 
-    if test_arithmetic_expr(expr):
-        if test_zero_arithmetic_expr(expr):
-            return Integer0
-        if test_positive_arithmetic_expr(expr):
-            return Integer1
-        if test_negative_arithmetic_expr(expr):
-            return IntegerM1
+def eval_Sign_number(n: Number) -> Number:
+    """
+    Evals the absolute value of a number.
+    """
+    if n.is_zero:
+        return Integer0
+    if isinstance(n, (Integer, Rational, Real)):
+        return Integer1 if n.value > 0 else IntegerM1
+    if isinstance(n, Complex):
+        abs_sq = eval_add_numbers(
+            *(eval_multiply_numbers(x, x) for x in (n.real, n.imag))
+        )
+        criteria = eval_add_numbers(abs_sq, IntegerM1)
+        if test_zero_arithmetic_expr(criteria):
+            return n
+        if n.is_inexact():
+            return eval_multiply_numbers(n, eval_Power_number(abs_sq, RealM0p5))
+        if test_zero_arithmetic_expr(criteria, numeric=True):
+            return n
+        return eval_multiply_numbers(n, eval_Power_number(abs_sq, RationalMOneHalf))
 
 
 def eval_mpmath_function(
@@ -390,6 +371,31 @@ def eval_mpmath_function(
             return call_mpmath(mpmath_function, tuple(mpmath_args), prec)
 
 
+def eval_Exponential(exp: BaseElement) -> BaseElement:
+    """
+    Eval E^exp
+    """
+    # If both base and exponent are exact quantities,
+    # use sympy.
+
+    if not exp.is_inexact():
+        exp_sp = exp.to_sympy()
+        if exp_sp is None:
+            return None
+        return from_sympy(sympy.Exp(exp_sp))
+
+    prec = exp.get_precision()
+    if prec is not None:
+        if exp.is_machine_precision():
+            number = mpmath.exp(exp.to_mpmath())
+            result = from_mpmath(number)
+            return result
+        else:
+            with mpmath.workprec(prec):
+                number = mpmath.exp(exp.to_mpmath())
+                return from_mpmath(number, prec)
+
+
 def eval_Plus(*items: BaseElement) -> BaseElement:
     "evaluate Plus for general elements"
     numbers, items_tuple = segregate_numbers_from_sorted_list(*items)
@@ -457,6 +463,13 @@ def append_last():
         elements_properties=ElementsProperties(False, False, True),
     )
 
+    elements.sort()
+    return Expression(
+        SymbolPlus,
+        *elements,
+        elements_properties=ElementsProperties(False, False, True),
+    )
+
 
 def eval_Power_number(base: Number, exp: Number) -> Optional[Number]:
     """
@@ -688,8 +701,88 @@ def eval_Times(*items: BaseElement) -> BaseElement:
     )
 
 
+# Here I used the convention of calling eval_* to functions that can produce a new expression, or None
+# if the result can not be evaluated, or is trivial. For example, if we call eval_Power_number(Integer2, RationalOneHalf)
+# it returns ``None`` instead of ``Expression(SymbolPower, Integer2, RationalOneHalf)``.
+# The reason is that these functions are written to be part of replacement rules, to be applied during the evaluation process.
+# In that process, a rule is considered applied if produces an expression that is different from the original one, or
+# if the replacement function returns (Python's) ``None``.
+#
+# For example, when the expression ``Power[4, 1/2]`` is evaluated, a (Builtin) rule  ``Power[base_, exp_]->eval_repl_rule(base, expr)``
+# is applied. If the rule matches, `repl_rule` is called with arguments ``(4, 1/2)`` and produces `2`. As `Integer2.sameQ(Power[4, 1/2])`
+# is False, then no new rules for `Power` are checked, and a new round of evaluation is atempted.
+#
+# On the other hand, if ``Power[3, 1/2]``, ``repl_rule`` can do two possible things: one is return ``Power[3, 1/2]``. If it does,
+# the rule is considered applied. Then, the evaluation method checks if `Power[3, 1/2].sameQ(Power[3, 1/2])`. In this case it is true,
+# and then the expression is kept as it is.
+# The other possibility is to return  (Python's) `None`. In that case, the evaluator considers that the rule failed to be applied,
+# and look for another rule associated to ``Power``. To return ``None`` produces then a faster evaluation, since no ``sameQ`` call is needed,
+# and do not prevent that other rules are attempted.
+#
+# The bad part of using ``None`` as a return is that I would expect that ``eval`` produces always a valid Expression, so if at some point of
+# the code I call ``eval_Power_number(Integer3, RationalOneHalf)`` I get ``Expression(SymbolPower, Integer3, RationalOneHalf)``.
+#
+# From my point of view, it would make more sense to use the following convention:
+#  * if the method has signature ``eval_method(...)->BaseElement:`` then use the prefix ``eval_``
+#  * if the method has the siguature ``apply_method(...)->Optional[BaseElement]`` use the prefix ``apply_`` or maybe ``repl_``.
+#
+# In any case, let's keep the current convention.
+#
+#
+
+
+def associate_powers(expr: BaseElement, power: BaseElement = Integer1) -> BaseElement:
+    """
+    base^a^b^c^...^power -> base^(a*b*c*...power)
+    provided one of the following cases
+    * `a`, `b`, ... `power` are all integer numbers
+    * `a`, `b`,... are Rational/Real number with absolute value <=1,
+      and the other powers are not integer numbers.
+    * `a` is not a Rational/Real number, and b, c, ... power are all
+      integer numbers.
+    """
+    powers = []
+    base = expr
+    if power is not Integer1:
+        powers.append(power)
+
+    while base.has_form("Power", 2):
+        previous_base, outer_power = base, power
+        base, power = base.elements
+        if len(powers) == 0:
+            if power is not Integer1:
+                powers.append(power)
+            continue
+        if power is IntegerM1:
+            powers.append(power)
+            continue
+        if isinstance(power, (Rational, Real)):
+            if abs(power.value) < 1:
+                powers.append(power)
+                continue
+        # power is not rational/real and outer_power is integer,
+        elif isinstance(outer_power, Integer):
+            if power is not Integer1:
+                powers.append(power)
+            if isinstance(power, Integer):
+                continue
+            else:
+                break
+        # in any other case, use the previous base and
+        # exit the loop
+        base = previous_base
+        break
+
+    if len(powers) == 0:
+        return base
+    elif len(powers) == 1:
+        return Expression(SymbolPower, base, powers[0])
+    result = Expression(SymbolPower, base, Expression(SymbolTimes, *powers))
+    return result
+
+
 def eval_add_numbers(
-    *numbers: Number,
+    *numbers: List[Number],
 ) -> BaseElement:
     """
     Add the elements in ``numbers``.
@@ -736,7 +829,7 @@ def eval_inverse_number(n: Number) -> Number:
     return eval_Power_number(n, IntegerM1)
 
 
-def eval_multiply_numbers(*numbers: Number) -> Number:
+def eval_multiply_numbers(*numbers: Number) -> BaseElement:
     """
     Multiply the elements in ``numbers``.
     """

test/builtin/arithmetic/test_basic.py

@@ -153,8 +153,8 @@ def test_multiply(str_expr, str_expected, msg):
         ("a  b  DirectedInfinity[q]", "a b (q Infinity)", ""),
         # Failing tests
         # Problem with formatting. Parenthezise are missing...
-        #        ("a  b  DirectedInfinity[-I]", "a b (-I Infinity)",  ""),
-        #        ("a  b  DirectedInfinity[-3]", "a b (-Infinity)",  ""),
+        ("a  b  DirectedInfinity[-I]", "a b (-I Infinity)", ""),
+        ("a  b  DirectedInfinity[-3]", "a b (-Infinity)", ""),
     ],
 )
 def test_directed_infinity_precedence(str_expr, str_expected, msg):
@@ -197,7 +197,7 @@ def test_directed_infinity_precedence(str_expr, str_expected, msg):
         ("I^(2/3)", "(-1) ^ (1 / 3)", None),
         # In WMA, the next test would return ``-(-I)^(2/3)``
         # which is less compact and elegant...
-        #        ("(-I)^(2/3)", "(-1) ^ (-1 / 3)", None),
+        ("(-I)^(2/3)", "(-1) ^ (-1 / 3)", None),
         ("(2+3I)^3", "-46 + 9 I", None),
         ("(1.+3. I)^.6", "1.46069 + 1.35921 I", None),
         ("3^(1+2 I)", "3 ^ (1 + 2 I)", None),
@@ -208,15 +208,15 @@ def test_directed_infinity_precedence(str_expr, str_expected, msg):
         # sympy, which produces the result
         ("(3/Pi)^(-I)", "(3 / Pi) ^ (-I)", None),
         # Association rules
-        #        ('(a^"w")^2', 'a^(2 "w")', "Integer power of a power with string exponent"),
+        ('(a^"w")^2', 'a^(2 "w")', "Integer power of a power with string exponent"),
         ('(a^2)^"w"', '(a ^ 2) ^ "w"', None),
         ('(a^2)^"w"', '(a ^ 2) ^ "w"', None),
         ("(a^2)^(1/2)", "Sqrt[a ^ 2]", None),
         ("(a^(1/2))^2", "a", None),
         ("(a^(1/2))^2", "a", None),
         ("(a^(3/2))^3.", "(a ^ (3 / 2)) ^ 3.", None),
-        #        ("(a^(1/2))^3.", "a ^ 1.5", "Power associativity rational, real"),
-        #        ("(a^(.3))^3.", "a ^ 0.9", "Power associativity for real powers"),
+        ("(a^(1/2))^3.", "a ^ 1.5", "Power associativity rational, real"),
+        ("(a^(.3))^3.", "a ^ 0.9", "Power associativity for real powers"),
         ("(a^(1.3))^3.", "(a ^ 1.3) ^ 3.", None),
         # Exponentials involving expressions
         ("(a^(p-2 q))^3", "a ^ (3 p - 6 q)", None),