diff --git a/mathics/builtin/numbers/calculus.py b/mathics/builtin/numbers/calculus.py index ad5c38a1d..556c8e5db 100644 --- a/mathics/builtin/numbers/calculus.py +++ b/mathics/builtin/numbers/calculus.py @@ -59,16 +59,19 @@ from mathics.core.systemsymbols import ( SymbolAnd, SymbolAutomatic, + SymbolComplex, SymbolConditionalExpression, SymbolD, SymbolDerivative, SymbolInfinity, SymbolInfix, + SymbolInteger, SymbolIntegrate, SymbolLeft, SymbolLog, SymbolNIntegrate, SymbolO, + SymbolReal, SymbolRule, SymbolSequence, SymbolSeries, @@ -76,6 +79,7 @@ SymbolSimplify, SymbolUndefined, ) +from mathics.eval.calculus import solve_sympy from mathics.eval.makeboxes import format_element from mathics.eval.nevaluator import eval_N from mathics.eval.numbers.calculus.integrators import ( @@ -330,9 +334,9 @@ def summand(element, index): Expression( SymbolDerivative, *( - [Integer0] * (index) - + [Integer1] - + [Integer0] * (len(f.elements) - index - 1) + [Integer0] * (index) + + [Integer1] + + [Integer0] * (len(f.elements) - index - 1) ), ), f.head, @@ -664,8 +668,8 @@ def eval(self, f, x, x0, evaluation: Evaluation, options: dict): # Determine the "jacobian"s if ( - method in ("Newton", "Automatic") - and options["System`Jacobian"] is SymbolAutomatic + method in ("Newton", "Automatic") and + options["System`Jacobian"] is SymbolAutomatic ): def diff(evaluation): @@ -1323,16 +1327,16 @@ class NIntegrate(Builtin): messages = { "bdmtd": "The Method option should be a built-in method name.", "inumr": ( - "The integrand `1` has evaluated to non-numerical " - + "values for all sampling points in the region " - + "with boundaries `2`" + "The integrand `1` has evaluated to non-numerical " + + "values for all sampling points in the region " + + "with boundaries `2`" ), "nlim": "`1` = `2` is not a valid limit of integration.", "ilim": "Invalid integration variable or limit(s) in `1`.", "mtdfail": ( - "The specified method failed to return a " - + "number. Falling back into the internal " - + "evaluator." + "The specified method failed to return a " + + "number. Falling back into the internal " + + "evaluator." ), "cmpint": ("Integration over a complex domain is not " + "implemented yet"), } @@ -1375,10 +1379,10 @@ class NIntegrate(Builtin): messages.update( { - "bdmtd": "The Method option should be a " - + "built-in method name in {`" - + "`, `".join(list(methods)) - + "`}. Using `Automatic`" + "bdmtd": "The Method option should be a " + + "built-in method name in {`" + + "`, `".join(list(methods)) + + "`}. Using `Automatic`" } ) @@ -1398,7 +1402,7 @@ def eval_with_func_domain( elif isinstance(method, Symbol): method = method.get_name() # strip context - method = method[method.rindex("`") + 1 :] + method = method[method.rindex("`") + 1:] else: evaluation.message("NIntegrate", "bdmtd", method) return @@ -2210,165 +2214,121 @@ class Solve(Builtin): messages = { "eqf": "`1` is not a well-formed equation.", "svars": 'Equations may not give solutions for all "solve" variables.', + "fulldim": "The solution set contains a full-dimensional component; use Reduce for complete solution information.", } - # FIXME: the problem with removing the domain parameter from the outside - # is that the we can't make use of this information inside - # the evaluation method where it is may be needed. rules = { - "Solve[eqs_, vars_, Complexes]": "Solve[eqs, vars]", - "Solve[eqs_, vars_, Reals]": ( - "Cases[Solve[eqs, vars], {Rule[x_,y_?RealValuedNumberQ]}]" - ), - "Solve[eqs_, vars_, Integers]": ( - "Cases[Solve[eqs, vars], {Rule[x_,y_Integer]}]" - ), + "Solve[eqs_, vars_]": "Solve[eqs, vars, Complexes]" } summary_text = "find generic solutions for variables" - def eval(self, eqs, vars, evaluation: Evaluation): - "Solve[eqs_, vars_]" + def eval(self, eqs, vars, domain, evaluation: Evaluation): + "Solve[eqs_, vars_, domain_]" - vars_original = vars - head_name = vars.get_head_name() + variables = vars + head_name = variables.get_head_name() if head_name == "System`List": - vars = vars.elements + variables = variables.elements else: - vars = [vars] - for var in vars: + variables = [variables] + for var in variables: if ( - (isinstance(var, Atom) and not isinstance(var, Symbol)) - or head_name in ("System`Plus", "System`Times", "System`Power") # noqa - or A_CONSTANT & var.get_attributes(evaluation.definitions) + (isinstance(var, Atom) and not isinstance(var, Symbol)) or + head_name in ("System`Plus", "System`Times", "System`Power") or # noqa + A_CONSTANT & var.get_attributes(evaluation.definitions) ): - evaluation.message("Solve", "ivar", vars_original) - return - if eqs.get_head_name() in ("System`List", "System`And"): - eq_list = eqs.elements - else: - eq_list = [eqs] - sympy_eqs = [] - sympy_denoms = [] - for eq in eq_list: - if eq is SymbolTrue: - pass - elif eq is SymbolFalse: - return ListExpression() - elif not eq.has_form("Equal", 2): - evaluation.message("Solve", "eqf", eqs) + evaluation.message("Solve", "ivar", vars) return - else: - left, right = eq.elements - left = left.to_sympy() - right = right.to_sympy() - if left is None or right is None: - return - eq = left - right - eq = sympy.together(eq) - eq = sympy.cancel(eq) - sympy_eqs.append(eq) - numer, denom = eq.as_numer_denom() - sympy_denoms.append(denom) - - vars_sympy = [var.to_sympy() for var in vars] - if None in vars_sympy: - return - # delete unused variables to avoid SymPy's - # PolynomialError: Not a zero-dimensional system - # in e.g. Solve[x^2==1&&z^2==-1,{x,y,z}] - all_vars = vars[:] - all_vars_sympy = vars_sympy[:] - vars = [] - vars_sympy = [] - for var, var_sympy in zip(all_vars, all_vars_sympy): - pattern = Pattern.create(var) - if not eqs.is_free(pattern, evaluation): - vars.append(var) - vars_sympy.append(var_sympy) - - def transform_dict(sols): - if not sols: - yield sols - for var, sol in sols.items(): - rest = sols.copy() - del rest[var] - rest = transform_dict(rest) - if not isinstance(sol, (tuple, list)): - sol = [sol] - if not sol: - for r in rest: - yield r - else: - for r in rest: - for item in sol: - new_sols = r.copy() - new_sols[var] = item - yield new_sols - break - - def transform_solution(sol): - if not isinstance(sol, dict): - if not isinstance(sol, (list, tuple)): - sol = [sol] - sol = dict(list(zip(vars_sympy, sol))) - return transform_dict(sol) - - if not sympy_eqs: - sympy_eqs = True - elif len(sympy_eqs) == 1: - sympy_eqs = sympy_eqs[0] - - try: - if isinstance(sympy_eqs, bool): - result = sympy_eqs + sympy_variables = [var.to_sympy() for var in variables] + if None in sympy_variables: + evaluation.message("Solve", "ivar") + return + variable_tuples = list(zip(variables, sympy_variables)) + + def solve_recur(expression: Expression): + '''solve And, Or and List within the scope of sympy, + but including the translation from Mathics to sympy + + returns: + solutions: a list of sympy solution dictionarys + conditions: a sympy condition object + + note: + for And and List, should always return either (solutions, None) or ([], conditions) + for Or, all combinations are possible. if Or is root, this should be handled outside''' + head = expression.get_head_name() + if head in ('System`And', 'System`List'): + solutions = [] + equations: list[Expression] = [] + inequations = [] + for child in expression.elements: + if child.has_form("Equal", 2): + equations.append(child) + elif child.get_head_name() in ('System`And', 'System`Or'): + sub_solution, sub_condition = solve_recur(child) + solutions.extend(sub_solution) + if sub_condition is not None: + inequations.append(sub_condition) + else: + inequations.append(child.to_sympy()) + solutions.extend(solve_sympy(evaluation, equations, variables, domain)) + conditions = sympy.And(*inequations) + result = [sol for sol in solutions if conditions.subs(sol)] + return result, None if solutions else conditions + else: # assume should be System`Or + assert head == 'System`Or' + solutions = [] + conditions = [] + for child in expression.elements: + if child.has_form("Equal", 2): + solutions.extend(solve_sympy(evaluation, child, variables, domain)) + elif child.get_head_name() in ('System`And', 'System`Or'): # I don't believe List would be in here + sub_solution, sub_condition = solve_recur(child) + solutions.extend(sub_solution) + if sub_condition is not None: + conditions.append(sub_condition) + else: + # SymbolTrue and SymbolFalse are allowed here since it's subtree context + # FIXME: None is not allowed, not sure what to do here + conditions.append(child.to_sympy()) + conditions = sympy.Or(*conditions) + return solutions, conditions + + if eqs.get_head_name() in ("System`List", "System`And", "System`Or"): + solutions, conditions = solve_recur(eqs) + # non True conditions are only accepted in subtrees only, not root + if conditions is not None: + evaluation.message("Solve", "fulldim") + else: + if eqs.get_head_name() == "System`Equal": + solutions = solve_sympy(evaluation, eqs, variables, domain) else: - result = sympy.solve(sympy_eqs, vars_sympy) - if not isinstance(result, list): - result = [result] - if isinstance(result, list) and len(result) == 1 and result[0] is True: + evaluation.message("Solve", "fulldim") return ListExpression(ListExpression()) - if result == [None]: - return ListExpression() - results = [] - for sol in result: - results.extend(transform_solution(sol)) - result = results - if any( - sol and any(var not in sol for var in all_vars_sympy) for sol in result - ): - evaluation.message("Solve", "svars") - # Filter out results for which denominator is 0 - # (SymPy should actually do that itself, but it doesn't!) - result = [ - sol - for sol in result - if all(sympy.simplify(denom.subs(sol)) != 0 for denom in sympy_denoms) - ] + if solutions is None: + evaluation.message("Solve", "ivars") + return ListExpression(ListExpression()) - return ListExpression( - *( - ListExpression( - *( - Expression(SymbolRule, var, from_sympy(sol[var_sympy])) - for var, var_sympy in zip(vars, vars_sympy) - if var_sympy in sol - ), - ) - for sol in result - ), - ) - except sympy.PolynomialError: - # raised for e.g. Solve[x^2==1&&z^2==-1,{x,y,z}] when not deleting - # unused variables beforehand - pass - except NotImplementedError: - pass - except TypeError as exc: - if str(exc).startswith("expected Symbol, Function or Derivative"): - evaluation.message("Solve", "ivar", vars_original) + if any( + sol and any(var not in sol for var in sympy_variables) for sol in solutions + ): + evaluation.message("Solve", "svars") + + return ListExpression( + *( + ListExpression( + *( + Expression(SymbolRule, var, from_sympy(sol[var_sympy])) + for var, var_sympy in variable_tuples + if var_sympy in sol + ), + ) + for sol in solutions + ), + ) # Auxiliary routines. Maybe should be moved to another module. diff --git a/mathics/core/atoms.py b/mathics/core/atoms.py index eadc66ad0..64249c2dd 100644 --- a/mathics/core/atoms.py +++ b/mathics/core/atoms.py @@ -774,9 +774,9 @@ def get_sort_key(self, pattern_sort=False) -> tuple: def sameQ(self, other) -> bool: """Mathics SameQ""" return ( - isinstance(other, Complex) - and self.real == other.real - and self.imag == other.imag + isinstance(other, Complex) and + self.real == other.real and + self.imag == other.imag ) def round(self, d=None) -> "Complex": diff --git a/mathics/core/convert/sympy.py b/mathics/core/convert/sympy.py index 843d679f3..fa333b906 100644 --- a/mathics/core/convert/sympy.py +++ b/mathics/core/convert/sympy.py @@ -5,6 +5,7 @@ Conversion to SymPy is handled directly in BaseElement descendants. """ +from collections.abc import Iterable from typing import Optional, Type, Union import sympy @@ -13,9 +14,6 @@ # Import the singleton class from sympy.core.numbers import S -BasicSympy = sympy.Expr - - from mathics.core.atoms import ( MATHICS3_COMPLEX_I, Complex, @@ -40,6 +38,7 @@ ) from mathics.core.list import ListExpression from mathics.core.number import FP_MANTISA_BINARY_DIGITS +from mathics.core.rules import Pattern from mathics.core.symbols import ( Symbol, SymbolFalse, @@ -62,16 +61,21 @@ SymbolGreater, SymbolGreaterEqual, SymbolIndeterminate, + SymbolIntegers, SymbolLess, SymbolLessEqual, SymbolMatrixPower, SymbolO, SymbolPi, SymbolPiecewise, + SymbolReals, SymbolSlot, SymbolUnequal, ) +BasicSympy = sympy.Expr + + SymbolPrime = Symbol("Prime") SymbolRoot = Symbol("Root") SymbolRootSum = Symbol("RootSum") @@ -130,6 +134,39 @@ def to_sympy_matrix(data, **kwargs) -> Optional[sympy.MutableDenseMatrix]: return None +def apply_domain_to_symbols(symbols: Iterable[sympy.Symbol], domain) -> dict[sympy.Symbol, sympy.Symbol]: + """Create new sympy symbols with domain applied. + Return a dict maps old to new. + """ + # FIXME: this substitute solution would break when Solve[Abs[x]==3, x],where x=-3 and x=3. + # However, substituting symbol prior to actual solving would cause sympy to have biased assumption, + # it would refuse to solve Abs() when symbol is in Complexes + result = {} + for symbol in symbols: + if domain == SymbolReals: + new_symbol = sympy.Symbol(repr(symbol), real=True) + elif domain == SymbolIntegers: + new_symbol = sympy.Symbol(repr(symbol), integer=True) + else: + new_symbol = symbol + result[symbol] = new_symbol + return result + + +def cut_dimension(evaluation, expressions: Union[Expression, list[Expression]], symbols: Iterable[sympy.Symbol]) -> set[sympy.Symbol]: + '''delete unused variables to avoid SymPy's PolynomialError + : Not a zero-dimensional system in e.g. Solve[x^2==1&&z^2==-1,{x,y,z}]''' + if not isinstance(expressions, list): + expressions = [expressions] + subset = set() + for symbol in symbols: + pattern = Pattern.create(symbol) + for equation in expressions: + if not equation.is_free(pattern, evaluation): + subset.add(symbol) + return subset + + class SympyExpression(BasicSympy): is_Function = True nargs = None @@ -363,9 +400,9 @@ def old_from_sympy(expr) -> BaseElement: if is_Cn_expr(name): return Expression(SymbolC, Integer(int(name[1:]))) if name.startswith(sympy_symbol_prefix): - name = name[len(sympy_symbol_prefix) :] + name = name[len(sympy_symbol_prefix):] if name.startswith(sympy_slot_prefix): - index = name[len(sympy_slot_prefix) :] + index = name[len(sympy_slot_prefix):] return Expression(SymbolSlot, Integer(int(index))) elif expr.is_NumberSymbol: name = str(expr) @@ -517,7 +554,7 @@ def old_from_sympy(expr) -> BaseElement: *[from_sympy(arg) for arg in expr.args] ) if name.startswith(sympy_symbol_prefix): - name = name[len(sympy_symbol_prefix) :] + name = name[len(sympy_symbol_prefix):] args = [from_sympy(arg) for arg in expr.args] builtin = sympy_to_mathics.get(name) if builtin is not None: diff --git a/mathics/core/systemsymbols.py b/mathics/core/systemsymbols.py index fa388a92d..fb1ea6d79 100644 --- a/mathics/core/systemsymbols.py +++ b/mathics/core/systemsymbols.py @@ -56,6 +56,7 @@ SymbolCompile = Symbol("System`Compile") SymbolCompiledFunction = Symbol("System`CompiledFunction") SymbolComplex = Symbol("System`Complex") +SymbolComplexes = Symbol("System`Complexes") SymbolComplexInfinity = Symbol("System`ComplexInfinity") SymbolCondition = Symbol("System`Condition") SymbolConditionalExpression = Symbol("System`ConditionalExpression") @@ -124,6 +125,7 @@ SymbolInfix = Symbol("System`Infix") SymbolInputForm = Symbol("System`InputForm") SymbolInteger = Symbol("System`Integer") +SymbolIntegers = Symbol("System`Integers") SymbolIntegrate = Symbol("System`Integrate") SymbolLeft = Symbol("System`Left") SymbolLength = Symbol("System`Length") @@ -200,6 +202,7 @@ SymbolRational = Symbol("System`Rational") SymbolRe = Symbol("System`Re") SymbolReal = Symbol("System`Real") +SymbolReals = Symbol("System`Reals") SymbolRealAbs = Symbol("System`RealAbs") SymbolRealDigits = Symbol("System`RealDigits") SymbolRealSign = Symbol("System`RealSign") diff --git a/test/builtin/calculus/test_solve.py b/test/builtin/calculus/test_solve.py index 7b1c79f07..234a641dc 100644 --- a/test/builtin/calculus/test_solve.py +++ b/test/builtin/calculus/test_solve.py @@ -42,6 +42,36 @@ def test_solve(): "{x->1.51213}", "Issue #1235", ), + ( + "Solve[Abs[-2/3*(lambda + 2) + 8/3 + 4] == 4, lambda,Reals]", + "{{lambda -> 2}, {lambda -> 14}}", + "abs()", + ), + ( + "Solve[q^3 == (20-12)/(4-3), q,Reals]", + "{{q -> 2}}", + "domain check", + ), + ( + "Solve[x + Pi/3 == 2k*Pi + Pi/6 || x + Pi/3 == 2k*Pi + 5Pi/6, x,Reals]", + "{{x -> -Pi / 6 + 2 k Pi}, {x -> Pi / 2 + 2 k Pi}}", + "logics involved", + ), + ( + "Solve[m - 1 == 0 && -(m + 1) != 0, m,Reals]", + "{{m -> 1}}", + "logics and constraints", + ), + ( + "Solve[(lambda + 1)/6 == 1/(mu - 1) == lambda/4, {lambda, mu},Reals]", + "{{lambda -> 2, mu -> 3}}", + "chained equations", + ), + ( + "Solve[2*x0*Log[x0] + x0 - 2*a*x0 == -1 && x0^2*Log[x0] - a*x0^2 + b == b - x0, {x0, a, b},Reals]", + "{{x0 -> 1, a -> 1}}", + "excess variable b", + ), ): session.evaluate("Clear[h]; Clear[g]; Clear[f];") check_evaluation(str_expr, str_expected, message)