Special Functions

from scipy import special
import numpy as np

# Common direct imports
from scipy.special import (
    # Gamma & related
    gamma, gammaln, gammasgn, loggamma,
    rgamma, multigammaln, gammainccinv,
    gammainc, gammaincc, gammaincinv,
    psi, digamma, polygamma, beta, betaln, betainc, betaincinv,
    # Error functions
    erf, erfc, erfinv, erfcinv, erfi,
    # Bessel functions
    jv, yv, iv, kv, jn, yn,
    j0, j1, y0, y1, i0, i1, k0, k1,
    jvp, yvp, ivp, kvp,
    jn_zeros, jnp_zeros,
    # Airy functions
    airy, airye, ai_zeros, bi_zeros,
    # Orthogonal polynomials
    eval_legendre, eval_chebyt, eval_chebyu,
    eval_hermite, eval_laguerre, eval_genlaguerre,
    # Hypergeometric
    hyp1f1, hyp2f1, hyp0f1,
    # Elliptic integrals
    ellipk, ellipe, ellipkm1, ellipj,
    # Combinatorics
    comb, perm, factorial, factorial2,
    # Information theory
    entr, rel_entr, xlogy, xlog1py,
    # Statistical
    ndtr, ndtri, chdtr, stdtr,
    # Misc
    expit, logit, log1p, expm1,
    sinc, sici, shichi,
    zeta, zetac,
    lambertw,
    # Spherical harmonics
    sph_harm_y,
)

from scipy.special import gamma, gammaln, loggamma, rgamma, multigammaln

# Gamma function Γ(z)
gamma(5)          # 24.0 = 4!
gamma(0.5)        # sqrt(pi)
gamma(np.array([1, 2, 3, 4]))  # vectorized

# Log-gamma (avoids overflow)
gammaln(100)      # log(Γ(100))
loggamma(100)     # same, returns complex for negative non-integers

# Reciprocal gamma 1/Γ(z)
rgamma(5)         # 1/24

# Multivariate log-gamma ln(Γ_p(a))
multigammaln(5, d=3)   # d = dimension

from scipy.special import gammainc, gammaincc, gammaincinv, gammainccinv

# Regularized lower incomplete gamma P(a, x) = γ(a,x)/Γ(a)
gammainc(2, 1)          # (a, x) -> P(2, 1) ≈ 0.264
gammaincc(2, 1)         # Q(a,x) = 1 - P(a,x) (upper)

# Inverse: x such that P(a,x) = p
gammaincinv(2, 0.264)   # (a, p)
gammainccinv(2, 0.736)  # inverse of upper

from scipy.special import digamma, psi, polygamma

psi(1)         # -0.5772... = -Euler-Mascheroni constant
digamma(1)     # same as psi

polygamma(0, 1)    # (n, x) — digamma
polygamma(1, 1)    # trigamma
polygamma(2, 1)    # tetragamma

from scipy.special import beta, betaln, betainc, betaincinv

beta(2, 3)         # B(2,3) = Γ(2)Γ(3)/Γ(5) = 1/12
betaln(2, 3)       # log(B(2,3)), more stable

# Regularized incomplete beta I_x(a,b)
betainc(2, 3, 0.5)         # (a, b, x)
betaincinv(2, 3, 0.6875)   # (a, b, p) — inverse

from scipy.special import erf, erfc, erfinv, erfcinv, erfi, erfcx

# erf(x) = (2/sqrt(pi)) * int_0^x exp(-t^2) dt
erf(1.0)      # ≈ 0.8427
erfc(1.0)     # 1 - erf(x) ≈ 0.1573 (more accurate for large x)
erfinv(0.8)   # inverse erf
erfcinv(0.2)  # inverse erfc

# Imaginary error function erf(ix)/i
erfi(1.0)

# Scaled complementary: erfcx(x) = exp(x^2) * erfc(x) (no overflow)
erfcx(10.0)   # avoids overflow where erfc underflows

from scipy.special import jv, yv, iv, kv, jvp, yvp, ivp, kvp

# J_v(z) — Bessel of the first kind; args are (v, z): order, then argument
jv(0, 1.0)
jv(1.5, np.linspace(0, 10, 100))

# Y_v(z) — Bessel of the second kind (Neumann)
yv(0, 2.0)

# I_v(z) — modified Bessel of the first kind
iv(0, 1.0)

