From 17e6ed1bf0e98bb67e29c96ce7a9ea4b7f18e177 Mon Sep 17 00:00:00 2001 From: mmatera Date: Thu, 25 May 2023 22:52:17 -0300 Subject: [PATCH 01/14] test_zero_arithmetic_expr --- mathics/builtin/arithmetic.py | 22 +- .../numerical_properties.py | 5 + mathics/core/systemsymbols.py | 6 + mathics/eval/arithmetic.py | 232 +++++++++++++++++- test/builtin/arithmetic/test_abs.py | 3 +- 5 files changed, 258 insertions(+), 10 deletions(-) diff --git a/mathics/builtin/arithmetic.py b/mathics/builtin/arithmetic.py index 5a85bbe7e..152e16d9f 100644 --- a/mathics/builtin/arithmetic.py +++ b/mathics/builtin/arithmetic.py @@ -629,10 +629,16 @@ class DirectedInfinity(SympyFunction): rules = { "DirectedInfinity[args___] ^ -1": "0", # Special arguments: - "DirectedInfinity[DirectedInfinity[args___]]": "DirectedInfinity[args]", "DirectedInfinity[Indeterminate]": "Indeterminate", "DirectedInfinity[Alternatives[0, 0.]]": "DirectedInfinity[]", # Plus + "DirectedInfinity[DirectedInfinity[args___]]": "DirectedInfinity[args]", + # "DirectedInfinity[a_?NumericQ] /; N[Abs[a]] != 1": "DirectedInfinity[a / Abs[a]]", + # "DirectedInfinity[a_] * DirectedInfinity[b_]": "DirectedInfinity[a*b]", + # "DirectedInfinity[] * DirectedInfinity[args___]": "DirectedInfinity[]", + # Rules already implemented in Times.eval + # "z_?NumberQ * DirectedInfinity[]": "DirectedInfinity[]", + # "z_?NumberQ * DirectedInfinity[a_]": "DirectedInfinity[z * a]", "DirectedInfinity[a_] + DirectedInfinity[b_] /; b == -a": ( "Message[Infinity::indet," " Unevaluated[DirectedInfinity[a] + DirectedInfinity[b]]];" @@ -644,23 +650,27 @@ class DirectedInfinity(SympyFunction): "Indeterminate" ), "DirectedInfinity[args___] + _?NumberQ": "DirectedInfinity[args]", - # Times. See if can be reinstalled in eval_Times + "DirectedInfinity[Alternatives[0, 0.]]": "DirectedInfinity[]", + "DirectedInfinity[0.]": ( + "Message[Infinity::indet," + " Unevaluated[DirectedInfinity[0.]]];" + "Indeterminate" + ), "Alternatives[0, 0.] DirectedInfinity[z___]": ( "Message[Infinity::indet," " Unevaluated[0 DirectedInfinity[z]]];" "Indeterminate" ), - "a_?NumericQ * DirectedInfinity[b_]": "DirectedInfinity[a * b]", "a_ DirectedInfinity[]": "DirectedInfinity[]", "DirectedInfinity[a_] * DirectedInfinity[b_]": "DirectedInfinity[a * b]", + "a_?NumericQ * DirectedInfinity[b_]": "DirectedInfinity[a * b]", } formats = { "DirectedInfinity[1]": "HoldForm[Infinity]", - "DirectedInfinity[-1]": "HoldForm[-Infinity]", + "DirectedInfinity[-1]": "PrecedenceForm[-HoldForm[Infinity], 390]", "DirectedInfinity[]": "HoldForm[ComplexInfinity]", - "DirectedInfinity[DirectedInfinity[z_]]": "DirectedInfinity[z]", - "DirectedInfinity[z_?NumericQ]": "HoldForm[z Infinity]", + "DirectedInfinity[z_]": "PrecedenceForm[z HoldForm[Infinity], 390]", } def eval_complex_infinity(self, evaluation: Evaluation): diff --git a/mathics/builtin/testing_expressions/numerical_properties.py b/mathics/builtin/testing_expressions/numerical_properties.py index 8de45e8b1..5f5822f9e 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 BooleanType, 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 @@ -460,6 +461,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/systemsymbols.py b/mathics/core/systemsymbols.py index 52659e775..f463e5f91 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") @@ -84,6 +85,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") @@ -109,6 +111,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") @@ -126,6 +129,7 @@ SymbolLess = Symbol("System`Less") SymbolLessEqual = Symbol("System`LessEqual") SymbolKey = Symbol("System`Key") +SymbolKhinchin = Symbol("System`Khinchin") SymbolLetterCharacter = Symbol("System`LetterCharacter") SymbolLine = Symbol("System`Line") SymbolLog = Symbol("System`Log") @@ -179,6 +183,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") @@ -242,6 +247,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 index 0b9564067..399d75b86 100644 --- a/mathics/eval/arithmetic.py +++ b/mathics/eval/arithmetic.py @@ -30,10 +30,30 @@ from mathics.core.element import BaseElement, ElementsProperties from mathics.core.expression import Expression from mathics.core.number import FP_MANTISA_BINARY_DIGITS, SpecialValueError, min_prec -from mathics.core.symbols import Symbol, SymbolPlus, SymbolPower, SymbolTimes -from mathics.core.systemsymbols import SymbolComplexInfinity, SymbolIndeterminate +from mathics.core.rules import Rule +from mathics.core.symbols import Atom, Symbol, SymbolPlus, SymbolPower, SymbolTimes +from mathics.core.systemsymbols import ( + SymbolComplexInfinity, + SymbolE, + SymbolEulerGamma, + SymbolI, + SymbolIndeterminate, + SymbolKhinchin, + SymbolLog, + SymbolPi, +) + +# FIXME: replace by numpy constants: +NUMERICAL_CONSTANTS = { + SymbolE: 2.718281828, + SymbolEulerGamma: 0.5772156649, + SymbolKhinchin: 2.685452001, + SymbolPi: 3.141592654, +} + RationalMOneHalf = Rational(-1, 2) +RealOne = Real(1.0) # This cache might not be used that much. @@ -374,3 +394,211 @@ def segregate_numbers_from_sorted_list( if not isinstance(element, Number): return list(elements[:pos]), list(elements[pos:]) return list(elements), [] + + +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 + if isinstance(expr, Symbol): + return False + + head, elements = expr.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/test/builtin/arithmetic/test_abs.py b/test/builtin/arithmetic/test_abs.py index c523bf1d3..11b3da92b 100644 --- a/test/builtin/arithmetic/test_abs.py +++ b/test/builtin/arithmetic/test_abs.py @@ -42,8 +42,7 @@ 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), - # FIXME: add rules to handle this kind of case - # ("Sign[I^(2 Pi)]", "I^(2 Pi)", None), + ("Sign[I^(2 Pi)]", "I^(2 Pi)", None), # ("Sign[4^(2 Pi I)]", "1", None), ], ) From 511073947fe301be62deab6d06194e9f4f3bd382 Mon Sep 17 00:00:00 2001 From: mmatera Date: Fri, 26 May 2023 10:54:20 -0300 Subject: [PATCH 02/14] improve positivity test for powers. Handle constants --- mathics/builtin/base.py | 6 ++- mathics/builtin/numbers/constants.py | 80 +++++++++++++++++----------- mathics/core/atoms.py | 10 ++++ mathics/core/number.py | 4 +- mathics/eval/arithmetic.py | 40 ++++++-------- 5 files changed, 83 insertions(+), 57 deletions(-) diff --git a/mathics/builtin/base.py b/mathics/builtin/base.py index 06a4337b0..8eaccab79 100644 --- a/mathics/builtin/base.py +++ b/mathics/builtin/base.py @@ -799,10 +799,14 @@ def get_operator_display(self) -> Optional[str]: class Predefined(Builtin): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self.symbol = Symbol(self.get_name()) + def get_functions(self, prefix="eval", is_pymodule=False) -> List[Callable]: functions = list(super().get_functions(prefix)) if prefix == "eval" and hasattr(self, "evaluate"): - functions.append((Symbol(self.get_name()), self.evaluate)) + functions.append((self.symbol, self.evaluate)) return functions diff --git a/mathics/builtin/numbers/constants.py b/mathics/builtin/numbers/constants.py index 787826869..a29aefe81 100644 --- a/mathics/builtin/numbers/constants.py +++ b/mathics/builtin/numbers/constants.py @@ -9,25 +9,19 @@ # This tells documentation how to sort this module sort_order = "mathics.builtin.mathematical-constants" - import math +from typing import Optional import mpmath import numpy import sympy from mathics.builtin.base import Builtin, Predefined, SympyObject -from mathics.core.atoms import MachineReal, PrecisionReal +from mathics.core.atoms import NUMERICAL_CONSTANTS, MachineReal, PrecisionReal from mathics.core.attributes import A_CONSTANT, A_PROTECTED, A_READ_PROTECTED +from mathics.core.element import BaseElement from mathics.core.evaluation import Evaluation -from mathics.core.number import ( - MACHINE_DIGITS, - MAX_MACHINE_NUMBER, - MIN_MACHINE_NUMBER, - PrecisionValueError, - get_precision, - prec, -) +from mathics.core.number import MACHINE_DIGITS, PrecisionValueError, get_precision, prec from mathics.core.symbols import Atom, Symbol, strip_context from mathics.core.systemsymbols import SymbolIndeterminate @@ -89,28 +83,33 @@ def eval_N(self, precision, evaluation): def is_constant(self) -> bool: return True - def get_constant(self, precision, evaluation): + def get_constant( + self, + precision: Optional[BaseElement] = None, + evaluation: Optional[Evaluation] = None, + ): # first, determine the precision d = None - if precision: - try: - d = get_precision(precision, evaluation) - except PrecisionValueError: - pass + preference = None + if evaluation: + if precision: + try: + d = get_precision(precision, evaluation) + except PrecisionValueError: + pass + + preflist = evaluation._preferred_n_method.copy() + while preflist: + pref_method = preflist.pop() + if pref_method in ("numpy", "mpmath", "sympy"): + preference = pref_method + break if d is None: d = MACHINE_DIGITS # If preference not especified, determine it # from the precision. - preference = None - preflist = evaluation._preferred_n_method.copy() - while preflist: - pref_method = preflist.pop() - if pref_method in ("numpy", "mpmath", "sympy"): - preference = pref_method - break - if preference is None: if d <= MACHINE_DIGITS: preference = "numpy" @@ -131,10 +130,16 @@ def get_constant(self, precision, evaluation): preference = "mpmath" else: preference = "" + if preference == "numpy": - value = numpy_constant(self.numpy_name) if d == MACHINE_DIGITS: - return MachineReal(value) + try: + return NUMERICAL_CONSTANTS[self.symbol] + except KeyError: + value = MachineReal(numpy_constant(self.numpy_name)) + NUMERICAL_CONSTANTS[self.symbol] = value + return value + value = numpy_constant(self.numpy_name) if preference == "sympy": value = sympy_constant(self.sympy_name, d + 2) if preference == "mpmath": @@ -177,13 +182,16 @@ class _NumpyConstant(_Constant_Common): def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) if self.numpy_name is None: - self.numpy_name = strip_context(self.get_name()).lower() + self.numpy_name = strip_context(self.symbol.name).lower() self.mathics_to_numpy[self.__class__.__name__] = self.numpy_name + try: + value_float = numpy_constant(self.numpy_name) + except AttributeError: + value_float = self.to_numpy(self.symbol) + NUMERICAL_CONSTANTS[self.symbol] = MachineReal(value_float) def to_numpy(self, args): - if self.numpy_name is None or len(args) != 0: - return None - return self.get_constant() + return NUMERICAL_CONSTANTS[self.symbol] class _SympyConstant(_Constant_Common, SympyObject): @@ -608,7 +616,7 @@ class MaxMachineNumber(Predefined): summary_text = "largest normalized positive machine number" def evaluate(self, evaluation: Evaluation) -> MachineReal: - return MachineReal(MAX_MACHINE_NUMBER) + return NUMERICAL_CONSTANTS[self.symbol] class MinMachineNumber(Predefined): @@ -635,10 +643,11 @@ class MinMachineNumber(Predefined): summary_text = "smallest normalized positive machine number" def evaluate(self, evaluation: Evaluation) -> MachineReal: - return MachineReal(MIN_MACHINE_NUMBER) + return NUMERICAL_CONSTANTS[self.symbol] class Pi(_MPMathConstant, _SympyConstant): + # class Pi(_MPMathConstant, _NumpyConstant, _SympyConstant): """ :Pi, \u03c0: https://en.wikipedia.org/wiki/Pi ( @@ -740,3 +749,10 @@ class Underflow(Builtin): "Underflow[] * x_Real": "0.", } summary_text = "underflow in numeric evaluation" + + +# Constants that are not numpy constants, +for cls in (Catalan, Degree, Glaisher, GoldenRatio, Khinchin, Pi): + instance = cls(expression=False) + val = instance.get_constant() + NUMERICAL_CONSTANTS[instance.symbol] = MachineReal(val.value) diff --git a/mathics/core/atoms.py b/mathics/core/atoms.py index f2b5c91a1..eadc66ad0 100644 --- a/mathics/core/atoms.py +++ b/mathics/core/atoms.py @@ -13,6 +13,8 @@ from mathics.core.number import ( FP_MANTISA_BINARY_DIGITS, MACHINE_PRECISION_VALUE, + MAX_MACHINE_NUMBER, + MIN_MACHINE_NUMBER, dps, min_prec, prec, @@ -946,6 +948,14 @@ def is_zero(self) -> bool: MATHICS3_COMPLEX_I = Complex(Integer0, Integer1) MATHICS3_COMPLEX_I_NEG = Complex(Integer0, IntegerM1) +# Numerical constants +# These constants are populated by the `Predefined` +# classes. See `mathics.builtin.numbers.constants` +NUMERICAL_CONSTANTS = { + Symbol("System`$MaxMachineNumber"): MachineReal(MAX_MACHINE_NUMBER), + Symbol("System`$MinMachineNumber"): MachineReal(MIN_MACHINE_NUMBER), +} + class String(Atom, BoxElementMixin): value: str diff --git a/mathics/core/number.py b/mathics/core/number.py index f30de3b92..01673f0ea 100644 --- a/mathics/core/number.py +++ b/mathics/core/number.py @@ -68,7 +68,9 @@ def _get_float_inf(value, evaluation) -> Optional[float]: return value.round_to_float(evaluation) -def get_precision(value, evaluation, show_messages=True) -> Optional[float]: +def get_precision( + value: BaseElement, evaluation, show_messages: bool = True +) -> Optional[float]: """ Returns the ``float`` in the interval [``$MinPrecision``, ``$MaxPrecision``] closest to ``value``. If ``value`` does not belongs to that interval, and ``show_messages`` is True, a Message warning is shown. diff --git a/mathics/eval/arithmetic.py b/mathics/eval/arithmetic.py index 399d75b86..c1cf6dd5d 100644 --- a/mathics/eval/arithmetic.py +++ b/mathics/eval/arithmetic.py @@ -14,12 +14,14 @@ import sympy from mathics.core.atoms import ( + NUMERICAL_CONSTANTS, Complex, Integer, Integer0, Integer1, Integer2, IntegerM1, + MachineReal, Number, Rational, RationalOneHalf, @@ -29,29 +31,22 @@ from mathics.core.convert.sympy import from_sympy from mathics.core.element import BaseElement, ElementsProperties from mathics.core.expression import Expression -from mathics.core.number import FP_MANTISA_BINARY_DIGITS, SpecialValueError, min_prec +from mathics.core.number import ( + FP_MANTISA_BINARY_DIGITS, + MAX_MACHINE_NUMBER, + MIN_MACHINE_NUMBER, + SpecialValueError, + min_prec, +) from mathics.core.rules import Rule from mathics.core.symbols import Atom, Symbol, SymbolPlus, SymbolPower, SymbolTimes from mathics.core.systemsymbols import ( SymbolComplexInfinity, - SymbolE, - SymbolEulerGamma, SymbolI, SymbolIndeterminate, - SymbolKhinchin, SymbolLog, - SymbolPi, ) -# FIXME: replace by numpy constants: -NUMERICAL_CONSTANTS = { - SymbolE: 2.718281828, - SymbolEulerGamma: 0.5772156649, - SymbolKhinchin: 2.685452001, - SymbolPi: 3.141592654, -} - - RationalMOneHalf = Rational(-1, 2) RealOne = Real(1.0) @@ -406,7 +401,7 @@ def test_arithmetic_expr(expr: BaseElement, only_real: bool = True) -> bool: return True if expr in NUMERICAL_CONSTANTS: return True - if isinstance(expr, Complex): + if isinstance(expr, Complex) or expr is SymbolI: return not only_real if isinstance(expr, Symbol): return False @@ -431,9 +426,6 @@ def test_arithmetic_expr(expr: BaseElement, only_real: bool = True) -> bool: 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 @@ -525,7 +517,7 @@ def test_positive_arithmetic_expr(expr: BaseElement) -> bool: 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) + return test_arithmetic_expr(exponent, only_real=True) if head is SymbolLog: if len(elements) > 2: return False @@ -573,9 +565,10 @@ def is_numeric_zero(z: Number): 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 test_zero_arithmetic_expr(base, numeric): + return test_nonnegative_arithmetic_expr(exp) + if base.has_form("DirectedInfinity", None): + return test_positive_arithmetic_expr(exp) if expr.has_form("Plus", None): result = eval_add_numbers(*expr.elements) if numeric: @@ -600,5 +593,6 @@ def to_inexact_value(expr: BaseElement) -> BaseElement: if isinstance(expr, Expression): for const, value in NUMERICAL_CONSTANTS.items(): - expr, success = expr.do_apply_rule(Rule(const, Real(value))) + expr, success = expr.do_apply_rule(Rule(const, value)) + return eval_multiply_numbers(RealOne, expr) From 8e7e7743e2f01c02313842e3b898330383c54f84 Mon Sep 17 00:00:00 2001 From: mmatera Date: Fri, 26 May 2023 11:13:29 -0300 Subject: [PATCH 03/14] adding tests. fix typo improving tests more fixes more tests --- mathics/eval/arithmetic.py | 18 +++- mathics/session.py | 6 ++ .../arithmetic/test_lowlevel_properties.py | 86 +++++++++++++++++++ 3 files changed, 107 insertions(+), 3 deletions(-) create mode 100644 test/builtin/arithmetic/test_lowlevel_properties.py diff --git a/mathics/eval/arithmetic.py b/mathics/eval/arithmetic.py index c1cf6dd5d..33d76c6e8 100644 --- a/mathics/eval/arithmetic.py +++ b/mathics/eval/arithmetic.py @@ -446,6 +446,9 @@ def test_nonnegative_arithmetic_expr(expr: BaseElement) -> bool: Check if the expression is an arithmetic expression representing a nonnegative number """ + if not test_arithmetic_expr(expr): + return False + if test_zero_arithmetic_expr(expr) or test_positive_arithmetic_expr(expr): return True @@ -455,6 +458,9 @@ def test_nonpositive_arithetic_expr(expr: BaseElement) -> bool: Check if the expression is an arithmetic expression representing a nonnegative number """ + if not test_arithmetic_expr(expr): + return False + if test_zero_arithmetic_expr(expr) or test_negative_arithmetic_expr(expr): return True return False @@ -473,6 +479,7 @@ def test_positive_arithmetic_expr(expr: BaseElement) -> bool: return False head, elements = expr.get_head(), expr.elements + print("check specifics", expr) if head is SymbolPlus: positive_nonpositive_terms = {True: [], False: []} for term in elements: @@ -517,7 +524,9 @@ def test_positive_arithmetic_expr(expr: BaseElement) -> bool: 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, only_real=True) + return test_arithmetic_expr(expr.elements[0], only_real=True) + if expr.has_form("Sqrt", 1): + return test_positive_arithmetic_expr(expr.elements[0]) if head is SymbolLog: if len(elements) > 2: return False @@ -526,8 +535,11 @@ def test_positive_arithmetic_expr(expr: BaseElement) -> bool: return False arg = elements[-1] return test_positive_arithmetic_expr(eval_add_numbers(arg, IntegerM1)) - if head.has_form("Abs", 1): - return True + if expr.has_form("Abs", 1): + arg = elements[0] + return test_arithmetic_expr( + arg, only_real=False + ) and not test_zero_arithmetic_expr(arg) if head.has_form("DirectedInfinity", 1): return test_positive_arithmetic_expr(elements[0]) diff --git a/mathics/session.py b/mathics/session.py index 410d9c2b1..043106b34 100644 --- a/mathics/session.py +++ b/mathics/session.py @@ -100,3 +100,9 @@ def format_result(self, str_expression=None, timeout=None, form=None): if form is None: form = self.form return res.do_format(self.evaluation, form) + + def parse(self, str_expression): + """ + Just parse the expression + """ + return parse(self.definitions, MathicsSingleLineFeeder(str_expression)) diff --git a/test/builtin/arithmetic/test_lowlevel_properties.py b/test/builtin/arithmetic/test_lowlevel_properties.py new file mode 100644 index 000000000..10cebfd06 --- /dev/null +++ b/test/builtin/arithmetic/test_lowlevel_properties.py @@ -0,0 +1,86 @@ +# -*- coding: utf-8 -*- +""" +Unit tests for mathics.eval.arithmetic low level positivity tests +""" +from test.helper import session + +import pytest + +from mathics.eval.arithmetic import ( + test_arithmetic_expr as check_arithmetic, + test_positive_arithmetic_expr as check_positive, + test_zero_arithmetic_expr as check_zero, +) + + +@pytest.mark.parametrize( + ("str_expr", "expected", "msg"), + [ + ("I", False, None), + ("0", False, None), + ("1", True, None), + ("Pi", True, None), + ("a", False, None), + ("-Pi", False, None), + ("(-1)^2", True, None), + ("(-1)^3", False, None), + ("Sqrt[2]", True, None), + ("Sqrt[-2]", False, None), + ("(-2)^(1/2)", False, None), + ("(2)^(1/2)", True, None), + ("Exp[a]", False, None), + ("Exp[2.3]", True, None), + ("Log[1/2]", False, None), + ("Exp[I]", False, None), + ("Log[3]", True, None), + ("Log[I]", False, None), + ("Abs[a]", False, None), + ("Abs[0]", False, None), + ("Abs[1+3 I]", True, None), + ("Sin[Pi]", False, None), + ], +) +def test_positivity(str_expr, expected, msg): + expr = session.parse(str_expr) + if msg: + assert check_positive(expr) == expected, msg + else: + assert check_positive(expr) == expected + + +@pytest.mark.parametrize( + ("str_expr", "expected", "msg"), + [ + ("I", False, None), + ("0", True, None), + ("1", False, None), + ("Pi", False, None), + ("a", False, None), + ("a-a", True, None), + ("3-3.", True, None), + ("2-Sqrt[4]", True, None), + ("-Pi", False, None), + ("(-1)^2", False, None), + ("(-1)^3", False, None), + ("Sqrt[2]", False, None), + ("Sqrt[-2]", False, None), + ("(-2)^(1/2)", False, None), + ("(2)^(1/2)", False, None), + ("Exp[a]", False, None), + ("Exp[2.3]", False, None), + ("Log[1/2]", False, None), + ("Exp[I]", False, None), + ("Log[3]", False, None), + ("Log[I]", False, None), + ("Abs[a]", False, None), + ("Abs[0]", False, None), + ("Abs[1+3 I]", False, None), + # ("Sin[Pi]", False, None), + ], +) +def test_zero(str_expr, expected, msg): + expr = session.parse(str_expr) + if msg: + assert check_zero(expr) == expected, msg + else: + assert check_zero(expr) == expected From 773aee0b09b423b0e20e97f142527ff8f1692b93 Mon Sep 17 00:00:00 2001 From: mmatera Date: Fri, 26 May 2023 14:07:37 -0300 Subject: [PATCH 04/14] Pi is a Numpy constant --- mathics/builtin/numbers/constants.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/mathics/builtin/numbers/constants.py b/mathics/builtin/numbers/constants.py index a29aefe81..73779ff6a 100644 --- a/mathics/builtin/numbers/constants.py +++ b/mathics/builtin/numbers/constants.py @@ -646,8 +646,7 @@ def evaluate(self, evaluation: Evaluation) -> MachineReal: return NUMERICAL_CONSTANTS[self.symbol] -class Pi(_MPMathConstant, _SympyConstant): - # class Pi(_MPMathConstant, _NumpyConstant, _SympyConstant): +class Pi(_MPMathConstant, _NumpyConstant, _SympyConstant): """ :Pi, \u03c0: https://en.wikipedia.org/wiki/Pi ( @@ -752,7 +751,7 @@ class Underflow(Builtin): # Constants that are not numpy constants, -for cls in (Catalan, Degree, Glaisher, GoldenRatio, Khinchin, Pi): +for cls in (Catalan, Degree, Glaisher, GoldenRatio, Khinchin): instance = cls(expression=False) val = instance.get_constant() NUMERICAL_CONSTANTS[instance.symbol] = MachineReal(val.value) From e8b90f82e482ce51d8e74597c1ebea7fffeb63a5 Mon Sep 17 00:00:00 2001 From: mmatera Date: Fri, 26 May 2023 15:11:22 -0300 Subject: [PATCH 05/14] arithmetic refactor --- mathics/builtin/arithfns/basic.py | 82 +- mathics/builtin/arithmetic.py | 22 +- mathics/builtin/base.py | 6 +- mathics/builtin/numbers/constants.py | 81 +- .../numerical_properties.py | 5 + mathics/core/atoms.py | 10 + mathics/core/number.py | 4 +- mathics/core/systemsymbols.py | 6 + mathics/eval/arithmetic.py | 720 +++++++++++++++++- mathics/eval/parts.py | 8 +- mathics/eval/testing_expressions.py | 1 - mathics/session.py | 6 + test/builtin/arithmetic/test_abs.py | 3 +- test/builtin/arithmetic/test_basic.py | 12 +- .../arithmetic/test_lowlevel_properties.py | 86 +++ test/format/test_format.py | 43 +- 16 files changed, 983 insertions(+), 112 deletions(-) create mode 100644 test/builtin/arithmetic/test_lowlevel_properties.py diff --git a/mathics/builtin/arithfns/basic.py b/mathics/builtin/arithfns/basic.py index 6fb5a8873..63752d4ab 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): diff --git a/mathics/builtin/arithmetic.py b/mathics/builtin/arithmetic.py index 5a85bbe7e..152e16d9f 100644 --- a/mathics/builtin/arithmetic.py +++ b/mathics/builtin/arithmetic.py @@ -629,10 +629,16 @@ class DirectedInfinity(SympyFunction): rules = { "DirectedInfinity[args___] ^ -1": "0", # Special arguments: - "DirectedInfinity[DirectedInfinity[args___]]": "DirectedInfinity[args]", "DirectedInfinity[Indeterminate]": "Indeterminate", "DirectedInfinity[Alternatives[0, 0.]]": "DirectedInfinity[]", # Plus + "DirectedInfinity[DirectedInfinity[args___]]": "DirectedInfinity[args]", + # "DirectedInfinity[a_?NumericQ] /; N[Abs[a]] != 1": "DirectedInfinity[a / Abs[a]]", + # "DirectedInfinity[a_] * DirectedInfinity[b_]": "DirectedInfinity[a*b]", + # "DirectedInfinity[] * DirectedInfinity[args___]": "DirectedInfinity[]", + # Rules already implemented in Times.eval + # "z_?NumberQ * DirectedInfinity[]": "DirectedInfinity[]", + # "z_?NumberQ * DirectedInfinity[a_]": "DirectedInfinity[z * a]", "DirectedInfinity[a_] + DirectedInfinity[b_] /; b == -a": ( "Message[Infinity::indet," " Unevaluated[DirectedInfinity[a] + DirectedInfinity[b]]];" @@ -644,23 +650,27 @@ class DirectedInfinity(SympyFunction): "Indeterminate" ), "DirectedInfinity[args___] + _?NumberQ": "DirectedInfinity[args]", - # Times. See if can be reinstalled in eval_Times + "DirectedInfinity[Alternatives[0, 0.]]": "DirectedInfinity[]", + "DirectedInfinity[0.]": ( + "Message[Infinity::indet," + " Unevaluated[DirectedInfinity[0.]]];" + "Indeterminate" + ), "Alternatives[0, 0.] DirectedInfinity[z___]": ( "Message[Infinity::indet," " Unevaluated[0 DirectedInfinity[z]]];" "Indeterminate" ), - "a_?NumericQ * DirectedInfinity[b_]": "DirectedInfinity[a * b]", "a_ DirectedInfinity[]": "DirectedInfinity[]", "DirectedInfinity[a_] * DirectedInfinity[b_]": "DirectedInfinity[a * b]", + "a_?NumericQ * DirectedInfinity[b_]": "DirectedInfinity[a * b]", } formats = { "DirectedInfinity[1]": "HoldForm[Infinity]", - "DirectedInfinity[-1]": "HoldForm[-Infinity]", + "DirectedInfinity[-1]": "PrecedenceForm[-HoldForm[Infinity], 390]", "DirectedInfinity[]": "HoldForm[ComplexInfinity]", - "DirectedInfinity[DirectedInfinity[z_]]": "DirectedInfinity[z]", - "DirectedInfinity[z_?NumericQ]": "HoldForm[z Infinity]", + "DirectedInfinity[z_]": "PrecedenceForm[z HoldForm[Infinity], 390]", } def eval_complex_infinity(self, evaluation: Evaluation): diff --git a/mathics/builtin/base.py b/mathics/builtin/base.py index 06a4337b0..8eaccab79 100644 --- a/mathics/builtin/base.py +++ b/mathics/builtin/base.py @@ -799,10 +799,14 @@ def get_operator_display(self) -> Optional[str]: class Predefined(Builtin): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self.symbol = Symbol(self.get_name()) + def get_functions(self, prefix="eval", is_pymodule=False) -> List[Callable]: functions = list(super().get_functions(prefix)) if prefix == "eval" and hasattr(self, "evaluate"): - functions.append((Symbol(self.get_name()), self.evaluate)) + functions.append((self.symbol, self.evaluate)) return functions diff --git a/mathics/builtin/numbers/constants.py b/mathics/builtin/numbers/constants.py index 787826869..73779ff6a 100644 --- a/mathics/builtin/numbers/constants.py +++ b/mathics/builtin/numbers/constants.py @@ -9,25 +9,19 @@ # This tells documentation how to sort this module sort_order = "mathics.builtin.mathematical-constants" - import math +from typing import Optional import mpmath import numpy import sympy from mathics.builtin.base import Builtin, Predefined, SympyObject -from mathics.core.atoms import MachineReal, PrecisionReal +from mathics.core.atoms import NUMERICAL_CONSTANTS, MachineReal, PrecisionReal from mathics.core.attributes import A_CONSTANT, A_PROTECTED, A_READ_PROTECTED +from mathics.core.element import BaseElement from mathics.core.evaluation import Evaluation -from mathics.core.number import ( - MACHINE_DIGITS, - MAX_MACHINE_NUMBER, - MIN_MACHINE_NUMBER, - PrecisionValueError, - get_precision, - prec, -) +from mathics.core.number import MACHINE_DIGITS, PrecisionValueError, get_precision, prec from mathics.core.symbols import Atom, Symbol, strip_context from mathics.core.systemsymbols import SymbolIndeterminate @@ -89,28 +83,33 @@ def eval_N(self, precision, evaluation): def is_constant(self) -> bool: return True - def get_constant(self, precision, evaluation): + def get_constant( + self, + precision: Optional[BaseElement] = None, + evaluation: Optional[Evaluation] = None, + ): # first, determine the precision d = None - if precision: - try: - d = get_precision(precision, evaluation) - except PrecisionValueError: - pass + preference = None + if evaluation: + if precision: + try: + d = get_precision(precision, evaluation) + except PrecisionValueError: + pass + + preflist = evaluation._preferred_n_method.copy() + while preflist: + pref_method = preflist.pop() + if pref_method in ("numpy", "mpmath", "sympy"): + preference = pref_method + break if d is None: d = MACHINE_DIGITS # If preference not especified, determine it # from the precision. - preference = None - preflist = evaluation._preferred_n_method.copy() - while preflist: - pref_method = preflist.pop() - if pref_method in ("numpy", "mpmath", "sympy"): - preference = pref_method - break - if preference is None: if d <= MACHINE_DIGITS: preference = "numpy" @@ -131,10 +130,16 @@ def get_constant(self, precision, evaluation): preference = "mpmath" else: preference = "" + if preference == "numpy": - value = numpy_constant(self.numpy_name) if d == MACHINE_DIGITS: - return MachineReal(value) + try: + return NUMERICAL_CONSTANTS[self.symbol] + except KeyError: + value = MachineReal(numpy_constant(self.numpy_name)) + NUMERICAL_CONSTANTS[self.symbol] = value + return value + value = numpy_constant(self.numpy_name) if preference == "sympy": value = sympy_constant(self.sympy_name, d + 2) if preference == "mpmath": @@ -177,13 +182,16 @@ class _NumpyConstant(_Constant_Common): def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) if self.numpy_name is None: - self.numpy_name = strip_context(self.get_name()).lower() + self.numpy_name = strip_context(self.symbol.name).lower() self.mathics_to_numpy[self.__class__.__name__] = self.numpy_name + try: + value_float = numpy_constant(self.numpy_name) + except AttributeError: + value_float = self.to_numpy(self.symbol) + NUMERICAL_CONSTANTS[self.symbol] = MachineReal(value_float) def to_numpy(self, args): - if self.numpy_name is None or len(args) != 0: - return None - return self.get_constant() + return NUMERICAL_CONSTANTS[self.symbol] class _SympyConstant(_Constant_Common, SympyObject): @@ -608,7 +616,7 @@ class MaxMachineNumber(Predefined): summary_text = "largest normalized positive machine number" def evaluate(self, evaluation: Evaluation) -> MachineReal: - return MachineReal(MAX_MACHINE_NUMBER) + return NUMERICAL_CONSTANTS[self.symbol] class MinMachineNumber(Predefined): @@ -635,10 +643,10 @@ class MinMachineNumber(Predefined): summary_text = "smallest normalized positive machine number" def evaluate(self, evaluation: Evaluation) -> MachineReal: - return MachineReal(MIN_MACHINE_NUMBER) + return NUMERICAL_CONSTANTS[self.symbol] -class Pi(_MPMathConstant, _SympyConstant): +class Pi(_MPMathConstant, _NumpyConstant, _SympyConstant): """ :Pi, \u03c0: https://en.wikipedia.org/wiki/Pi ( @@ -740,3 +748,10 @@ class Underflow(Builtin): "Underflow[] * x_Real": "0.", } summary_text = "underflow in numeric evaluation" + + +# Constants that are not numpy constants, +for cls in (Catalan, Degree, Glaisher, GoldenRatio, Khinchin): + instance = cls(expression=False) + val = instance.get_constant() + NUMERICAL_CONSTANTS[instance.symbol] = MachineReal(val.value) diff --git a/mathics/builtin/testing_expressions/numerical_properties.py b/mathics/builtin/testing_expressions/numerical_properties.py index 8de45e8b1..5f5822f9e 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 BooleanType, 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 @@ -460,6 +461,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/atoms.py b/mathics/core/atoms.py index f2b5c91a1..eadc66ad0 100644 --- a/mathics/core/atoms.py +++ b/mathics/core/atoms.py @@ -13,6 +13,8 @@ from mathics.core.number import ( FP_MANTISA_BINARY_DIGITS, MACHINE_PRECISION_VALUE, + MAX_MACHINE_NUMBER, + MIN_MACHINE_NUMBER, dps, min_prec, prec, @@ -946,6 +948,14 @@ def is_zero(self) -> bool: MATHICS3_COMPLEX_I = Complex(Integer0, Integer1) MATHICS3_COMPLEX_I_NEG = Complex(Integer0, IntegerM1) +# Numerical constants +# These constants are populated by the `Predefined` +# classes. See `mathics.builtin.numbers.constants` +NUMERICAL_CONSTANTS = { + Symbol("System`$MaxMachineNumber"): MachineReal(MAX_MACHINE_NUMBER), + Symbol("System`$MinMachineNumber"): MachineReal(MIN_MACHINE_NUMBER), +} + class String(Atom, BoxElementMixin): value: str diff --git a/mathics/core/number.py b/mathics/core/number.py index f30de3b92..01673f0ea 100644 --- a/mathics/core/number.py +++ b/mathics/core/number.py @@ -68,7 +68,9 @@ def _get_float_inf(value, evaluation) -> Optional[float]: return value.round_to_float(evaluation) -def get_precision(value, evaluation, show_messages=True) -> Optional[float]: +def get_precision( + value: BaseElement, evaluation, show_messages: bool = True +) -> Optional[float]: """ Returns the ``float`` in the interval [``$MinPrecision``, ``$MaxPrecision``] closest to ``value``. If ``value`` does not belongs to that interval, and ``show_messages`` is True, a Message warning is shown. diff --git a/mathics/core/systemsymbols.py b/mathics/core/systemsymbols.py index 52659e775..f463e5f91 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") @@ -84,6 +85,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") @@ -109,6 +111,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") @@ -126,6 +129,7 @@ SymbolLess = Symbol("System`Less") SymbolLessEqual = Symbol("System`LessEqual") SymbolKey = Symbol("System`Key") +SymbolKhinchin = Symbol("System`Khinchin") SymbolLetterCharacter = Symbol("System`LetterCharacter") SymbolLine = Symbol("System`Line") SymbolLog = Symbol("System`Log") @@ -179,6 +183,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") @@ -242,6 +247,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 index 0b9564067..2e7dad161 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. @@ -14,6 +16,7 @@ import sympy from mathics.core.atoms import ( + NUMERICAL_CONSTANTS, Complex, Integer, Integer0, @@ -30,10 +33,21 @@ from mathics.core.element import BaseElement, ElementsProperties from mathics.core.expression import Expression from mathics.core.number import FP_MANTISA_BINARY_DIGITS, SpecialValueError, min_prec -from mathics.core.symbols import Symbol, SymbolPlus, SymbolPower, SymbolTimes -from mathics.core.systemsymbols import SymbolComplexInfinity, SymbolIndeterminate +from mathics.core.rules import Rule +from mathics.core.symbols import Atom, Symbol, SymbolPlus, SymbolPower, SymbolTimes +from mathics.core.systemsymbols import ( + SymbolAbs, + SymbolComplexInfinity, + SymbolExp, + SymbolI, + SymbolIndeterminate, + SymbolLog, + SymbolSign, +) RationalMOneHalf = Rational(-1, 2) +RealM0p5 = Real(-0.5) +RealOne = Real(1.0) # This cache might not be used that much. @@ -70,31 +84,110 @@ def eval_Abs(expr: BaseElement) -> Optional[BaseElement]: """ if expr is a number, return the absolute value. """ - if isinstance(expr, (Integer, Rational, Real)): - if expr.value >= 0: - return expr - return eval_multiply_numbers(*[IntegerM1, expr]) - if isinstance(expr, Complex): - re, im = expr.real, expr.imag - sqabs = eval_add_numbers( - eval_multiply_numbers(re, re), eval_multiply_numbers(im, im) - ) - return Expression(SymbolPower, sqabs, RationalOneHalf) + + 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_Sign(expr: BaseElement) -> Optional[BaseElement]: """ if expr is a number, return its sign. """ - if isinstance(expr, (Integer, Rational, Real)): - if expr.value > 0: + 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 - elif expr.value == 0: - return Integer0 - else: - return IntegerM1 - + 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 if isinstance(expr, Complex): re, im = expr.real, expr.imag sqabs = eval_add_numbers(eval_Times(re, re), eval_Times(im, im)) @@ -106,6 +199,28 @@ def eval_Sign(expr: BaseElement) -> Optional[BaseElement]: 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)) + + def eval_mpmath_function( mpmath_function: Callable, *args: Number, prec: Optional[int] = None ) -> Optional[Number]: @@ -133,6 +248,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) @@ -200,6 +340,157 @@ 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]: + """ + 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, prec) + def eval_Times(*items: BaseElement) -> BaseElement: elements = [] @@ -287,9 +578,112 @@ 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 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_add_numbers( - *numbers: Number, -) -> BaseElement: + *numbers: List[Number], +) -> Number: """ Add the elements in ``numbers``. """ @@ -315,7 +709,37 @@ def eval_add_numbers( return from_sympy(sum(item.to_sympy() for item in numbers)) -def eval_multiply_numbers(*numbers: Number) -> BaseElement: +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_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``. """ @@ -340,6 +764,37 @@ def eval_multiply_numbers(*numbers: Number) -> BaseElement: 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 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 segregate_numbers( *elements: BaseElement, ) -> Tuple[List[Number], List[BaseElement]]: @@ -374,3 +829,222 @@ def segregate_numbers_from_sorted_list( if not isinstance(element, Number): return list(elements[:pos]), list(elements[pos:]) return list(elements), [] + + +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) or expr is SymbolI: + return not only_real + if isinstance(expr, Symbol): + return False + + head, elements = expr.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) + 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 not test_arithmetic_expr(expr): + return False + + 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 not test_arithmetic_expr(expr): + return False + + 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(expr.elements[0], only_real=True) + if expr.has_form("Sqrt", 1): + return test_positive_arithmetic_expr(expr.elements[0]) + 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 expr.has_form("Abs", 1): + arg = elements[0] + return test_arithmetic_expr( + arg, only_real=False + ) and not test_zero_arithmetic_expr(arg) + 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 + + if test_zero_arithmetic_expr(base, numeric): + return test_nonnegative_arithmetic_expr(exp) + if base.has_form("DirectedInfinity", None): + return test_positive_arithmetic_expr(exp) + 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, value)) + + return eval_multiply_numbers(RealOne, expr) diff --git a/mathics/eval/parts.py b/mathics/eval/parts.py index 3ecf356d1..61a1adf33 100644 --- a/mathics/eval/parts.py +++ b/mathics/eval/parts.py @@ -6,7 +6,7 @@ from typing import List -from mathics.core.atoms import Integer, Integer1 +from mathics.core.atoms import Integer from mathics.core.convert.expression import make_expression from mathics.core.element import BaseElement, BoxElementMixin from mathics.core.exceptions import ( @@ -20,11 +20,7 @@ from mathics.core.list import ListExpression from mathics.core.subexpression import SubExpression from mathics.core.symbols import Atom, Symbol, SymbolList -from mathics.core.systemsymbols import ( - SymbolDirectedInfinity, - SymbolInfinity, - SymbolNothing, -) +from mathics.core.systemsymbols import SymbolInfinity, SymbolNothing from mathics.eval.patterns import Matcher diff --git a/mathics/eval/testing_expressions.py b/mathics/eval/testing_expressions.py index c15471f48..cab4a7d4b 100644 --- a/mathics/eval/testing_expressions.py +++ b/mathics/eval/testing_expressions.py @@ -13,7 +13,6 @@ def cmp(a, b) -> int: def do_cmp(x1, x2) -> Optional[int]: - # don't attempt to compare complex numbers for x in (x1, x2): # TODO: Send message General::nord diff --git a/mathics/session.py b/mathics/session.py index 410d9c2b1..043106b34 100644 --- a/mathics/session.py +++ b/mathics/session.py @@ -100,3 +100,9 @@ def format_result(self, str_expression=None, timeout=None, form=None): if form is None: form = self.form return res.do_format(self.evaluation, form) + + def parse(self, str_expression): + """ + Just parse the expression + """ + return parse(self.definitions, MathicsSingleLineFeeder(str_expression)) diff --git a/test/builtin/arithmetic/test_abs.py b/test/builtin/arithmetic/test_abs.py index c523bf1d3..11b3da92b 100644 --- a/test/builtin/arithmetic/test_abs.py +++ b/test/builtin/arithmetic/test_abs.py @@ -42,8 +42,7 @@ 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), - # FIXME: add rules to handle this kind of case - # ("Sign[I^(2 Pi)]", "I^(2 Pi)", None), + ("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 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/builtin/arithmetic/test_lowlevel_properties.py b/test/builtin/arithmetic/test_lowlevel_properties.py new file mode 100644 index 000000000..10cebfd06 --- /dev/null +++ b/test/builtin/arithmetic/test_lowlevel_properties.py @@ -0,0 +1,86 @@ +# -*- coding: utf-8 -*- +""" +Unit tests for mathics.eval.arithmetic low level positivity tests +""" +from test.helper import session + +import pytest + +from mathics.eval.arithmetic import ( + test_arithmetic_expr as check_arithmetic, + test_positive_arithmetic_expr as check_positive, + test_zero_arithmetic_expr as check_zero, +) + + +@pytest.mark.parametrize( + ("str_expr", "expected", "msg"), + [ + ("I", False, None), + ("0", False, None), + ("1", True, None), + ("Pi", True, None), + ("a", False, None), + ("-Pi", False, None), + ("(-1)^2", True, None), + ("(-1)^3", False, None), + ("Sqrt[2]", True, None), + ("Sqrt[-2]", False, None), + ("(-2)^(1/2)", False, None), + ("(2)^(1/2)", True, None), + ("Exp[a]", False, None), + ("Exp[2.3]", True, None), + ("Log[1/2]", False, None), + ("Exp[I]", False, None), + ("Log[3]", True, None), + ("Log[I]", False, None), + ("Abs[a]", False, None), + ("Abs[0]", False, None), + ("Abs[1+3 I]", True, None), + ("Sin[Pi]", False, None), + ], +) +def test_positivity(str_expr, expected, msg): + expr = session.parse(str_expr) + if msg: + assert check_positive(expr) == expected, msg + else: + assert check_positive(expr) == expected + + +@pytest.mark.parametrize( + ("str_expr", "expected", "msg"), + [ + ("I", False, None), + ("0", True, None), + ("1", False, None), + ("Pi", False, None), + ("a", False, None), + ("a-a", True, None), + ("3-3.", True, None), + ("2-Sqrt[4]", True, None), + ("-Pi", False, None), + ("(-1)^2", False, None), + ("(-1)^3", False, None), + ("Sqrt[2]", False, None), + ("Sqrt[-2]", False, None), + ("(-2)^(1/2)", False, None), + ("(2)^(1/2)", False, None), + ("Exp[a]", False, None), + ("Exp[2.3]", False, None), + ("Log[1/2]", False, None), + ("Exp[I]", False, None), + ("Log[3]", False, None), + ("Log[I]", False, None), + ("Abs[a]", False, None), + ("Abs[0]", False, None), + ("Abs[1+3 I]", False, None), + # ("Sin[Pi]", False, None), + ], +) +def test_zero(str_expr, expected, msg): + expr = session.parse(str_expr) + if msg: + assert check_zero(expr) == expected, msg + else: + assert check_zero(expr) == expected 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 From 757011f51a32235548ff691ea0bd51a01fb00e80 Mon Sep 17 00:00:00 2001 From: mmatera Date: Fri, 26 May 2023 15:24:20 -0300 Subject: [PATCH 06/14] flake8 mathics.eval.parts --- mathics/eval/parts.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/mathics/eval/parts.py b/mathics/eval/parts.py index 3ecf356d1..61a1adf33 100644 --- a/mathics/eval/parts.py +++ b/mathics/eval/parts.py @@ -6,7 +6,7 @@ from typing import List -from mathics.core.atoms import Integer, Integer1 +from mathics.core.atoms import Integer from mathics.core.convert.expression import make_expression from mathics.core.element import BaseElement, BoxElementMixin from mathics.core.exceptions import ( @@ -20,11 +20,7 @@ from mathics.core.list import ListExpression from mathics.core.subexpression import SubExpression from mathics.core.symbols import Atom, Symbol, SymbolList -from mathics.core.systemsymbols import ( - SymbolDirectedInfinity, - SymbolInfinity, - SymbolNothing, -) +from mathics.core.systemsymbols import SymbolInfinity, SymbolNothing from mathics.eval.patterns import Matcher From 7b6f6508f83d2bffc1e73668ca1d5c71e0d64a42 Mon Sep 17 00:00:00 2001 From: mmatera Date: Fri, 26 May 2023 14:07:37 -0300 Subject: [PATCH 07/14] Pi is a Numpy constant --- mathics/builtin/numbers/constants.py | 5 ++--- mathics/eval/arithmetic.py | 2 +- mathics/eval/parts.py | 8 ++------ 3 files changed, 5 insertions(+), 10 deletions(-) diff --git a/mathics/builtin/numbers/constants.py b/mathics/builtin/numbers/constants.py index a29aefe81..73779ff6a 100644 --- a/mathics/builtin/numbers/constants.py +++ b/mathics/builtin/numbers/constants.py @@ -646,8 +646,7 @@ def evaluate(self, evaluation: Evaluation) -> MachineReal: return NUMERICAL_CONSTANTS[self.symbol] -class Pi(_MPMathConstant, _SympyConstant): - # class Pi(_MPMathConstant, _NumpyConstant, _SympyConstant): +class Pi(_MPMathConstant, _NumpyConstant, _SympyConstant): """ :Pi, \u03c0: https://en.wikipedia.org/wiki/Pi ( @@ -752,7 +751,7 @@ class Underflow(Builtin): # Constants that are not numpy constants, -for cls in (Catalan, Degree, Glaisher, GoldenRatio, Khinchin, Pi): +for cls in (Catalan, Degree, Glaisher, GoldenRatio, Khinchin): instance = cls(expression=False) val = instance.get_constant() NUMERICAL_CONSTANTS[instance.symbol] = MachineReal(val.value) diff --git a/mathics/eval/arithmetic.py b/mathics/eval/arithmetic.py index 33d76c6e8..00976ac7b 100644 --- a/mathics/eval/arithmetic.py +++ b/mathics/eval/arithmetic.py @@ -479,7 +479,7 @@ def test_positive_arithmetic_expr(expr: BaseElement) -> bool: return False head, elements = expr.get_head(), expr.elements - print("check specifics", expr) + if head is SymbolPlus: positive_nonpositive_terms = {True: [], False: []} for term in elements: diff --git a/mathics/eval/parts.py b/mathics/eval/parts.py index 3ecf356d1..61a1adf33 100644 --- a/mathics/eval/parts.py +++ b/mathics/eval/parts.py @@ -6,7 +6,7 @@ from typing import List -from mathics.core.atoms import Integer, Integer1 +from mathics.core.atoms import Integer from mathics.core.convert.expression import make_expression from mathics.core.element import BaseElement, BoxElementMixin from mathics.core.exceptions import ( @@ -20,11 +20,7 @@ from mathics.core.list import ListExpression from mathics.core.subexpression import SubExpression from mathics.core.symbols import Atom, Symbol, SymbolList -from mathics.core.systemsymbols import ( - SymbolDirectedInfinity, - SymbolInfinity, - SymbolNothing, -) +from mathics.core.systemsymbols import SymbolInfinity, SymbolNothing from mathics.eval.patterns import Matcher From f44dc779f7b2bf3e67f40a5841a80bb3ef905aac Mon Sep 17 00:00:00 2001 From: mmatera Date: Fri, 26 May 2023 15:35:12 -0300 Subject: [PATCH 08/14] adding session.parse method, which is handy for low-level tests. removing unused imports. adding symbols to mathics.core.systemsymbols --- mathics/core/number.py | 4 +++- mathics/core/systemsymbols.py | 6 ++++++ mathics/eval/parts.py | 8 ++------ mathics/session.py | 6 ++++++ 4 files changed, 17 insertions(+), 7 deletions(-) diff --git a/mathics/core/number.py b/mathics/core/number.py index f30de3b92..01673f0ea 100644 --- a/mathics/core/number.py +++ b/mathics/core/number.py @@ -68,7 +68,9 @@ def _get_float_inf(value, evaluation) -> Optional[float]: return value.round_to_float(evaluation) -def get_precision(value, evaluation, show_messages=True) -> Optional[float]: +def get_precision( + value: BaseElement, evaluation, show_messages: bool = True +) -> Optional[float]: """ Returns the ``float`` in the interval [``$MinPrecision``, ``$MaxPrecision``] closest to ``value``. If ``value`` does not belongs to that interval, and ``show_messages`` is True, a Message warning is shown. diff --git a/mathics/core/systemsymbols.py b/mathics/core/systemsymbols.py index 52659e775..f463e5f91 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") @@ -84,6 +85,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") @@ -109,6 +111,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") @@ -126,6 +129,7 @@ SymbolLess = Symbol("System`Less") SymbolLessEqual = Symbol("System`LessEqual") SymbolKey = Symbol("System`Key") +SymbolKhinchin = Symbol("System`Khinchin") SymbolLetterCharacter = Symbol("System`LetterCharacter") SymbolLine = Symbol("System`Line") SymbolLog = Symbol("System`Log") @@ -179,6 +183,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") @@ -242,6 +247,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/parts.py b/mathics/eval/parts.py index 3ecf356d1..61a1adf33 100644 --- a/mathics/eval/parts.py +++ b/mathics/eval/parts.py @@ -6,7 +6,7 @@ from typing import List -from mathics.core.atoms import Integer, Integer1 +from mathics.core.atoms import Integer from mathics.core.convert.expression import make_expression from mathics.core.element import BaseElement, BoxElementMixin from mathics.core.exceptions import ( @@ -20,11 +20,7 @@ from mathics.core.list import ListExpression from mathics.core.subexpression import SubExpression from mathics.core.symbols import Atom, Symbol, SymbolList -from mathics.core.systemsymbols import ( - SymbolDirectedInfinity, - SymbolInfinity, - SymbolNothing, -) +from mathics.core.systemsymbols import SymbolInfinity, SymbolNothing from mathics.eval.patterns import Matcher diff --git a/mathics/session.py b/mathics/session.py index 410d9c2b1..043106b34 100644 --- a/mathics/session.py +++ b/mathics/session.py @@ -100,3 +100,9 @@ def format_result(self, str_expression=None, timeout=None, form=None): if form is None: form = self.form return res.do_format(self.evaluation, form) + + def parse(self, str_expression): + """ + Just parse the expression + """ + return parse(self.definitions, MathicsSingleLineFeeder(str_expression)) From 2481118bc4902388ac6213fc86e453f0d42a9d73 Mon Sep 17 00:00:00 2001 From: mmatera Date: Fri, 26 May 2023 14:07:37 -0300 Subject: [PATCH 09/14] Pi is a Numpy constant --- mathics/builtin/arithmetic.py | 22 ++++++---------------- mathics/builtin/numbers/constants.py | 5 ++--- mathics/eval/arithmetic.py | 2 +- mathics/eval/parts.py | 8 ++------ 4 files changed, 11 insertions(+), 26 deletions(-) diff --git a/mathics/builtin/arithmetic.py b/mathics/builtin/arithmetic.py index 152e16d9f..5a85bbe7e 100644 --- a/mathics/builtin/arithmetic.py +++ b/mathics/builtin/arithmetic.py @@ -629,16 +629,10 @@ class DirectedInfinity(SympyFunction): rules = { "DirectedInfinity[args___] ^ -1": "0", # Special arguments: + "DirectedInfinity[DirectedInfinity[args___]]": "DirectedInfinity[args]", "DirectedInfinity[Indeterminate]": "Indeterminate", "DirectedInfinity[Alternatives[0, 0.]]": "DirectedInfinity[]", # Plus - "DirectedInfinity[DirectedInfinity[args___]]": "DirectedInfinity[args]", - # "DirectedInfinity[a_?NumericQ] /; N[Abs[a]] != 1": "DirectedInfinity[a / Abs[a]]", - # "DirectedInfinity[a_] * DirectedInfinity[b_]": "DirectedInfinity[a*b]", - # "DirectedInfinity[] * DirectedInfinity[args___]": "DirectedInfinity[]", - # Rules already implemented in Times.eval - # "z_?NumberQ * DirectedInfinity[]": "DirectedInfinity[]", - # "z_?NumberQ * DirectedInfinity[a_]": "DirectedInfinity[z * a]", "DirectedInfinity[a_] + DirectedInfinity[b_] /; b == -a": ( "Message[Infinity::indet," " Unevaluated[DirectedInfinity[a] + DirectedInfinity[b]]];" @@ -650,27 +644,23 @@ class DirectedInfinity(SympyFunction): "Indeterminate" ), "DirectedInfinity[args___] + _?NumberQ": "DirectedInfinity[args]", - "DirectedInfinity[Alternatives[0, 0.]]": "DirectedInfinity[]", - "DirectedInfinity[0.]": ( - "Message[Infinity::indet," - " Unevaluated[DirectedInfinity[0.]]];" - "Indeterminate" - ), + # Times. See if can be reinstalled in eval_Times "Alternatives[0, 0.] DirectedInfinity[z___]": ( "Message[Infinity::indet," " Unevaluated[0 DirectedInfinity[z]]];" "Indeterminate" ), + "a_?NumericQ * DirectedInfinity[b_]": "DirectedInfinity[a * b]", "a_ DirectedInfinity[]": "DirectedInfinity[]", "DirectedInfinity[a_] * DirectedInfinity[b_]": "DirectedInfinity[a * b]", - "a_?NumericQ * DirectedInfinity[b_]": "DirectedInfinity[a * b]", } formats = { "DirectedInfinity[1]": "HoldForm[Infinity]", - "DirectedInfinity[-1]": "PrecedenceForm[-HoldForm[Infinity], 390]", + "DirectedInfinity[-1]": "HoldForm[-Infinity]", "DirectedInfinity[]": "HoldForm[ComplexInfinity]", - "DirectedInfinity[z_]": "PrecedenceForm[z HoldForm[Infinity], 390]", + "DirectedInfinity[DirectedInfinity[z_]]": "DirectedInfinity[z]", + "DirectedInfinity[z_?NumericQ]": "HoldForm[z Infinity]", } def eval_complex_infinity(self, evaluation: Evaluation): diff --git a/mathics/builtin/numbers/constants.py b/mathics/builtin/numbers/constants.py index a29aefe81..73779ff6a 100644 --- a/mathics/builtin/numbers/constants.py +++ b/mathics/builtin/numbers/constants.py @@ -646,8 +646,7 @@ def evaluate(self, evaluation: Evaluation) -> MachineReal: return NUMERICAL_CONSTANTS[self.symbol] -class Pi(_MPMathConstant, _SympyConstant): - # class Pi(_MPMathConstant, _NumpyConstant, _SympyConstant): +class Pi(_MPMathConstant, _NumpyConstant, _SympyConstant): """ :Pi, \u03c0: https://en.wikipedia.org/wiki/Pi ( @@ -752,7 +751,7 @@ class Underflow(Builtin): # Constants that are not numpy constants, -for cls in (Catalan, Degree, Glaisher, GoldenRatio, Khinchin, Pi): +for cls in (Catalan, Degree, Glaisher, GoldenRatio, Khinchin): instance = cls(expression=False) val = instance.get_constant() NUMERICAL_CONSTANTS[instance.symbol] = MachineReal(val.value) diff --git a/mathics/eval/arithmetic.py b/mathics/eval/arithmetic.py index 33d76c6e8..00976ac7b 100644 --- a/mathics/eval/arithmetic.py +++ b/mathics/eval/arithmetic.py @@ -479,7 +479,7 @@ def test_positive_arithmetic_expr(expr: BaseElement) -> bool: return False head, elements = expr.get_head(), expr.elements - print("check specifics", expr) + if head is SymbolPlus: positive_nonpositive_terms = {True: [], False: []} for term in elements: diff --git a/mathics/eval/parts.py b/mathics/eval/parts.py index 3ecf356d1..61a1adf33 100644 --- a/mathics/eval/parts.py +++ b/mathics/eval/parts.py @@ -6,7 +6,7 @@ from typing import List -from mathics.core.atoms import Integer, Integer1 +from mathics.core.atoms import Integer from mathics.core.convert.expression import make_expression from mathics.core.element import BaseElement, BoxElementMixin from mathics.core.exceptions import ( @@ -20,11 +20,7 @@ from mathics.core.list import ListExpression from mathics.core.subexpression import SubExpression from mathics.core.symbols import Atom, Symbol, SymbolList -from mathics.core.systemsymbols import ( - SymbolDirectedInfinity, - SymbolInfinity, - SymbolNothing, -) +from mathics.core.systemsymbols import SymbolInfinity, SymbolNothing from mathics.eval.patterns import Matcher From a0b53105ae75a1b913305038e4513fddc02daffe Mon Sep 17 00:00:00 2001 From: mmatera Date: Sat, 27 May 2023 13:16:09 -0300 Subject: [PATCH 10/14] flake8 --- mathics/eval/arithmetic.py | 1 - mathics/eval/testing_expressions.py | 1 + 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/mathics/eval/arithmetic.py b/mathics/eval/arithmetic.py index fb511271c..d196a20b1 100644 --- a/mathics/eval/arithmetic.py +++ b/mathics/eval/arithmetic.py @@ -23,7 +23,6 @@ Integer1, Integer2, IntegerM1, - MachineReal, Number, Rational, RationalOneHalf, diff --git a/mathics/eval/testing_expressions.py b/mathics/eval/testing_expressions.py index cab4a7d4b..c15471f48 100644 --- a/mathics/eval/testing_expressions.py +++ b/mathics/eval/testing_expressions.py @@ -13,6 +13,7 @@ def cmp(a, b) -> int: def do_cmp(x1, x2) -> Optional[int]: + # don't attempt to compare complex numbers for x in (x1, x2): # TODO: Send message General::nord From a2ada0b5598d387dd36637dc40b741b2c3b41d69 Mon Sep 17 00:00:00 2001 From: rocky Date: Sun, 16 Jul 2023 16:35:19 -0400 Subject: [PATCH 11/14] Correct typo in error message pattern --- .gitignore | 5 +++-- mathics/builtin/tensors.py | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/.gitignore b/.gitignore index ccd7736f5..2d93eb961 100644 --- a/.gitignore +++ b/.gitignore @@ -14,9 +14,10 @@ .settings .vscode /.cache +/.gdbinit /.python-version -/Mathics3.egg-info /Mathics.egg-info +/Mathics3.egg-info ChangeLog Documents/ Homepage/ @@ -33,9 +34,9 @@ mathics/doc/tex/logo-heptatom.pdf mathics/doc/tex/logo-text-nodrop.pdf mathics/doc/tex/mathics-*.asy mathics/doc/tex/mathics-*.dvi +mathics/doc/tex/mathics-*.dvi mathics/doc/tex/mathics-*.eps mathics/doc/tex/mathics-*.pdf -mathics/doc/tex/mathics-*.dvi mathics/doc/tex/mathics-*.tex mathics/doc/tex/mathics.aux mathics/doc/tex/mathics.dvi diff --git a/mathics/builtin/tensors.py b/mathics/builtin/tensors.py index 57e8198ab..c76f74486 100644 --- a/mathics/builtin/tensors.py +++ b/mathics/builtin/tensors.py @@ -216,7 +216,7 @@ class Inner(Builtin): messages = { "incom": ( "Length `1` of dimension `2` in `3` is incommensurate with " - "length `4` of dimension 1 in `5." + "length `4` of dimension 1 in `5`." ), } From 87316dcd5ef84d1843fbf71720dd7f4f372776f9 Mon Sep 17 00:00:00 2001 From: mmatera Date: Mon, 29 May 2023 10:17:14 -0300 Subject: [PATCH 12/14] RealAbs and RealSign --- CHANGES.rst | 4 +- SYMBOLS_MANIFEST.txt | 2 + mathics/builtin/arithmetic.py | 82 ++- .../equality_inequality.py | 2 +- mathics/core/systemsymbols.py | 4 +- mathics/eval/arithmetic.py | 664 +++++++++++++----- mathics/eval/testing_expressions.py | 4 +- .../arithmetic/test_lowlevel_properties.py | 6 +- 8 files changed, 597 insertions(+), 171 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 76238a823..6bfe2990d 100644 --- a/CHANGES.rst +++ b/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 ---- diff --git a/SYMBOLS_MANIFEST.txt b/SYMBOLS_MANIFEST.txt index 1bf18d1ea..d3ecf4204 100644 --- a/SYMBOLS_MANIFEST.txt +++ b/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 diff --git a/mathics/builtin/arithmetic.py b/mathics/builtin/arithmetic.py index d54f8f7bd..51980d633 100644 --- a/mathics/builtin/arithmetic.py +++ b/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 @@ -1207,6 +1213,44 @@ class Real_(Builtin): name = "Real" +class RealAbs(Builtin): + """ + :WMA link:https://reference.wolfram.com/language/ref/RealAbs.html + +
+
'RealAbs[$x$]' +
returns the absolute value of a real number $x$. +
+ '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,6 +1281,42 @@ def test(self, expr) -> bool: return isinstance(expr, (Integer, Rational, Real)) +class RealSign(Builtin): + """ + :WMA link:https://reference.wolfram.com/language/ref/RealAbs.html + +
+
'RealSign[$x$]' +
returns $-1$, $0$ or $1$ depending on whether $x$ is negative, + zero or positive. +
+ '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): """ :WMA link:https://reference.wolfram.com/language/ref/Sign.html diff --git a/mathics/builtin/testing_expressions/equality_inequality.py b/mathics/builtin/testing_expressions/equality_inequality.py index f95289f17..d5e0ce6e5 100644 --- a/mathics/builtin/testing_expressions/equality_inequality.py +++ b/mathics/builtin/testing_expressions/equality_inequality.py @@ -70,7 +70,7 @@ def numerify_args(items, evaluation) -> list: 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: diff --git a/mathics/core/systemsymbols.py b/mathics/core/systemsymbols.py index 31d945385..fa388a92d 100644 --- a/mathics/core/systemsymbols.py +++ b/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") diff --git a/mathics/eval/arithmetic.py b/mathics/eval/arithmetic.py index 2e893b014..abd4e51a8 100644 --- a/mathics/eval/arithmetic.py +++ b/mathics/eval/arithmetic.py @@ -31,16 +31,21 @@ from mathics.core.element import BaseElement, ElementsProperties from mathics.core.expression import Expression from mathics.core.number import FP_MANTISA_BINARY_DIGITS, SpecialValueError, min_prec -from mathics.core.rules import Rule from mathics.core.symbols import Atom, Symbol, SymbolPlus, SymbolPower, SymbolTimes from mathics.core.systemsymbols import ( + SymbolAbs, SymbolComplexInfinity, + SymbolExp, SymbolI, SymbolIndeterminate, SymbolLog, + SymbolRealSign, + SymbolSign, + SymbolSqrt, ) RationalMOneHalf = Rational(-1, 2) +RealM0p5 = Real(-0.5) RealOne = Real(1.0) @@ -78,40 +83,284 @@ def eval_Abs(expr: BaseElement) -> Optional[BaseElement]: """ if expr is a number, return the absolute value. """ - if isinstance(expr, (Integer, Rational, Real)): - if expr.value >= 0: + + 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.has_form("Exp", 1): + exp = expr.elements[0] + if isinstance(exp, (Integer, Real, Rational)): return expr - return eval_multiply_numbers(*[IntegerM1, expr]) - if isinstance(expr, Complex): - re, im = expr.real, expr.imag - sqabs = eval_add_numbers( - eval_multiply_numbers(re, re), eval_multiply_numbers(im, im) - ) - return Expression(SymbolPower, sqabs, RationalOneHalf) + if isinstance(exp, Complex): + return Expression(SymbolExp, exp.real) + if expr.get_head() is SymbolTimes: + factors = [] + rest = [] + for x in expr.elements: + factor = eval_Abs(x) + if factor: + factors.append(factor) + else: + rest.append(x) + if factors: + return Expression(SymbolTimes, eval_multiply_numbers(*factors), *rest) + if test_nonnegative_arithmetic_expr(expr): + return expr + if test_negative_arithmetic_expr(expr): + return eval_multiply_numbers(IntegerM1, expr) return None -def eval_Sign(expr: BaseElement) -> Optional[BaseElement]: +def eval_Abs_number(n: Number) -> Number: """ - if expr is a number, return its sign. + Evals the absolute value of a number + """ + if isinstance(n, (Real, Integer)): + n_val = n.value + if n_val >= 0: + return n + return eval_negate_number(n) + 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, 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_Exp(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_RealSign(expr: BaseElement) -> Optional[Integer]: """ + If the argument is a real algebraic expression, + return the sign of the expression. + """ + if expr.is_zero: + return Integer0 if isinstance(expr, (Integer, Rational, Real)): - if expr.value > 0: + return Integer1 if expr.value > 0 else IntegerM1 + if expr in NUMERICAL_CONSTANTS: + return Integer1 + if expr.has_form("Abs", 1): + arg = expr.elements[0] + arg_sign = eval_Sign(arg) + if arg_sign is None: + return None + if arg_sign.is_zero: + return Integer0 + if isinstance(arg_sign, Number): return Integer1 - elif expr.value == 0: + return None + if expr.has_form("Sqrt", 1): + return Integer1 if eval_Sign(expr.elements[0]) is Integer1 else None + if expr.has_form("Exp", 1): + return Integer1 if test_arithmetic_expr(expr.elements[0]) else None + if expr.has_form("Log", 1) or expr.has_form("DirectedInfinity", 1): + return eval_RealSign(eval_add_numbers(expr.elements[0], IntegerM1)) + if expr.has_form("Times", None): + sign = 1 + for factor in expr.elements: + factor_sign = eval_RealSign(factor) + if factor_sign in (None, Integer0): + return factor_sign + if factor_sign is IntegerM1: + sign = -sign + return Integer1 if sign == 1 else IntegerM1 + if expr.has_form("Power", 2): + base, exp = expr.elements + base_sign = eval_RealSign(base) + if base_sign is None: + return None + if base_sign is Integer0: + if eval_RealSign(exp) in (IntegerM1, Integer0, None): + return None return Integer0 - else: + # The exponent must represent a real number to continue: + if not test_arithmetic_expr(exp): + return None + # At this point, the exponent is a real number, so if the base + # is 1, does not matter its value: + if base_sign is Integer1: + return Integer1 + if base_sign is IntegerM1: + if not isinstance(base, Integer): + return None + if isinstance(exp, Integer): + return base_sign if (exp.value % 2 == 1) else Integer1 + return None + if expr.has_form("Plus", None): + signed = {Integer1: [], IntegerM1: []} + for term in expr.elements: + rsign = eval_RealSign(term) + if rsign is Integer0: + continue + elif rsign is None: + return None + signed[rsign].append(term) + if len(signed[IntegerM1]) == 0: + return Integer0 if len(signed[Integer1]) == 0 else Integer1 + if len(signed[Integer1]) == 0: return IntegerM1 + # Try to explicitly add the numbers: + try_add = eval_add_numbers(*(term for term in expr.elements)) + if try_add is not None and not try_add.sameQ(expr): + return eval_RealSign(try_add) + # Now, try to convert to inexact values: + try_add = eval_add_numbers(*(to_inexact_value(term) for term in expr.elements)) + if try_add is not None and try_add is not expr: + return eval_RealSign(try_add) - if isinstance(expr, Complex): - re, im = expr.real, expr.imag - sqabs = eval_add_numbers(eval_Times(re, re), eval_Times(im, im)) - norm = Expression(SymbolPower, sqabs, RationalMOneHalf) - result = eval_Times(expr, norm) - if result is None: - return Expression(SymbolTimes, expr, norm) - return result - return None + +def eval_Sign(expr: BaseElement) -> Optional[BaseElement]: + """ + if expr is a number, return its sign. + """ + + def eval_complex_sign(n: BaseElement) -> Optional[BaseElement]: + 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)) + if isinstance(n, Atom): + return None + if n.has_form("Abs", 1): + inner_sign = eval_Sign(n.elements[0]) + if inner_sign is Integer0: + return Integer0 + if isinstance(inner_sign, Number): + return Integer1 + + if n.has_form("Exp", 1): + exponent = n.elements[0] + if isinstance(exponent, Complex): + return Expression(SymbolExp, exponent.imag) + return None + if n.has_form("DirectedInfinity", 1): + return eval_Sign(n.elements[0]) + if n.has_form("Power", 2): + base, exponent = expr.elements + base_rsign = eval_RealSign(base) + if exponent.is_zero: + return SymbolIndeterminate if base_rsign is Integer0 else Integer1 + if test_arithmetic_expr(exponent): + base_sign = eval_Sign(base) or Expression(SymbolSign, base) + return eval_Power_number(base_sign, exponent) + if isinstance(exponent, Complex): + if base_rsign is Integer1: + exp_im = exponent.imag + return eval_Power_number(base, Complex(Integer0, exp_im)) + + if test_arithmetic_expr(base): + eval_Power_number(base_sign, exponent) + base_sign = eval_Sign(base) + return eval_Power_number(base_sign, exponent) + if n.head is SymbolTimes: + signs = [] + for factor in expr.elements: + factor_sign = eval_Sign(factor) + if factor_sign in (None, Integer0): + return factor_sign + if factor_sign is not Integer1: + signs.append(factor_sign) + return Integer1 if len(signs) == 0 else eval_multiply_numbers(*signs) + + try_inexact = to_inexact_value(n) + if try_inexact: + return eval_Sign(try_inexact) + return None + + 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_mpmath_function( @@ -209,6 +458,150 @@ def append_last(): ) +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, prec) + + def eval_Times(*items: BaseElement) -> BaseElement: elements = [] numbers = [] @@ -323,7 +716,27 @@ def eval_add_numbers( return from_sympy(sum(item.to_sympy() for item in numbers)) -def eval_multiply_numbers(*numbers: Number) -> BaseElement: +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: Number) -> Number: """ Multiply the elements in ``numbers``. """ @@ -348,6 +761,19 @@ def eval_multiply_numbers(*numbers: Number) -> BaseElement: 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 segregate_numbers( *elements: BaseElement, ) -> Tuple[List[Number], List[BaseElement]]: @@ -398,11 +824,19 @@ def test_arithmetic_expr(expr: BaseElement, only_real: bool = True) -> bool: return not only_real if isinstance(expr, Symbol): return False + if isinstance(expr, Atom): + return False head, elements = expr.head, expr.elements if head in (SymbolPlus, SymbolTimes): return all(test_arithmetic_expr(term, only_real) for term in elements) + 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 expr.has_form("Exp", 1): return test_arithmetic_expr(elements[0], only_real) if head is SymbolLog: @@ -410,15 +844,16 @@ def test_arithmetic_expr(expr: BaseElement, only_real: bool = True) -> bool: return False if len(elements) == 2: base = elements[0] - if not test_positive_arithmetic_expr(base): + if only_real and eval_RealSign(base) is not Integer1: + return False + elif not test_arithmetic_expr(base): return False return test_arithmetic_expr(elements[-1], only_real) - if expr.has_form("Power", 2): - base, exponent = elements + if expr.has_form("Sqrt", 1): + radicand = elements[0] if only_real: - if isinstance(exponent, Integer): - return test_arithmetic_expr(base) - return all(test_arithmetic_expr(item, only_real) for item in elements) + return eval_RealSign(radicand) in (Integer0, Integer1) + return test_arithmetic_expr(radicand, only_real) return False @@ -427,11 +862,7 @@ 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) + return eval_RealSign(expr) is IntegerM1 def test_nonnegative_arithmetic_expr(expr: BaseElement) -> bool: @@ -439,11 +870,7 @@ def test_nonnegative_arithmetic_expr(expr: BaseElement) -> bool: Check if the expression is an arithmetic expression representing a nonnegative number """ - if not test_arithmetic_expr(expr): - return False - - if test_zero_arithmetic_expr(expr) or test_positive_arithmetic_expr(expr): - return True + return eval_RealSign(expr) in (Integer0, Integer1) def test_nonpositive_arithetic_expr(expr: BaseElement) -> bool: @@ -451,12 +878,7 @@ def test_nonpositive_arithetic_expr(expr: BaseElement) -> bool: Check if the expression is an arithmetic expression representing a nonnegative number """ - if not test_arithmetic_expr(expr): - return False - - if test_zero_arithmetic_expr(expr) or test_negative_arithmetic_expr(expr): - return True - return False + return eval_RealSign(expr) in (Integer0, IntegerM1) def test_positive_arithmetic_expr(expr: BaseElement) -> bool: @@ -464,79 +886,7 @@ 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(expr.elements[0], only_real=True) - if expr.has_form("Sqrt", 1): - return test_positive_arithmetic_expr(expr.elements[0]) - 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 expr.has_form("Abs", 1): - arg = elements[0] - return test_arithmetic_expr( - arg, only_real=False - ) and not test_zero_arithmetic_expr(arg) - if head.has_form("DirectedInfinity", 1): - return test_positive_arithmetic_expr(elements[0]) - - return False + return eval_RealSign(expr) is Integer1 def test_zero_arithmetic_expr(expr: BaseElement, numeric: bool = False) -> bool: @@ -544,47 +894,28 @@ 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: + if numeric: + if isinstance(expr, Complex): + if abs(expr.real.value) + abs(expr.imag.value) < 2.0e-10: return True - if isinstance(z, Number): - if abs(z.value) < 1e-10: + if isinstance(expr, Number): + if abs(expr.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 - if test_zero_arithmetic_expr(base, numeric): - return test_nonnegative_arithmetic_expr(exp) - if base.has_form("DirectedInfinity", None): - return test_positive_arithmetic_expr(exp) - 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 + + return eval_RealSign(expr) is Integer0 + + +EVAL_TO_INEXACT_DISPATCH = { + SymbolPlus: eval_add_numbers, + SymbolTimes: eval_multiply_numbers, + SymbolPower: eval_Power_number, + SymbolExp: eval_Exp, + SymbolSqrt: (lambda x: eval_Power_number(x, RationalOneHalf)), + SymbolAbs: eval_Abs, + SymbolSign: eval_Sign, + SymbolRealSign: eval_RealSign, +} def to_inexact_value(expr: BaseElement) -> BaseElement: @@ -595,9 +926,18 @@ def to_inexact_value(expr: BaseElement) -> BaseElement: """ if expr.is_inexact(): return expr + if isinstance(expr, Number): + return expr.round() + if expr is SymbolI: + return Complex(Integer0, RealOne) + if isinstance(expr, Symbol): + return NUMERICAL_CONSTANTS.get(expr, None) if isinstance(expr, Expression): - for const, value in NUMERICAL_CONSTANTS.items(): - expr, _ = expr.do_apply_rule(Rule(const, value)) - - return eval_multiply_numbers(RealOne, expr) + try: + head = expr.head + elements = tuple(to_inexact_value(element) for element in expr.elements) + return EVAL_TO_INEXACT_DISPATCH[head](*elements) + except Exception: + pass + return None diff --git a/mathics/eval/testing_expressions.py b/mathics/eval/testing_expressions.py index c15471f48..b0783ceb9 100644 --- a/mathics/eval/testing_expressions.py +++ b/mathics/eval/testing_expressions.py @@ -6,10 +6,12 @@ from mathics.core.expression import Expression from mathics.core.systemsymbols import SymbolDirectedInfinity - +# TODO: Remove me. The following function is not used anywhere +""" def cmp(a, b) -> int: "Returns 0 if a == b, -1 if a < b and 1 if a > b" return (a > b) - (a < b) +""" def do_cmp(x1, x2) -> Optional[int]: diff --git a/test/builtin/arithmetic/test_lowlevel_properties.py b/test/builtin/arithmetic/test_lowlevel_properties.py index 10cebfd06..14f3836ff 100644 --- a/test/builtin/arithmetic/test_lowlevel_properties.py +++ b/test/builtin/arithmetic/test_lowlevel_properties.py @@ -56,7 +56,7 @@ def test_positivity(str_expr, expected, msg): ("1", False, None), ("Pi", False, None), ("a", False, None), - ("a-a", True, None), + ("a-a", False, "the low-level check does not try to evaluate the input"), ("3-3.", True, None), ("2-Sqrt[4]", True, None), ("-Pi", False, None), @@ -73,9 +73,9 @@ def test_positivity(str_expr, expected, msg): ("Log[3]", False, None), ("Log[I]", False, None), ("Abs[a]", False, None), - ("Abs[0]", False, None), + ("Abs[0]", True, None), ("Abs[1+3 I]", False, None), - # ("Sin[Pi]", False, None), + ("Sin[Pi]", False, None), ], ) def test_zero(str_expr, expected, msg): From 4d8f5b18bdb2a7f9ab5a89c3b8dbd2ac48d9421b Mon Sep 17 00:00:00 2001 From: mmatera Date: Thu, 20 Jul 2023 22:03:25 -0300 Subject: [PATCH 13/14] tiny tweaks in documentation --- mathics/builtin/arithmetic.py | 26 +++++++++++++++---- .../equality_inequality.py | 1 - mathics/eval/testing_expressions.py | 7 ----- 3 files changed, 21 insertions(+), 13 deletions(-) diff --git a/mathics/builtin/arithmetic.py b/mathics/builtin/arithmetic.py index 51980d633..3ed71b068 100644 --- a/mathics/builtin/arithmetic.py +++ b/mathics/builtin/arithmetic.py @@ -829,7 +829,9 @@ def evaluate(self, evaluation: Evaluation): class Im(SympyFunction): """ - :WMA link:https://reference.wolfram.com/language/ref/Im.html + + :WMA link: + https://reference.wolfram.com/language/ref/Im.html
'Im[$z$]' @@ -1215,7 +1217,12 @@ class Real_(Builtin): class RealAbs(Builtin): """ - :WMA link:https://reference.wolfram.com/language/ref/RealAbs.html + + :Abs (Real): + https://en.wikipedia.org/wiki/Absolute_value ( + :WMA link: + https://reference.wolfram.com/language/ref/RealAbs.html + )
'RealAbs[$x$]' @@ -1283,8 +1290,12 @@ def test(self, expr) -> bool: class RealSign(Builtin): """ - :WMA link:https://reference.wolfram.com/language/ref/RealAbs.html - + + :Signum: + https://en.wikipedia.org/wiki/Sign_function ( + :WMA link: + https://reference.wolfram.com/language/ref/RealSign.html + )
'RealSign[$x$]'
returns $-1$, $0$ or $1$ depending on whether $x$ is negative, @@ -1319,7 +1330,12 @@ def eval(self, x: Number, evaluation: Evaluation) -> Optional[Integer]: class Sign(SympyFunction): """ - :WMA link:https://reference.wolfram.com/language/ref/Sign.html + + :Sign: + https://en.wikipedia.org/wiki/Sign_function ( + :WMA link: + https://reference.wolfram.com/language/ref/Sign.html + )
'Sign[$x$]' diff --git a/mathics/builtin/testing_expressions/equality_inequality.py b/mathics/builtin/testing_expressions/equality_inequality.py index d5e0ce6e5..59152e31e 100644 --- a/mathics/builtin/testing_expressions/equality_inequality.py +++ b/mathics/builtin/testing_expressions/equality_inequality.py @@ -69,7 +69,6 @@ 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, SymbolMaxExtraPrecision) n_items.append(item) items = n_items diff --git a/mathics/eval/testing_expressions.py b/mathics/eval/testing_expressions.py index b0783ceb9..4046d0c8c 100644 --- a/mathics/eval/testing_expressions.py +++ b/mathics/eval/testing_expressions.py @@ -6,13 +6,6 @@ from mathics.core.expression import Expression from mathics.core.systemsymbols import SymbolDirectedInfinity -# TODO: Remove me. The following function is not used anywhere -""" -def cmp(a, b) -> int: - "Returns 0 if a == b, -1 if a < b and 1 if a > b" - return (a > b) - (a < b) -""" - def do_cmp(x1, x2) -> Optional[int]: From c16e3bcf2f008f502d130d8c42f562fa1e22b9b8 Mon Sep 17 00:00:00 2001 From: mmatera Date: Thu, 20 Jul 2023 23:11:06 -0300 Subject: [PATCH 14/14] remove duplicated --- mathics/eval/arithmetic.py | 295 ------------------------------------- 1 file changed, 295 deletions(-) diff --git a/mathics/eval/arithmetic.py b/mathics/eval/arithmetic.py index 93501af17..8c25a9dc7 100644 --- a/mathics/eval/arithmetic.py +++ b/mathics/eval/arithmetic.py @@ -463,301 +463,6 @@ 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]: - """ - 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, prec) - - -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, prec) - def eval_Power_number(base: Number, exp: Number) -> Optional[Number]: """