From c46f005c2588d75fa92c88fe8e3e1e206ac98c8b Mon Sep 17 00:00:00 2001
From: rocky <rb@dustyfeet.com>
Date: Sat, 26 Jul 2025 17:24:00 -0400
Subject: [PATCH 01/17] First cut at Hypergeometric2F1

---
 CHANGES.rst                             |  1 +
 mathics/builtin/specialfns/hypergeom.py | 69 ++++++++++++++-------
 mathics/eval/specialfns/hypergeom.py    | 82 +++++++++++++++++++++++++
 3 files changed, 131 insertions(+), 21 deletions(-)
 create mode 100644 mathics/eval/specialfns/hypergeom.py

diff --git a/CHANGES.rst b/CHANGES.rst
index 780c19762..fc4af20b8 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -9,6 +9,7 @@ New Builtins
 
 * ``$SessionID``
 * ``BinaryReadList``
+* ``Hypergeometric2F1``
 
 Internals
 ---------
diff --git a/mathics/builtin/specialfns/hypergeom.py b/mathics/builtin/specialfns/hypergeom.py
index c83bd576f..4513e5423 100644
--- a/mathics/builtin/specialfns/hypergeom.py
+++ b/mathics/builtin/specialfns/hypergeom.py
@@ -12,7 +12,6 @@
 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,
@@ -21,8 +20,12 @@
 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.systemsymbols import SymbolComplexInfinity, SymbolMachinePrecision
+from mathics.eval.specialfns.hypergeom import (
+    eval_Hypergeometric2F1,
+    eval_HypergeometricPQF,
+    eval_N_HypergeometricPQF,
+)
 
 
 class HypergeometricPFQ(MPMathFunction):
@@ -70,27 +73,12 @@ class HypergeometricPFQ(MPMathFunction):
 
     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
+        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
+        # FIXME: prec is not used. Why?
+        return eval_N_HypergeometricPQF(a, b, z)
 
     def eval_numeric(self, a, b, z, evaluation: Evaluation):
         "HypergeometricPFQ[a:{__?NumericQ}, b:{__?NumericQ}, z_?MachineNumberQ]"
@@ -124,7 +112,7 @@ class Hypergeometric1F1(MPMathFunction):
     """
 
     attributes = A_LISTABLE | A_NUMERIC_FUNCTION | A_PROTECTED
-    mpmath_name = ""
+    mpmath_name = "hymp1f1"
     nargs = {3}
     rules = {
         "Hypergeometric1F1[a_, b_, z_]": "HypergeometricPFQ[{a},{b},z]",
@@ -133,6 +121,45 @@ class Hypergeometric1F1(MPMathFunction):
     sympy_name = ""
 
 
+class Hypergeometric2F1(MPMathFunction):
+    """
+    <url>
+    :Hypergeometric function: https://en.wikipedia.org/wiki/Hypergeometric_function</url> (<url>
+    :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#mpmath.hyp2f1</url>, <url>
+    :WMA: https://reference.wolfram.com/language/ref/Hypergeometric2F1.html</url>)
+    <dl>
+      <dt>'Hypergeometric2F1'[$a$, $b$, $c$, $z$]
+      <dd>returns ${}_2 F_1(a; b; c; z)$.
+    </dl>
+
+    >> 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 MeijerG(MPMathFunction):
     """
     <url>
diff --git a/mathics/eval/specialfns/hypergeom.py b/mathics/eval/specialfns/hypergeom.py
new file mode 100644
index 000000000..3499e3160
--- /dev/null
+++ b/mathics/eval/specialfns/hypergeom.py
@@ -0,0 +1,82 @@
+"""
+Evaluation functions for built-in Hypergeometric functions
+"""
+
+from typing import Sequence, cast
+
+import mpmath
+import sympy
+
+import mathics.eval.tracing as tracing
+from mathics.core.atoms import 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_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 = True
+    for arg in args:
+        if isinstance(arg, Number):
+            sympy_arg = arg.value
+        else:
+            sympy_arg = arg.to_sympy()
+            all_numeric = False
+
+        sympy_args.append(sympy_arg)
+
+    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:
+        sympy_result = tracing.run_sympy(
+            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)
+    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

From 0c9c25ecce7bcfb9473bce72c849d9831b80d52b Mon Sep 17 00:00:00 2001
From: rocky <rb@dustyfeet.com>
Date: Sat, 26 Jul 2025 20:11:58 -0400
Subject: [PATCH 02/17] Move more fns to mpmath.eval.specialfns...

Also some results were corrected by making more use of SymPy's
hyperexpand function.
---
 mathics/builtin/specialfns/hypergeom.py   | 92 +++++++++++++++++------
 mathics/eval/specialfns/hypergeom.py      | 76 +++++++++++++++++--
 test/builtin/specialfns/test_hypergeom.py |  6 +-
 3 files changed, 139 insertions(+), 35 deletions(-)

diff --git a/mathics/builtin/specialfns/hypergeom.py b/mathics/builtin/specialfns/hypergeom.py
index 4513e5423..2d1433f9d 100644
--- a/mathics/builtin/specialfns/hypergeom.py
+++ b/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"
 
@@ -89,37 +102,66 @@ class Hypergeometric1F1(MPMathFunction):
     """
     <url>
     :Kummer confluent hypergeometric function: https://en.wikipedia.org/wiki/Confluent_hypergeometric_function</url> (<url>
-    :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#hyper</url>, <url>
+    :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#hyp1f1</url>, <url>
     :WMA: https://reference.wolfram.com/language/ref/Hypergeometric1F1.html</url>)
     <dl>
       <dt>'Hypergeometric1F1'[$a$, $b$, $z$]
       <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
+
+    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 = "hymp1f1"
+    mpmath_name = "hyp1f1"
     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):
     """
diff --git a/mathics/eval/specialfns/hypergeom.py b/mathics/eval/specialfns/hypergeom.py
index 3499e3160..dd33a17c1 100644
--- a/mathics/eval/specialfns/hypergeom.py
+++ b/mathics/eval/specialfns/hypergeom.py
@@ -8,7 +8,7 @@
 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
@@ -16,6 +16,61 @@
 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,10 +109,9 @@ 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):
@@ -62,8 +119,7 @@ def eval_HypergeometricPQF(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)
diff --git a/test/builtin/specialfns/test_hypergeom.py b/test/builtin/specialfns/test_hypergeom.py
index 7bdbc594c..9f20b6329 100644
--- a/test/builtin/specialfns/test_hypergeom.py
+++ b/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),
     ],

From 2a3fb4fe44c78bd6f7c433bc3467d7e06803a5c7 Mon Sep 17 00:00:00 2001
From: rocky <rb@dustyfeet.com>
Date: Sun, 27 Jul 2025 09:23:12 -0400
Subject: [PATCH 03/17] Move more functions to eval. Simplify

---
 mathics/builtin/specialfns/hypergeom.py | 95 +++++++++++--------------
 mathics/eval/specialfns/hypergeom.py    | 94 +++++++++++++++---------
 2 files changed, 105 insertions(+), 84 deletions(-)

diff --git a/mathics/builtin/specialfns/hypergeom.py b/mathics/builtin/specialfns/hypergeom.py
index 2d1433f9d..01a577497 100644
--- a/mathics/builtin/specialfns/hypergeom.py
+++ b/mathics/builtin/specialfns/hypergeom.py
@@ -7,7 +7,6 @@
 """
 
 import mpmath
