Add scipy-gated special functions

.github/workflows/python-app.yml

@@ -20,7 +20,7 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        pip install flake8 pytest google-re2 numpy
+        pip install flake8 pytest google-re2 numpy scipy
         if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
     - name: Lint with flake8
       run: |

.github/workflows/test-coverage.yml

@@ -21,7 +21,7 @@ jobs:
 
       - name: Install Dependencies
         run: |
-          pip install pytest pytest-cov google-re2 numpy
+          pip install pytest pytest-cov google-re2 numpy scipy
 
       - name: Run pytest, then run again with Coverage
         run: |

CHANGELOG.md

@@ -29,7 +29,8 @@ 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²`), `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.
+- **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: `JacobiTheta`, `DedekindEta` - `scipy.special` (see below) doesn't implement these either, so they'd need a from-scratch q-series with the same "subtly wrong across the domain" risk as any other naive truncated series.
+- **`BesselJ`/`BesselY`/`BesselI`/`BesselK`, `AiryAi`/`AiryBi`/`AiryAiPrime`/`AiryBiPrime`, `Zeta`, `GammaRegularized`, `BetaRegularized`** - thin wrappers around `scipy.special`, behind a new optional `special-functions` extra (`pip install mathjson-solver[special-functions]`), same pattern as `integration`/`regex`. `GammaRegularized(a, z)` is the *upper* regularized incomplete gamma (`scipy`'s own `gammainc` is the lower form; this uses `gammaincc`); `BetaRegularized(x, a, b)` puts `x` first, unlike `scipy.special.betainc`'s `(a, b, x)`. Without the extra installed, these raise a clear `ImportError` naming the missing package.
 - **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).

README.md

