diff --git a/CHANGES.rst b/CHANGES.rst index 780c19762..e20db6344 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -9,6 +9,7 @@ New Builtins * ``$SessionID`` * ``BinaryReadList`` +* ``Hypergeometric2F1`` Internals --------- @@ -17,6 +18,19 @@ Mathics scanner exceptions of class TranslateError are incompatible with previous versions, and now store error parameters, "name", "tag", and "args". +WMA Compatibility +----------------- + +Hypergeometric functions have been revised to conform better to WMA behavior by +expanding hypergeometric results. + +Bugs +---- + +* Fixed #1187 + + + 8.0.1 ----- diff --git a/mathics/builtin/files_io/files.py b/mathics/builtin/files_io/files.py index a14da432d..1b198afe7 100644 --- a/mathics/builtin/files_io/files.py +++ b/mathics/builtin/files_io/files.py @@ -49,6 +49,7 @@ read_name_and_stream, ) from mathics.eval.makeboxes import do_format, format_element +from mathics.eval.stackframe import get_eval_Expression class Input_(Predefined): @@ -567,7 +568,7 @@ class Put(InfixOperator): summary_text = "write an expression to a file" - def eval(self, exprs, filename, evaluation): + def eval(self, exprs, filename, evaluation: Evaluation): "Put[exprs___, filename_String]" instream = to_expression("OpenWrite", filename).evaluate(evaluation) if len(instream.elements) == 2: @@ -581,7 +582,7 @@ def eval(self, exprs, filename, evaluation): close_stream(py_instream, instream_number) return result - def eval_input(self, exprs, name, n, evaluation): + def eval_input(self, exprs, name, n, evaluation: Evaluation): "Put[exprs___, OutputStream[name_, n_]]" stream = stream_manager.lookup_stream(n.get_int_value()) @@ -608,11 +609,10 @@ def do_format_output(expr, evaluation): return SymbolNull - def eval_default(self, exprs, filename, evaluation): + def eval_default(self, exprs, filename, evaluation: Evaluation): "Put[exprs___, filename_]" - expr = to_expression("Put", exprs, filename) evaluation.message("Put", "stream", filename) - return expr + return to_expression("Put", exprs, filename) class PutAppend(InfixOperator): @@ -661,7 +661,7 @@ class PutAppend(InfixOperator): summary_text = "append an expression to a file" - def eval(self, exprs, filename, evaluation): + def eval(self, exprs, filename: String, evaluation): "PutAppend[exprs___, filename_String]" instream = to_expression("OpenAppend", filename).evaluate(evaluation) if len(instream.elements) == 2: @@ -672,12 +672,12 @@ def eval(self, exprs, filename, evaluation): to_expression("Close", instream).evaluate(evaluation) return result - def eval_input(self, exprs, name, n, evaluation): + def eval_input(self, exprs, name, n, evaluation: Evaluation): "PutAppend[exprs___, OutputStream[name_, n_]]" stream = stream_manager.lookup_stream(n.get_int_value()) if stream is None or stream.io.closed: - evaluation.message("Put", "openx", to_expression("OutputSteam", name, n)) + evaluation.message("Put", "openx", get_eval_Expression()) return text = [ @@ -691,11 +691,10 @@ def eval_input(self, exprs, name, n, evaluation): return SymbolNull - def eval_default(self, exprs, filename, evaluation): + def eval_default(self, exprs, filename, evaluation: Evaluation): "PutAppend[exprs___, filename_]" - expr = to_expression("PutAppend", exprs, filename) evaluation.message("PutAppend", "stream", filename) - return expr + return get_eval_Expression() def validate_read_type(name: str, typ, evaluation: Evaluation): @@ -1100,11 +1099,9 @@ def eval_n(self, file, types, n: Integer, evaluation: Evaluation, options: dict) # token_words = py_options['TokenWords'] # word_separators = py_options['WordSeparators'] - py_n = n.get_int_value() + py_n = n.value if py_n < 0: - evaluation.message( - "ReadList", "intnm", to_expression("ReadList", file, types, n) - ) + evaluation.message("ReadList", "intnm", get_eval_Expression()) return result = [] @@ -1215,9 +1212,7 @@ def eval_input(self, name, n, m, evaluation): seekpos = m.to_python() if not (isinstance(seekpos, int) or seekpos == float("inf")): - evaluation.message( - "SetStreamPosition", "stmrng", to_expression("InputStream", name, n), m - ) + evaluation.message("SetStreamPosition", "stmrng", get_eval_Expression(), m) return try: @@ -1312,7 +1307,7 @@ def eval(self, name, n, typ, m, evaluation: Evaluation, options: dict): evaluation.message( "Skip", "intm", - to_expression("Skip", to_expression("InputStream", name, n), typ, m), + to_expression("Skip", stream, typ, m), ) return for i in range(py_m): @@ -1478,7 +1473,7 @@ def eval(self, evaluation): "Streams[]" return self.eval_name(None, evaluation) - def eval_name(self, name, evaluation): + def eval_name(self, name: String, evaluation): "Streams[name_String]" result = [] for stream in stream_manager.STREAMS.values(): diff --git a/mathics/builtin/messages.py b/mathics/builtin/messages.py index 09b82b98e..8168d3d33 100644 --- a/mathics/builtin/messages.py +++ b/mathics/builtin/messages.py @@ -192,6 +192,7 @@ class General(Builtin): "fnsym": ( "First argument in `1` is not a symbol " "or a string naming a symbol." ), + "hdiv": "`1` does not exist. Arguments are not consistent.", "heads": "Heads `1` and `2` are expected to be the same.", "ilsnn": ( "Single or list of non-negative integers expected at " "position `1`." diff --git a/mathics/builtin/specialfns/hypergeom.py b/mathics/builtin/specialfns/hypergeom.py index c83bd576f..397ff4792 100644 --- a/mathics/builtin/specialfns/hypergeom.py +++ b/mathics/builtin/specialfns/hypergeom.py @@ -7,22 +7,138 @@ """ import mpmath -import sympy import mathics.eval.tracing as tracing from mathics.core.attributes import ( A_LISTABLE, - A_N_HOLD_FIRST, A_NUMERIC_FUNCTION, A_PROTECTED, A_READ_PROTECTED, ) from mathics.core.builtin import MPMathFunction from mathics.core.convert.mpmath import from_mpmath -from mathics.core.convert.sympy import from_sympy from mathics.core.evaluation import Evaluation -from mathics.core.number import FP_MANTISA_BINARY_DIGITS +from mathics.core.list import ListExpression from mathics.core.systemsymbols import SymbolComplexInfinity, SymbolMachinePrecision +from mathics.eval.specialfns.hypergeom import ( + eval_Hypergeometric1F1, + eval_Hypergeometric2F1, + eval_HypergeometricPQF, + eval_MeijerG, + eval_N_HypergeometricPQF, +) +from mathics.eval.stackframe import get_eval_Expression + + +class Hypergeometric1F1(MPMathFunction): + """ + + :Kummer confluent hypergeometric function: https://en.wikipedia.org/wiki/Confluent_hypergeometric_function ( + :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#hyp1f1, + :WMA: https://reference.wolfram.com/language/ref/Hypergeometric1F1.html) +
+
'Hypergeometric1F1'[$a$, $b$, $z$] +
returns ${}_1 F_1(a; b; z)$. +
+ + Numeric evaluation: + >> Hypergeometric1F1[1, 2, 3.0] + = 6.36185 + + Plot over a subset of reals: + >> Plot[Hypergeometric1F1[1, 2, x], {x, -5, 5}] + = -Graphics- + + >> Plot[{Hypergeometric1F1[1/2, Sqrt[2], x], Hypergeometric1F1[1/2, Sqrt[3], x], Hypergeometric1F1[1/2, Sqrt[5], x]}, {x, -4, 4}] + = -Graphics- + + >> Plot[{Hypergeometric1F1[Sqrt[3], Sqrt[2], z], -0.01}, {z, -10, -2}] + = -Graphics- + + >> Plot[{Hypergeometric1F1[Sqrt[2], b, 1], Hypergeometric1F1[Sqrt[5], b, 1], Hypergeometric1F1[Sqrt[7], b, 1]}, {b, -3, 3}] + = -Graphics- + + Compute the elementwise values of an array: + >> Hypergeometric1F1[1, 1, {{1, 0}, {0, 1}}] + = {{E, 1}, {1, E}} + + >> Hypergeometric1F1[1/2, 1, x] + = BesselI[0, x / 2] E ^ (x / 2) + + Evaluate using complex arguments: + >> Hypergeometric1F1[2 + I, 2, 0.5] + = 1.61833 + 0.379258 I + + Large numbers are supported: + >> Hypergeometric1F1[3, 4, 10^10] + = -3 / 500000000000000000000000000000 + 149999999970000000003 E ^ 10000000000 / 500000000000000000000000000000 + + 'Hypergeometric1F1' evaluates to simpler functions for certain parameters: + >> Hypergeometric1F1[1/2, 1, x] + = BesselI[0, x / 2] E ^ (x / 2) + + >> Hypergeometric1F1[2, 1, x] + = (1 + x) E ^ x + + >> Hypergeometric1F1[1, 1/2, x] + = -Sqrt[x] (-E ^ (-x) / Sqrt[x] - Sqrt[Pi] Erf[Sqrt[x]]) E ^ x + """ + + attributes = A_LISTABLE | A_NUMERIC_FUNCTION | A_PROTECTED + mpmath_name = "hyp1f1" + nargs = {3} + + rules = { + "Hypergeometric1F1[0, c_, z_?MachineNumberQ]": "1.0", + "Hypergeometric1F1[0, c_, z_]": "1", + } + + summary_text = "compute Kummer confluent hypergeometric function" + sympy_name = "" + + def eval(self, a, b, z, evaluation: Evaluation): + "Hypergeometric1F1[a_, b_, z_]" + + return eval_Hypergeometric1F1(a, b, z) + + +class Hypergeometric2F1(MPMathFunction): + """ + + :Hypergeometric function: https://en.wikipedia.org/wiki/Hypergeometric_function ( + :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#mpmath.hyp2f1, + :WMA: https://reference.wolfram.com/language/ref/Hypergeometric2F1.html) +
+
'Hypergeometric2F1'[$a$, $b$, $c$, $z$] +
returns ${}_2 F_1(a; b; c; z)$. +
+ + >> Hypergeometric2F1[2., 3., 4., 5.0] + = 0.156542 + 0.150796 I + + Evaluate symbolically: + >> Hypergeometric2F1[2, 3, 4, x] + = 6 Log[1 - x] / x ^ 3 + (-6 + 3 x) / (-x ^ 2 + x ^ 3) + + Evaluate using complex arguments: + >> Hypergeometric2F1[2 + I, -I, 3/4, 0.5 - 0.5 I] + = -0.972167 - 0.181659 I + + Plot over a subset of the reals: + >> Plot[Hypergeometric2F1[1/3, 1/3, 2/3, x], {x, -1, 1}] + = -Graphics- + """ + + attributes = A_LISTABLE | A_NUMERIC_FUNCTION | A_PROTECTED + mpmath_name = "hyp2f1" + nargs = {4} + summary_text = "compute Gauss hypergeometric function function" + sympy_name = "hyper" + + def eval(self, a, b, c, z, evaluation: Evaluation): + "Hypergeometric2F1[a_, b_, c_, z_]" + + return eval_Hypergeometric2F1(a, b, c, z) class HypergeometricPFQ(MPMathFunction): @@ -30,7 +146,7 @@ class HypergeometricPFQ(MPMathFunction): :Generalized hypergeometric function: https://en.wikipedia.org/wiki/Generalized_hypergeometric_function ( :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#hyper, - :sympy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.hyper, + :SymPy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.hyper, :WMA: https://reference.wolfram.com/language/ref/HypergeometricPFQ.html)
'HypergeometricPFQ'[${a_1, ..., a_p}, {b_1, ..., b_q}, z$] @@ -41,95 +157,99 @@ class HypergeometricPFQ(MPMathFunction): Result is symbollicaly simplified by default: >> HypergeometricPFQ[{3}, {2}, 1] - = HypergeometricPFQ[{3}, {2}, 1] + = 3 E / 2 + unless a numerical evaluation is explicitly requested: >> HypergeometricPFQ[{3}, {2}, 1] // N = 4.07742 + >> HypergeometricPFQ[{3}, {2}, 1.] = 4.07742 + >> Plot[HypergeometricPFQ[{1, 1}, {3, 3, 3}, x], {x, -30, 30}] + = -Graphics- + + >> HypergeometricPFQ[{1, 1, 2}, {3, 3}, z] + = -4 PolyLog[2, z] / z ^ 2 + 4 Log[1 - z] / z ^ 2 - 4 Log[1 - z] / z + 8 / z + The following special cases are handled: >> HypergeometricPFQ[{}, {}, z] - = 1 + = E ^ z >> HypergeometricPFQ[{0}, {b}, z] = 1 - >> Hypergeometric1F1[b, b, z] - = E ^ z + + >> HypergeometricPFQ[{1, 1, 3}, {2, 2}, x] + = -Log[1 - x] / (2 x) - 1 / (-2 + 2 x) + + 'HypergeometricPFQ' evaluates to a polynomial if any of the parameters $a_k$ is a non-positive integer: + >> HypergeometricPFQ[{-2, a}, {b}, x] + = (-2 a x (1 + b) + a x ^ 2 (1 + a) + b (1 + b)) / (b (1 + b)) + + Value at origin: + >> HypergeometricPFQ[{a1, b2, a3}, {b1, b2, b3}, 0] + = 1 """ attributes = A_NUMERIC_FUNCTION | A_PROTECTED | A_READ_PROTECTED mpmath_name = "hyper" nargs = {3} - rules = { - "HypergeometricPFQ[{}, {}, z_]": "1", - "HypergeometricPFQ[{0}, b_, z_]": "1", - "HypergeometricPFQ[b_, b_, z_]": "Exp[z]", - } + summary_text = "compute the generalized hypergeometric function" sympy_name = "hyper" def eval(self, a, b, z, evaluation: Evaluation): "HypergeometricPFQ[a_, b_, z_]" - try: - a_sympy = [e.to_sympy() for e in a] - b_sympy = [e.to_sympy() for e in b] - result_sympy = tracing.run_sympy( - sympy.hyper, a_sympy, b_sympy, z.to_sympy() - ) - return from_sympy(result_sympy) - except Exception: - pass + + # FIXME: a lot more checking could be done here. + if not (isinstance(a, ListExpression)): + evaluation.message("HypergeometricPFQ", "hdiv", get_eval_Expression()) + + return eval_HypergeometricPQF(a, b, z) def eval_N(self, a, b, z, prec, evaluation: Evaluation): "N[HypergeometricPFQ[a_, b_, z_], prec_]" - try: - result_mpmath = tracing.run_mpmath( - mpmath.hyper, a.to_python(), b.to_python(), z.to_python() - ) - return from_mpmath(result_mpmath) - except ZeroDivisionError: - return SymbolComplexInfinity - except Exception as ex: - pass - - def eval_numeric(self, a, b, z, evaluation: Evaluation): - "HypergeometricPFQ[a:{__?NumericQ}, b:{__?NumericQ}, z_?MachineNumberQ]" - return self.eval_N(a, b, z, SymbolMachinePrecision, evaluation) + # FIXME: prec is not used. It should be though. + return eval_N_HypergeometricPQF(a, b, z) -class Hypergeometric1F1(MPMathFunction): +class HypergeometricU(MPMathFunction): """ - :Kummer confluent hypergeometric function: https://en.wikipedia.org/wiki/Confluent_hypergeometric_function ( - :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#hyper, - :WMA: https://reference.wolfram.com/language/ref/Hypergeometric1F1.html) + :Confluent hypergeometric function: https://en.wikipedia.org/wiki/Confluent_hypergeometric_function ( + :mpmath: https://mpmath.org/doc/current/functions/bessel.html#mpmath.hyperu, + :WMA: https://reference.wolfram.com/language/ref/HypergeometricU.html)
-
'Hypergeometric1F1'[$a$, $b$, $z$] -
returns ${}_1 F_1(a; b; z)$. +
'HypergeometricU'[$a$, $b$, $z$] +
returns $U(a, b, z)$.
- - Result is symbollicaly simplified by default: - >> Hypergeometric1F1[3, 2, 1] - = HypergeometricPFQ[{3}, {2}, 1] + Result is symbollicaly simplified, where possible: + >> HypergeometricU[3, 2, 1] + = MeijerG[{{-2}, {}}, {{0, -1}, {}}, 1] / 2 + >> HypergeometricU[1,4,8] + = HypergeometricU[1, 4, 8] unless a numerical evaluation is explicitly requested: - >> Hypergeometric1F1[3, 2, 1] // N - = 4.07742 - >> Hypergeometric1F1[3, 2, 1.] - = 4.07742 + >> HypergeometricU[3, 2, 1] // N + = 0.105479 + >> HypergeometricU[3, 2, 1.] + = 0.105479 - Plot 'M'[3, 2, x] from 0 to 2 in steps of 0.5: - >> Plot[Hypergeometric1F1[3, 2, x], {x, 0.5, 2}] + Plot 'U'[3, 2, x] from 0 to 10 in steps of 0.5: + >> Plot[HypergeometricU[3, 2, x], {x, 0.5, 10}] = -Graphics- - Here, plot explicitly requests a numerical evaluation. + + We handle this special case: + >> HypergeometricU[0, b, z] + = 1 """ - attributes = A_LISTABLE | A_NUMERIC_FUNCTION | A_PROTECTED - mpmath_name = "" + attributes = A_LISTABLE | A_NUMERIC_FUNCTION | A_PROTECTED | A_READ_PROTECTED + mpmath_name = "hyperu" nargs = {3} rules = { - "Hypergeometric1F1[a_, b_, z_]": "HypergeometricPFQ[{a},{b},z]", + "HypergeometricU[0, c_, z_]": "1", + "HypergeometricU[a_, b_, z_] /; (a > 0) && (a-b+1 > 0)": "MeijerG[{{1-a},{}},{{0,1-b},{}},z]/Gamma[a]/Gamma[a-b+1]", } - summary_text = "compute Kummer confluent hypergeometric function" + summary_text = "compute the Tricomi confluent hypergeometric function" sympy_name = "" @@ -138,7 +258,7 @@ class MeijerG(MPMathFunction): :Meijer G-function: https://en.wikipedia.org/wiki/Meijer_G-function ( :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#meijerg, - :sympy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.meijerg, + :SymPy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.meijerg, :WMA: https://reference.wolfram.com/language/ref/MeijerG.html)
'MeijerG'[${{a_1, ..., a_n}, {a_{n+1}, ..., a_p}}, {{b_1, ..., b_m}, {b_{m+1}, ..., a_q}}, z$] @@ -163,16 +283,10 @@ class MeijerG(MPMathFunction): def eval(self, a, b, z, evaluation: Evaluation): "MeijerG[a_, b_, z_]" - try: - a_sympy = [[e2.to_sympy() for e2 in e1] for e1 in a] - b_sympy = [[e2.to_sympy() for e2 in e1] for e1 in b] - result_sympy = tracing.run_sympy( - sympy.meijerg, a_sympy, b_sympy, z.to_sympy() - ) - return from_sympy(result_sympy) - except Exception: - pass + return eval_MeijerG(a, b, z) + # FIXME: Moving this causes test/doc/test_common.py to fail! It + # reports that Hypergeometric has been included more than once Weird! def eval_N(self, a, b, z, prec, evaluation: Evaluation): "N[MeijerG[a_, b_, z_], prec_]" try: @@ -188,44 +302,3 @@ def eval_N(self, a, b, z, prec, evaluation: Evaluation): def eval_numeric(self, a, b, z, evaluation: Evaluation): "MeijerG[a:{___List?(AllTrue[#, NumericQ, Infinity]&)}, b:{___List?(AllTrue[#, NumericQ, Infinity]&)}, z_?MachineNumberQ]" return self.eval_N(a, b, z, SymbolMachinePrecision, evaluation) - - -class HypergeometricU(MPMathFunction): - """ - - :Confluent hypergeometric function: https://en.wikipedia.org/wiki/Confluent_hypergeometric_function ( - :mpmath: https://mpmath.org/doc/current/functions/bessel.html#mpmath.hyperu, - :WMA: https://reference.wolfram.com/language/ref/HypergeometricU.html) -
-
'HypergeometricU'[$a$, $b$, $z$] -
returns $U(a, b, z)$. -
- Result is symbollicaly simplified, where possible: - >> HypergeometricU[3, 2, 1] - = MeijerG[{{-2}, {}}, {{0, -1}, {}}, 1] / 2 - >> HypergeometricU[1,4,8] - = HypergeometricU[1, 4, 8] - unless a numerical evaluation is explicitly requested: - >> HypergeometricU[3, 2, 1] // N - = 0.105479 - >> HypergeometricU[3, 2, 1.] - = 0.105479 - - Plot 'U'[3, 2, x] from 0 to 10 in steps of 0.5: - >> Plot[HypergeometricU[3, 2, x], {x, 0.5, 10}] - = -Graphics- - - We handle this special case: - >> HypergeometricU[0, b, z] - = 1 - """ - - attributes = A_LISTABLE | A_NUMERIC_FUNCTION | A_PROTECTED | A_READ_PROTECTED - mpmath_name = "hyperu" - nargs = {3} - rules = { - "HypergeometricU[0, c_, z_]": "1", - "HypergeometricU[a_, b_, z_] /; (a > 0) && (a-b+1 > 0)": "MeijerG[{{1-a},{}},{{0,1-b},{}},z]/Gamma[a]/Gamma[a-b+1]", - } - summary_text = "compute the Tricomi confluent hypergeometric function" - sympy_name = "" diff --git a/mathics/builtin/trace.py b/mathics/builtin/trace.py index f8b3954a3..dc27fde52 100644 --- a/mathics/builtin/trace.py +++ b/mathics/builtin/trace.py @@ -357,11 +357,22 @@ class TraceEvaluation(Builtin): ## :trace native symbol:
-
'TraceEvaluation'[$expr$] +
'TraceEvaluation'[$expr$, $options$]
Evaluate $expr$ and print each step of the evaluation.
- The 'ShowTimeBySteps' option prints the elapsed time before an evaluation occurs. + Options adjust output and filtering behavior +
+
'ShowTimeBySteps' +
Print the elapsed time before an evaluation occurs. \ + default is 'False'. +
'ShowEvaluation' +
Show evaluation calls and returns. The default is 'True'. +
'ShowRewrite' +
Show the effect of rewrite rules. The default is 'True'. +
+ + Note: It does not make sense to set both 'ShowRewrite' and 'ShowEvaluation' to 'False'. >> TraceEvaluation[(x + x)^2] | ... @@ -370,43 +381,72 @@ class TraceEvaluation(Builtin): >> TraceEvaluation[(x + x)^2, ShowTimeBySteps->True] | ... = ... + + Now consider this function which consists of a function call that involves a rewrite rule: + >> TraceEvaluation[BesselK[0, 0]] + | ... + = ... + + Sometimes, 'TraceEvaluation' traces can get quite large. To reduce the size, it may be helpful \ + to filter on either the evaluations or the replacement rules. + + To see just the evaluations and return values, but not rewrite that occurs: + >> TraceEvaluation[BesselK[0, 0], ShowRewrite-> False] + | ... + = ... + + To see just the rewrite that occurs, which tends to summarizes even more: + >> TraceEvaluation[BesselK[0, 0], ShowEvaluation-> False] + | ... + = ... + """ attributes = A_HOLD_ALL | A_PROTECTED options = { "System`ShowTimeBySteps": "False", + "System`ShowRewrite": "True", # Do we want to see rewrite rules? + "System`ShowEvaluation": "True", # Do we want to see Evaluate and Returns? } summary_text = "trace expression evaluation" def eval(self, expr, evaluation: Evaluation, options: dict): "TraceEvaluation[expr_, OptionsPattern[]]" - curr_trace_evaluation = evaluation.definitions.trace_evaluation - curr_time_by_steps = evaluation.definitions.timing_trace_evaluation - + # Save various trace settings before changing them. old_evaluation_call_hook = mathics.eval.tracing.trace_evaluate_on_call old_evaluation_return_hook = mathics.eval.tracing.trace_evaluate_on_return + old_time_by_steps = evaluation.definitions.timing_trace_evaluation + old_trace_evaluation = evaluation.definitions.trace_evaluation + old_trace_show_rewrite = evaluation.definitions.trace_show_rewrite + # Adjust trace settings based on the options given. + evaluation.definitions.timing_trace_evaluation = ( + options["System`ShowTimeBySteps"] is SymbolTrue + ) + evaluation.definitions.trace_evaluation = ( + options["System`ShowEvaluation"] is SymbolTrue + ) + evaluation.definitions.trace_show_rewrite = ( + options["System`ShowRewrite"] is SymbolTrue + ) mathics.eval.tracing.trace_evaluate_on_call = ( mathics.eval.tracing.print_evaluate ) - mathics.eval.tracing.trace_evaluate_on_return = ( mathics.eval.tracing.print_evaluate ) - evaluation.definitions.trace_evaluation = True - evaluation.definitions.timing_trace_evaluation = ( - options["System`ShowTimeBySteps"] is SymbolTrue - ) + # Now perform the evaluation... try: return expr.evaluate(evaluation) except Exception: raise finally: - evaluation.definitions.trace_evaluation = curr_trace_evaluation - evaluation.definitions.timing_trace_evaluation = curr_time_by_steps - + # Restore settings to the way the were before the TraceEvaluation. + evaluation.definitions.trace_evaluation = old_trace_evaluation + evaluation.definitions.timing_trace_evaluation = old_time_by_steps + evaluation.definitions.trace_show_rewrite = old_trace_show_rewrite mathics.eval.tracing.trace_evaluate_on_call = old_evaluation_call_hook mathics.eval.tracing.trace_evaluate_on_return = old_evaluation_return_hook diff --git a/mathics/core/atoms.py b/mathics/core/atoms.py index a5cadae53..61b075fe4 100644 --- a/mathics/core/atoms.py +++ b/mathics/core/atoms.py @@ -532,6 +532,7 @@ def to_sympy(self, *args, **kwargs): MachineReal0 = MachineReal(0) +MachineReal1 = MachineReal(1) class PrecisionReal(Real[sympy.Float]): diff --git a/mathics/core/definitions.py b/mathics/core/definitions.py index fbe01b36a..96e866f28 100644 --- a/mathics/core/definitions.py +++ b/mathics/core/definitions.py @@ -162,7 +162,10 @@ def __init__( ) self.inputfile = "" + # These are used by TraceEvaluation to``` + # whether what information to show. self.trace_evaluation = False + self.trace_show_rewrite = False self.timing_trace_evaluation = False # Importing "mathics.format" populates the Symbol of the diff --git a/mathics/eval/specialfns/hypergeom.py b/mathics/eval/specialfns/hypergeom.py new file mode 100644 index 000000000..eff373b0c --- /dev/null +++ b/mathics/eval/specialfns/hypergeom.py @@ -0,0 +1,204 @@ +""" +Evaluation functions for built-in Hypergeometric functions +""" + +from typing import Sequence, Tuple, cast + +import mpmath +import sympy + +import mathics.eval.tracing as tracing +from mathics.core.atoms import Complex, Integer1, MachineReal1, Number +from mathics.core.convert.mpmath import from_mpmath +from mathics.core.convert.sympy import from_sympy +from mathics.core.number import RECONSTRUCT_MACHINE_PRECISION_DIGITS, dps, min_prec +from mathics.core.systemsymbols import SymbolComplexInfinity +from mathics.eval.arithmetic import eval_mpmath_function + + +def eval_Hypergeometric1F1(a, b, z): + """Hypergeometric2F1[a_, b_, z_] + + We use SymPy's hyper and expand the symbolic result. + But if that fails to expand and all arugments are numeric, use mpmath hyp1f1. + + We prefer SymPy because it preserves constants like E whereas mpmath will + convert E to a precisioned number. + """ + + # SymPy returns E ^ z for Hypergeometric1F1[0,0,z], but + # WMA gives 1. Therefore, we add the below code to give the WMA + # behavior. If SymPy switches, this code be eliminated. + if hasattr(a, "is_zero") and a.is_zero: + return ( + MachineReal1 + if a.is_machine_precision() + or hasattr(z, "machine_precision") + and z.is_machine_precision() + else Integer1 + ) + + args = (a, b, z) + sympy_args, all_numeric = to_sympy_with_classification(args) + sympy_result = tracing.run_sympy( + sympy.hyper, [sympy_args[0]], [sympy_args[1]], sympy_args[2] + ) + expanded_result = sympy.hyperexpand(sympy_result) + + # Oddly, SymPy expansion sometimes doesn't work when complex arguments are given. + # However mpmath can handle this. + # I imagine at some point in the future this will be fixed. + if isinstance(expanded_result, sympy.hyper) and all_numeric: + args = cast(Sequence[Number], args) + + if any(arg.is_machine_precision() for arg in args): + prec = None + else: + prec = min_prec(*args) + if prec is None: + prec = RECONSTRUCT_MACHINE_PRECISION_DIGITS + d = dps(prec) + args = tuple([arg.round(d) for arg in args]) + + return eval_mpmath_function( + mpmath.hyp1f1, *cast(Sequence[Number], args), prec=prec + ) + else: + return from_sympy(expanded_result) + + +def eval_Hypergeometric2F1(a, b, c, z): + """Hypergeometric2F1[a_, b_, c_, z_] + + Uses mpmath hyp2f1 if all args are numeric. Otherwise, + we use SymPy's hyper and expand the symbolic result. + """ + + args = (a, b, c, z) + sympy_args, all_numeric = to_sympy_with_classification(args) + if all_numeric: + args = cast(Sequence[Number], args) + + if any(arg.is_machine_precision() for arg in args): + prec = None + else: + prec = min_prec(*args) + if prec is None: + prec = RECONSTRUCT_MACHINE_PRECISION_DIGITS + d = dps(prec) + args = tuple([arg.round(d) for arg in args]) + + return eval_mpmath_function( + mpmath.hyp2f1, *cast(Sequence[Number], args), prec=prec + ) + else: + return run_sympy_hyper_and_expand( + [sympy_args[0], sympy_args[1]], [sympy_args[2]], sympy_args[3] + ) + + +def eval_HypergeometricPQF(a, b, z): + "HypergeometricPFQ[a_, b_, z_]" + + # SymPy returns E for HypergeometricPFQ[{0},{0},Number], but + # WMA gives 1. Therefore, we add the below code to give the WMA + # behavior. If SymPy switches, this code be eliminated. + if ( + len(a.elements) > 0 + and hasattr(a[0], "is_zero") + and a[0].is_zero + and isinstance(z, Number) + ): + return MachineReal1 if a[0].is_machine_precision() else Integer1 + + # FIXME: This isn't complete. If parameters "a" or "b" contain MachineReal + # numbers then the results should be MachineReal as well. + if z.is_machine_precision(): + return eval_N_HypergeometricPQF(a, b, z) + + try: + a_sympy = [e.to_sympy() for e in a] + b_sympy = [e.to_sympy() for e in b] + return run_sympy_hyper_and_expand(a_sympy, b_sympy, z.to_sympy()) + except Exception: + return None + + +# FIXME, this should take a "prec" parameter? +def eval_N_HypergeometricPQF(a, b, z): + "N[HypergeometricPFQ[a_, b_, z_]]" + try: + result_mpmath = tracing.run_mpmath( + mpmath.hyper, a.to_python(), b.to_python(), z.to_python() + ) + return from_mpmath(result_mpmath) + except ZeroDivisionError: + return SymbolComplexInfinity + except Exception: + return None + + +def eval_MeijerG(a, b, z): + """ + Use sympy.meijerg to compute MeijerG(a, b, z) + """ + try: + a_sympy = [[e2.to_sympy() for e2 in e1] for e1 in a] + b_sympy = [[e2.to_sympy() for e2 in e1] for e1 in b] + result_sympy = tracing.run_sympy(sympy.meijerg, a_sympy, b_sympy, z.to_sympy()) + # For now, we don't allow simplification and conversion + # to other functions like Bessel, because this can introduce + # SymPy's exp_polar() function for which we don't have a + # Mathics3 equivalent for yet. + # return from_sympy(sympy.hyperexpand(result_sympy)) + return from_sympy(result_sympy) + except Exception: + return None + + +# FIXME: ADDING THIS causes test/doc/test_common.py to fail! It reports that Hypergeometric has been included more than once +# Weird! + +# def eval_N_MeijerG(a, b, z): +# "N[MeijerG[a_, b_, z_], prec_]" +# try: +# result_mpmath = tracing.run_mpmath( +# mpmath.meijerg, a.to_python(), b.to_python(), z.to_python() +# ) +# return from_mpmath(result_mpmath) +# except ZeroDivisionError: +# return SymbolComplexInfinity +# except Exception: +# return None + + +def run_sympy_hyper_and_expand(a, b, z): + """Ruy SymPy's hyper function and expand the result.""" + sympy_result = tracing.run_sympy(sympy.hyper, a, b, z) + result = sympy.hyperexpand(sympy_result) + return from_sympy(result) + + +def to_sympy_with_classification(args: tuple) -> Tuple[list, bool]: + """Converts `args` to its corresponding SymPy form. However, if + all elements of args are numeric, then we detect and report that. + We do this so that the caller can decide whether to use mpmath if + SymPy fails. One might expect SymPy to do this automatically, but + it doesn't catch all opportunites. + """ + sympy_args = [] + all_numeric = True + for arg in args: + if isinstance(arg, Number): + # FIXME: in the future, .value + # should work on Complex numbers. + if isinstance(arg, Complex): + sympy_arg = arg.to_python() + else: + sympy_arg = arg.value + else: + sympy_arg = arg.to_sympy() + all_numeric = False + + sympy_args.append(sympy_arg) + return sympy_args, all_numeric diff --git a/mathics/eval/stackframe.py b/mathics/eval/stackframe.py new file mode 100644 index 000000000..f65e1465b --- /dev/null +++ b/mathics/eval/stackframe.py @@ -0,0 +1,82 @@ +""" +Mathics3 Introspection routines to get Mathics3-oriented Python frames. +""" + +import inspect +from types import FrameType +from typing import Optional + +from mathics.core.builtin import Builtin +from mathics.core.expression import Expression + + +def find_Mathics3_evaluation_method(frame: Optional[FrameType]) -> Optional[FrameType]: + """ + Returns the most recent Mathics3 evaluation frame given "frame". + """ + if not inspect.isframe(frame): + return None + + while frame is not None: + if is_Mathics3_eval_method(frame): + return frame + frame = frame.f_back + return None + + +def get_self(frame: FrameType) -> Optional[FrameType]: + """Returns the `self` object in a Python frame if it exists, or None if + it does not exist. + """ + return frame.f_locals.get("self", None) + + +def get_eval_Expression() -> Optional[Expression]: + """Returns the Expression used in a Mathics3 evaluation() or + eval() Builtin method. This is often needed in an error message to + report what Expression was getting evaluated. None is returned if + not found. + + The method is fragile in that it relies on the Mathics3 implementation + having a Expression.rewrite_apply_eval_step() method. It walks to + call stack to find that Expression. + + None is returned if we can't find this. + """ + frame = inspect.currentframe() + if frame is None: + return None + + frame = frame.f_back + while True: + if frame is None: + return None + method_code = frame.f_code + if method_code.co_name == "rewrite_apply_eval_step": + if (self_obj := get_self(frame)) is not None and isinstance( + self_obj, Expression + ): + return self_obj + + frame = frame.f_back + + +def is_Mathics3_eval_method(frame) -> bool: + """ + Returns True if frame is Python frame object for a Mathics3 Builtin, and + returns False otherwise. + """ + if not inspect.isframe(frame): + return False + + method_code = frame.f_code + if not method_code.co_name.startswith("eval"): + return False + if frame.f_locals.get("evaluation", None) is None: + return False + + if (self_obj := get_self(frame)) is None or frame.f_locals.get( + "evaluation", None + ) is None: + return False + return isinstance(self_obj, Builtin) diff --git a/mathics/eval/tracing.py b/mathics/eval/tracing.py index d458f41f4..11fee0593 100644 --- a/mathics/eval/tracing.py +++ b/mathics/eval/tracing.py @@ -102,7 +102,7 @@ def print_evaluate(expr, evaluation, status: str, fn: Callable, orig_expr=None): assert isinstance(expr, tuple) if orig_expr != expr[0]: if status == "Returning": - if expr[1]: + if expr[1] and evaluation.definitions.trace_show_rewrites: status = "Rewriting" arrow = " -> " else: @@ -116,11 +116,17 @@ def print_evaluate(expr, evaluation, status: str, fn: Callable, orig_expr=None): ) else: if status == "Returning" and isinstance(expr, tuple): - status = "Evaluating/Replacing" + if not evaluation.definitions.trace_show_rewrite: + return + status = "Replacing" expr = expr[0] + elif not evaluation.definitions.trace_evaluation: + return evaluation.print_out(f"{indents}{status}: {orig_expr} = " + str(expr)) elif not is_performing_rewrite(fn): + if not evaluation.definitions.trace_evaluation: + return evaluation.print_out(f"{indents}{status}: {expr}") return diff --git a/test/builtin/specialfns/test_hypergeom.py b/test/builtin/specialfns/test_hypergeom.py index 7bdbc594c..8e408fea3 100644 --- a/test/builtin/specialfns/test_hypergeom.py +++ b/test/builtin/specialfns/test_hypergeom.py @@ -7,6 +7,50 @@ import pytest +@pytest.mark.parametrize( + ("str_expr", "msgs", "str_expected", "fail_msg"), + [ + ( + "Hypergeometric1F1[{3,5},{7,1},{4,2}]", + None, + "{-1545 / 128 + 45 E ^ 4 / 128, 27 E ^ 2}", + None, + ), + ("N[Hypergeometric1F1[{3,5},{7,1},{4,2}]]", None, "{7.12435, 199.505}", None), + ("Hypergeometric1F1[b,b,z]", None, "E ^ z", None), + ( + "Hypergeometric1F1[0,0,z]", + None, + "1", + "Special case: Integer 1 (SymPy may work differently)", + ), + ( + "Hypergeometric1F1[0.0,0,z]", + None, + "1.", + "Special case: MachineReal a parameter (SymPy may work differently)", + ), + ( + "Hypergeometric1F1[0,0,1.]", + None, + "1.", + "Special case: MachineReal z parameter (SymPy may work differently)", + ), + ], +) +def test_Hypergeometric1F1(str_expr, msgs, str_expected, fail_msg): + """ """ + check_evaluation( + str_expr, + str_expected, + to_string_expr=True, + to_string_expected=True, + hold_expected=True, + failure_message=fail_msg, + expected_messages=msgs, + ) + + @pytest.mark.parametrize( ("str_expr", "msgs", "str_expected", "fail_msg"), [ @@ -16,37 +60,83 @@ "ComplexInfinity", None, ), + ("HypergeometricPFQ[{6},{1},2]", None, "719 E ^ 2 / 15", None), + ("N[HypergeometricPFQ[{6},{1},2]]", None, "354.182", None), + ("HypergeometricPFQ[{},{},z]", None, "E ^ z", None), + ("HypergeometricPFQ[{0},{c1,c2},z]", None, "1", None), + ("HypergeometricPFQ[{c1,c2},{c1,c2},z]", None, "E ^ z", None), ( - "MeijerG[{{},{}},{{0,0},{0,0}},100^2]", - None, - "MeijerG[{{}, {}}, {{0, 0}, {0, 0}}, 10000]", + "HypergeometricPFQ[{0},{0},2]", None, + "1", + "Special case: using Integer 'z' parameter", ), - ("N[MeijerG[{{},{}},{{0,0},{0,0}},100^2]]", None, "0.000893912", None), ( - "HypergeometricU[{3,1},{2,4},{7,8}]", + "HypergeometricPFQ[{0},{0},3.0]", None, - "{MeijerG[{{-2}, {}}, {{0, -1}, {}}, 7] / 2, HypergeometricU[1, 4, 8]}", + "1", + "Special case: using MachineInteger 'z' parameter", + ), + ( + "HypergeometricPFQ[{0.0},{0},3.0]", None, + "1.", + "Special case: using MachineInteger 'a' parameter", ), + ], +) +def test_HypergeometricPFQ(str_expr, msgs, str_expected, fail_msg): + """ """ + check_evaluation( + str_expr, + str_expected, + to_string_expr=True, + to_string_expected=True, + hold_expected=True, + failure_message=fail_msg, + expected_messages=msgs, + ) + + +@pytest.mark.parametrize( + ("str_expr", "msgs", "str_expected", "fail_msg"), + [ ("N[HypergeometricU[{3,1},{2,4},{7,8}]]", None, "{0.00154364, 0.160156}", None), ("HypergeometricU[0,c,z]", None, "1", None), + ], +) +def test_HypergeometricU(str_expr, msgs, str_expected, fail_msg): + """ """ + check_evaluation( + str_expr, + str_expected, + to_string_expr=True, + to_string_expected=True, + hold_expected=True, + failure_message=fail_msg, + expected_messages=msgs, + ) + + +@pytest.mark.parametrize( + ("str_expr", "msgs", "str_expected", "fail_msg"), + [ ( - "Hypergeometric1F1[{3,5},{7,1},{4,2}]", + "MeijerG[{{},{}},{{0,0},{0,0}},100^2]", None, - "{HypergeometricPFQ[{3}, {7}, 4], HypergeometricPFQ[{5}, {1}, 2]}", + "MeijerG[{{}, {}}, {{0, 0}, {0, 0}}, 10000]", + None, + ), + ("N[MeijerG[{{},{}},{{0,0},{0,0}},100^2]]", None, "0.000893912", None), + ( + "HypergeometricU[{3,1},{2,4},{7,8}]", + None, + "{MeijerG[{{-2}, {}}, {{0, -1}, {}}, 7] / 2, HypergeometricU[1, 4, 8]}", None, ), - ("N[Hypergeometric1F1[{3,5},{7,1},{4,2}]]", None, "{7.12435, 199.505}", None), - ("Hypergeometric1F1[b,b,z]", None, "E ^ z", None), - ("HypergeometricPFQ[{6},{1},2]", None, "HypergeometricPFQ[{6}, {1}, 2]", None), - ("N[HypergeometricPFQ[{6},{1},2]]", None, "354.182", None), - ("HypergeometricPFQ[{},{},z]", None, "1", None), - ("HypergeometricPFQ[{0},{c1,c2},z]", None, "1", None), - ("HypergeometricPFQ[{c1,c2},{c1,c2},z]", None, "E ^ z", None), ], ) -def test_private_hypergeom(str_expr, msgs, str_expected, fail_msg): +def test_MeijerG(str_expr, msgs, str_expected, fail_msg): """ """ check_evaluation( str_expr, diff --git a/test/builtin/test_trace.py b/test/builtin/test_trace.py index a5c159371..e1a085975 100644 --- a/test/builtin/test_trace.py +++ b/test/builtin/test_trace.py @@ -119,9 +119,9 @@ def capture_print(s: str): assert [ " Evaluating: System`Plus[System`Times[2, 3], 4]", " Evaluating: System`Times[2, 3]", - " Evaluating/Replacing: System`Times[2, 3] = 6", + " Replacing: System`Times[2, 3] = 6", " Returning: System`Times[2, 3] = 6", - " Evaluating/Replacing: System`Plus[System`Times[2, 3], 4] = 10", + " Replacing: System`Plus[System`Times[2, 3], 4] = 10", " Returning: System`Plus[System`Times[2, 3], 4] = 10", ] == event_queue # print()