-import sympy
 
 import mathics.eval.tracing as tracing
 from mathics.core.attributes import (
@@ -18,13 +17,13 @@
 )
 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.systemsymbols import SymbolComplexInfinity, SymbolMachinePrecision
 from mathics.eval.specialfns.hypergeom import (
     eval_Hypergeometric1F1,
     eval_Hypergeometric2F1,
     eval_HypergeometricPQF,
+    eval_MeijerG,
     eval_N_HypergeometricPQF,
 )
 
@@ -202,6 +201,47 @@ def eval(self, a, b, c, z, evaluation: Evaluation):
         return eval_Hypergeometric2F1(a, b, c, z)
 
 
+class HypergeometricU(MPMathFunction):
+    """
+    <url>
+    :Confluent hypergeometric function: https://en.wikipedia.org/wiki/Confluent_hypergeometric_function</url> (<url>
+    :mpmath: https://mpmath.org/doc/current/functions/bessel.html#mpmath.hyperu</url>, <url>
+    :WMA: https://reference.wolfram.com/language/ref/HypergeometricU.html</url>)
+    <dl>
+      <dt>'HypergeometricU'[$a$, $b$, $z$]
+      <dd>returns $U(a, b, z)$.
+    </dl>
+    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 = ""
+
+
 class MeijerG(MPMathFunction):
     """
     <url>
@@ -232,15 +272,7 @@ 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)
 
     def eval_N(self, a, b, z, prec, evaluation: Evaluation):
         "N[MeijerG[a_, b_, z_], prec_]"
@@ -257,44 +289,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):
-    """
-    <url>
-    :Confluent hypergeometric function: https://en.wikipedia.org/wiki/Confluent_hypergeometric_function</url> (<url>
-    :mpmath: https://mpmath.org/doc/current/functions/bessel.html#mpmath.hyperu</url>, <url>
-    :WMA: https://reference.wolfram.com/language/ref/HypergeometricU.html</url>)
-    <dl>
-      <dt>'HypergeometricU'[$a$, $b$, $z$]
-      <dd>returns $U(a, b, z)$.
-    </dl>
-    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/eval/specialfns/hypergeom.py b/mathics/eval/specialfns/hypergeom.py
index dd33a17c1..86b911066 100644
--- a/mathics/eval/specialfns/hypergeom.py
+++ b/mathics/eval/specialfns/hypergeom.py
@@ -2,7 +2,7 @@
 Evaluation functions for built-in Hypergeometric functions
 """
 
-from typing import Sequence, cast
+from typing import Sequence, Tuple, cast
 
 import mpmath
 import sympy
@@ -27,23 +27,7 @@ def eval_Hypergeometric1F1(a, b, z):
     """
 
     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_args, all_numeric = to_sympy_with_classification(args)
     sympy_result = tracing.run_sympy(
         sympy.hyper, [sympy_args[0]], [sympy_args[1]], sympy_args[2]
     )
@@ -79,20 +63,7 @@ def eval_Hypergeometric2F1(a, b, c, z):
     """
 
     args = (a, b, c, z)
-    sympy_args = []
-    all_numeric = True
-    for arg in args:
-        if isinstance(arg, Number):
-            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_args, all_numeric = to_sympy_with_classification(args)
     if all_numeric:
         args = cast(Sequence[Number], args)
 
@@ -138,7 +109,66 @@ def eval_N_HypergeometricPQF(a, b, z):
         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(a, b, z):
     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

From 1c8422825678ca1bc396b72fbf374cb7eef449d2 Mon Sep 17 00:00:00 2001
From: rocky <rb@dustyfeet.com>
Date: Sun, 27 Jul 2025 10:28:15 -0400
Subject: [PATCH 04/17] Add docstring to a new function and update CHANGES

---
 CHANGES.rst                          | 13 +++++++++++++
 mathics/eval/specialfns/hypergeom.py | 13 +++++++------
 2 files changed, 20 insertions(+), 6 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index fc4af20b8..e20db6344 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -18,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/eval/specialfns/hypergeom.py b/mathics/eval/specialfns/hypergeom.py
index 86b911066..bc83a68ea 100644
--- a/mathics/eval/specialfns/hypergeom.py
+++ b/mathics/eval/specialfns/hypergeom.py
@@ -33,7 +33,7 @@ def eval_Hypergeometric1F1(a, b, z):
     )
     expanded_result = sympy.hyperexpand(sympy_result)
 
-    # Oddly, expansion sometimes doesn't work for when complex arguments are given.
+    # 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:
@@ -58,8 +58,8 @@ def eval_Hypergeometric1F1(a, b, z):
 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.