@@ -34,6 +34,7 @@ pip install mathjson-solver
 **Optional:**
 - `pip install mathjson-solver[integration]` (installs numpy; only required for `TrapezoidalIntegrate`)
 - `pip install mathjson-solver[regex]` (installs [google-re2](https://pypi.org/project/google-re2/); only required for `RegExp`/`IsMatch`/`StringMatch`/`StringMatchAll`)
+- `pip install mathjson-solver[special-functions]` (installs [scipy](https://scipy.org/); only required for `BesselJ`/`BesselY`/`BesselI`/`BesselK`, `AiryAi`/`AiryBi`/`AiryAiPrime`/`AiryBiPrime`, `Zeta`, `GammaRegularized`, `BetaRegularized`)
 
 ## Quick Start
 
@@ -128,7 +129,8 @@ The library supports a comprehensive set of mathematical operations:
 * **Pattern Matching (requires the optional `regex` extra, `pip install mathjson-solver[regex]`):** RegExp, IsMatch, StringMatch, StringMatchAll — backed by [RE2](https://github.com/google/re2) rather than Python's `re`, for a hard guarantee against catastrophic backtracking (no backreferences/lookaround, as a deliberate trade-off for that guarantee)
 * **Date/Time:** Strptime, Strftime, Today, Now, TimeDelta functions (Weeks, Days, Hours, Minutes)
 * **Number Theory:** Chop, Mod, Clamp, GCD, LCM, Factorial, Binomial, IsPrime, Erf, Erfc, Rational, Numerator, Denominator, MachineEpsilon, CatalanConstant, EulerGamma, PowerMod, ModularInverse, IntegerSqrt, FactorInteger, PrimeFactors, PrimeNu, PrimeOmega, Radical, IsSquareFree, Divisors, Sigma0, Sigma1, SigmaMinus1, DivisorSigma, Divides, Totient, IsPerfectPower, NthPrime, NextPrime, PrimePi, ExtendedGCD, ChineseRemainder, CarmichaelLambda, JacobiSymbol, LegendreSymbol, MultiplicativeOrder, PrimitiveRoot, LucasL, CatalanNumber, BernoulliB, ContinuedFraction, FromContinuedFraction, IntegerDigits, DigitCount, DigitSum, FromDigits, IsSquare, IsTriangular, IsPentagonal, IsOctahedral, IsCenteredSquare, IsPerfect, IsAbundant, IsHappy
-* **Special Functions:** Gamma, GammaLn, Beta, Factorial2, ErfInv, LambertW, AGM, EllipticK, EllipticE
+* **Special Functions:** Gamma, GammaLn, Beta, Factorial2, ErfInv, LambertW, AGM, EllipticK, EllipticE, Hypergeometric1F1, Hypergeometric2F1
+* **Special Functions (requires the optional `special-functions` extra, `pip install mathjson-solver[special-functions]`):** BesselJ, BesselY, BesselI, BesselK, AiryAi, AiryBi, AiryAiPrime, AiryBiPrime, Zeta, GammaRegularized, BetaRegularized — backed by [scipy](https://scipy.org/)
 * **Combinatorics:** Choose, Fibonacci, Multinomial, Subfactorial, BellNumber, PowerSet, Permutations, Combinations, CartesianProduct (the last four have no output-size limit - see [Restricting Available Functions](#restricting-available-functions))
 * **Core (structural introspection):** Head, Tail, Hold, Identity, Type, IsSame, Same
 * **Integration (requires the optional `integration` extra, `pip install mathjson-solver[integration]`):** TrapezoidalIntegrate. Also in this group but with no extra dependency: Interp, FindIntervalIndex, Variable references

docs/README.md

@@ -460,15 +460,55 @@ defined for `x ≥ -1/e`.
 ["Round", ["EllipticE", 0.5], 5]   # 1.35064
 ```
 
+### Hypergeometric1F1, Hypergeometric2F1
+Confluent (`1F1`) and Gauss (`2F1`) hypergeometric functions, via their
+defining power series - verified against `mpmath` and CortexJS's own
+documented examples before shipping. `Hypergeometric1F1` is capped at
+`|z| ≤ 500` (the naive series overflows a float beyond that) and uses
+Kummer's transformation for `z < 0` to avoid catastrophic cancellation
+that would otherwise silently corrupt the result well within that
+range. `Hypergeometric2F1` only supports `|z| < 1`, its actual
+mathematical radius of convergence.
+
+```python
+["Round", ["Hypergeometric1F1", 1, 2, 2], 5]        # 3.19453
+["Round", ["Hypergeometric2F1", 1, 1, 2, 0.5], 5]   # 1.38629
+```
+
+### BesselJ, BesselY, BesselI, BesselK, AiryAi, AiryBi, AiryAiPrime, AiryBiPrime, Zeta, GammaRegularized, BetaRegularized
+Require the optional `special-functions` extra
+(`pip install mathjson-solver[special-functions]`), which installs
+[scipy](https://scipy.org/) - these need real numerical algorithms
+(stable series/asymptotic-expansion switching depending on argument
+regime) that this project doesn't reimplement from scratch, unlike the
+verified-in-house algorithms above. Thin wrappers around
+`scipy.special`, verified against CortexJS's own documented examples.
+Two argument-convention gotchas worth knowing, since they don't match
+`scipy`'s own conventions directly:
+- `GammaRegularized(a, z)` is the *upper* regularized incomplete gamma,
+  `Q(a, z) = Γ(a, z) / Γ(a)` - `scipy.special.gammainc` computes the
+  *lower* form (`P`); this uses `gammaincc` instead.
+- `BetaRegularized(x, a, b)` puts `x` first, while
+  `scipy.special.betainc` takes `(a, b, x)`.
+
+```python
+["Round", ["BesselJ", 0, 1], 4]              # 0.7652
+["Round", ["AiryAi", 0], 4]                  # 0.355
+["Round", ["Zeta", 2], 6]                    # 1.644934  (π²/6)
+["Round", ["GammaRegularized", 3, 5], 5]     # 0.12465
+["Round", ["BetaRegularized", 0.5, 2, 3], 4] # 0.6875
+```
+
+Without the extra installed, these constructs raise a clear
+`ImportError` (naming the missing package) rather than failing with an
+unrelated `ModuleNotFoundError` deep in the call stack.
+
 ### Deliberately excluded
-`BesselJ`/`BesselY`/`BesselI`/`BesselK`, `AiryAi`/`AiryBi` (and their
-derivatives), `JacobiTheta`, `DedekindEta`, `Zeta`, `GammaRegularized`,
-`BetaRegularized` - these need real numerical algorithms (stable
-series/asymptotic-expansion switching depending on argument regime) to
-not be subtly wrong across their domain. Doing them properly means
-depending on `scipy`, a materially heavier dependency than `numpy`/
-`google-re2` - a separate dependency conversation, not a default
-inclusion.
+`JacobiTheta`, `DedekindEta` - `scipy.special` doesn't implement these
+either (no theta-function support), so the `special-functions` extra
+above doesn't unlock them. A from-scratch q-series implementation would
+carry the same "naive truncated series can be subtly wrong" risk noted
+elsewhere in this section - a separate decision, not bundled in here.
 
 ---
 
@@ -2090,6 +2130,8 @@ Evaluating a `["Function", ...]` expression outside of such a context (i.e. not
 ### Special Functions
 - [Gamma, GammaLn, Beta, Factorial2](#gamma-gammaln-beta-factorial2) - Gamma-based functions
 - [ErfInv, LambertW, AGM, EllipticK, EllipticE](#erfinv-lambertw-agm-elliptick-elliptice) - Verified iterative algorithms
+- [Hypergeometric1F1, Hypergeometric2F1](#hypergeometric1f1-hypergeometric2f1) - Confluent and Gauss hypergeometric functions
+- [BesselJ, BesselY, BesselI, BesselK, AiryAi, AiryBi, AiryAiPrime, AiryBiPrime, Zeta, GammaRegularized, BetaRegularized](#besselj-bessely-besseli-besselk-airyai-airybi-airyaiprime-airybiprime-zeta-gammaregularized-betaregularized) - Requires the optional `special-functions` extra (scipy)
 
 ### Combinatorics
 - [Choose, Fibonacci, Multinomial, Subfactorial, BellNumber](#combinatorics) - Counting functions with no explosion risk

pyproject.toml

@@ -26,6 +26,11 @@ integration = ["numpy"]
 # solver that evaluates untrusted expressions with no execution budget
 # of its own. Optional, same as `integration` above.
 regex = ["google-re2"]
+# Bessel/Airy/Zeta/GammaRegularized/BetaRegularized need scipy.special -
+# real numerical algorithms (stable series/asymptotic-expansion
+# switching depending on argument regime), not something to reimplement
+# from scratch. Optional, same as `integration`/`regex` above.
+special-functions = ["scipy"]
 
 [project.urls]
 Homepage = "https://github.com/LongenesisLtd/mathjson-solver"

src/mathjson_solver/__main__.py

@@ -243,6 +243,17 @@ def translate_v1_mathjson(expr: MathJSONExpression) -> MathJSONExpression:
     EllipticE,
     Hypergeometric1F1,
     Hypergeometric2F1,
+    BesselJ,
+    BesselY,
+    BesselI,
+    BesselK,
+    AiryAi,
+    AiryBi,
+    AiryAiPrime,
+    AiryBiPrime,
+    Zeta,
+    GammaRegularized,
+    BetaRegularized,
 )
 from ._statistics_constructs import (
     Variance,
@@ -613,6 +624,17 @@ def If(f, c, solver_parameters, s):
     "EllipticE": EllipticE,
     "Hypergeometric1F1": Hypergeometric1F1,
     "Hypergeometric2F1": Hypergeometric2F1,
+    "BesselJ": BesselJ,
+    "BesselY": BesselY,
+    "BesselI": BesselI,
+    "BesselK": BesselK,
+    "AiryAi": AiryAi,
+    "AiryBi": AiryBi,
+    "AiryAiPrime": AiryAiPrime,
+    "AiryBiPrime": AiryBiPrime,
+    "Zeta": Zeta,
+    "GammaRegularized": GammaRegularized,
+    "BetaRegularized": BetaRegularized,
     # --- Combinatorics (trivial subset) ---
     "Choose": lambda f, c, solver_parameters, s: math.comb(
         int(f(s[1], c)), int(f(s[2], c))
@@ -1049,6 +1071,17 @@ def extract_variables(
         "EllipticE",
         "Hypergeometric1F1",
         "Hypergeometric2F1",
+        "BesselJ",
+        "BesselY",
+        "BesselI",
+        "BesselK",
+        "AiryAi",
+        "AiryBi",
+        "AiryAiPrime",
+        "AiryBiPrime",
+        "Zeta",
+        "GammaRegularized",
+        "BetaRegularized",
         "Choose",
         "Fibonacci",
         "Multinomial",

src/mathjson_solver/_special_function_constructs.py

@@ -1,6 +1,6 @@
 """
 Special-function constructs, backed by the iterative algorithms in
-_special_functions.py.
+_special_functions.py, plus (where available) scipy.special.
 """
 
 import math
@@ -13,6 +13,24 @@
     _lambert_w,
 )
 
+SCIPY_AVAILABLE = False
+try:
+    from scipy import special as _scipy_special
+
+    SCIPY_AVAILABLE = True
+except ImportError:
+    pass
+
+
+def _require_scipy(construct_name):
+    if not SCIPY_AVAILABLE:
+        raise ImportError(
+            f"'{construct_name}' requires 'scipy'. Install with "
+            f"'pip install scipy' (or 'pip install "
+            f"mathjson-solver[special-functions]')."
+        )
+
+
 # --- Special functions ---
 
 
@@ -58,3 +76,91 @@ def Hypergeometric1F1(f, c, solver_parameters, s):
 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))
+
+
+# --- scipy-gated: Bessel, Airy, Zeta, regularized gamma/beta ---
+#
+# Thin wrappers around scipy.special - unlike ErfInv/LambertW/AGM/
+# EllipticK/EllipticE/Hypergeometric1F1/Hypergeometric2F1 above, there's
+# no algorithm of this project's own to verify here, just an argument-
+# convention mapping (verified against CortexJS's own documented
+# examples below, each matching to full precision).
+
+
+def BesselJ(f, c, solver_parameters, s):
+    """["BesselJ", n, x] - Bessel function of the first kind."""
+    _require_scipy("BesselJ")
+    return _scipy_special.jv(f(s[1], c), f(s[2], c))
+
+
+def BesselY(f, c, solver_parameters, s):
+    """["BesselY", n, x] - Bessel function of the second kind."""
+    _require_scipy("BesselY")
+    return _scipy_special.yv(f(s[1], c), f(s[2], c))
+
+
+def BesselI(f, c, solver_parameters, s):
+    """["BesselI", n, x] - modified Bessel function of the first kind."""
+    _require_scipy("BesselI")
+    return _scipy_special.iv(f(s[1], c), f(s[2], c))
+
+
+def BesselK(f, c, solver_parameters, s):
+    """["BesselK", n, x] - modified Bessel function of the second kind."""
+    _require_scipy("BesselK")
+    return _scipy_special.kv(f(s[1], c), f(s[2], c))
+
+
+def AiryAi(f, c, solver_parameters, s):
+    """["AiryAi", x] - Airy function Ai(x)."""
+    _require_scipy("AiryAi")
+    return _scipy_special.airy(f(s[1], c))[0]
+
+
+def AiryBi(f, c, solver_parameters, s):
+    """["AiryBi", x] - Airy function Bi(x)."""
+    _require_scipy("AiryBi")
+    return _scipy_special.airy(f(s[1], c))[2]
+
+
+def AiryAiPrime(f, c, solver_parameters, s):
+    """["AiryAiPrime", x] - derivative of the Airy function Ai(x)."""
+    _require_scipy("AiryAiPrime")
+    return _scipy_special.airy(f(s[1], c))[1]
+
+
+def AiryBiPrime(f, c, solver_parameters, s):
+    """["AiryBiPrime", x] - derivative of the Airy function Bi(x)."""
+    _require_scipy("AiryBiPrime")
+    return _scipy_special.airy(f(s[1], c))[3]
+
+
+def Zeta(f, c, solver_parameters, s):
+    """["Zeta", s] - the Riemann zeta function."""
+    _require_scipy("Zeta")
+    return _scipy_special.zeta(f(s[1], c), 1)
+
+
+def GammaRegularized(f, c, solver_parameters, s):
+    """
+    ["GammaRegularized", a, z] - the *upper* regularized incomplete
+    gamma function, Q(a, z) = Γ(a, z) / Γ(a) - CortexJS's own
+    definition (not the lower one, P(a, z), that a function named
+    "regularized incomplete gamma" more often defaults to elsewhere;
+    verified against CortexJS's own documented example,
+    `GammaRegularized(3, 5)` = 0.12465201948308115 - `scipy`'s
+    `gammainc` computes the *lower* form, `gammaincc` the upper one
+    needed here).
+    """
+    _require_scipy("GammaRegularized")
+    return _scipy_special.gammaincc(f(s[1], c), f(s[2], c))
+
+
+def BetaRegularized(f, c, solver_parameters, s):
+    """
+    ["BetaRegularized", x, a, b] - the regularized incomplete beta
+    function, I_x(a, b). Note the argument order: CortexJS puts `x`
+    first, while `scipy.special.betainc` takes `(a, b, x)`.
+    """
+    _require_scipy("BetaRegularized")
+    return _scipy_special.betainc(f(s[2], c), f(s[3], c), f(s[1], c))