# K_v(z) — modified Bessel of the second kind
kv(0, 1.0)

# Derivatives
jvp(0, 1.0, n=1)   # (v, z, n): d/dz J_0(z) at z=1
ivp(0, 1.0, n=2)   # second derivative of I_0

from scipy.special import j0, j1, y0, y1, i0, i1, i0e, i1e, k0, k1, k0e, k1e

j0(x), j1(x)     # J_0, J_1 (fast C implementations)
y0(x), y1(x)     # Y_0, Y_1
i0(x), i1(x)     # I_0, I_1
i0e(x), i1e(x)   # exp(-abs(x)) * I_0(x): exponentially scaled (avoids overflow)
k0(x), k1(x)     # K_0, K_1
k0e(x), k1e(x)   # exp(x) * K_0(x): exponentially scaled

from scipy.special import jn_zeros, jnp_zeros, yn_zeros, ynp_zeros

# First nt zeros of J_n(x)
jn_zeros(0, 5)    # [2.4048, 5.5201, 8.6537, 11.7915, 14.9309]
jn_zeros(1, 3)

# First nt zeros of J_n'(x)
jnp_zeros(0, 5)

from scipy.special import spherical_jn, spherical_yn, spherical_in, spherical_kn

spherical_jn(0, 1.0)
spherical_jn(1, 1.0, derivative=True)  # derivative
spherical_yn(0, 1.0)

from scipy.special import airy, airye, ai_zeros, bi_zeros

# Ai(z), Ai'(z), Bi(z), Bi'(z)
Ai, Aip, Bi, Bip = airy(1.0)

# Exponentially scaled (avoids overflow for large z)
Ai_e, Aip_e, Bi_e, Bip_e = airye(10.0)

# Zeros
a, ap, ai_ap, aip_a = ai_zeros(5)   # zeros of Ai, zeros of Ai', values Ai(a'), values Ai'(a)

from scipy.special import (
    eval_legendre, eval_chebyt, eval_chebyu,
    eval_hermite, eval_hermitenorm,
    eval_laguerre, eval_genlaguerre,
    eval_jacobi
)

# Legendre P_n(x), x in [-1, 1] — args are (n, x)
eval_legendre(3, 0.5)

# Chebyshev T_n(x) (first kind), U_n(x) (second kind)
eval_chebyt(4, 0.5)
eval_chebyu(4, 0.5)

# Hermite H_n(x) (physicists') and He_n(x) (probabilists')
eval_hermite(3, 1.0)
eval_hermitenorm(3, 1.0)

# Laguerre L_n(x), generalized L_n^alpha(x): (n, alpha, x)
eval_laguerre(3, 1.0)
eval_genlaguerre(3, 0.5, 1.0)

# Jacobi P_n^(alpha,beta)(x): (n, alpha, beta, x)
eval_jacobi(3, 0.5, 0.5, 0.5)

from scipy.special import roots_legendre, roots_chebyt, roots_hermite, roots_laguerre

# Gauss-Legendre quadrature nodes and weights
x, w = roots_legendre(10)   # on [-1, 1]

# Gauss-Chebyshev
x, w = roots_chebyt(10)

# Gauss-Hermite (physicists')
x, w = roots_hermite(10)

# Gauss-Laguerre
x, w = roots_laguerre(10)   # on [0, ∞)

from scipy.special import hyp1f1, hyp2f1, hyp0f1

# Confluent hypergeometric: M(a, b, z) = 1F1(a; b; z) — args (a, b, z)
hyp1f1(1, 2, 1.0)

# Gauss hypergeometric: 2F1(a, b; c; z), |z| < 1 — args (a, b, c, z)
hyp2f1(0.5, 0.5, 1.0, 0.5)

# 0F1(; b; z) — Bessel-related — args (b, z)
hyp0f1(1.0, 1.0)

from scipy.special import ellipk, ellipe, ellipkm1, ellipj

# Complete elliptic integral of the first kind K(m), second kind E(m)
# Note: argument is m = k^2, not k
K = ellipk(0.5)    # K(k^2 = 0.5)
E = ellipe(0.5)    # E(k^2 = 0.5)

# K(1-m) for near m=1 (avoids cancellation)
K1m = ellipkm1(0.001)   # = ellipk(1 - 0.001) but more accurate