+    Uses mpmath hyp2f1 if all args are numeric. Otherwise,
+    we use SymPy's hyper and expand the symbolic result.
     """
 
     args = (a, b, c, z)
@@ -80,7 +80,7 @@ def eval_Hypergeometric2F1(a, b, c, z):
             mpmath.hyp2f1, *cast(Sequence[Number], args), prec=prec
         )
     else:
-        return run_sympy_hyper(
+        return run_sympy_hyper_and_expand(
             [sympy_args[0], sympy_args[1]], [sympy_args[2]], sympy_args[3]
         )
 
@@ -90,7 +90,7 @@ def eval_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(a_sympy, b_sympy, z.to_sympy())
+        return run_sympy_hyper_and_expand(a_sympy, b_sympy, z.to_sympy())
     except Exception:
         return None
 
@@ -143,7 +143,8 @@ def eval_MeijerG(a, b, z):
 #         return None
 
 
-def run_sympy_hyper(a, b, z):
+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)

From 049d2262ce2cf1f7c509b111c2c6057bac12aa9e Mon Sep 17 00:00:00 2001
From: rocky <rb@dustyfeet.com>
Date: Mon, 28 Jul 2025 11:57:11 -0400
Subject: [PATCH 05/17] :sympy: -> :SymPy: in URL tag

---
 mathics/builtin/specialfns/hypergeom.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/mathics/builtin/specialfns/hypergeom.py b/mathics/builtin/specialfns/hypergeom.py
index 01a577497..4002f9fc6 100644
--- a/mathics/builtin/specialfns/hypergeom.py
+++ b/mathics/builtin/specialfns/hypergeom.py
@@ -33,7 +33,7 @@ class HypergeometricPFQ(MPMathFunction):
     <url>
     :Generalized hypergeometric function: https://en.wikipedia.org/wiki/Generalized_hypergeometric_function</url> (<url>
     :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#hyper</url>, <url>
-    :sympy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.hyper</url>, <url>
+    :SymPy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.hyper</url>, <url>
     :WMA: https://reference.wolfram.com/language/ref/HypergeometricPFQ.html</url>)
     <dl>
       <dt>'HypergeometricPFQ'[${a_1, ..., a_p}, {b_1, ..., b_q}, z$]
@@ -247,7 +247,7 @@ class MeijerG(MPMathFunction):
     <url>
     :Meijer G-function: https://en.wikipedia.org/wiki/Meijer_G-function</url> (<url>
     :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#meijerg</url>, <url>
-    :sympy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.meijerg</url>, <url>
+    :SymPy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.meijerg</url>, <url>
     :WMA: https://reference.wolfram.com/language/ref/MeijerG.html</url>)
     <dl>
       <dt>'MeijerG'[${{a_1, ..., a_n}, {a_{n+1}, ..., a_p}}, {{b_1, ..., b_m}, {b_{m+1}, ..., a_q}}, z$]

From 101ee0a5b92ef052c42724514a36627ae05ee68d Mon Sep 17 00:00:00 2001
From: rocky <rb@dustyfeet.com>
Date: Tue, 29 Jul 2025 08:43:43 -0400
Subject: [PATCH 06/17] Add special cases for hypergeometric fns...

based on Aravindh's review
---
 mathics/builtin/messages.py               |   1 +
 mathics/builtin/specialfns/hypergeom.py   | 183 ++++++++++++++--------
 mathics/core/atoms.py                     |   1 +
 test/builtin/specialfns/test_hypergeom.py | 120 ++++++++++++--
 4 files changed, 221 insertions(+), 84 deletions(-)

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 4002f9fc6..b33c0cdc6 100644
--- a/mathics/builtin/specialfns/hypergeom.py
+++ b/mathics/builtin/specialfns/hypergeom.py
@@ -9,6 +9,7 @@
 import mpmath
 
 import mathics.eval.tracing as tracing
+from mathics.core.atoms import Integer1, MachineReal1, Number
 from mathics.core.attributes import (
     A_LISTABLE,
     A_NUMERIC_FUNCTION,
@@ -18,6 +19,9 @@
 from mathics.core.builtin import MPMathFunction
 from mathics.core.convert.mpmath import from_mpmath
 from mathics.core.evaluation import Evaluation
+from mathics.core.expression import Expression
+from mathics.core.list import ListExpression
+from mathics.core.symbols import Symbol
 from mathics.core.systemsymbols import SymbolComplexInfinity, SymbolMachinePrecision
 from mathics.eval.specialfns.hypergeom import (
     eval_Hypergeometric1F1,
@@ -28,75 +32,6 @@
 )
 
 
-class HypergeometricPFQ(MPMathFunction):
-    """
-    <url>
-    :Generalized hypergeometric function: https://en.wikipedia.org/wiki/Generalized_hypergeometric_function</url> (<url>
-    :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#hyper</url>, <url>
-    :SymPy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.hyper</url>, <url>
-    :WMA: https://reference.wolfram.com/language/ref/HypergeometricPFQ.html</url>)
-    <dl>
-      <dt>'HypergeometricPFQ'[${a_1, ..., a_p}, {b_1, ..., b_q}, z$]
-      <dd>returns ${}_p F_q({a_1, ..., a_p}; {b_1, ..., b_q}; z)$.
-    </dl>
-    >> HypergeometricPFQ[{2}, {2}, 1]
-     = E
-
-    Result is symbollicaly simplified by default:
-    >> 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]
-     = E ^ z
-    >> HypergeometricPFQ[{0}, {b}, z]
-     = 1
-
-     >> 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}
-    summary_text = "compute the generalized hypergeometric function"
-    sympy_name = "hyper"
-
-    def eval(self, a, b, z, evaluation: Evaluation):
-        "HypergeometricPFQ[a_, b_, z_]"
-        return eval_HypergeometricPQF(a, b, z)
-
-    def eval_N(self, a, b, z, prec, evaluation: Evaluation):
-        "N[HypergeometricPFQ[a_, b_, z_], prec_]"
-        # FIXME: prec is not used. Why?
-        return eval_N_HypergeometricPQF(a, b, z)
-
-    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)
-
-
 class Hypergeometric1F1(MPMathFunction):
     """
     <url>
@@ -154,11 +89,30 @@ class Hypergeometric1F1(MPMathFunction):
     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_]"
+
+        # 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
+            )
+
         return eval_Hypergeometric1F1(a, b, z)
 
 
@@ -201,6 +155,97 @@ def eval(self, a, b, c, z, evaluation: Evaluation):
         return eval_Hypergeometric2F1(a, b, c, z)
 
 
+class HypergeometricPFQ(MPMathFunction):
+    """
+    <url>
+    :Generalized hypergeometric function: https://en.wikipedia.org/wiki/Generalized_hypergeometric_function</url> (<url>
+    :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#hyper</url>, <url>
+    :SymPy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.hyper</url>, <url>
+    :WMA: https://reference.wolfram.com/language/ref/HypergeometricPFQ.html</url>)
+    <dl>
+      <dt>'HypergeometricPFQ'[${a_1, ..., a_p}, {b_1, ..., b_q}, z$]
+      <dd>returns ${}_p F_q({a_1, ..., a_p}; {b_1, ..., b_q}; z)$.
+    </dl>
+    >> HypergeometricPFQ[{2}, {2}, 1]
+     = E
+
+    Result is symbollicaly simplified by default:
+    >> 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]
+     = E ^ z
+    >> HypergeometricPFQ[{0}, {b}, z]
+     = 1
+
+     >> 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}
+
+    summary_text = "compute the generalized hypergeometric function"
+    sympy_name = "hyper"
+
+    def eval(self, a, b, z, evaluation: Evaluation):
+        "HypergeometricPFQ[a_, b_, z_]"
+
+        # FIXME: a lot more checking could be done here.
+        if not (isinstance(a, ListExpression)):
+            evaluation.message(
+                "HypergeometricPQF",
+                "hdiv",
+                Expression(Symbol("Hypergeometric"), 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)
+
+        return eval_HypergeometricPQF(a, b, z)
+
+    def eval_N(self, a, b, z, prec, evaluation: Evaluation):
+        "N[HypergeometricPFQ[a_, b_, z_], prec_]"
+        # FIXME: prec is not used. It should be though.
+        return eval_N_HypergeometricPQF(a, b, z)
+
+
 class HypergeometricU(MPMathFunction):
     """
     <url>
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/test/builtin/specialfns/test_hypergeom.py b/test/builtin/specialfns/test_hypergeom.py
index 9f20b6329..8e408fea3 100644
--- a/test/builtin/specialfns/test_hypergeom.py
+++ b/test/builtin/specialfns/test_hypergeom.py
@@ -11,42 +11,132 @@
     ("str_expr", "msgs", "str_expected", "fail_msg"),
     [
         (
-            "N[HypergeometricPFQ[{4},{},1]]",
+            "Hypergeometric1F1[{3,5},{7,1},{4,2}]",
             None,
-            "ComplexInfinity",
+            "{-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),
         (
-            "MeijerG[{{},{}},{{0,0},{0,0}},100^2]",
-            None,
-            "MeijerG[{{}, {}}, {{0, 0}, {0, 0}}, 10000]",
+            "Hypergeometric1F1[0,0,z]",
             None,
+            "1",
+            "Special case: Integer 1 (SymPy may work differently)",
         ),
-        ("N[MeijerG[{{},{}},{{0,0},{0,0}},100^2]]", None, "0.000893912", None),
         (
-            "HypergeometricU[{3,1},{2,4},{7,8}]",
+            "Hypergeometric1F1[0.0,0,z]",
             None,
-            "{MeijerG[{{-2}, {}}, {{0, -1}, {}}, 7] / 2, HypergeometricU[1, 4, 8]}",
+            "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)",
         ),
-        ("N[HypergeometricU[{3,1},{2,4},{7,8}]]", None, "{0.00154364, 0.160156}", None),
-        ("HypergeometricU[0,c,z]", None, "1", None),
+    ],
+)
+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"),
+    [
         (
-            "Hypergeometric1F1[{3,5},{7,1},{4,2}]",
+            "N[HypergeometricPFQ[{4},{},1]]",
             None,
-            "{-1545 / 128 + 45 E ^ 4 / 128, 27 E ^ 2}",
+            "ComplexInfinity",
             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, "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),
+        (
+            "HypergeometricPFQ[{0},{0},2]",
+            None,
+            "1",
+            "Special case: using Integer 'z' parameter",
+        ),
+        (
+            "HypergeometricPFQ[{0},{0},3.0]",
+            None,
+            "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"),
+    [
+        (
+            "MeijerG[{{},{}},{{0,0},{0,0}},100^2]",
+            None,
+            "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,
+        ),
     ],
 )
-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,

From 40c310de52f3de970abe5467f8095b870476067b Mon Sep 17 00:00:00 2001
From: rocky <rb@dustyfeet.com>
Date: Tue, 29 Jul 2025 10:40:54 -0400
Subject: [PATCH 07/17] Move more eval code mathics.builtin to mathics.eval

---
 mathics/builtin/specialfns/hypergeom.py | 31 ++-----------------------
 mathics/eval/specialfns/hypergeom.py    | 31 ++++++++++++++++++++++++-
 2 files changed, 32 insertions(+), 30 deletions(-)

diff --git a/mathics/builtin/specialfns/hypergeom.py b/mathics/builtin/specialfns/hypergeom.py
index b33c0cdc6..3317d77c5 100644
--- a/mathics/builtin/specialfns/hypergeom.py
+++ b/mathics/builtin/specialfns/hypergeom.py
@@ -9,7 +9,6 @@
 import mpmath
 
 import mathics.eval.tracing as tracing
-from mathics.core.atoms import Integer1, MachineReal1, Number
 from mathics.core.attributes import (
     A_LISTABLE,
     A_NUMERIC_FUNCTION,
@@ -101,18 +100,6 @@ class Hypergeometric1F1(MPMathFunction):
     def eval(self, a, b, z, evaluation: Evaluation):
         "Hypergeometric1F1[a_, b_, z_]"
 
-        # 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
-            )
-
         return eval_Hypergeometric1F1(a, b, z)
 
 
@@ -222,22 +209,6 @@ def eval(self, a, b, z, evaluation: Evaluation):
                 Expression(Symbol("Hypergeometric"), 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)
-
         return eval_HypergeometricPQF(a, b, z)
 
     def eval_N(self, a, b, z, prec, evaluation: Evaluation):
@@ -319,6 +290,8 @@ def eval(self, a, b, z, evaluation: Evaluation):
         "MeijerG[a_, b_, z_]"
         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:
diff --git a/mathics/eval/specialfns/hypergeom.py b/mathics/eval/specialfns/hypergeom.py
index bc83a68ea..eff373b0c 100644
--- a/mathics/eval/specialfns/hypergeom.py
+++ b/mathics/eval/specialfns/hypergeom.py
@@ -8,7 +8,7 @@
 import sympy
 
 import mathics.eval.tracing as tracing
-from mathics.core.atoms import Complex, Number
+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
@@ -26,6 +26,18 @@ def eval_Hypergeometric1F1(a, b, z):
     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(
@@ -87,6 +99,23 @@ def eval_Hypergeometric2F1(a, b, c, z):
 
 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]

From 2d17cf7d48d96c356df59868bd5cd78a7e194ae6 Mon Sep 17 00:00:00 2001
From: "R. Bernstein" <rocky@users.noreply.github.com>
Date: Wed, 30 Jul 2025 18:02:56 -0400
Subject: [PATCH 08/17] TraceEvaluation can filter rewrites or evaluations...
 (#1436)

New options booleans ShowRewrites and ShowEvaluation can be added to
allow filtering TraceEvalation output to filter other either rewrites
rules or evaluation expressions. Presumably you don't want to filter
both.
---
 mathics/builtin/trace.py    | 66 +++++++++++++++++++++++++++++--------
 mathics/core/definitions.py |  3 ++
 mathics/eval/tracing.py     | 10 ++++--
 test/builtin/test_trace.py  |  4 +--
 4 files changed, 66 insertions(+), 17 deletions(-)

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):
     ## <url>:trace native symbol:</url>
 
     <dl>
-      <dt>'TraceEvaluation'[$expr$]
+      <dt>'TraceEvaluation'[$expr$, $options$]
       <dd>Evaluate $expr$ and print each step of the evaluation.
     </dl>
 
-    The 'ShowTimeBySteps' option prints the elapsed time before an evaluation occurs.
+    Options adjust output and filtering behavior
+    <dl>
+      <dt>'ShowTimeBySteps'
+      <dd>Print the elapsed time before an evaluation occurs. \
+      default is 'False'.
+      <dt>'ShowEvaluation'
+      <dd>Show evaluation calls and returns. The default is 'True'.
+      <dt>'ShowRewrite'
+      <dd>Show the effect of rewrite rules. The default is 'True'.
+    </dl>
+
+    <i>Note:</i> It does not make sense to set <i>both</i> '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/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/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/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()

From 13a2ccdd8197ae8f073df75e064b5c6941b68dd6 Mon Sep 17 00:00:00 2001
From: "R. Bernstein" <rocky@users.noreply.github.com>
Date: Thu, 31 Jul 2025 21:03:04 -0400
Subject: [PATCH 09/17] Add get_eval_Expression() (#1443)

Introduce `get_eval_Expression()`

`get_eval_Expression()` can be used in Mathics3 Builtin evaluation
methods to find the Expression that caused the evaluation to get
triggered.

mathics.file_io.files module converted to use this function to test
worthiness.
---
 mathics/builtin/files_io/files.py | 35 ++++++-------
 mathics/eval/stackframe.py        | 82 +++++++++++++++++++++++++++++++
 2 files changed, 97 insertions(+), 20 deletions(-)
 create mode 100644 mathics/eval/stackframe.py

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/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)

From c202f56cc1321f8ad5c4ba66831fad491c49ee8c Mon Sep 17 00:00:00 2001
From: rocky <rb@dustyfeet.com>
Date: Sat, 26 Jul 2025 17:24:00 -0400
Subject: [PATCH 10/17] First cut at Hypergeometric2F1

---
 CHANGES.rst                             |  1 +
 mathics/builtin/specialfns/hypergeom.py | 69 ++++++++++++++-------
 mathics/eval/specialfns/hypergeom.py    | 82 +++++++++++++++++++++++++
 3 files changed, 131 insertions(+), 21 deletions(-)
 create mode 100644 mathics/eval/specialfns/hypergeom.py

diff --git a/CHANGES.rst b/CHANGES.rst
index 780c19762..fc4af20b8 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -9,6 +9,7 @@ New Builtins
 
 * ``$SessionID``
 * ``BinaryReadList``
+* ``Hypergeometric2F1``
 
 Internals
 ---------
diff --git a/mathics/builtin/specialfns/hypergeom.py b/mathics/builtin/specialfns/hypergeom.py
index c83bd576f..4513e5423 100644
--- a/mathics/builtin/specialfns/hypergeom.py
+++ b/mathics/builtin/specialfns/hypergeom.py
@@ -12,7 +12,6 @@
 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,
@@ -21,8 +20,12 @@
 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.systemsymbols import SymbolComplexInfinity, SymbolMachinePrecision
+from mathics.eval.specialfns.hypergeom import (
+    eval_Hypergeometric2F1,
+    eval_HypergeometricPQF,
+    eval_N_HypergeometricPQF,
+)
 
 
 class HypergeometricPFQ(MPMathFunction):
@@ -70,27 +73,12 @@ class HypergeometricPFQ(MPMathFunction):
 
     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
+        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
+        # FIXME: prec is not used. Why?
+        return eval_N_HypergeometricPQF(a, b, z)
 
     def eval_numeric(self, a, b, z, evaluation: Evaluation):
         "HypergeometricPFQ[a:{__?NumericQ}, b:{__?NumericQ}, z_?MachineNumberQ]"
@@ -124,7 +112,7 @@ class Hypergeometric1F1(MPMathFunction):
     """
 
     attributes = A_LISTABLE | A_NUMERIC_FUNCTION | A_PROTECTED
-    mpmath_name = ""
+    mpmath_name = "hymp1f1"
     nargs = {3}
     rules = {
         "Hypergeometric1F1[a_, b_, z_]": "HypergeometricPFQ[{a},{b},z]",
@@ -133,6 +121,45 @@ class Hypergeometric1F1(MPMathFunction):
     sympy_name = ""
 
 
+class Hypergeometric2F1(MPMathFunction):
+    """
+    <url>
+    :Hypergeometric function: https://en.wikipedia.org/wiki/Hypergeometric_function</url> (<url>
+    :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#mpmath.hyp2f1</url>, <url>
+    :WMA: https://reference.wolfram.com/language/ref/Hypergeometric2F1.html</url>)
+    <dl>
+      <dt>'Hypergeometric2F1'[$a$, $b$, $c$, $z$]
+      <dd>returns ${}_2 F_1(a; b; c; z)$.
+    </dl>
+
+    >> 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 MeijerG(MPMathFunction):
     """
     <url>
diff --git a/mathics/eval/specialfns/hypergeom.py b/mathics/eval/specialfns/hypergeom.py
new file mode 100644
index 000000000..3499e3160
--- /dev/null
+++ b/mathics/eval/specialfns/hypergeom.py
@@ -0,0 +1,82 @@
+"""
+Evaluation functions for built-in Hypergeometric functions
+"""
+
+from typing import Sequence, cast
+
+import mpmath
+import sympy
+
+import mathics.eval.tracing as tracing
+from mathics.core.atoms import 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_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 = True
+    for arg in args:
+        if isinstance(arg, Number):
+            sympy_arg = arg.value
+        else:
+            sympy_arg = arg.to_sympy()
+            all_numeric = False
+
+        sympy_args.append(sympy_arg)
+
+    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:
+        sympy_result = tracing.run_sympy(
+            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)
+    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

From 35f07fe648934d12d58366bc5e3896740b8048af Mon Sep 17 00:00:00 2001
From: rocky <rb@dustyfeet.com>
Date: Sat, 26 Jul 2025 20:11:58 -0400
Subject: [PATCH 11/17] Move more fns to mpmath.eval.specialfns...

Also some results were corrected by making more use of SymPy's
hyperexpand function.
---
 mathics/builtin/specialfns/hypergeom.py   | 92 +++++++++++++++++------
 mathics/eval/specialfns/hypergeom.py      | 76 +++++++++++++++++--
 test/builtin/specialfns/test_hypergeom.py |  6 +-
 3 files changed, 139 insertions(+), 35 deletions(-)

diff --git a/mathics/builtin/specialfns/hypergeom.py b/mathics/builtin/specialfns/hypergeom.py
index 4513e5423..2d1433f9d 100644
--- a/mathics/builtin/specialfns/hypergeom.py
+++ b/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"
 
@@ -89,37 +102,66 @@ class Hypergeometric1F1(MPMathFunction):
     """
     <url>
     :Kummer confluent hypergeometric function: https://en.wikipedia.org/wiki/Confluent_hypergeometric_function</url> (<url>
-    :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#hyper</url>, <url>
+    :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#hyp1f1</url>, <url>
     :WMA: https://reference.wolfram.com/language/ref/Hypergeometric1F1.html</url>)
     <dl>
       <dt>'Hypergeometric1F1'[$a$, $b$, $z$]
       <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
+
+    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 = "hymp1f1"
+    mpmath_name = "hyp1f1"
     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):
     """
diff --git a/mathics/eval/specialfns/hypergeom.py b/mathics/eval/specialfns/hypergeom.py
index 3499e3160..dd33a17c1 100644
--- a/mathics/eval/specialfns/hypergeom.py
+++ b/mathics/eval/specialfns/hypergeom.py
@@ -8,7 +8,7 @@
 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
@@ -16,6 +16,61 @@
 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,10 +109,9 @@ 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):
@@ -62,8 +119,7 @@ def eval_HypergeometricPQF(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)
diff --git a/test/builtin/specialfns/test_hypergeom.py b/test/builtin/specialfns/test_hypergeom.py
index 7bdbc594c..9f20b6329 100644
--- a/test/builtin/specialfns/test_hypergeom.py
+++ b/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),
     ],

