diff --git a/CHANGES.rst b/CHANGES.rst
index 774387872..249d1fa6b 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -76,7 +76,9 @@ Internals
#. ``_SetOperator.assign_elementary`` was renamed as ``_SetOperator.assign``. All the special cases are not handled by the ``_SetOperator.special_cases`` dict.
#. ``isort`` run over all Python files. More type annotations and docstrings on functions added.
#. caching on immutable atoms like, ``String``, ``Integer``, ``Real``, etc. was improved; the ``__hash__()`` function was sped up. There is a small speedup overall from this at the expense of increased memory.
-#. more type annotations added to functions, especially builtin functions
+#. more type annotations added to functions, especially builtin functions.
+#. Improved implementation for ``Accuracy`` and ``Precision``.
+#. Improved implementation for ``Plus``, ``Times`` and ``Power``, ``Exp``, ``Abs`` and ``Sign`` for numerical elements. Also, the new versions improves the compatibility with WMA in the results of arithmetic symbolic evaluation.
Bugs
@@ -93,6 +95,8 @@ Bugs
#. Units and Quantities were sometimes failing. Also they were omitted from documentation.
#. Better handling of ``Infinite`` quantities.
#. Fix ``Precision`` compatibility with WMA.
+#. Prevent that `DirectedInfinity[dir]` be transformed into a `ComplexInfinity`. Rules for `Infinity`, `ComplexInfinity` and `DirectedInfinity wre improved`.
+
PyPI Package requirements
diff --git a/mathics/builtin/arithfns/basic.py b/mathics/builtin/arithfns/basic.py
index f09dd735c..6260c68ec 100644
--- a/mathics/builtin/arithfns/basic.py
+++ b/mathics/builtin/arithfns/basic.py
@@ -6,8 +6,6 @@
"""
-
-import mpmath
import sympy
from mathics.builtin.arithmetic import _MPMathFunction, create_infix
@@ -37,22 +35,19 @@
A_READ_PROTECTED,
)
from mathics.core.convert.expression import to_expression
-from mathics.core.convert.mpmath import from_mpmath
from mathics.core.convert.sympy import from_sympy
+from mathics.core.element import ElementsProperties
from mathics.core.expression import Expression
from mathics.core.list import ListExpression
-from mathics.core.number import dps, min_prec
from mathics.core.symbols import (
Symbol,
SymbolDivide,
SymbolHoldForm,
- SymbolNull,
SymbolPlus,
SymbolPower,
SymbolTimes,
)
from mathics.core.systemsymbols import (
- SymbolAccuracy,
SymbolBlank,
SymbolComplexInfinity,
SymbolDirectedInfinity,
@@ -61,8 +56,17 @@
SymbolInfix,
SymbolLeft,
SymbolMinus,
+ SymbolOverflow,
SymbolPattern,
- SymbolSequence,
+)
+from mathics.eval.arithmetic import (
+ associate_powers,
+ eval_add_numbers,
+ eval_exponential,
+ eval_multiply_numbers,
+ eval_power_inexact,
+ eval_power_number,
+ segregate_numbers_from_sorted_list as segregate_numbers,
)
from mathics.eval.nevaluator import eval_N
from mathics.eval.numerify import numerify
@@ -396,13 +400,13 @@ def eval(self, items, evaluation):
"Plus[items___]"
items_tuple = numerify(items, evaluation).get_sequence()
+ numbers, items_tuple = segregate_numbers(*items_tuple)
elements = []
last_item = last_count = None
+ number = eval_add_numbers(numbers)
- prec = min_prec(*items_tuple)
- is_machine_precision = any(item.is_machine_precision() for item in items_tuple)
- numbers = []
-
+ # This reduces common factors
+ # TODO: Check if it possible to avoid the conversions back and forward to sympy.
def append_last():
if last_item is not None:
if last_count == 1:
@@ -420,65 +424,46 @@ def append_last():
)
for item in items_tuple:
- if isinstance(item, Number):
- numbers.append(item)
+ count = rest = None
+ if item.has_form("Times", None):
+ for element in item.elements:
+ if isinstance(element, Number):
+ count = element.to_sympy()
+ rest = item.get_mutable_elements()
+ rest.remove(element)
+ if len(rest) == 1:
+ rest = rest[0]
+ else:
+ rest.sort()
+ rest = Expression(SymbolTimes, *rest)
+ break
+ if count is None:
+ count = sympy.Integer(1)
+ rest = item
+ if last_item is not None and last_item == rest:
+ last_count = last_count + count
else:
- count = rest = None
- if item.has_form("Times", None):
- for element in item.elements:
- if isinstance(element, Number):
- count = element.to_sympy()
- rest = item.get_mutable_elements()
- rest.remove(element)
- if len(rest) == 1:
- rest = rest[0]
- else:
- rest.sort()
- rest = Expression(SymbolTimes, *rest)
- break
- if count is None:
- count = sympy.Integer(1)
- rest = item
- if last_item is not None and last_item == rest:
- last_count = last_count + count
- else:
- append_last()
- last_item = rest
- last_count = count
+ append_last()
+ last_item = rest
+ last_count = count
append_last()
- if numbers:
- if prec is not None:
- if is_machine_precision:
- numbers = [item.to_mpmath() for item in numbers]
- number = mpmath.fsum(numbers)
- number = from_mpmath(number)
- else:
- # For a sum, what is relevant is the minimum accuracy of the terms
- acc = (
- Expression(SymbolAccuracy, ListExpression(items))
- .evaluate(evaluation)
- .to_python()
- )
- with mpmath.workprec(prec):
- numbers = [item.to_mpmath() for item in numbers]
- number = mpmath.fsum(numbers)
- number = from_mpmath(number, acc=acc)
- else:
- number = from_sympy(sum(item.to_sympy() for item in numbers))
- else:
- number = Integer0
+ # now elements contains the symbolic terms which can not be simplified.
+ # by collecting common symbolic factors.
+ if not elements:
+ return number
- if not number.sameQ(Integer0):
+ if number is not Integer0:
elements.insert(0, number)
-
- if not elements:
- return Integer0
elif len(elements) == 1:
return elements[0]
- else:
- elements.sort()
- return Expression(SymbolPlus, *elements)
+
+ elements.sort()
+ return Expression(
+ SymbolPlus,
+ *elements,
+ elements_properties=ElementsProperties(False, False, True),
+ )
class Power(BinaryOperator, _MPMathFunction):
@@ -614,6 +599,8 @@ class Power(BinaryOperator, _MPMathFunction):
rules = {
"Power[]": "1",
"Power[x_]": "x",
+ "Power[I,-1]": "-I",
+ "Power[-1, 1/2]": "I",
}
summary_text = "exponentiate"
@@ -622,15 +609,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:
@@ -644,17 +631,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)
- )
+ return SymbolComplexInfinity
- result = self.eval(Expression(SymbolSequence, x, y), evaluation)
- if result is None or result != SymbolNull:
- return result
+ # 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):
@@ -932,90 +949,112 @@ def format_outputform(self, items, evaluation):
def eval(self, items, evaluation):
"Times[items___]"
- items = numerify(items, evaluation).get_sequence()
elements = []
numbers = []
infinity_factor = False
- prec = min_prec(*items)
- is_machine_precision = any(item.is_machine_precision() for item in items)
-
# find numbers and simplify Times -> Power
- for item in items:
- if isinstance(item, Number):
- numbers.append(item)
- elif elements and item == elements[-1]:
- elements[-1] = Expression(SymbolPower, elements[-1], Integer2)
- elif (
- elements
- and item.has_form("Power", 2)
- and elements[-1].has_form("Power", 2)
- and item.elements[0].sameQ(elements[-1].elements[0])
- ):
- elements[-1] = Expression(
- SymbolPower,
- elements[-1].elements[0],
- Expression(SymbolPlus, item.elements[1], elements[-1].elements[1]),
- )
- elif (
- elements
- and item.has_form("Power", 2)
- and item.elements[0].sameQ(elements[-1])
- ):
- elements[-1] = Expression(
- SymbolPower,
- elements[-1],
- Expression(SymbolPlus, item.elements[1], Integer1),
- )
- elif (
- elements
- and elements[-1].has_form("Power", 2)
- and elements[-1].elements[0].sameQ(item)
- ):
- elements[-1] = Expression(
- SymbolPower,
- item,
- Expression(SymbolPlus, Integer1, elements[-1].elements[1]),
- )
- elif item.get_head().sameQ(SymbolDirectedInfinity):
+
+ numbers, symbolic_items = segregate_numbers(
+ *(numerify(items, evaluation).get_sequence())
+ )
+
+ # This loop handles factors representing infinite quantities,
+ # and factors which are powers of the same basis.
+ for item in symbolic_items:
+ # Process Infinity factors
+ if item.get_head() is SymbolDirectedInfinity:
+ if len(item.elements) == 0:
+ return (
+ SymbolIndeterminate
+ if any(n.is_zero for n in numbers)
+ else SymbolComplexInfinity
+ )
infinity_factor = True
- if len(item.elements) > 0:
- direction = item.elements[0]
- if isinstance(direction, Number):
- numbers.append(direction)
- else:
- elements.append(direction)
- elif item.sameQ(SymbolInfinity) or item.sameQ(SymbolComplexInfinity):
+ direction = item.elements[0]
+ if isinstance(direction, Number):
+ numbers.append(direction)
+ continue
+ else:
+ item = direction
+ elif item is SymbolInfinity:
infinity_factor = True
+ continue
+ elif item is SymbolComplexInfinity:
+ return (
+ SymbolIndeterminate
+ if any(n.is_zero for n in numbers)
+ else SymbolComplexInfinity
+ )
+ elif item is SymbolIndeterminate:
+ return item
+ # Process powers
+ if elements:
+ previous_elem = elements[-1]
+ if item == previous_elem:
+ elements[-1] = Expression(SymbolPower, previous_elem, Integer2)
+ continue
+ elif item.has_form("Power", 2):
+ base, exp = item.elements
+ if previous_elem.has_form("Power", 2) and base.sameQ(
+ previous_elem.elements[0]
+ ):
+ exp = Expression(
+ SymbolPlus, exp, previous_elem.elements[1]
+ ).evaluate(evaluation)
+ elements[-1] = Expression(
+ SymbolPower,
+ base,
+ exp,
+ )
+ continue
+ if base.sameQ(previous_elem):
+ exp = Expression(SymbolPlus, Integer1, exp).evaluate(evaluation)
+ elements[-1] = Expression(
+ SymbolPower,
+ base,
+ exp,
+ )
+ continue
+ elif previous_elem.has_form("Power", 2) and previous_elem.elements[
+ 0
+ ].sameQ(item):
+ exp = Expression(
+ SymbolPlus, Integer1, previous_elem.elements[1]
+ ).evaluate(evaluation)
+ elements[-1] = Expression(
+ SymbolPower,
+ item,
+ exp,
+ )
+ continue
else:
- elements.append(item)
-
- if numbers:
- if prec is not None:
- if is_machine_precision:
- numbers = [item.to_mpmath() for item in numbers]
- number = mpmath.fprod(numbers)
- number = from_mpmath(number)
+ item = item
+ # Otherwise, just append the element...
+ elements.append(item)
+
+ number = eval_multiply_numbers(numbers)
+
+ if infinity_factor:
+ numeric = [number]
+ non_numeric = []
+ for element in elements:
+ if element.is_numeric(evaluation):
+ numeric.append(element)
else:
- with mpmath.workprec(prec):
- numbers = [item.to_mpmath() for item in numbers]
- number = mpmath.fprod(numbers)
- number = from_mpmath(number, dps(prec))
- else:
- number = sympy.Mul(*[item.to_sympy() for item in numbers])
- number = from_sympy(number)
- else:
+ non_numeric.append(element)
+ elements = non_numeric
+ direction = (
+ numeric[0]
+ if len(numeric) == 1
+ else Times(*numeric).evaluate(evaluation)
+ )
+ elements.append(Expression(SymbolDirectedInfinity, direction))
number = Integer1
-
- if number.sameQ(Integer1):
- number = None
- elif number.is_zero:
- if infinity_factor:
- return SymbolIndeterminate
+ elif len(elements) == 0 or number is Integer0:
return number
- elif (
- number.sameQ(IntegerM1) and elements and elements[0].has_form("Plus", None)
- ):
+
+ if number is IntegerM1 and elements and elements[0].has_form("Plus", None):
elements[0] = Expression(
elements[0].get_head(),
*[
@@ -1023,21 +1062,22 @@ def eval(self, items, evaluation):
for element in elements[0].elements
],
)
- number = None
+ number = Integer1
- if number is not None:
+ if number is not Integer1:
elements.insert(0, number)
- if not elements:
- if infinity_factor:
- return SymbolComplexInfinity
- return Integer1
-
if len(elements) == 1:
- ret = elements[0]
- else:
- ret = Expression(SymbolTimes, *elements)
- if infinity_factor:
- return Expression(SymbolDirectedInfinity, ret)
- else:
- return ret
+ return elements[0]
+
+ elements = sorted(elements)
+ items_elements = items.elements
+ if len(elements) == len(items_elements) and all(
+ elem.sameQ(item) for elem, item in zip(elements, items_elements)
+ ):
+ return None
+ return Expression(
+ SymbolTimes,
+ *elements,
+ elements_properties=ElementsProperties(False, False, True),
+ )
diff --git a/mathics/builtin/arithmetic.py b/mathics/builtin/arithmetic.py
index 28c279d7d..f5f4bfd72 100644
--- a/mathics/builtin/arithmetic.py
+++ b/mathics/builtin/arithmetic.py
@@ -7,6 +7,7 @@
Basic arithmetic functions, including complex number arithmetic.
"""
+from mathics.eval.arithmetic import eval_abs, eval_sign
from mathics.eval.numerify import numerify
# This tells documentation how to sort this module
@@ -56,17 +57,14 @@
from mathics.core.symbols import (
Atom,
Symbol,
- SymbolAbs,
SymbolFalse,
SymbolList,
SymbolPlus,
- SymbolPower,
SymbolTimes,
SymbolTrue,
)
from mathics.core.systemsymbols import (
SymbolAnd,
- SymbolComplexInfinity,
SymbolDirectedInfinity,
SymbolIndeterminate,
SymbolInfix,
@@ -78,6 +76,8 @@
)
from mathics.eval.nevaluator import eval_N
+ExpressionComplexInfinity = Expression(SymbolDirectedInfinity)
+
@lru_cache(maxsize=4096)
def call_mpmath(mpmath_function, mpmath_args):
@@ -86,7 +86,7 @@ def call_mpmath(mpmath_function, mpmath_args):
except ValueError as exc:
text = str(exc)
if text == "gamma function pole":
- return SymbolComplexInfinity
+ return ExpressionComplexInfinity
else:
raise
except ZeroDivisionError:
@@ -130,6 +130,8 @@ def eval(self, z, evaluation):
result = self.prepare_mathics(result)
result = from_sympy(result)
# evaluate elements to convert e.g. Plus[2, I] -> Complex[2, 1]
+ if result is None:
+ return result
return result.evaluate_elements(evaluation)
elif mpmath_function is None:
return
@@ -150,7 +152,7 @@ def eval(self, z, evaluation):
if isinstance(result, (mpmath.mpc, mpmath.mpf)):
if mpmath.isinf(result) and isinstance(result, mpmath.mpc):
- result = SymbolComplexInfinity
+ result = ExpressionComplexInfinity
elif mpmath.isinf(result) and result > 0:
result = Expression(SymbolDirectedInfinity, Integer1)
elif mpmath.isinf(result) and result < 0:
@@ -262,6 +264,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):
"""
@@ -636,7 +645,6 @@ class DirectedInfinity(SympyFunction):
= ComplexInfinity
>> DirectedInfinity[1 + I]
= (1 / 2 + I / 2) Sqrt[2] Infinity
-
>> 1 / DirectedInfinity[1 + I]
= 0
>> DirectedInfinity[1] + DirectedInfinity[-1]
@@ -648,7 +656,7 @@ class DirectedInfinity(SympyFunction):
= Indeterminate
#> DirectedInfinity[1+I]+DirectedInfinity[2+I]
- = (2 / 5 + I / 5) Sqrt[5] Infinity + (1 / 2 + I / 2) Sqrt[2] Infinity
+ = (2 / 5 + I / 5) Sqrt[5] Infinity + (1 / 2 + I / 2) Sqrt[2] Infinity
#> DirectedInfinity[Sqrt[3]]
= Infinity
@@ -658,7 +666,8 @@ 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]",
+ # "0 * DirectedInfinity[args___]": "Message[Infinity::indet, Unevaluated[0 DirectedInfinity[args]]]; Indeterminate",
# "DirectedInfinity[a_?NumericQ] /; N[Abs[a]] != 1": "DirectedInfinity[a / Abs[a]]",
# "DirectedInfinity[a_] * DirectedInfinity[b_]": "DirectedInfinity[a*b]",
# "DirectedInfinity[] * DirectedInfinity[args___]": "DirectedInfinity[]",
@@ -691,24 +700,27 @@ class DirectedInfinity(SympyFunction):
formats = {
"DirectedInfinity[1]": "HoldForm[Infinity]",
- "DirectedInfinity[-1]": "HoldForm[-Infinity]",
+ "DirectedInfinity[-1]": "-HoldForm[Infinity]",
"DirectedInfinity[]": "HoldForm[ComplexInfinity]",
"DirectedInfinity[DirectedInfinity[z_]]": "DirectedInfinity[z]",
"DirectedInfinity[z_?NumericQ]": "HoldForm[z Infinity]",
}
- 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:
@@ -1203,18 +1215,15 @@ 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__]"
diff --git a/mathics/builtin/atomic/numbers.py b/mathics/builtin/atomic/numbers.py
index c5c4b405f..631102e01 100644
--- a/mathics/builtin/atomic/numbers.py
+++ b/mathics/builtin/atomic/numbers.py
@@ -24,12 +24,10 @@
from mathics.builtin.base import Builtin, Predefined, Test
from mathics.core.atoms import (
- Complex,
Integer,
Integer0,
Integer10,
MachineReal,
- MachineReal0,
Number,
Rational,
Real,
@@ -38,18 +36,20 @@
from mathics.core.convert.python import from_python
from mathics.core.expression import Expression
from mathics.core.list import ListExpression
-from mathics.core.number import dps, machine_epsilon, machine_precision
+from mathics.core.number import machine_epsilon, machine_precision
from mathics.core.symbols import Symbol, SymbolDivide
from mathics.core.systemsymbols import (
SymbolIndeterminate,
SymbolInfinity,
SymbolLog,
+ SymbolMachinePrecision,
SymbolN,
SymbolPrecision,
SymbolRealDigits,
SymbolRound,
)
from mathics.eval.nevaluator import eval_N
+from mathics.eval.numbers import MACHINE_PRECISION_VALUE, eval_accuracy, eval_precision
SymbolIntegerDigits = Symbol("IntegerDigits")
SymbolIntegerExponent = Symbol("IntegerExponent")
@@ -203,32 +203,9 @@ class Accuracy(Builtin):
def eval(self, z, evaluation):
"Accuracy[z_]"
- if isinstance(z, Real):
- if z.is_zero:
- return MachineReal(dps(z.get_precision()))
- z_f = z.to_python()
- log10_z = mpmath.log((-z_f if z_f < 0 else z_f), 10)
- return MachineReal(dps(z.get_precision()) - log10_z)
-
- if isinstance(z, Complex):
- acc_real = self.eval(z.real, evaluation)
- acc_imag = self.eval(z.imag, evaluation)
- if acc_real is SymbolInfinity:
- return acc_imag
- if acc_imag is SymbolInfinity:
- return acc_real
- return Real(min(acc_real.to_python(), acc_imag.to_python()))
-
- if isinstance(z, Expression):
- result = None
- for element in z.elements:
- candidate = self.eval(element, evaluation)
- if isinstance(candidate, Real):
- candidate_f = candidate.to_python()
- if result is None or candidate_f < result:
- result = candidate_f
- if result is not None:
- return Real(result)
+ acc = eval_accuracy(z)
+ if acc is not None:
+ return Real(acc)
return SymbolInfinity
@@ -954,37 +931,16 @@ class Precision(Builtin):
/doc/reference-of-built-in-symbols/atomic-elements-of-expressions/representation-of-numbers/accuracy/.
"""
- rules = {
- "Precision[z_?MachineNumberQ]": "MachinePrecision",
- }
-
summary_text = "find the precision of a number"
def eval(self, z, evaluation):
"Precision[z_]"
- if isinstance(z, Real):
- if z.is_zero:
- return MachineReal0
- return MachineReal(dps(z.get_precision()))
-
- if isinstance(z, Complex):
- prec_real = self.eval(z.real, evaluation)
- prec_imag = self.eval(z.imag, evaluation)
- if prec_real is SymbolInfinity:
- return prec_imag
- if prec_imag is SymbolInfinity:
- return prec_real
-
- return Real(min(prec_real.to_python(), prec_imag.to_python()))
-
- if isinstance(z, Expression):
- result = None
- for element in z.elements:
- candidate = self.eval(element, evaluation)
- if isinstance(candidate, Real):
- candidate_f = candidate.to_python()
- if result is None or candidate_f < result:
- result = candidate_f
- if result is not None:
- return Real(result)
- return SymbolInfinity
+ if isinstance(z, MachineReal):
+ return SymbolMachinePrecision
+ prec = eval_precision(z)
+ if prec is None:
+ return SymbolInfinity
+ if prec == MACHINE_PRECISION_VALUE:
+ return SymbolMachinePrecision
+
+ return MachineReal(prec)
diff --git a/mathics/builtin/testing_expressions/numerical_properties.py b/mathics/builtin/testing_expressions/numerical_properties.py
index 2a7af66b4..3c95f38ac 100644
--- a/mathics/builtin/testing_expressions/numerical_properties.py
+++ b/mathics/builtin/testing_expressions/numerical_properties.py
@@ -13,6 +13,7 @@
from mathics.core.expression import Expression
from mathics.core.symbols import 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
@@ -373,6 +374,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:
diff --git a/mathics/core/number.py b/mathics/core/number.py
index ba166ed61..9608a9c67 100644
--- a/mathics/core/number.py
+++ b/mathics/core/number.py
@@ -119,12 +119,10 @@ def prec(dps) -> int:
def min_prec(*args):
- result = None
- for arg in args:
- prec = arg.get_precision()
- if result is None or (prec is not None and prec < result):
- result = prec
- return result
+ args_prec = (arg.get_precision() for arg in args)
+ return min(
+ (arg_prec for arg_prec in args_prec if arg_prec is not None), default=None
+ )
def pickle_mp(value):
diff --git a/mathics/core/systemsymbols.py b/mathics/core/systemsymbols.py
index 05635182e..0c435f139 100644
--- a/mathics/core/systemsymbols.py
+++ b/mathics/core/systemsymbols.py
@@ -24,6 +24,7 @@
# This list is sorted in alphabetic order.
SymbolAborted = Symbol("System`$Aborted")
+SymbolAbs = Symbol("System`Abs")
SymbolAccuracy = Symbol("System`Accuracy")
SymbolAll = Symbol("System`All")
SymbolAlternatives = Symbol("System`Alternatives")
@@ -81,6 +82,7 @@
SymbolEquivalent = Symbol("System`Equivalent")
SymbolEulerGamma = Symbol("System`EulerGamma")
SymbolExactNumberQ = Symbol("System`ExactNumberQ")
+SymbolExp = Symbol("System`Exp")
SymbolExpandAll = Symbol("System`ExpandAll")
SymbolExport = Symbol("System`Export")
SymbolExportString = Symbol("System`ExportString")
@@ -105,6 +107,7 @@
SymbolHoldForm = Symbol("System`HoldForm")
SymbolHoldPattern = Symbol("System`HoldPattern")
SymbolHue = Symbol("System`Hue")
+SymbolI = Symbol("System`I")
SymbolIf = Symbol("System`If")
SymbolIm = Symbol("System`Im")
SymbolImage = Symbol("System`Image")
@@ -122,6 +125,7 @@
SymbolLess = Symbol("System`Less")
SymbolLessEqual = Symbol("System`LessEqual")
SymbolKey = Symbol("System`Key")
+SymbolKhinchin = Symbol("System`Khinchin")
SymbolLine = Symbol("System`Line")
SymbolLog = Symbol("System`Log")
SymbolLog10 = Symbol("System`Log10")
@@ -173,6 +177,7 @@
SymbolPi = Symbol("System`Pi")
SymbolPiecewise = Symbol("System`Piecewise")
SymbolPlot = Symbol("System`Plot")
+SymbolPlus = Symbol("System`Plus")
SymbolPoint = Symbol("System`Point")
SymbolPower = Symbol("System`Power")
SymbolPolygon = Symbol("System`Polygon")
@@ -233,6 +238,7 @@
SymbolTan = Symbol("System`Tan")
SymbolTanh = Symbol("System`Tanh")
SymbolTeXForm = Symbol("System`TeXForm")
+SymbolTimes = Symbol("System`Times")
SymbolThrow = Symbol("System`Throw")
SymbolThreshold = Symbol("System`Threshold")
SymbolToString = Symbol("System`ToString")
diff --git a/mathics/eval/arithmetic.py b/mathics/eval/arithmetic.py
new file mode 100644
index 000000000..fdc191d28
--- /dev/null
+++ b/mathics/eval/arithmetic.py
@@ -0,0 +1,815 @@
+# -*- coding: utf-8 -*-
+
+"""
+helper functions for arithmetic evaluation, which do not
+depends on the evaluation context. Conversions to Sympy are
+used just as a last resource.
+"""
+
+from typing import List, Optional, Tuple
+
+import mpmath
+import sympy
+
+from mathics.core.atoms import (
+ Complex,
+ Integer,
+ Integer0,
+ Integer1,
+ IntegerM1,
+ Number,
+ Rational,
+ RationalOneHalf,
+ Real,
+)
+from mathics.core.convert.mpmath import from_mpmath
+from mathics.core.convert.sympy import from_sympy
+from mathics.core.element import BaseElement
+from mathics.core.expression import Expression
+from mathics.core.number import dps, min_prec
+from mathics.core.rules import Rule
+from mathics.core.symbols import Atom
+from mathics.core.systemsymbols import (
+ SymbolAbs,
+ SymbolE,
+ SymbolEulerGamma,
+ SymbolExp,
+ SymbolI,
+ SymbolKhinchin,
+ SymbolLog,
+ SymbolPi,
+ SymbolPlus,
+ SymbolPower,
+ SymbolSign,
+ SymbolTimes,
+)
+from mathics.eval.numbers import eval_accuracy
+
+RationalMOneHalf = Rational(-1, 2)
+RealM0p5 = Real(-0.5)
+RealOne = Real(1.0)
+
+
+# Here we store numerical constants with their approximate
+# numerical values.
+
+NUMERICAL_CONSTANTS = {
+ SymbolE: 2.718281828,
+ SymbolEulerGamma: 0.5772156649,
+ SymbolKhinchin: 2.685452001,
+ SymbolPi: 3.141592654,
+}
+
+
+# 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 distribute_factor(expr: BaseElement, factor: BaseElement) -> BaseElement:
+ """
+ q * (a + b + c) -> (q a + q b + q c)
+ """
+ if not expr.has_form("Plus", None):
+ return expr
+ terms = (Expression(SymbolTimes, factor, term) for term in expr.elements)
+ return Expression(SymbolPlus, *terms)
+
+
+def distribute_powers(expr: BaseElement) -> BaseElement:
+ """
+ (a b c)^p -> (a^p b^p c^p)
+ """
+ if not expr.has_form("Power", 2):
+ return expr
+ base, exp = expr.elements
+ if not base.has_form("Times", None):
+ return expr
+ factors = (Expression(SymbolPower, factor, exp) for factor in base.elements)
+ return Expression(SymbolTimes, *factors)
+
+
+def eval_abs(expr: BaseElement) -> Optional[BaseElement]:
+ """
+ evaluates Abs[expr]
+ """
+ if isinstance(expr, Number):
+ return eval_abs_number(expr)
+ if expr.has_form("Power", 2):
+ base, exp = expr.elements
+ if exp.is_zero:
+ return Integer1
+ if test_arithmetic_expr(expr):
+ abs_base = eval_abs(base)
+ if abs_base is None:
+ abs_base = Expression(SymbolAbs, base)
+ return Expression(SymbolPower, abs_base, exp)
+ if expr.get_head() is SymbolTimes:
+ return eval_multiply_numbers([eval_abs(x) for x in expr.elements])
+ if test_nonnegative_arithmetic_expr(expr):
+ return expr
+ if test_negative_arithmetic_expr(expr):
+ return eval_multiply_numbers([IntegerM1, expr])
+ return None
+
+
+def eval_abs_number(n: Number) -> Number:
+ """
+ Evals the absolute value of a number
+ """
+ if isinstance(n, Integer):
+ n_val = n.value
+ if n_val >= 0:
+ return n
+ return Integer(-n_val)
+ if isinstance(n, Rational):
+ n_num, n_den = n.value.as_numer_denom()
+ if n_num >= 0:
+ return n
+ return Rational(-n_num, n_den)
+ if isinstance(n, Real):
+ n_val = n.value
+ if n_val >= 0:
+ return n
+ return eval_multiply_numbers([IntegerM1, n])
+ if isinstance(n, Complex):
+ if n.real.is_zero:
+ return eval_abs_number(n.imag)
+ sq_comp = tuple((eval_multiply_numbers([x, x]) for x in (n.real, n.imag)))
+ sq_abs = eval_add_numbers(sq_comp)
+ result = eval_power_number(sq_abs, RationalOneHalf) or Expression(
+ SymbolPower, sq_abs, RationalOneHalf
+ )
+ return result
+
+
+def eval_add_numbers(
+ numbers: List[Number],
+) -> Number:
+ """
+ Add the elements in ``numbers``.
+ """
+ if len(numbers) == 0:
+ return Integer0
+
+ prec = min_prec(*numbers)
+ if prec is not None:
+ is_machine_precision = any(number.is_machine_precision() for number in numbers)
+ if is_machine_precision:
+ numbers = (item.to_mpmath() for item in numbers)
+ number = mpmath.fsum(numbers)
+ return from_mpmath(number)
+ else:
+ # For a sum, what is relevant is the minimum accuracy of the terms
+ item_accuracies = (eval_accuracy(item) for item in numbers)
+ acc = min(
+ (item_acc for item_acc in item_accuracies if item_acc is not None),
+ default=None,
+ )
+ with mpmath.workprec(prec):
+ numbers = (item.to_mpmath() for item in numbers)
+ number = mpmath.fsum(numbers)
+ return from_mpmath(number, acc=acc)
+ else:
+ return from_sympy(sum(item.to_sympy() for item in numbers))
+
+
+def eval_complex_conjugate(z: Number) -> Number:
+ """
+ Evaluates the complex conjugate of z.
+ """
+ if isinstance(z, Complex):
+ re, im = z.real, z.imag
+ return Complex(re, eval_negate_number(im))
+ return z
+
+
+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, dps(prec))
+
+
+def eval_inverse_number(n: Number) -> Number:
+ """
+ Eval 1/n
+ """
+ if isinstance(n, Integer):
+ n_value = n.value
+ if n_value == 1 or n_value == -1:
+ return n
+ return Rational(-1, -n_value) if n_value < 0 else Rational(1, n_value)
+ if isinstance(n, Rational):
+ n_num, n_den = n.value.as_numer_denom()
+ if n_num < 0:
+ n_num, n_den = -n_num, -n_den
+ if n_num == 1:
+ return Integer(n_den)
+ return Rational(n_den, n_num)
+ # Otherwise, use power....
+ return eval_power_number(n, IntegerM1)
+
+
+def eval_multiply_numbers(numbers: List[Number]) -> Number:
+ """
+ Multiply the elements in ``numbers``.
+ """
+ if len(numbers) == 0:
+ return Integer1
+
+ prec = min_prec(*numbers)
+ if prec is not None:
+ is_machine_precision = any(number.is_machine_precision() for number in numbers)
+ if is_machine_precision:
+ numbers = (item.to_mpmath() for item in numbers)
+ number = mpmath.fprod(numbers)
+ return from_mpmath(number)
+ else:
+ with mpmath.workprec(prec):
+ numbers = (item.to_mpmath() for item in numbers)
+ number = mpmath.fprod(numbers)
+ return from_mpmath(number, dps(prec))
+ else:
+ return from_sympy(sympy.Mul(*(item.to_sympy() for item in numbers)))
+
+
+def eval_negate_number(n: Number) -> Number:
+ """
+ Changes the sign of n
+ """
+ if isinstance(n, Integer):
+ return Integer(-n.value)
+ if isinstance(n, Rational):
+ n_num, n_den = n.value.as_numer_denom()
+ return Rational(-n_num, n_den)
+ # Otherwise, multiply by -1:
+ return eval_multiply_numbers([IntegerM1, n])
+
+
+def eval_power_number(base: Number, exp: Number) -> Optional[Number]:
+ """
+ Eval base^exp for `base` and `exp` two numbers. If the expression
+ remains the same, return None.
+ """
+ # If both base and exponent are exact quantities,
+ # use sympy.
+ # If base or exp are inexact quantities, use
+ # the inexact routine.
+ if base.is_inexact() or exp.is_inexact():
+ return eval_power_inexact(base, exp)
+
+ # Trivial special cases
+ if exp is Integer1:
+ return base
+ if exp is Integer0:
+ return Integer1
+ if base is Integer1:
+ return Integer1
+
+ def eval_power_sympy() -> Optional[Number]:
+ """
+ Tries to compute x^p using sympy rules.
+ If the answer is again x^p, return None.
+ """
+ # This function is called just if useful native rules
+ # are available.
+ result = from_sympy(sympy.Pow(base.to_sympy(), exp.to_sympy()))
+ if result.has_form("Power", 2):
+ # If the expression didn´t change, return None
+ if result.elements[0].sameQ(base):
+ return None
+ return result
+
+ # Rational exponent
+ if isinstance(exp, Rational):
+ exp_p, exp_q = exp.value.as_numer_denom()
+ if abs(exp_p) > exp_q:
+ exp_int, exp_num = divmod(exp_p, exp_q)
+ exp_rem = Rational(exp_num, exp_q)
+ factor_1 = eval_power_number(base, Integer(exp_int))
+ factor_2 = eval_power_number(base, exp_rem) or Expression(
+ SymbolPower, base, exp_rem
+ )
+ if factor_1 is Integer1:
+ return factor_2
+ return Expression(SymbolTimes, factor_1, factor_2)
+
+ # Integer base
+ if isinstance(base, Integer):
+ base_value = base.value
+ if base_value == -1:
+ if isinstance(exp, Rational):
+ if exp.sameQ(RationalOneHalf):
+ return SymbolI
+ return None
+ return eval_power_sympy()
+ elif base_value < 0:
+ neg_base = eval_negate_number(base)
+ candidate = eval_power_number(neg_base, exp)
+ if candidate is None:
+ return None
+ sign_factor = eval_power_number(IntegerM1, exp)
+ if candidate is Integer1:
+ return sign_factor
+ return Expression(SymbolTimes, candidate, sign_factor)
+
+ # Rational base
+ if isinstance(base, Rational):
+ # If the exponent is an Integer or Rational negative value
+ # restate as a positive power
+ if (
+ isinstance(exp, Integer)
+ and exp.value < 0
+ or isinstance(exp, Rational)
+ and exp.value.p < 0
+ ):
+ base, exp = eval_inverse_number(base), eval_negate_number(exp)
+ return eval_power_number(base, exp) or Expression(SymbolPower, base, exp)
+
+ p, q = (Integer(u) for u in base.value.as_numer_denom())
+ p_eval, q_eval = (eval_power_number(u, exp) for u in (p, q))
+ # If neither p^exp or q^exp produced a new result,
+ # leave it alone
+ if q_eval is None and p_eval is None:
+ return None
+ # if q^exp == 1: return p_eval
+ # (should not happen)
+ if q_eval is Integer1:
+ return p_eval
+ if isinstance(q_eval, Integer):
+ if isinstance(p_eval, Integer):
+ return Rational(p_eval.value, q_eval.value)
+
+ if p_eval is None:
+ p_eval = Expression(SymbolPower, p, exp)
+
+ if q_eval is None:
+ q_eval = Expression(SymbolPower, q, exp)
+ return Expression(
+ SymbolTimes, p_eval, Expression(SymbolPower, q_eval, IntegerM1)
+ )
+ # Pure imaginary base case
+ elif isinstance(base, Complex) and base.real.is_zero:
+ base = base.imag
+ if base.value < 0:
+ base = eval_negate_number(base)
+ phase = Expression(
+ SymbolPower,
+ IntegerM1,
+ eval_multiply_numbers([IntegerM1, RationalOneHalf, exp]),
+ )
+ else:
+ phase = Expression(
+ SymbolPower, IntegerM1, eval_multiply_numbers([RationalOneHalf, exp])
+ )
+ real_factor = eval_power_number(base, exp)
+
+ if real_factor is None:
+ return None
+ return Expression(SymbolTimes, real_factor, phase)
+
+ # Generic case
+ return eval_power_sympy()
+
+
+def eval_power_inexact(base: Number, exp: Number) -> BaseElement:
+ """
+ Eval base^exp for `base` and `exp` inexact numbers
+ """
+ # If both base and exponent are exact quantities,
+ # use sympy.
+ prec = min_prec(base, exp)
+ if prec is not None:
+ is_machine_precision = base.is_machine_precision() or exp.is_machine_precision()
+ if is_machine_precision:
+ number = mpmath.power(base.to_mpmath(), exp.to_mpmath())
+ return from_mpmath(number)
+ else:
+ with mpmath.workprec(prec):
+ number = mpmath.power(base.to_mpmath(), exp.to_mpmath())
+ return from_mpmath(number, dps(prec))
+
+
+def eval_sign(expr: BaseElement) -> Optional[BaseElement]:
+ """
+ evaluates Sign[expr]
+ """
+ if isinstance(expr, Atom):
+ return eval_sign_number(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.has_form("Exp", 1):
+ exp = expr.elements[0]
+ if isinstance(exp, (Integer, Real, Rational)):
+ return Integer1
+ if isinstance(exp, Complex):
+ return Expression(SymbolExp, exp.imag)
+ 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_nonnegative_arithmetic_expr(expr):
+ return Integer1
+ if test_negative_arithmetic_expr(expr):
+ return IntegerM1
+ if test_zero_arithmetic_expr:
+ return Integer0
+ return None
+
+
+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)])
+
+
+# TODO: Annotate me
+def segregate_numbers(
+ *elements: BaseElement,
+) -> Tuple[List[Number], List[BaseElement]]:
+ """
+ From a list of elements, produce two lists, one with the numeric items
+ and the other with the remaining
+ """
+ items = {True: [], False: []}
+ for element in elements:
+ items[isinstance(element, Number)].append(element)
+ return items[True], items[False]
+
+
+# TODO: Annotate me
+def segregate_numbers_from_sorted_list(
+ *elements: BaseElement,
+) -> Tuple[List[Number], List[BaseElement]]:
+ """
+ From a list of elements, produce two lists, one with the numeric items
+ and the other with the remaining. Differently than `segregate_numbers`,
+ this function assumes that elements are sorted with the numbers at
+ the begining.
+ """
+ for pos, element in enumerate(elements):
+ if not isinstance(element, Number):
+ return list(elements[:pos]), list(elements[pos:])
+ return list(elements), []
+
+
+def flat_arithmetic_operators(expr: Expression) -> Expression:
+ """
+ operation[a_number, b, operation[c_number, d], e]-> operation[a, c, b, c, d, e]
+ """
+ # items is a dict with two keys: True and False.
+ # In True we store numeric items, and in False the symbolic ones.
+ items = {True: [], False: []}
+ head = expr.get_head()
+ for element in expr.elements:
+ # If the element is also head[elements],
+ # take its elements, and append to the main expression.
+ if element.get_head() is head:
+ for item in flat_arithmetic_operators(element).elements:
+ item[isinstance(item, Number)].append(item)
+ item[isinstance(item, Number)].append(item)
+ return Expression(head, *items[True], *items[False])
+
+
+def test_arithmetic_expr(expr: BaseElement, only_real: bool = True) -> bool:
+ """
+ Check if an expression `expr` is an arithmetic expression
+ composed only by numbers and arithmetic operations.
+ If only_real is set to True, then `I` is not considered a number.
+ """
+ if isinstance(expr, (Integer, Rational, Real)):
+ return True
+ if expr in NUMERICAL_CONSTANTS:
+ return True
+ if isinstance(expr, Complex):
+ return not only_real
+
+ head, elements = expr.get_head(), expr.elements
+
+ if head in (SymbolPlus, SymbolTimes):
+ return all(test_arithmetic_expr(term, only_real) for term in elements)
+ if expr.has_form("Exp", 1):
+ return test_arithmetic_expr(elements[0], only_real)
+ if head is SymbolLog:
+ if len(elements) > 2:
+ return False
+ if len(elements) == 2:
+ base = elements[0]
+ if not test_positive_arithmetic_expr(base):
+ return False
+ return test_arithmetic_expr(elements[-1], only_real)
+ if expr.has_form("Power", 2):
+ base, exponent = elements
+ if only_real:
+ if isinstance(exponent, Integer):
+ return test_arithmetic_expr(base)
+ return all(test_arithmetic_expr(item, only_real) for item in elements)
+ if not only_real:
+ if expr is SymbolI or isinstance(expr, Complex):
+ return True
+ return False
+
+
+def test_negative_arithmetic_expr(expr: BaseElement) -> bool:
+ """
+ Check if the expression is an arithmetic expression
+ representing a negative value.
+ """
+ if isinstance(expr, (Integer, Rational, Real)):
+ return expr.value < 0
+
+ expr = eval_multiply_numbers([IntegerM1, expr])
+ return test_positive_arithmetic_expr(expr)
+
+
+def test_nonnegative_arithmetic_expr(expr: BaseElement) -> bool:
+ """
+ Check if the expression is an arithmetic expression
+ representing a nonnegative number
+ """
+ if test_zero_arithmetic_expr(expr) or test_positive_arithmetic_expr(expr):
+ return True
+
+
+def test_nonpositive_arithetic_expr(expr: BaseElement) -> bool:
+ """
+ Check if the expression is an arithmetic expression
+ representing a nonnegative number
+ """
+ if test_zero_arithmetic_expr(expr) or test_negative_arithmetic_expr(expr):
+ return True
+ return False
+
+
+def test_positive_arithmetic_expr(expr: BaseElement) -> bool:
+ """
+ Check if the expression is an arithmetic expression
+ representing a positive value.
+ """
+ if isinstance(expr, (Integer, Rational, Real)):
+ return expr.value > 0
+ if expr in NUMERICAL_CONSTANTS:
+ return True
+ if isinstance(expr, Atom):
+ return False
+
+ head, elements = expr.get_head(), expr.elements
+ if head is SymbolPlus:
+ positive_nonpositive_terms = {True: [], False: []}
+ for term in elements:
+ positive_nonpositive_terms[test_positive_arithmetic_expr(term)].append(term)
+
+ if len(positive_nonpositive_terms[False]) == 0:
+ return True
+ if len(positive_nonpositive_terms[True]) == 0:
+ return False
+
+ pos, neg = (
+ eval_add_numbers(items) for items in positive_nonpositive_terms.values()
+ )
+ if neg.is_zero:
+ return True
+ if not test_arithmetic_expr(neg):
+ return False
+
+ total = eval_add_numbers([pos, neg])
+ # Check positivity of the evaluated expression
+ if isinstance(total, (Integer, Rational, Real)):
+ return total.value > 0
+ if isinstance(total, Complex):
+ return False
+ if total.sameQ(expr):
+ return False
+ return test_positive_arithmetic_expr(total)
+
+ if head is SymbolTimes:
+ nonpositive_factors = tuple(
+ (item for item in elements if not test_positive_arithmetic_expr(item))
+ )
+ if len(nonpositive_factors) == 0:
+ return True
+ evaluated_expr = eval_multiply_numbers(nonpositive_factors)
+ if evaluated_expr.sameQ(expr):
+ return False
+ return test_positive_arithmetic_expr(evaluated_expr)
+ if expr.has_form("Power", 2):
+ base, exponent = elements
+ if isinstance(exponent, Integer) and exponent.value % 2 == 0:
+ return test_arithmetic_expr(base)
+ return test_arithmetic_expr(exponent) and test_positive_arithmetic_expr(base)
+ if expr.has_form("Exp", 1):
+ return test_arithmetic_expr(exponent)
+ if head is SymbolLog:
+ if len(elements) > 2:
+ return False
+ if len(elements) == 2:
+ if not test_positive_arithmetic_expr(elements[0]):
+ return False
+ arg = elements[-1]
+ return test_positive_arithmetic_expr(eval_add_numbers([arg, IntegerM1]))
+ if head.has_form("Abs", 1):
+ return True
+ if head.has_form("DirectedInfinity", 1):
+ return test_positive_arithmetic_expr(elements[0])
+
+ return False
+
+
+def test_zero_arithmetic_expr(expr: BaseElement, numeric: bool = False) -> bool:
+ """
+ return True if expr evaluates to a number compatible
+ with 0
+ """
+
+ def is_numeric_zero(z: Number):
+ if isinstance(z, Complex):
+ if abs(z.real.value) + abs(z.imag.value) < 2.0e-10:
+ return True
+ if isinstance(z, Number):
+ if abs(z.value) < 1e-10:
+ return True
+ return False
+
+ if expr.is_zero:
+ return True
+ if numeric:
+ if is_numeric_zero(expr):
+ return True
+ expr = to_inexact_value(expr)
+ if expr.has_form("Times", None):
+ if any(
+ test_zero_arithmetic_expr(element, numeric=numeric)
+ for element in expr.elements
+ ) and not any(
+ element.has_form("DirectedInfinity", None) for element in expr.elements
+ ):
+ return True
+ if expr.has_form("Power", 2):
+ base, exp = expr.elements
+ return test_zero_arithmetic_expr(base, numeric) and test_arithmetic_expr(
+ exp, numeric
+ )
+ if expr.has_form("Plus", None):
+ result = eval_add_numbers(expr.elements)
+ if numeric:
+ if isinstance(result, complex):
+ if abs(result.real.value) + abs(result.imag.value) < 2.0e-10:
+ return True
+ if isinstance(result, Number):
+ if abs(result.value) < 1e-10:
+ return True
+ return result.is_zero
+ return False
+
+
+def to_inexact_value(expr: BaseElement) -> BaseElement:
+ """
+ Converts an expression into an inexact expression.
+ Replaces numerical constants by their numerical approximation,
+ and then multiplies the expression by Real(1.)
+ """
+ if expr.is_inexact():
+ return expr
+
+ if isinstance(expr, Expression):
+ for const, value in NUMERICAL_CONSTANTS.items():
+ expr, success = expr.do_apply_rule(Rule(const, Real(value)))
+ return eval_multiply_numbers([RealOne, expr])
diff --git a/mathics/eval/numbers.py b/mathics/eval/numbers.py
index 4e1db1dc3..48beeff26 100644
--- a/mathics/eval/numbers.py
+++ b/mathics/eval/numbers.py
@@ -1,9 +1,102 @@
+from typing import Optional
+
+import mpmath
import sympy
+from mathics.core.atoms import Complex, MachineReal, PrecisionReal, Real
from mathics.core.convert.sympy import from_sympy
+from mathics.core.element import BaseElement
from mathics.core.expression import Expression
+from mathics.core.number import dps, machine_precision
from mathics.core.symbols import SymbolPlus
+MACHINE_PRECISION_VALUE = machine_precision * mpmath.log(2, 10.0)
+
+
+def eval_accuracy(z: BaseElement) -> Optional[float]:
+ """
+ Determine the accuracy of an expression expr.
+ If z is a Real value, returns the difference between
+ the number of significant decimal figures (Precision) and
+ log_10(z).
+
+ For example,
+ ```
+ 12.345`2
+ ```
+ which is equivalent to 12.`2 has an accuracy of
+ ```
+ 0.908509 == 2. - log(10, 12.345)
+ ```
+
+ If the expression contains Real values, returns
+ the minimal accuracy of all the numbers in the expression.
+
+ Otherwise returns None, representing infinite accuracy.
+ """
+ if isinstance(z, Real):
+ if z.is_zero:
+ # WMA returns 323.607 ~ $MachinePrecision - Log[10, 2.225073*^-308]
+ # i.e. the accuracy for the smallest positive double precision value.
+ return dps(z.get_precision())
+ z_f = z.to_python()
+ log10_z = mpmath.log((-z_f if z_f < 0 else z_f), 10)
+ return dps(z.get_precision()) - log10_z
+
+ if isinstance(z, Complex):
+ acc_real = eval_accuracy(z.real)
+ acc_imag = eval_accuracy(z.imag)
+ if acc_real is None:
+ return acc_imag
+ if acc_imag is None:
+ return acc_real
+ return min(acc_real, acc_imag)
+
+ if isinstance(z, Expression):
+ elem_accuracies = (eval_accuracy(z_elem) for z_elem in z.elements)
+ return min((acc for acc in elem_accuracies if acc is not None), default=None)
+ return None
+
+
+def eval_precision(z: BaseElement) -> Optional[float]:
+ """
+ Determine the precision of an expression expr.
+ If z is a Real value, returns the number of significant
+ decimal figures of z. For example,
+ ```
+ 12.345`2
+ ```
+ which is equivalent to 12.`2 has a precision of 2.
+
+ If the expression contains Real values, returns
+ the minimal accuracy of all the numbers in the expression.
+
+ Otherwise returns None, representing infinite precision.
+ """
+
+ if isinstance(z, MachineReal):
+ return MACHINE_PRECISION_VALUE
+
+ if isinstance(z, PrecisionReal):
+ if z.is_zero:
+ return MACHINE_PRECISION_VALUE
+ return float(dps(z.get_precision()))
+
+ if isinstance(z, Complex):
+ prec_real = eval_precision(z.real)
+ prec_imag = eval_precision(z.imag)
+ if prec_real is None:
+ return prec_imag
+ if prec_imag is None:
+ return prec_real
+ return min(prec_real, prec_imag)
+
+ if isinstance(z, Expression):
+ elem_prec = (eval_precision(z_elem) for z_elem in z.elements)
+ return min((prec for prec in elem_prec if prec is not None), default=None)
+
+ return None
+
def cancel(expr):
if expr.has_form("Plus", None):
diff --git a/test/builtin/arithmetic/test_abs.py b/test/builtin/arithmetic/test_abs.py
index f17bb24ea..890e80d6f 100644
--- a/test/builtin/arithmetic/test_abs.py
+++ b/test/builtin/arithmetic/test_abs.py
@@ -39,8 +39,8 @@ def test_abs(str_expr, str_expected, msg):
("Sign[2+3 I]", "(2 + 3 I)/(13^(1/2))", None),
("Sign[2.+3 I]", "0.5547 + 0.83205 I", None),
("Sign[4^(2 Pi)]", "1", None),
+ ("Sign[I^(2 Pi)]", "I^(2 Pi)", None),
# FixME: add rules to handle this kind of case
- # ("Sign[I^(2 Pi)]", "I^(2 Pi)", None),
# ("Sign[4^(2 Pi I)]", "1", None),
],
)
diff --git a/test/builtin/arithmetic/test_basic.py b/test/builtin/arithmetic/test_basic.py
index cc771eb00..0708b02c7 100644
--- a/test/builtin/arithmetic/test_basic.py
+++ b/test/builtin/arithmetic/test_basic.py
@@ -129,7 +129,6 @@ def test_exponential(str_expr, str_expected):
# ("a b DirectedInfinity[-3]", "a b (-Infinity)", ""),
],
)
-@pytest.mark.xfail
def test_multiply(str_expr, str_expected, msg):
check_evaluation(str_expr, str_expected, failure_message=msg, hold_expected=True)
@@ -190,6 +189,5 @@ def test_multiply(str_expr, str_expected, msg):
("(a^(p-2 q))^3.", "(a ^ (p - 2 q)) ^ 3.", None),
],
)
-@pytest.mark.xfail
def test_power(str_expr, str_expected, msg):
check_evaluation(str_expr, str_expected, failure_message=msg)
diff --git a/test/builtin/atomic/test_numbers.py b/test/builtin/atomic/test_numbers.py
index 9039abc63..0bdd83ee0 100644
--- a/test/builtin/atomic/test_numbers.py
+++ b/test/builtin/atomic/test_numbers.py
@@ -205,9 +205,9 @@ def test_precision():
for str_expr, str_expected in (
("0`4", "0"),
("Precision[0.0]", "MachinePrecision"),
- ("Precision[0.000000000000000000000000000000000000]", "0."),
+ ("Precision[0.000000000000000000000000000000000000]", "MachinePrecision"),
("Precision[-0.0]", "MachinePrecision"),
- ("Precision[-0.000000000000000000000000000000000000]", "0."),
+ ("Precision[-0.000000000000000000000000000000000000]", "MachinePrecision"),
("1.0000000000000000 // Precision", "MachinePrecision"),
("1.00000000000000000 // Precision", "17."),
(" 0.4 + 2.4 I // Precision", "MachinePrecision"),
diff --git a/test/format/test_format.py b/test/format/test_format.py
index eaca6fec2..2130d1ed6 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