diff --git a/mathics/builtin/arithfns/basic.py b/mathics/builtin/arithfns/basic.py
index 6fb5a8873..e13c7341a 100644
--- a/mathics/builtin/arithfns/basic.py
+++ b/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()
diff --git a/mathics/eval/arithmetic.py b/mathics/eval/arithmetic.py
index abd4e51a8..1a97f7991 100644
--- a/mathics/eval/arithmetic.py
+++ b/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``.
"""
diff --git a/test/builtin/arithmetic/test_basic.py b/test/builtin/arithmetic/test_basic.py
index d99b0b9dc..9a05aac28 100644
--- a/test/builtin/arithmetic/test_basic.py
+++ b/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),
diff --git a/test/format/test_format.py b/test/format/test_format.py
index 161ebc5df..ee81add3c 100644
--- a/test/format/test_format.py
+++ b/test/format/test_format.py
@@ -456,34 +456,53 @@
"Sqrt[1/(1+1/(1+1/a))]": {
"msg": "SqrtBox",
"text": {
- "System`StandardForm": "Sqrt[1 / (1+1 / (1+1 / a))]",
- "System`TraditionalForm": "Sqrt[1 / (1+1 / (1+1 / a))]",
- "System`InputForm": "Sqrt[1 / (1 + 1 / (1 + 1 / a))]",
- "System`OutputForm": "Sqrt[1 / (1 + 1 / (1 + 1 / a))]",
+ "System`StandardForm": "1 / Sqrt[1+1 / (1+1 / a)]",
+ "System`TraditionalForm": "1 / Sqrt[1+1 / (1+1 / a)]",
+ "System`InputForm": "1 / Sqrt[1 + 1 / (1 + 1 / a)]",
+ "System`OutputForm": "1 / Sqrt[1 + 1 / (1 + 1 / a)]",
},
"mathml": {
"System`StandardForm": (
- " 1 1 + 1 1 + 1 a ",
+ (
+ r"1 1 + 1 "
+ r"1 + 1 a "
+ r""
+ ),
"Fragile!",
),
"System`TraditionalForm": (
- " 1 1 + 1 1 + 1 a ",
+ (
+ r"1 1 + 1 "
+ r"1 + 1 a "
+ r""
+ ),
"Fragile!",
),
"System`InputForm": (
- "Sqrt [ 1 / ( 1 + 1 / ( 1 + 1 / a ) ) ]",
+ (
+ r"1 / Sqrt [ "
+ r"1 + 1 / "
+ r"( 1 + 1 "
+ r" / a ) ]"
+ ),
"Fragile!",
),
"System`OutputForm": (
- "Sqrt [ 1 / ( 1 + 1 / ( 1 + 1 / a ) ) ]",
+ (
+ r"1 / Sqrt ["
+ r" 1 + 1 "
+ r" / ( 1 "
+ r" + 1 / "
+ r"a ) ]"
+ ),
"Fragile!",
),
},
"latex": {
- "System`StandardForm": "\\sqrt{\\frac{1}{1+\\frac{1}{1+\\frac{1}{a}}}}",
- "System`TraditionalForm": "\\sqrt{\\frac{1}{1+\\frac{1}{1+\\frac{1}{a}}}}",
- "System`InputForm": "\\text{Sqrt}\\left[1\\text{ / }\\left(1\\text{ + }1\\text{ / }\\left(1\\text{ + }1\\text{ / }a\\right)\\right)\\right]",
- "System`OutputForm": "\\text{Sqrt}\\left[1\\text{ / }\\left(1\\text{ + }1\\text{ / }\\left(1\\text{ + }1\\text{ / }a\\right)\\right)\\right]",
+ "System`StandardForm": "\\frac{1}{\\sqrt{1+\\frac{1}{1+\\frac{1}{a}}}}",
+ "System`TraditionalForm": "\\frac{1}{\\sqrt{1+\\frac{1}{1+\\frac{1}{a}}}}",
+ "System`InputForm": r"1\text{ / }\text{Sqrt}\left[1\text{ + }1\text{ / }\left(1\text{ + }1\text{ / }a\right)\right]",
+ "System`OutputForm": r"1\text{ / }\text{Sqrt}\left[1\text{ + }1\text{ / }\left(1\text{ + }1\text{ / }a\right)\right]",
},
},
# Grids, arrays and matrices