From d893f788cf8a962830eda1a82e1f5884aa243691 Mon Sep 17 00:00:00 2001
From: rocky <rb@dustyfeet.com>
Date: Sun, 27 Jul 2025 09:23:12 -0400
Subject: [PATCH 12/17] Move more functions to eval. Simplify

---
 mathics/builtin/specialfns/hypergeom.py | 95 +++++++++++--------------
 mathics/eval/specialfns/hypergeom.py    | 94 +++++++++++++++---------
 2 files changed, 105 insertions(+), 84 deletions(-)

diff --git a/mathics/builtin/specialfns/hypergeom.py b/mathics/builtin/specialfns/hypergeom.py
index 2d1433f9d..01a577497 100644
--- a/mathics/builtin/specialfns/hypergeom.py
+++ b/mathics/builtin/specialfns/hypergeom.py
@@ -7,7 +7,6 @@
 """
 
 import mpmath
-import sympy
 
 import mathics.eval.tracing as tracing
 from mathics.core.attributes import (
@@ -18,13 +17,13 @@
 )
 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.systemsymbols import SymbolComplexInfinity, SymbolMachinePrecision
 from mathics.eval.specialfns.hypergeom import (
     eval_Hypergeometric1F1,
     eval_Hypergeometric2F1,
     eval_HypergeometricPQF,
+    eval_MeijerG,
     eval_N_HypergeometricPQF,
 )
 
@@ -202,6 +201,47 @@ def eval(self, a, b, c, z, evaluation: Evaluation):
         return eval_Hypergeometric2F1(a, b, c, z)
 
 
+class HypergeometricU(MPMathFunction):
+    """
+    <url>
+    :Confluent hypergeometric function: https://en.wikipedia.org/wiki/Confluent_hypergeometric_function</url> (<url>
+    :mpmath: https://mpmath.org/doc/current/functions/bessel.html#mpmath.hyperu</url>, <url>
+    :WMA: https://reference.wolfram.com/language/ref/HypergeometricU.html</url>)
+    <dl>
+      <dt>'HypergeometricU'[$a$, $b$, $z$]
+      <dd>returns $U(a, b, z)$.
+    </dl>
+    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 = ""
+
+
 class MeijerG(MPMathFunction):
     """
     <url>
@@ -232,15 +272,7 @@ 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)
 
     def eval_N(self, a, b, z, prec, evaluation: Evaluation):
         "N[MeijerG[a_, b_, z_], prec_]"
@@ -257,44 +289,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):
-    """
-    <url>
-    :Confluent hypergeometric function: https://en.wikipedia.org/wiki/Confluent_hypergeometric_function</url> (<url>
-    :mpmath: https://mpmath.org/doc/current/functions/bessel.html#mpmath.hyperu</url>, <url>
-    :WMA: https://reference.wolfram.com/language/ref/HypergeometricU.html</url>)
-    <dl>
-      <dt>'HypergeometricU'[$a$, $b$, $z$]
-      <dd>returns $U(a, b, z)$.
-    </dl>
-    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/eval/specialfns/hypergeom.py b/mathics/eval/specialfns/hypergeom.py
index dd33a17c1..86b911066 100644
--- a/mathics/eval/specialfns/hypergeom.py
+++ b/mathics/eval/specialfns/hypergeom.py
@@ -2,7 +2,7 @@
 Evaluation functions for built-in Hypergeometric functions
 """
 
-from typing import Sequence, cast
+from typing import Sequence, Tuple, cast
 
 import mpmath
 import sympy
@@ -27,23 +27,7 @@ def eval_Hypergeometric1F1(a, b, z):
     """
 
     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_args, all_numeric = to_sympy_with_classification(args)
     sympy_result = tracing.run_sympy(
         sympy.hyper, [sympy_args[0]], [sympy_args[1]], sympy_args[2]
     )
