Add Hypergeometric1F1 and Hypergeometric2F1

CHANGELOG.md

@@ -29,7 +29,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/).
   - Figurate numbers and predicates: `IsSquare`, `IsTriangular`, `IsPentagonal`, `IsOctahedral`, `IsCenteredSquare`, `IsPerfect`, `IsAbundant`, `IsHappy`
 
   Not included: `RandomPrime` (nondeterminism, same reasoning as `Random`/`RandomChoice`/`RandomSample`).
-- **Special functions:** `Gamma`, `GammaLn` (`math.gamma`/`math.lgamma`), `Beta`, `Factorial2`, and verified iterative algorithms `ErfInv`, `LambertW`, `AGM`, `EllipticK`, `EllipticE` (the "parameter" convention `K(m)`/`E(m)`, `m = k²`). Not included: `BesselJ`/`BesselY`/`BesselI`/`BesselK`, `AiryAi`/`AiryBi` (+ derivatives), `JacobiTheta`, `DedekindEta`, `Zeta`, `GammaRegularized`, `BetaRegularized` - these need `scipy`-grade numerics to not be subtly wrong across their domain.
+- **Special functions:** `Gamma`, `GammaLn` (`math.gamma`/`math.lgamma`), `Beta`, `Factorial2`, and verified iterative algorithms `ErfInv`, `LambertW`, `AGM`, `EllipticK`, `EllipticE` (the "parameter" convention `K(m)`/`E(m)`, `m = k²`), `Hypergeometric1F1`, `Hypergeometric2F1` (defining power series; `Hypergeometric1F1` capped at `|z| <= 500` and uses Kummer's transformation for `z < 0` to avoid catastrophic cancellation, `Hypergeometric2F1` restricted to `|z| < 1`, its actual radius of convergence). Not included: `BesselJ`/`BesselY`/`BesselI`/`BesselK`, `AiryAi`/`AiryBi` (+ derivatives), `JacobiTheta`, `DedekindEta`, `Zeta`, `GammaRegularized`, `BetaRegularized` - these need `scipy`-grade numerics to not be subtly wrong across their domain.
 - **Combinatorics:** `Choose` (alias for `Binomial`), `Fibonacci`, `Multinomial`, `Subfactorial`, `BellNumber`, and the enumeration functions `PowerSet`, `Permutations`, `Combinations`, `CartesianProduct` - the latter four have **no built-in output-size limit** (a set of `n` elements has `2^n` subsets, etc.); see the new `blacklist` parameter below for disabling them in deployments that accept untrusted expressions.
 - **Core (structural introspection):** `Head`/`Tail` (operator name / argument list of a compound expression, read from the raw unevaluated tree), `Hold` (returns an argument unevaluated), `Identity`, `Type` (runtime kind of an evaluated value), `IsSame`/`Same` (structural, pre-evaluation equality - distinct from `Equal`/`StrictEqual`/`IdenticallyEqual`, which all compare evaluated values). Not included: CAS functions (`Evaluate`, `Expand`, `Simplify`, `Solve`, ...), mutable-state functions (`Declare`, `Assign`, `Assume`, ...), and LaTeX serialization (`Parse`, `Latex`, `Subscript`, ...) - same reasoning as elsewhere in this solver (no CAS, no mutable state, no rendering surface). `Error`/`IsError` also not included - this solver already has an established, different error model (`MathJSONException`).
 - **`create_solver(parameters, blacklist=[...])`**: disables the named constructs for that solver instance, raising `MathJSONException` if an expression tries to use one, rather than evaluating it or silently ignoring it. An access-policy control, not a resource limiter - it doesn't make an enabled construct safe against pathological input (in particular, it's the recommended way to disable the four uncapped enumeration functions above for untrusted expression sources). A blacklisted name that isn't an actual construct has no effect (not validated against the known-construct list).

src/mathjson_solver/__main__.py

@@ -241,6 +241,8 @@ def translate_v1_mathjson(expr: MathJSONExpression) -> MathJSONExpression:
     AGM,
     EllipticK,
     EllipticE,
