Move more fns to mpmath.eval.specialfns...

Also some results were corrected by making more use of SymPy's
hyperexpand function.

Author: rocky

mathics/builtin/specialfns/hypergeom.py

@@ -22,6 +22,7 @@
 from mathics.core.evaluation import Evaluation
 from mathics.core.systemsymbols import SymbolComplexInfinity, SymbolMachinePrecision
 from mathics.eval.specialfns.hypergeom import (
+    eval_Hypergeometric1F1,
     eval_Hypergeometric2F1,
     eval_HypergeometricPQF,
     eval_N_HypergeometricPQF,
@@ -44,30 +45,42 @@ 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"
 
@@ -96,30 +109,55 @@ class Hypergeometric1F1(MPMathFunction):
       <dd>returns ${}_1 F_1(a; b; z)$.
     </dl>
 
-    Result is symbollicaly simplified by default:
-    >> Hypergeometric1F1[3, 2, 1]
-     = HypergeometricPFQ[{3}, {2}, 1]
-    unless a numerical evaluation is explicitly requested:
-    >> Hypergeometric1F1[3, 2, 1] // N
-     = 4.07742
-    >> Hypergeometric1F1[3, 2, 1.]
-     = 4.07742
+    Numeric evaluation:
+    >> Hypergeometric1F1[1, 2, 3.0]
+     = 6.36185
 
-    Plot 'M'[3, 2, x] from 0 to 2 in steps of 0.5:
-    >> Plot[Hypergeometric1F1[3, 2, x], {x, 0.5, 2}]
+    Plot over a subset of reals:
+    >> Plot[Hypergeometric1F1[1, 2, x], {x, -5, 5}]
      = -Graphics-
-    Here, plot explicitly requests a numerical evaluation.
+
+    >> 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
+
+    '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 = "hymp1f1"
     nargs = {3}
-    rules = {
-        "Hypergeometric1F1[a_, b_, z_]": "HypergeometricPFQ[{a},{b},z]",
-    }
     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):
     """

mathics/eval/specialfns/hypergeom.py

@@ -8,14 +8,69 @@
 import sympy
 
 import mathics.eval.tracing as tracing
-from mathics.core.atoms import Number
+from mathics.core.atoms import Complex, 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.
+    """
+
+    args = (a, b, z)
+
+    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)
+
+    sympy_result = tracing.run_sympy(
+        sympy.hyper, [sympy_args[0]], [sympy_args[1]], sympy_args[2]
+    )
+    expanded_result = sympy.hyperexpand(sympy_result)
+
+    # Oddly, expansion sometimes doesn't work for 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_]
 
@@ -28,7 +83,10 @@ def eval_Hypergeometric2F1(a, b, c, z):
     all_numeric = True
     for arg in args:
         if isinstance(arg, Number):
-            sympy_arg = arg.value
+            if isinstance(arg, Complex):
+                sympy_arg = arg.to_python()
+            else:
+                sympy_arg = arg.value
         else:
             sympy_arg = arg.to_sympy()
             all_numeric = False
@@ -51,19 +109,17 @@ def eval_Hypergeometric2F1(a, b, c, z):
             mpmath.hyp2f1, *cast(Sequence[Number], args), prec=prec
         )
     else:
-        sympy_result = tracing.run_sympy(
-            sympy.hyper, [sympy_args[0], sympy_args[1]], [sympy_args[2]], sympy_args[3]
+        return run_sympy_hyper(
+            [sympy_args[0], sympy_args[1]], [sympy_args[2]], sympy_args[3]
         )
-        return from_sympy(sympy.hyperexpand(sympy_result))
 
 
 def eval_HypergeometricPQF(a, b, z):
     "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)
+        return run_sympy_hyper(a_sympy, b_sympy, z.to_sympy())
     except Exception:
         return None
 
@@ -80,3 +136,9 @@ def eval_N_HypergeometricPQF(a, b, z):
         return SymbolComplexInfinity
     except Exception:
         return None
+
+
+def run_sympy_hyper(a, b, z):
+    sympy_result = tracing.run_sympy(sympy.hyper, a, b, z)
+    result = sympy.hyperexpand(sympy_result)
+    return from_sympy(result)

test/builtin/specialfns/test_hypergeom.py

@@ -34,14 +34,14 @@
         (
             "Hypergeometric1F1[{3,5},{7,1},{4,2}]",
             None,
-            "{HypergeometricPFQ[{3}, {7}, 4], HypergeometricPFQ[{5}, {1}, 2]}",
+            "{-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),
-        ("HypergeometricPFQ[{6},{1},2]", None, "HypergeometricPFQ[{6}, {1}, 2]", None),
+        ("HypergeometricPFQ[{6},{1},2]", None, "719 E ^ 2 / 15", None),
         ("N[HypergeometricPFQ[{6},{1},2]]", None, "354.182", None),
-        ("HypergeometricPFQ[{},{},z]", None, "1", 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),
     ],