@@ -79,20 +63,7 @@ def eval_Hypergeometric2F1(a, b, c, z):
     """
 
     args = (a, b, c, z)
-    sympy_args = []
-    all_numeric = True
-    for arg in args:
-        if isinstance(arg, Number):
-            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_args, all_numeric = to_sympy_with_classification(args)
     if all_numeric:
         args = cast(Sequence[Number], args)
 
@@ -138,7 +109,66 @@ def eval_N_HypergeometricPQF(a, b, z):
         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(a, b, z):
     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

From ac9d61a654d4ceaeeaff97cf7420a11ff079dd69 Mon Sep 17 00:00:00 2001
From: rocky <rb@dustyfeet.com>
Date: Sun, 27 Jul 2025 10:28:15 -0400
Subject: [PATCH 13/17] Add docstring to a new function and update CHANGES

---
 CHANGES.rst                          | 13 +++++++++++++
 mathics/eval/specialfns/hypergeom.py | 13 +++++++------
 2 files changed, 20 insertions(+), 6 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index fc4af20b8..e20db6344 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -18,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/eval/specialfns/hypergeom.py b/mathics/eval/specialfns/hypergeom.py
index 86b911066..bc83a68ea 100644
--- a/mathics/eval/specialfns/hypergeom.py
+++ b/mathics/eval/specialfns/hypergeom.py
@@ -33,7 +33,7 @@ def eval_Hypergeometric1F1(a, b, z):
     )
     expanded_result = sympy.hyperexpand(sympy_result)
 
-    # Oddly, expansion sometimes doesn't work for when complex arguments are given.
+    # 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:
@@ -58,8 +58,8 @@ def eval_Hypergeometric1F1(a, b, z):
 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.
+    Uses mpmath hyp2f1 if all args are numeric. Otherwise,
+    we use SymPy's hyper and expand the symbolic result.
     """
 
     args = (a, b, c, z)
