diff --git a/CHANGELOG.md b/CHANGELOG.md index 1388518..8993daf 100644 --- a/CHANGELOG.md +++ b/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). diff --git a/src/mathjson_solver/__main__.py b/src/mathjson_solver/__main__.py index f09d919..f1c15a2 100644 --- a/src/mathjson_solver/__main__.py +++ b/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", diff --git a/src/mathjson_solver/_special_function_constructs.py b/src/mathjson_solver/_special_function_constructs.py index d3cc988..600ea36 100644 --- a/src/mathjson_solver/_special_function_constructs.py +++ b/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)) diff --git a/src/mathjson_solver/_special_functions.py b/src/mathjson_solver/_special_functions.py index c254520..3f391ba 100644 --- a/src/mathjson_solver/_special_functions.py +++ b/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.""" diff --git a/tests/test_cortexjs_special_functions.py b/tests/test_cortexjs_special_functions.py index b00d849..e7a0261 100644 --- a/tests/test_cortexjs_special_functions.py +++ b/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])