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..3ed71b068 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 @@ -823,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$]' @@ -1207,6 +1215,49 @@ class Real_(Builtin): name = "Real" +class RealAbs(Builtin): + """ + + :Abs (Real): + https://en.wikipedia.org/wiki/Absolute_value ( + :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,9 +1288,54 @@ def test(self, expr) -> bool: return isinstance(expr, (Integer, Rational, Real)) +class RealSign(Builtin): + """ + + :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, + 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 + + :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 f95289f17..59152e31e 100644 --- a/mathics/builtin/testing_expressions/equality_inequality.py +++ b/mathics/builtin/testing_expressions/equality_inequality.py @@ -69,8 +69,7 @@ def numerify_args(items, evaluation) -> list: n_items = [] for item in items: if not isinstance(item, Number): - # TODO: use $MaxExtraPrecision insterad of hard-coded 50 - item = eval_N(item, evaluation, SymbolMaxPrecision) + item = eval_N(item, evaluation, SymbolMaxExtraPrecision) n_items.append(item) items = n_items else: 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..4046d0c8c 100644 --- a/mathics/eval/testing_expressions.py +++ b/mathics/eval/testing_expressions.py @@ -7,11 +7,6 @@ from mathics.core.systemsymbols import SymbolDirectedInfinity -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]: # don't attempt to compare complex numbers 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):