@@ -80,7 +80,7 @@ def eval_Hypergeometric2F1(a, b, c, z):
             mpmath.hyp2f1, *cast(Sequence[Number], args), prec=prec
         )
     else:
-        return run_sympy_hyper(
+        return run_sympy_hyper_and_expand(
             [sympy_args[0], sympy_args[1]], [sympy_args[2]], sympy_args[3]
         )
 
@@ -90,7 +90,7 @@ def eval_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(a_sympy, b_sympy, z.to_sympy())
+        return run_sympy_hyper_and_expand(a_sympy, b_sympy, z.to_sympy())
     except Exception:
         return None
 
@@ -143,7 +143,8 @@ def eval_MeijerG(a, b, z):
 #         return None
 
 
-def run_sympy_hyper(a, b, z):
+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)

From 968fa187b3f9606fa0e99a6b073d0e2529e77a25 Mon Sep 17 00:00:00 2001
From: rocky <rb@dustyfeet.com>
Date: Mon, 28 Jul 2025 11:57:11 -0400
Subject: [PATCH 14/17] :sympy: -> :SymPy: in URL tag

---
 mathics/builtin/specialfns/hypergeom.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/mathics/builtin/specialfns/hypergeom.py b/mathics/builtin/specialfns/hypergeom.py
index 01a577497..4002f9fc6 100644
--- a/mathics/builtin/specialfns/hypergeom.py
+++ b/mathics/builtin/specialfns/hypergeom.py
@@ -33,7 +33,7 @@ class HypergeometricPFQ(MPMathFunction):
     <url>
     :Generalized hypergeometric function: https://en.wikipedia.org/wiki/Generalized_hypergeometric_function</url> (<url>
     :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#hyper</url>, <url>
