Bessel, Airy, zeta and polylogarithm functions

Bessel function of the first kind J_nu(x) (DLMF 10.2; SymPy's besselj).

A half-integer order (Fraction(1, 2), pool.rational(3, 2)) with expand=True (the default) gives the elementary form — bessel_j(1/2, x) = sqrt(2)*sin(x)/(sqrt(pi)*sqrt(x)). Otherwise the node bessel_j(nu, x), canonicalised: bessel_j(0, x) is bessel_j0(x), bessel_j(1, x) is bessel_j1(x), J_{-n} = (-1)^n J_n, J_n(-x) = (-1)^n J_n(x), J_0(0) = 1 and J_nu(0) = 0 for nu > 0.

The node differentiates in x (d/dx J_nu = (J_{nu-1} - J_{nu+1})/2; a derivative in nu is refused with E-DIFF-001), evaluates in f64, in rigorous balls (interval_eval) and at complex points (evaluate(e, {x: 1+2j})), through Arb.

Raises SpecialFunctionError (E-SPECFN-004) for an expansion past order 200 — keep it symbolic with expand=False.

>>> import alkahest as ak
>>> from fractions import Fraction
>>> p = ak.ExprPool(); x = p.symbol("x")
>>> ak.bessel_j(Fraction(1, 2), x)     # sqrt(2)*pi^(-1/2)*x^(-1/2)*sin(x)

Bessel function of the second kind Y_nu(x) (DLMF 10.2; SymPy's bessely). Half-integer orders expand (Y_{1/2}(x) = -sqrt(2)*cos(x)/(sqrt(pi)*sqrt(x))); Y_{-n} = (-1)^n Y_n`. Real only for `x > 0`: `f64` evaluation at `x <= 0 declines.

Modified Bessel function of the first kind I_nu(x) (DLMF 10.25; SymPy's besseli). Half-integer orders expand into sinh/cosh; I_{-n} = I_n; d/dx I_nu = (I_{nu-1} + I_{nu+1})/2.

Modified Bessel function of the second kind K_nu(x) (DLMF 10.25; SymPy's besselk). Half-integer orders expand (K_{1/2}(x) = sqrt(pi/(2x))*exp(-x)); K_{-nu} = K_nu; d/dx K_nu = -(K_{nu-1} + K_{nu+1})/2.

Spherical Bessel function of the first kind j_n(x) = sqrt(pi/(2x)) J_{n+1/2}(x) (DLMF 10.47; SymPy's jn).

An integer n expands (j_1(x) = sin(x)/x**2 - cos(x)/x); a symbolic n keeps the node, with d/dx j_n = -j_{n+1} + (n/x) j_n. A non-integer literal order raises E-SPECFN-007.

Spherical Bessel function of the second kind y_n(x) = sqrt(pi/(2x)) Y_{n+1/2}(x) (DLMF 10.47; SymPy's yn). y_0(x) = -cos(x)/x.

Hankel function of the first kind H1_nu(x) = J_nu(x) + i*Y_nu(x) (DLMF 10.4; SymPy's hankel1, mpmath's hankel1).

Complex-valued. A half-integer order with expand=True (the default) gives the elementary form — hankel1(1/2, x) = -I*sqrt(2)*exp(I*x)/(sqrt(pi)*sqrt(x)) — and the reflection H1_{-nu} = exp(I*pi*nu)*H1_nu is applied exactly; H1_{-n} = (-1)**n H1_n for an integer order. The node differentiates in x ((H1_{nu-1} - H1_{nu+1})/2); a derivative in nu is refused (E-DIFF-001). There is no real f64 / ball value (the value is complex): evaluate(e, {x: 2.0}) returns the complex number, through Arb's principal branch.

Raises SpecialFunctionError (E-SPECFN-004) for an expansion past order 200 (use expand=False).

>>> import alkahest as ak
>>> from fractions import Fraction
>>> p = ak.ExprPool(); x = p.symbol("x")
>>> ak.hankel1(Fraction(1, 2), x)      # -I*sqrt(2)*pi^(-1/2)*x^(-1/2)*exp(I*x)

Hankel function of the second kind H2_nu(x) = J_nu(x) - i*Y_nu(x) (DLMF 10.4); see hankel1. For real x it is the complex conjugate of hankel1 at the same real order.

Airy function Ai(x) (DLMF 9.2; SymPy's airyai). d/dx Ai = airy_ai_prime`; `Ai(0) = 3**(-2/3)/gamma(2/3).

Airy function Bi(x) (DLMF 9.2; SymPy's airybi). d/dx Bi = airy_bi_prime`; `Bi(0) = 3**(-1/6)/gamma(2/3).

Ai'(x) (SymPy's airyaiprime). d/dx Ai'(x) = x*Ai(x); Ai'(0) = -3**(-1/3)/gamma(1/3).

Bi'(x) (SymPy's airybiprime). d/dx Bi'(x) = x*Bi(x); Bi'(0) = 3**(1/6)/gamma(1/3).

Riemann zeta function zeta(s) (DLMF 25.2).

Exact at even positive integers (zeta(2) = pi**2/6, zeta(4) = pi**4/90, via Bernoulli numbers) and at non-positive integers (zeta(0) = -1/2, zeta(-1) = -1/12, zeta(-2n) = 0); zeta(3) stays zeta(3) (no closed form is known). Evaluates in f64, balls and at complex s through Arb. diff(zeta(s), s) is zeta_prime(s); higher derivatives are zeta_deriv(n, s).

Raises SpecialFunctionError (E-SPECFN-006) at the pole s = 1.

Derivative of the Riemann zeta function, zeta'(s).

Exact at s = 0 (-log(2*pi)/2), s = -1 (1/12 - log(Glaisher)), s = 2 (pi**2/6 (EulerGamma + log(2 pi) - 12 log(Glaisher))) and the negative even integers (zeta'(-2) = -zeta(3)/(4 pi**2)); the node zeta_prime(s) elsewhere. This is what diff(zeta(s), s) returns; its derivative is zeta_deriv(2, s). Evaluates in f64, balls (interval_eval) and at complex s through Arb (acb_dirichlet_zeta_jet).

Raises SpecialFunctionError (E-SPECFN-006) at the pole s = 1.

>>> import alkahest as ak
>>> p = ak.ExprPool(); s = p.symbol("s")
>>> ak.diff(ak.zeta(s), s).value           # zeta_prime(s)

Hurwitz zeta function zeta(s, a) = sum_{k>=0} (k + a)**(-s) (DLMF 25.11; SymPy's two-argument zeta).

hurwitz_zeta(s, 1) = zeta(s), hurwitz_zeta(s, 1/2) = (2**s - 1)*zeta(s), hurwitz_zeta(-n, a) = -B_{n+1}(a)/(n+1) (a Bernoulli polynomial in any a); d/da zeta(s, a) = -s*zeta(s+1, a). E-SPECFN-006 at s = 1.

Polylogarithm Li_s(z) (DLMF 25.12; SymPy's polylog), principal branch with the cut [1, oo).

polylog(2, z) is dilog(z); polylog(1, z) = -log(1 - z); polylog(0, z) = z/(1 - z); polylog(-n, z) is the Eulerian rational function (Li_{-1}(z) = z/(1 - z)**2); polylog(s, 1) = zeta(s) for s > 1 and polylog(s, -1) = -(1 - 2**(1-s))*zeta(s). d/dz Li_s(z) = Li_{s-1}(z)/z; a derivative in s is refused.

Raises SpecialFunctionError (E-SPECFN-006) at z = 1 for s <= 1 (a pole).

Rewrite every Bessel node whose order is now a half-integer literal (and every spherical Bessel node with an integer order) into its elementary form, and re-canonicalise zeta / hurwitz_zeta / polylog nodes — e.g. after subs(bessel_j(nu, x), {nu: Fraction(1, 2)}). SymPy's expand_func for these families.

The Euler–Mascheroni constant gamma = 0.5772... as the exact symbol EulerGamma in pool (SymPy's EulerGamma).

Prints and parses as EulerGamma (LaTeX \gamma); every numeric evaluator resolves it (evaluate, interval_eval at any precision, compile_expr, emit_c), and no symbolic operation folds it into a float: simplify(digamma(1)) is -EulerGamma.

>>> import alkahest as ak
>>> p = ak.ExprPool()
>>> ak.euler_gamma(p)                      # EulerGamma

Catalan's constant G = sum (-1)**k/(2k+1)**2 = 0.9159... as the exact symbol Catalan (LaTeX G); trigamma(1/4) = pi**2 + 8*Catalan.

The Glaisher–Kinkelin constant A = 1.2824... as the exact symbol Glaisher (LaTeX A).

The Stieltjes constant gamma_n (n >= 0) as an exact symbol.

gamma_0 is EulerGamma; gamma_n for 1 <= n <= 30 is the symbol StieltjesGamma{n} (LaTeX \\gamma_{n}). They are the coefficients of the Laurent series of the Riemann zeta function, zeta(s) = 1/(s-1) + sum_n (-1)**n * gamma_n * (s-1)**n / n!, and appear in series(zeta(s), s, 1, order). No closed form exists, so they are never folded into floats; evaluate and interval_eval resolve them.

Raises ValueError (E-SPECFN-004) for n > 30.

>>> import alkahest as ak
>>> p = ak.ExprPool()
>>> ak.stieltjes_constant(p, 1)            # StieltjesGamma1  (-0.0728...)

The polygamma function polygamma(n, x) = psi^(n)(x), the (n+1)-th derivative of log(gamma(x)) (DLMF 5.15; SymPy's polygamma).

polygamma(0, x) is digamma(x) and polygamma(1, x) is trigamma(x). Exact values use EulerGamma, pi, log, zeta and Catalan — never floats: polygamma(0, 1) = -EulerGamma, polygamma(1, 1) = pi**2/6, polygamma(2, 1) = -2*zeta(3), polygamma(0, 1/2) = -EulerGamma - 2*log(2), at any integer shift. d/dx polygamma(n, x) = polygamma(n+1, x); a derivative in n is refused (E-DIFF-001). Evaluates in f64 and balls through Arb.

Raises SpecialFunctionError: E-SPECFN-007 for a non-integer order, E-SPECFN-003 for a negative one, E-SPECFN-006 at a pole (x = 0, -1, -2, ...).

>>> import alkahest as ak
>>> p = ak.ExprPool()
>>> ak.polygamma(2, p.integer(1))          # -2*zeta(3)

Log-gamma function loggamma(x) = log(Gamma(x)) (DLMF 5.4; SymPy's and mpmath's loggamma).

Real for x > 0; for complex x the principal branch (the analytic continuation from the positive axis, Arb's acb_lgamma), which is not log(gamma(x)) off the positive axis. Exact at 1 and 2 (0), at integers n <= 100 (log((n-1)!)) and positive half-integers (loggamma(1/2) = log(pi)/2); the node elsewhere. diff(loggamma(x), x) is digamma(x). Evaluates in f64, balls and Taylor models (x > 0) and at complex x through Arb. Its Stirling series comes from experimental.asymptotic_expand.

Raises SpecialFunctionError (E-SPECFN-006) at a pole (a non-positive integer literal).

>>> import alkahest as ak
>>> p = ak.ExprPool(); x = p.symbol("x")
>>> ak.loggamma(p.integer(5))              # log(24)
>>> ak.diff(ak.loggamma(x), x).value        # digamma(x)

Rewrite products and quotients of gamma with exact identities (SymPy's gammasimp):

• recurrence: gamma(x + n)/gamma(x) becomes the rising factorial x (x+1) ... (x+n-1) for an integer n (|n| <= 64);
• reflection: gamma(x) gamma(1 - x) = pi/sin(pi x) (and gamma(x) gamma(m - x) for an integer m through the recurrence);
• duplication: gamma(x) gamma(x + 1/2) = 2**(1 - 2x) sqrt(pi) gamma(2x).

Each rewrite is an identity of meromorphic functions. The recurrence result is the continuous extension at the poles x = 0, -1, ... where gamma(x + n)/gamma(x) itself is undefined, so add x not a non-positive integer as a side condition when that matters. Nothing else is changed.

>>> import alkahest as ak
>>> p = ak.ExprPool(); x = p.symbol("x")
>>> ak.gammasimp(ak.gamma(x) * ak.gamma(1 - x))      # pi/sin(pi*x)
>>> ak.gammasimp(ak.gamma(x + 3) / ak.gamma(x))      # x*(x+1)*(x+2)

The generalized hypergeometric function pFq(a; b; z) (DLMF 16.2; SymPy's hyper(a, b, z)).

a and b are sequences of parameters (Expr, int, Fraction). The node prints as hyper(p, q, a1, ..., ap, b1, ..., bq, z) and parses back; parameters common to both lists cancel, and an upper 0 or z = 0 gives 1. d/dz shifts every parameter by one (prod(a)/prod(b) * hyper(a+1, b+1, z)); a derivative in a parameter is refused. f64 and ball values come from Arb (real values, the principal branch for p = q + 1); evaluate(..., mode="complex") gives complex values, with 2F1 continued analytically past |z| = 1 (on the cut [1, oo) the limit from below, as in mpmath; other p = q + 1 series decline for |z| >= 1). Use expand_hyper for closed forms and hyper_transform for the Euler/Pfaff/connection/quadratic/Kummer transformations.

Raises SpecialFunctionError (E-SPECFN-006) for a lower parameter that is a non-positive integer (unless the series terminates first).

>>> import alkahest as ak
>>> p = ak.ExprPool(); z = p.symbol("z")
>>> ak.expand_hyper(ak.hyper([1, 1], [2], z))   # -log(1 - z)/z

Rewrite every hyper node into a polynomial (a terminating series) or a closed form where one is known — SymPy's hyperexpand: 0F0(;;z) = exp(z), 1F0(a;;z) = (1-z)**(-a), 2F1(1,1;2;z) = -log(1-z)/z, 2F1(1/2,1;3/2;z) = atanh(sqrt(z))/sqrt(z), 2F1(1/2,1/2;3/2;z) = asin(sqrt(z))/sqrt(z), Gauss's 2F1(a,b;c;1) = gamma(c)gamma(c-a-b)/(gamma(c-a)gamma(c-b)) when c-a-b is a positive literal, 1F1(1;2;z) = (exp(z)-1)/z, 1F1(1/2;3/2;-w) = sqrt(pi)*erf(sqrt(w))/(2 sqrt(w)), and 0F1 as cosh/sinh or Bessel J/I; the 2F1(a, a+1/2; c; z) families for c = 1/2, 3/2, 2a, 2a+1, 1F1(a; a+1; -w) = a*w**(-a)*lowergamma(a, w), 1F1(a; 2a; z) as Bessel I, and any node one Euler, Pfaff or Kummer transformation away from one of these (1F1(1; 3/2; z) becomes an erf). Other nodes are kept.

Raises SpecialFunctionError (E-SPECFN-004) for a terminating series of degree above 200.

Rewrite one hypergeometric node by a classical transformation and return the exact right-hand side (Gamma prefactors as gamma nodes, the transformed series as new hyper nodes).

kind is one of (non-regularized A&S 15.3 / DLMF 15.8, principal branches; each identity holds on the region given, which is part of the result's contract):

• "euler" — (1-z)**(c-a-b) * 2F1(c-a, c-b; c; z), z off [1, oo);
• "pfaff" (alias "z/(z-1)") — (1-z)**(-a) * 2F1(a, c-b; c; z/(z-1)); "pfaff_b" the same in b;
• "1-z", "1/z", "1/(1-z)", "1-1/z" — the two-term connection formulas with Gamma coefficients (A&S 15.3.6-15.3.9); 1-z and 1-1/z hold for z off (-oo, 0] U [1, oo), 1/z and 1/(1-z) for z off [0, oo); a symbolic c-a-b / a-b is taken to be a non-integer;
• "quadratic" — the quadratic transformation of the family the parameters are in: c = 2b (DLMF 15.8.13, z off [1, oo)), c = a-b+1 (15.8.15, |z| < 1), 2c = a+b+1 (15.8.18, Re z < 1/2);
• "kummer" (1F1) — exp(z) * 1F1(b-a; b; -z), every z.

Raises SpecialFunctionError (E-SPECFN-009) when expr is not a 2F1 (1F1 for "kummer") node, for the logarithmic case of a two-term formula (c-a-b or a-b an integer), for c a non-positive integer, and for parameters in no quadratic family; ValueError for an unknown kind.

>>> import alkahest as ak
>>> p = ak.ExprPool(); z = p.symbol("z")
>>> ak.hyper_transform(ak.hyper([1, 1], [3], z), "euler")   # (1-z) * 2F1(2, 2; 3; z)

The Euler Beta function beta(a, b) = gamma(a)*gamma(b)/gamma(a + b) (DLMF 5.12.1; SymPy's beta, mpmath's beta).

Args: a, b (Expr, int or Fraction; at least one an Expr). Returns the exact value when the gamma ratio is exact (beta(2, 3) = 1/12`, `beta(1/2, 1/2) = pi), the rational function (n-1)!/(a*(a+1)*...*(a+n-1)) when one argument is a positive integer n (beta(a, 1) = 1/a), and the node beta(a, b) otherwise. Differentiates through digamma; evaluates in f64 and balls via Arb.

Raises SpecialFunctionError (E-SPECFN-006) when a or b is a non-positive integer (a pole).

>>> import alkahest as ak
>>> p = ak.ExprPool(); a = p.symbol("a")
>>> ak.beta(a, p.integer(2))      # 1/(a*(a + 1))

The upper incomplete gamma function uppergamma(s, z) = integral(t**(s-1)*exp(-t), (t, z, oo)) (DLMF 8.2.2; SymPy's uppergamma, mpmath's gammainc(s, z)).

Args: s (order), z (lower limit of the integral). Returns the elementary form for a positive integer s (uppergamma(1, z) = exp(-z)), the erfc form for a half-integer s > 0 (uppergamma(1/2, z) = sqrt(pi)*erfc(sqrt(z))), gamma(s) at z = 0 for a literal s > 0, else the node. d/dz uppergamma(s, z) = -z**(s-1)*exp(-z); a derivative in s is refused (E-DIFF-001). Real values via Arb.

Raises SpecialFunctionError (E-SPECFN-006) for z = 0 with a literal s <= 0.

>>> import alkahest as ak
>>> p = ak.ExprPool(); z = p.symbol("z")
>>> ak.uppergamma(p.integer(2), z)     # (z + 1)*exp(-z)

The lower incomplete gamma function lowergamma(s, z) = integral(t**(s-1)*exp(-t), (t, 0, z)) (DLMF 8.2.1; SymPy's lowergamma, mpmath's gammainc(s, 0, z)).

Args: s (order), z (upper limit). Elementary for a positive integer s, erf form for a half-integer s > 0 (lowergamma(1/2, z) = sqrt(pi)*erf(sqrt(z))), 0 at z = 0 for a literal s > 0; else the node. uppergamma(s, z) + lowergamma(s, z) = gamma(s)`. `d/dz lowergamma(s, z) = z**(s-1)*exp(-z).

Raises SpecialFunctionError (E-SPECFN-006) for a non-positive integer s.

>>> import alkahest as ak
>>> p = ak.ExprPool(); z = p.symbol("z")
>>> ak.lowergamma(p.integer(1), z)     # 1 - exp(-z)

No docstring.

No docstring.

No docstring.

Principal argument (Arg ∈ (−π, π]; stable top-level export).

Only domain-safe literal simplifications are applied; branch-cut cases remain symbolic.

Symbolic real-part constructor (stable top-level export).

Symbolic imaginary-part constructor (stable top-level export).

Symbolic complex conjugation (stable top-level export).

Bessel function of the first kind, order 0: J₀(x).

Stable top-level export (also re-exported from alkahest.experimental).

Bessel function of the first kind, order 1: J₁(x).

Stable top-level export (also re-exported from alkahest.experimental).

Digamma ψ(x) = Γ′(x)/Γ(x).

Stable top-level export (also re-exported from alkahest.experimental).

Trigamma ψ₁(x) = ψ′(x) = Σ_{k≥0} (x+k)⁻².

This is what digamma differentiates to. trigamma itself does not differentiate: the polygamma ladder ψ₀ → ψ₁ → ψ₂ → … has no closed-form terminator without a binary polygamma(n, x), so diff(trigamma(x), x) raises E-DIFF-001 rather than returning something wrong.

Stable top-level export.

Dilogarithm Li₂(x) = Σ_{k≥1} xᵏ/k², principal branch.

The branch cut runs along [1, ∞) (DLMF §25.12(i), Lewin §1.1, Mathematica PolyLog[2, z]). Li₂ is real on (−∞, 1] — the endpoint included, Li₂(1) = π²/6 — and numeric evaluation declines for x > 1, where the principal value is complex.

d/dx Li₂(x) = −log(1 − x)/x.

Alkahest ships dilog rather than a general polylog(s, x): see the module documentation of alkahest-core/src/primitive/polylog.rs for the reasoning (∂Li_s/∂s has no closed form, and the rigorous Taylor tier is unary). Li₁(x) needs no primitive — it is -log(1 - x).

Stable top-level export.

Complete elliptic integral of the first kind, EllipticK(m).

Elliptic integral of the second kind.

elliptic_e(m) is the complete integral EllipticE(m). elliptic_e(phi, m) is the incomplete integral EllipticE(phi, m).

Incomplete elliptic integral of the first kind, EllipticF(phi, m).

Incomplete elliptic integral of the third kind, EllipticPi(n, phi, m).

Ei(x), the exponential integral (DLMF 6.2.5). d/dx Ei(x) = exp(x)/x.

li(x), the logarithmic integral (DLMF 6.2.8), a Cauchy principal value for x > 1. d/dx li(x) = 1/log(x).

Si(x), the sine integral (DLMF 6.2.9). d/dx Si(x) = sin(x)/x.

Ci(x), the cosine integral (DLMF 6.2.11). d/dx Ci(x) = cos(x)/x.

Shi(x), the hyperbolic sine integral (DLMF 6.2.15). d/dx Shi(x) = sinh(x)/x.

Chi(x), the hyperbolic cosine integral (DLMF 6.2.16). d/dx Chi(x) = cosh(x)/x.

Fresnel sine integral S(x) = ∫₀ˣ sin(π t²/2) dt.

Normalised (π/2) convention — DLMF §7.2(iii), Abramowitz & Stegun §7.3, SymPy fresnels, SciPy scipy.special.fresnel, Mathematica FresnelS. S(∞) = 1/2. The unnormalised ∫₀ˣ sin(t²) dt is a different function (it tends to √(π/8) ≈ 0.6267), and mixing the two is a silent factor-of-√(π/2) error.

d/dx S(x) = sin(π x²/2).

Stable top-level export.

Fresnel cosine integral C(x) = ∫₀ˣ cos(π t²/2) dt.

Same normalised (π/2) convention as fresnels(); C(∞) = 1/2 and d/dx C(x) = cos(π x²/2).

Stable top-level export.

Lambert W, the inverse of w*exp(w): W(x)*exp(W(x)) == x.

k = 0 (default) is the principal branch W₀, real on [-1/e, ∞) with values ≥ -1. k = -1 is the lower real branch W₋₁ (the same node as lambert_wm1()), real on [-1/e, 0) with values ≤ -1. Together they give every real solution of w*exp(w) = x: two for -1/e < x < 0. Any other k raises ValueError.

Stable top-level export (also re-exported from alkahest.experimental).

Lower real branch of Lambert W, W₋₁(x).

Defined for -1/e <= x < 0, where it decreases from W₋₁(-1/e) = -1 to -∞ as x → 0⁻; it solves w*exp(w) = x with w <= -1 (lambert_w gives the other root, w >= -1). Evaluates in f64 and in certified ball arithmetic; outside the domain it has no real value. d/dx W₋₁(x) = W₋₁(x)/(x*(1 + W₋₁(x))). Simplifies W₋₁(a*exp(a)) = a for rational a <= -1 and W₋₁(-log(2)/2) = -2*log(2).

Example:

>>> w = alkahest.lambert_wm1(pool.rational(-1, 4))
>>> alkahest.eval_expr(w, {})   # mpmath.lambertw(-0.25, -1)
-2.15329236411035