+    Hypergeometric1F1,
+    Hypergeometric2F1,
 )
 from ._statistics_constructs import (
     Variance,
@@ -609,6 +611,8 @@ def If(f, c, solver_parameters, s):
     "AGM": AGM,
     "EllipticK": EllipticK,
     "EllipticE": EllipticE,
+    "Hypergeometric1F1": Hypergeometric1F1,
+    "Hypergeometric2F1": Hypergeometric2F1,
     # --- Combinatorics (trivial subset) ---
     "Choose": lambda f, c, solver_parameters, s: math.comb(
         int(f(s[1], c)), int(f(s[2], c))
@@ -1043,6 +1047,8 @@ def extract_variables(
         "AGM",
         "EllipticK",
         "EllipticE",
+        "Hypergeometric1F1",
+        "Hypergeometric2F1",
         "Choose",
         "Fibonacci",
         "Multinomial",

src/mathjson_solver/_special_function_constructs.py

@@ -4,7 +4,14 @@
 """
 
 import math
-from ._special_functions import _agm, _elliptic_k_e, _erf_inv, _lambert_w
+from ._special_functions import (
+    _agm,
+    _elliptic_k_e,
+    _erf_inv,
+    _hyp1f1,
+    _hyp2f1,
+    _lambert_w,
+)
 
 # --- Special functions ---
 
@@ -41,3 +48,13 @@ def EllipticK(f, c, solver_parameters, s):
 
 def EllipticE(f, c, solver_parameters, s):
     return _elliptic_k_e(f(s[1], c))[1]
+
+
+def Hypergeometric1F1(f, c, solver_parameters, s):
+    """["Hypergeometric1F1", a, b, z] - confluent hypergeometric function 1F1(a;b;z)."""
+    return _hyp1f1(f(s[1], c), f(s[2], c), f(s[3], c))
+
+
+def Hypergeometric2F1(f, c, solver_parameters, s):
+    """["Hypergeometric2F1", a, b, c, z] - Gauss hypergeometric function 2F1(a,b;c;z)."""
+    return _hyp2f1(f(s[1], c), f(s[2], c), f(s[3], c), f(s[4], c))

src/mathjson_solver/_special_functions.py

@@ -1,12 +1,99 @@
 """
 Iterative special-function algorithms not available in the stdlib
-`math` module - backing ErfInv, LambertW, AGM, EllipticK, EllipticE.
-Each was verified against known reference values before shipping; see
-the docstrings below and the constructs' own docstrings in __main__.py.
+`math` module - backing ErfInv, LambertW, AGM, EllipticK, EllipticE,
+Hypergeometric1F1, Hypergeometric2F1. Each was verified against known
+reference values before shipping; see the docstrings below and the
+constructs' own docstrings in __main__.py.
 """
 
 import math
 
+# Hypergeometric1F1's defining series is entire (converges for any z),
+# but for large |z| the terms grow like e^|z| before the series
+# settles down, overflowing a float well before the true result would -
+# empirically, around |z| ~ 720. Capped well short of that with margin.
+_MAX_HYP1F1_Z = 500
+
+# Both hypergeometric series need more terms as their argument nears
+# the edge of where the naive term-by-term sum stays reliable (for
+# Hypergeometric2F1, its actual radius of convergence, |z| = 1) -
+# empirically, |z| = 0.99 needs ~2000 terms; a pathological or
+# malicious z closer still to that edge would need arbitrarily many,
+# so a cap is required for bounded latency, not just informational
+# accuracy - only a few ms even for the full 2000 (loop-bound, not
+# allocation-bound), verified empirically before choosing this number.
+_MAX_HYPERGEOMETRIC_TERMS = 2000
+
+
+def _hyp1f1_series(a, b, z, max_terms=_MAX_HYPERGEOMETRIC_TERMS, tol=1e-15):
+    total = 0.0
+    term = 1.0
+    for n in range(max_terms):
+        total += term
+        if abs(term) < tol * max(abs(total), 1):
+            return total
+        term *= (a + n) * z / ((b + n) * (n + 1))
+    raise ValueError(
+        f"'Hypergeometric1F1' did not converge within {max_terms} terms."
+    )
+
+
+def _hyp1f1(a, b, z):
+    """
+    Confluent hypergeometric function 1F1(a;b;z) (Kummer's function M),
+    via its defining power series. Verified against `mpmath.hyp1f1`
+    (temporary, verification-only install - never a project dependency)
+    across a range of a/b/z, matching to relative error < 1e-13, and
+    against CortexJS's own documented example
+    (`Hypergeometric1F1(1, 2, 2)` = 3.19452804946533).
+
+    For z < 0, Kummer's transformation
+    `1F1(a,b,z) = e^z * 1F1(b-a,b,-z)` is applied first - evaluating
+    the series directly at a negative z suffers catastrophic
+    cancellation (individually huge, alternating-sign terms summing to
+    something orders of magnitude smaller) that gets arbitrarily wrong
+    well within this function's other limits (e.g. a=10,b=1.5,z=-20 is
+    already off by a factor of ~1750 without the transformation) - the
+    transformed series always evaluates at a non-negative argument, so
+    it has no such cancellation (verified to match `mpmath` on exactly
+    the parameter combinations that broke the direct series).
+    """
+    if abs(z) > _MAX_HYP1F1_Z:
+        raise ValueError(
+            f"'Hypergeometric1F1' is limited to |z| <= {_MAX_HYP1F1_Z} "
+            f"(the naive series overflows a float, or needs excessive "
+            f"terms, beyond that)."
+        )
+    if z < 0:
+        return math.exp(z) * _hyp1f1_series(b - a, b, -z)
+    return _hyp1f1_series(a, b, z)
+
+
+def _hyp2f1(a, b, cc, z, max_terms=_MAX_HYPERGEOMETRIC_TERMS, tol=1e-15):
+    """
+    Gauss hypergeometric function 2F1(a,b;c;z), via its defining power
+    series - mathematically convergent only for |z| < 1. Verified
+    against `mpmath.hyp2f1` to relative error < 1e-10 even close to
+    that boundary (|z| = 0.99), and against CortexJS's own documented
+    example (`Hypergeometric2F1(1, 1, 2, 0.5)` = 1.38629436111989).
+    """
+    if not (-1 < z < 1):
+        raise ValueError(
+            "'Hypergeometric2F1' only supports |z| < 1 (the defining "
+            "power series' radius of convergence)."
+        )
+    total = 0.0
+    term = 1.0
+    for n in range(max_terms):
+        total += term
+        if abs(term) < tol * max(abs(total), 1):
+            return total
+        term *= (a + n) * (b + n) * z / ((cc + n) * (n + 1))
+    raise ValueError(
+        f"'Hypergeometric2F1' did not converge within {max_terms} terms "
+        f"(z is too close to the radius of convergence, |z| = 1)."
+    )
+
 
 def _erf_inv(x):
     """Inverse error function via Newton's method against math.erf."""

tests/test_cortexjs_special_functions.py

@@ -42,6 +42,30 @@
             ["Round", ["EllipticE", 0.5], 10],
             round(1.3506438810476755, 10),
         ),
+        # --- Hypergeometric: verified against CortexJS's own documented
+        # examples and against mpmath across several a/b/(c)/z combinations
+        # (see _special_functions.py's docstrings) ---
+        (
+            {},
+            ["Round", ["Hypergeometric1F1", 1, 2, 2], 10],
+            round(3.19452804946533, 10),
+        ),
+        (
+            {},
+            ["Round", ["Hypergeometric2F1", 1, 1, 2, 0.5], 10],
+            round(1.38629436111989, 10),
+        ),
+        # 1F1(a,a,z) = e^z (a closed-form identity, independent check)
+        ({}, ["Round", ["Hypergeometric1F1", 3, 3, 1], 10], round(math.e, 10)),
+        # Kummer's transformation path (z < 0) - without it, this
+        # particular a/b/z combination is off by a factor of ~1750.
+        (
+            {},
+            ["Round", ["Hypergeometric1F1", 10, 1.5, -20], 12],
+            round(-1.17307479e-06, 12),
+        ),
+        # 2F1(a,b;b;z) = (1-z)^-a (a closed-form identity, independent check)
+        ({}, ["Round", ["Hypergeometric2F1", 2, 3, 3, 0.3], 10], round((1 - 0.3) ** -2, 10)),
     ],
 )
 def test_special_functions(parameters, expression, expected_result):
@@ -69,3 +93,21 @@ def test_elliptic_functions_reject_out_of_domain():
         solver(["EllipticK", -0.5])
     with pytest.raises(Exception):
         solver(["EllipticK", 1.5])
+
+
+def test_hypergeometric1f1_rejects_out_of_domain():
+    solver = create_solver({})
+    with pytest.raises(Exception):
+        solver(["Hypergeometric1F1", 1, 2, 501])
+    with pytest.raises(Exception):
+        solver(["Hypergeometric1F1", 1, 2, -501])
+
+
+def test_hypergeometric2f1_rejects_out_of_domain():
+    solver = create_solver({})
+    with pytest.raises(Exception):
+        solver(["Hypergeometric2F1", 1, 1, 2, 1])
+    with pytest.raises(Exception):
+        solver(["Hypergeometric2F1", 1, 1, 2, -1])
+    with pytest.raises(Exception):
+        solver(["Hypergeometric2F1", 1, 1, 2, 1.5])