-    :sympy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.hyper</url>, <url>
+    :SymPy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.hyper</url>, <url>
     :WMA: https://reference.wolfram.com/language/ref/HypergeometricPFQ.html</url>)
     <dl>
       <dt>'HypergeometricPFQ'[${a_1, ..., a_p}, {b_1, ..., b_q}, z$]
@@ -247,7 +247,7 @@ class MeijerG(MPMathFunction):
     <url>
     :Meijer G-function: https://en.wikipedia.org/wiki/Meijer_G-function</url> (<url>
     :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#meijerg</url>, <url>
-    :sympy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.meijerg</url>, <url>
+    :SymPy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.meijerg</url>, <url>
     :WMA: https://reference.wolfram.com/language/ref/MeijerG.html</url>)
     <dl>
       <dt>'MeijerG'[${{a_1, ..., a_n}, {a_{n+1}, ..., a_p}}, {{b_1, ..., b_m}, {b_{m+1}, ..., a_q}}, z$]

From 7d5355c8906a77967229b263b01bc7c97f160fe2 Mon Sep 17 00:00:00 2001
From: rocky <rb@dustyfeet.com>
Date: Tue, 29 Jul 2025 08:43:43 -0400
Subject: [PATCH 15/17] Add special cases for hypergeometric fns...

based on Aravindh's review
---
 mathics/builtin/messages.py               |   1 +
 mathics/builtin/specialfns/hypergeom.py   | 183 ++++++++++++++--------
 mathics/core/atoms.py                     |   1 +
 test/builtin/specialfns/test_hypergeom.py | 120 ++++++++++++--
 4 files changed, 221 insertions(+), 84 deletions(-)

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 4002f9fc6..b33c0cdc6 100644
--- a/mathics/builtin/specialfns/hypergeom.py
+++ b/mathics/builtin/specialfns/hypergeom.py
@@ -9,6 +9,7 @@
 import mpmath
 
 import mathics.eval.tracing as tracing
+from mathics.core.atoms import Integer1, MachineReal1, Number
 from mathics.core.attributes import (
     A_LISTABLE,
     A_NUMERIC_FUNCTION,
@@ -18,6 +19,9 @@
 from mathics.core.builtin import MPMathFunction
 from mathics.core.convert.mpmath import from_mpmath
 from mathics.core.evaluation import Evaluation
+from mathics.core.expression import Expression
+from mathics.core.list import ListExpression
+from mathics.core.symbols import Symbol
 from mathics.core.systemsymbols import SymbolComplexInfinity, SymbolMachinePrecision
 from mathics.eval.specialfns.hypergeom import (
     eval_Hypergeometric1F1,
@@ -28,75 +32,6 @@
 )
 
 
-class HypergeometricPFQ(MPMathFunction):
-    """
-    <url>
-    :Generalized hypergeometric function: https://en.wikipedia.org/wiki/Generalized_hypergeometric_function</url> (<url>
-    :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#hyper</url>, <url>
-    :SymPy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.hyper</url>, <url>
-    :WMA: https://reference.wolfram.com/language/ref/HypergeometricPFQ.html</url>)
-    <dl>
-      <dt>'HypergeometricPFQ'[${a_1, ..., a_p}, {b_1, ..., b_q}, z$]
-      <dd>returns ${}_p F_q({a_1, ..., a_p}; {b_1, ..., b_q}; z)$.
-    </dl>
-    >> HypergeometricPFQ[{2}, {2}, 1]
-     = E
-
-    Result is symbollicaly simplified by default:
-    >> 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]
-     = E ^ z
-    >> HypergeometricPFQ[{0}, {b}, z]
-     = 1
-
-     >> 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}
-    summary_text = "compute the generalized hypergeometric function"
-    sympy_name = "hyper"
-
-    def eval(self, a, b, z, evaluation: Evaluation):
-        "HypergeometricPFQ[a_, b_, z_]"
-        return eval_HypergeometricPQF(a, b, z)
-
-    def eval_N(self, a, b, z, prec, evaluation: Evaluation):
-        "N[HypergeometricPFQ[a_, b_, z_], prec_]"
-        # FIXME: prec is not used. Why?
-        return eval_N_HypergeometricPQF(a, b, z)
-
-    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)
-
-
 class Hypergeometric1F1(MPMathFunction):
     """
     <url>
@@ -154,11 +89,30 @@ class Hypergeometric1F1(MPMathFunction):
     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_]"
+
+        # 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
+            )
+
         return eval_Hypergeometric1F1(a, b, z)
 
 
@@ -201,6 +155,97 @@ def eval(self, a, b, c, z, evaluation: Evaluation):
         return eval_Hypergeometric2F1(a, b, c, z)
 
 