# Jacobi elliptic functions sn(u|m), cn(u|m), dn(u|m), ph — args (u, m)
sn, cn, dn, ph = ellipj(1.0, 0.5)

from scipy.special import comb, perm, factorial, factorial2

# Binomial coefficient C(N, k)
comb(10, 3)                      # 120.0 (float)
comb(10, 3, exact=True)          # 120 (int, slower for large N)
comb(10, 3, repetition=True)     # with repetition = C(N+k-1, k)

# Permutations P(N, k)
perm(10, 3)                      # 720.0
perm(10, 3, exact=True)          # 720 (int)

# Factorial
factorial(10)                    # 3628800.0
factorial(10, exact=True)        # 3628800 (int)
factorial2(7)                    # 7!! = 7*5*3*1 = 105

from scipy.special import expit, logit, log_expit

# Sigmoid: expit(x) = 1 / (1 + exp(-x))
expit(0.0)   # 0.5
expit(np.array([-3, 0, 3]))

# Inverse sigmoid: logit(p) = log(p/(1-p))
logit(0.5)   # 0.0
logit(0.9)   # ≈ 2.197

# Numerically stable log(sigmoid): log_expit(x) = log(1/(1+exp(-x)))
log_expit(-100.0)   # no underflow

from scipy.special import entr, rel_entr, xlogy, xlog1py, kl_div

# entr(x) = -x*log(x) (0 if x=0, -inf if x<0)
entr(0.5)        # ≈ 0.347

# rel_entr(x, y) = x*log(x/y) (0 if x=0, inf if y=0 and x!=0)
rel_entr(0.3, 0.7)

# xlogy(x, y) = x*log(y) (0 if x=0, even if y=0)
xlogy(0, 0)   # 0 (not NaN)
xlogy(2, 3)   # 2*log(3)

# kl_div(x, y) = x*log(x/y) - x + y
kl_div(0.3, 0.7)

# KL divergence between two distributions
p = np.array([0.2, 0.5, 0.3])
q = np.array([0.1, 0.6, 0.3])
kl = np.sum(rel_entr(p, q))   # KL(P || Q)

from scipy.special import sinc, sici, shichi

# Normalized sinc: sin(pi*x) / (pi*x)
sinc(1.5)

# Sine and cosine integrals: Si(x), Ci(x)
Si, Ci = sici(np.pi)

# Hyperbolic sine and cosine integrals: Shi(x), Chi(x)
Shi, Chi = shichi(1.0)

from scipy.special import sph_harm_y   # SciPy 1.15+

# Y_n^m(theta, phi): args are (n, m, theta, phi)
# n = degree (l), m = order
# theta = POLAR angle (colatitude) in [0, pi]
# phi   = AZIMUTHAL angle in [0, 2*pi]
Y = sph_harm_y(2, 1, np.pi/4, np.pi/3)
# Returns complex value

# Grid of values
theta = np.linspace(0, np.pi, 50)
phi = np.linspace(0, 2*np.pi, 100)
T, P = np.meshgrid(theta, phi)
Y_grid = sph_harm_y(3, 2, T, P)

from scipy.special import zeta, zetac

# Riemann zeta: one-argument zeta(x) = sum_{k=1}^∞ 1/k^x
zeta(2)         # ≈ pi^2/6 = 1.6449...

# Hurwitz zeta: zeta(x, q) = sum_{k=0}^∞ 1/(k+q)^x
zeta(2, 1)      # q=1 reduces to the Riemann zeta

# zetac(x) = zeta(x) - 1 (accurate for large x, where zeta ≈ 1)
zetac(2)

from scipy.special import lambertw

# Principal branch W_0: W(z) * exp(W(z)) = z
lambertw(1.0)         # ≈ 0.5671 (Omega constant)
lambertw(0)           # 0
lambertw(-1/np.e)     # -1 (branch point)

# Secondary real branch W_{-1} for z in [-1/e, 0)
lambertw(-0.1, k=-1)  # k=-1 branch

# Complex input
lambertw(-1.0 + 0j)

from scipy.special import log1p, expm1

# Precision near x=0
np.log(1 + 1e-15)    # 0.0 (catastrophic cancellation)
log1p(1e-15)          # 9.99e-16 (correct)

np.exp(1e-15) - 1     # 0.0 (precision lost)
expm1(1e-15)          # 1e-15 (correct)