+class HypergeometricPFQ(MPMathFunction):
+    """
+    <url>
+    :Generalized hypergeometric function: https://en.wikipedia.org/wiki/Generalized_hypergeometric_function</url> (<url>
+    :mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#hyper</url>, <url>
+    :SymPy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.hyper</url>, <url>
+    :WMA: https://reference.wolfram.com/language/ref/HypergeometricPFQ.html</url>)
+    <dl>
+      <dt>'HypergeometricPFQ'[${a_1, ..., a_p}, {b_1, ..., b_q}, z$]
+      <dd>returns ${}_p F_q({a_1, ..., a_p}; {b_1, ..., b_q}; z)$.
+    </dl>
+    >> HypergeometricPFQ[{2}, {2}, 1]
+     = E
+
+    Result is symbollicaly simplified by default:
+    >> 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]
+     = E ^ z
+    >> HypergeometricPFQ[{0}, {b}, z]
+     = 1
+
+     >> 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}
+
+    summary_text = "compute the generalized hypergeometric function"
+    sympy_name = "hyper"
+
+    def eval(self, a, b, z, evaluation: Evaluation):
+        "HypergeometricPFQ[a_, b_, z_]"
+
+        # FIXME: a lot more checking could be done here.
+        if not (isinstance(a, ListExpression)):
+            evaluation.message(
+                "HypergeometricPQF",
+                "hdiv",
+                Expression(Symbol("Hypergeometric"), 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)
+
+        return eval_HypergeometricPQF(a, b, z)
+
+    def eval_N(self, a, b, z, prec, evaluation: Evaluation):
+        "N[HypergeometricPFQ[a_, b_, z_], prec_]"
+        # FIXME: prec is not used. It should be though.
+        return eval_N_HypergeometricPQF(a, b, z)
+
+
 class HypergeometricU(MPMathFunction):
     """
     <url>
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/test/builtin/specialfns/test_hypergeom.py b/test/builtin/specialfns/test_hypergeom.py
index 9f20b6329..8e408fea3 100644
--- a/test/builtin/specialfns/test_hypergeom.py
+++ b/test/builtin/specialfns/test_hypergeom.py
@@ -11,42 +11,132 @@
     ("str_expr", "msgs", "str_expected", "fail_msg"),
     [
         (
-            "N[HypergeometricPFQ[{4},{},1]]",
+            "Hypergeometric1F1[{3,5},{7,1},{4,2}]",
             None,
-            "ComplexInfinity",
+            "{-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),
         (
-            "MeijerG[{{},{}},{{0,0},{0,0}},100^2]",
-            None,
-            "MeijerG[{{}, {}}, {{0, 0}, {0, 0}}, 10000]",
+            "Hypergeometric1F1[0,0,z]",
             None,
+            "1",
+            "Special case: Integer 1 (SymPy may work differently)",
         ),
-        ("N[MeijerG[{{},{}},{{0,0},{0,0}},100^2]]", None, "0.000893912", None),
         (
-            "HypergeometricU[{3,1},{2,4},{7,8}]",
+            "Hypergeometric1F1[0.0,0,z]",
             None,
-            "{MeijerG[{{-2}, {}}, {{0, -1}, {}}, 7] / 2, HypergeometricU[1, 4, 8]}",
+            "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)",
         ),
-        ("N[HypergeometricU[{3,1},{2,4},{7,8}]]", None, "{0.00154364, 0.160156}", None),
-        ("HypergeometricU[0,c,z]", None, "1", None),
+    ],
+)
+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"),
+    [
         (
-            "Hypergeometric1F1[{3,5},{7,1},{4,2}]",
+            "N[HypergeometricPFQ[{4},{},1]]",
             None,
-            "{-1545 / 128 + 45 E ^ 4 / 128, 27 E ^ 2}",
+            "ComplexInfinity",
             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, "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),
+        (
+            "HypergeometricPFQ[{0},{0},2]",
+            None,
+            "1",
+            "Special case: using Integer 'z' parameter",
+        ),
+        (
+            "HypergeometricPFQ[{0},{0},3.0]",
+            None,
+            "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"),
+    [
+        (
+            "MeijerG[{{},{}},{{0,0},{0,0}},100^2]",
+            None,
+            "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,
+        ),
     ],
 )
-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,

From e0a1cb1b3400bf697fb12a56dea71713b49492f4 Mon Sep 17 00:00:00 2001
From: rocky <rb@dustyfeet.com>
Date: Tue, 29 Jul 2025 10:40:54 -0400
Subject: [PATCH 16/17] Move more eval code mathics.builtin to mathics.eval

---
 mathics/builtin/specialfns/hypergeom.py | 31 ++-----------------------
 mathics/eval/specialfns/hypergeom.py    | 31 ++++++++++++++++++++++++-
 2 files changed, 32 insertions(+), 30 deletions(-)

diff --git a/mathics/builtin/specialfns/hypergeom.py b/mathics/builtin/specialfns/hypergeom.py
index b33c0cdc6..3317d77c5 100644
--- a/mathics/builtin/specialfns/hypergeom.py
+++ b/mathics/builtin/specialfns/hypergeom.py
@@ -9,7 +9,6 @@
 import mpmath
 
 import mathics.eval.tracing as tracing
-from mathics.core.atoms import Integer1, MachineReal1, Number
 from mathics.core.attributes import (
     A_LISTABLE,
     A_NUMERIC_FUNCTION,
@@ -101,18 +100,6 @@ class Hypergeometric1F1(MPMathFunction):
     def eval(self, a, b, z, evaluation: Evaluation):
         "Hypergeometric1F1[a_, b_, z_]"
 
-        # 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
-            )
-
         return eval_Hypergeometric1F1(a, b, z)
 
 
@@ -222,22 +209,6 @@ def eval(self, a, b, z, evaluation: Evaluation):
                 Expression(Symbol("Hypergeometric"), 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)
-
         return eval_HypergeometricPQF(a, b, z)
 
     def eval_N(self, a, b, z, prec, evaluation: Evaluation):
@@ -319,6 +290,8 @@ def eval(self, a, b, z, evaluation: Evaluation):
         "MeijerG[a_, b_, z_]"
         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:
diff --git a/mathics/eval/specialfns/hypergeom.py b/mathics/eval/specialfns/hypergeom.py
index bc83a68ea..eff373b0c 100644
--- a/mathics/eval/specialfns/hypergeom.py
+++ b/mathics/eval/specialfns/hypergeom.py
@@ -8,7 +8,7 @@
 import sympy
 
 import mathics.eval.tracing as tracing
-from mathics.core.atoms import Complex, Number
+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
@@ -26,6 +26,18 @@ def eval_Hypergeometric1F1(a, b, z):
     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(
@@ -87,6 +99,23 @@ def eval_Hypergeometric2F1(a, b, c, z):
 
 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]

From 1273130bf0a9730309567b4a2bc49e56e5d48bde Mon Sep 17 00:00:00 2001
From: rocky <rb@dustyfeet.com>
Date: Thu, 31 Jul 2025 21:12:00 -0400
Subject: [PATCH 17/17] Use get_eval_Expression()

---
 mathics/builtin/specialfns/hypergeom.py | 9 ++-------
 1 file changed, 2 insertions(+), 7 deletions(-)

diff --git a/mathics/builtin/specialfns/hypergeom.py b/mathics/builtin/specialfns/hypergeom.py
index 3317d77c5..397ff4792 100644
--- a/mathics/builtin/specialfns/hypergeom.py
+++ b/mathics/builtin/specialfns/hypergeom.py
@@ -18,9 +18,7 @@
 from mathics.core.builtin import MPMathFunction
 from mathics.core.convert.mpmath import from_mpmath
 from mathics.core.evaluation import Evaluation
-from mathics.core.expression import Expression
 from mathics.core.list import ListExpression
-from mathics.core.symbols import Symbol
 from mathics.core.systemsymbols import SymbolComplexInfinity, SymbolMachinePrecision
 from mathics.eval.specialfns.hypergeom import (
     eval_Hypergeometric1F1,
@@ -29,6 +27,7 @@
     eval_MeijerG,
     eval_N_HypergeometricPQF,
 )
+from mathics.eval.stackframe import get_eval_Expression
 
 
 class Hypergeometric1F1(MPMathFunction):
@@ -203,11 +202,7 @@ def eval(self, a, b, z, evaluation: Evaluation):
 
         # FIXME: a lot more checking could be done here.
         if not (isinstance(a, ListExpression)):
-            evaluation.message(
-                "HypergeometricPQF",
-                "hdiv",
-                Expression(Symbol("Hypergeometric"), a, b, z),
-            )
+            evaluation.message("HypergeometricPFQ", "hdiv", get_eval_Expression())
 
         return eval_HypergeometricPQF(a, b, z)