mathematicskit.special_functions#
The gamma and beta functions (scipy.special.gamma/beta); Bessel
functions of the first/second kind (scipy.special.jv/yv);
orthogonal polynomial families (Legendre, Chebyshev, Hermite, Laguerre)
via numpy.polynomial; and the discrete Fourier transform via
numpy.fft, plus a from-scratch radix-2 FFT kept as a pedagogical
comparison against a naive DFT and numpy.fft.
mathematicskit.special_functions: special functions and transforms.
The gamma and beta functions (scipy.special.gamma/beta) with
Stirling’s asymptotic series; Bessel, Airy, and Mathieu functions;
elliptic integrals, the arithmetic-geometric mean, and Jacobi elliptic
functions; Gauss’s and Kummer’s hypergeometric functions; the error
function and Fresnel integrals; the Riemann zeta function and Euler’s
prime product; the Lambert W function; orthogonal polynomial families
(Legendre, Chebyshev, Hermite, Laguerre) via numpy.polynomial, with
orthogonality verified numerically in tests; and the discrete Fourier
transform via numpy.fft, plus a from-scratch radix-2 FFT kept
specifically as a pedagogical comparison against a naive \(O(n^2)\)
DFT and numpy.fft – the one sub-module in this domain where
hand-rolling is the point.
- class mathematicskit.special_functions.AiryResult(x, ai, ai_prime, bi, bi_prime)[source]#
Bases:
objectThe Airy functions
Ai,Biand their derivatives at the sample pointsx.
- class mathematicskit.special_functions.FFTComparisonResult(naive_result, radix2_result, numpy_result, naive_time, radix2_time, numpy_time, max_error_naive_vs_numpy, max_error_radix2_vs_numpy)[source]#
Bases:
objectContainer comparing naive DFT, hand-rolled radix-2 FFT, and
numpy.fft.- Parameters:
- class mathematicskit.special_functions.FresnelResult(t, s, c)[source]#
Bases:
objectThe Fresnel integrals \(S(t)\) and \(C(t)\) at the sample points
t.
- class mathematicskit.special_functions.JacobiEllipticResult(sn, cn, dn, amplitude)[source]#
Bases:
objectJacobi elliptic functions \(\operatorname{sn}\), \(\operatorname{cn}\), \(\operatorname{dn}\) at
(u, m).
- mathematicskit.special_functions.airy_functions(x)[source]#
The Airy functions \(\operatorname{Ai}, \operatorname{Ai}', \operatorname{Bi}, \operatorname{Bi}'\) at
x.\(\operatorname{Ai}(x) = \frac1\pi\int_0^\infty \cos(t^3/3 + xt)\,dt\) decays as \(x\to+\infty\); \(\operatorname{Bi}\) grows. Both satisfy \(y'' = xy\). Via
scipy.special.airy().Examples
>>> from scipy.special import gamma >>> r = airy_functions(0.0) >>> round(float(r.ai), 12) == round(1 / (3 ** (2 / 3) * float(gamma(2 / 3))), 12) True
- mathematicskit.special_functions.arithmetic_geometric_mean(a, b)[source]#
Gauss’s arithmetic-geometric mean \(\mathrm{AGM}(a, b)\).
Iterate \(a_{n+1} = (a_n + b_n)/2\), \(b_{n+1} = \sqrt{a_n b_n}\); both sequences converge quadratically to a common limit. Via
scipy.special.agm(). See C. F. Gauss, Werke, vol. 3 (1866), 361-403.- Parameters:
- Returns:
float or ndarray
Examples
>>> import math >>> # Gauss's constant 1/AGM(1, sqrt 2) = 0.8346268... >>> round(1 / float(arithmetic_geometric_mean(1.0, math.sqrt(2.0))), 7) 0.8346268
- mathematicskit.special_functions.bessel_first_kind(nu, x)[source]#
Bessel function of the first kind, \(J_\nu(x)\).
Solves Bessel’s differential equation \(x^2y'' + xy' + (x^2 - \nu^2)y = 0\), regular (finite) at \(x=0\). Via
scipy.special.jv(). See Abramowitz & Stegun, Handbook of Mathematical Functions, Sec. 9.1.Examples
>>> round(float(bessel_first_kind(0.0, 0.0)), 6) # J_0(0) = 1 1.0 >>> round(float(bessel_first_kind(1.0, 0.0)), 6) # J_nu(0) = 0 for nu > 0 0.0
- mathematicskit.special_functions.bessel_second_kind(nu, x)[source]#
Bessel function of the second kind (Weber/Neumann function), \(Y_\nu(x)\).
The second, linearly independent solution of Bessel’s equation, singular (diverging to \(-\infty\)) at \(x=0\). Via
scipy.special.yv(). See Abramowitz & Stegun, Handbook of Mathematical Functions, Sec. 9.1.Examples
>>> bool(bessel_second_kind(0.0, 0.1) < bessel_second_kind(0.0, 1.0)) # diverges toward x=0 True
- mathematicskit.special_functions.beta_function(a, b)[source]#
The beta function \(B(a,b) = \dfrac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)} = \int_0^1 t^{a-1}(1-t)^{b-1}\,dt\).
Via
scipy.special.beta(). See Abramowitz & Stegun, Handbook of Mathematical Functions, Sec. 6.2.- Parameters:
- Returns:
float or ndarray
Examples
>>> from scipy import special >>> round(float(beta_function(2.0, 3.0)), 6) == round(float(special.gamma(2.0) * special.gamma(3.0) / special.gamma(5.0)), 6) True
- mathematicskit.special_functions.chebyshev_polynomial(n, x)[source]#
The Chebyshev polynomial (first kind) \(T_n(x) = \cos(n\arccos x)\), orthogonal on \([-1,1]\) with weight \(1/\sqrt{1-x^2}\).
Via
numpy.polynomial.chebyshev.chebval(). The minimal-sup-norm property of (monic-scaled) Chebyshev polynomials is exactly why Chebyshev nodes (seemathematicskit.numerical_analysis.systems.chebyshev) avoid the Runge phenomenon. See Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., Ch. 18.4.- Parameters:
n (
int)x (float or array-like of float, in
[-1, 1])
- Returns:
float or ndarray
Examples
>>> round(float(chebyshev_polynomial(3, 0.5)), 6) # T_3(cos(pi/3)) = cos(pi) = -1 -1.0
- mathematicskit.special_functions.compare_fft_methods(n, seed=0)[source]#
Time and cross-check the naive DFT, radix-2 FFT, and
numpy.ffton the same random signal.Demonstrates the \(O(n^2)\) vs. \(O(n\log n)\) gap directly: the naive DFT’s time grows quadratically with n while the radix-2 FFT’s (and
numpy.fft’s) grows only log-linearly, becoming dramatically faster for even moderately large n.- Parameters:
n (
int) – Signal length; must be a power of 2 (forfft_radix2()).seed (
int)
- Return type:
- Returns:
FFTComparisonResult
Examples
>>> result = compare_fft_methods(256, seed=0) >>> result.max_error_naive_vs_numpy < 1e-8 True >>> result.max_error_radix2_vs_numpy < 1e-8 True
- mathematicskit.special_functions.complementary_error_function(x)[source]#
The complementary error function \(\operatorname{erfc}(x) = 1 - \operatorname{erf}(x)\).
Via
scipy.special.erfc(), which stays accurate in the far tail where computing1 - erf(x)directly would cancel to zero.Examples
>>> float(complementary_error_function(10.0)) > 0.0 # 1 - erf(10) is exactly 0.0 in float64 True
- mathematicskit.special_functions.complete_elliptic_integral_first_kind(m)[source]#
The complete elliptic integral of the first kind, \(K(m) = \int_0^{\pi/2} \frac{d\theta}{\sqrt{1 - m\sin^2\theta}}\).
Via
scipy.special.ellipk()(parameter \(m = k^2\)). See Abramowitz & Stegun, Sec. 17.3.Examples
>>> import math >>> abs(float(complete_elliptic_integral_first_kind(0.0)) - math.pi / 2) < 1e-15 True
- mathematicskit.special_functions.complete_elliptic_integral_second_kind(m)[source]#
The complete elliptic integral of the second kind, \(E(m) = \int_0^{\pi/2} \sqrt{1 - m\sin^2\theta}\,d\theta\).
\(4aE(e^2)\) is the perimeter of an ellipse with semi-major axis \(a\) and eccentricity \(e\), the problem that gave elliptic integrals their name. Via
scipy.special.ellipe(). See Abramowitz & Stegun, Sec. 17.3.Examples
>>> float(complete_elliptic_integral_second_kind(1.0)) # degenerate ellipse: a segment 1.0
- mathematicskit.special_functions.confluent_hypergeometric_1f1(a, b, z)[source]#
Kummer’s confluent hypergeometric function \(M(a, b, z) = {}_1F_1(a; b; z) = \sum_{n=0}^\infty \frac{(a)_n}{(b)_n}\frac{z^n}{n!}\).
It solves Kummer’s equation \(zw'' + (b - z)w' - aw = 0\) and arises as the confluent limit \({}_1F_1(a; b; z) = \lim_{c\to\infty} {}_2F_1(a, c; b; z/c)\). Kummer’s transformation reads \(M(a, b, z) = e^z M(b - a, b, -z)\). Via
scipy.special.hyp1f1(). See E. E. Kummer, “De integralibus quibusdam definitis et seriebus infinitis,” Journal für die reine und angewandte Mathematik 17 (1837), 228-242.Examples
>>> import math >>> round(float(confluent_hypergeometric_1f1(2.0, 2.0, 1.0)), 12) == round(math.e, 12) # 1F1(a; a; z) = e^z True
- mathematicskit.special_functions.dft_naive(x)[source]#
The discrete Fourier transform, by direct evaluation of its defining sum.
\(X_k = \sum_{n=0}^{N-1} x_n e^{-2\pi i kn/N}\), computed as a dense matrix-vector product – \(O(N^2)\) time. Kept hand-rolled purely to demonstrate the speedup
fft_radix2()achieves over it. See Cormen et al., Introduction to Algorithms, 3rd ed., Ch. 30.1.- Parameters:
x (
ndarray) – Complex (or real) input signal, any length.- Return type:
- Returns:
ndarray, complex
Examples
>>> import numpy as np >>> x = np.array([1.0, 2.0, 3.0, 4.0]) >>> np.allclose(dft_naive(x), np.fft.fft(x)) True
- mathematicskit.special_functions.error_function(x)[source]#
The error function \(\operatorname{erf}(x) = \frac{2}{\sqrt\pi}\int_0^x e^{-t^2}\,dt\).
The probability that a normal variable lies within \(x\sqrt2\) standard deviations of its mean. Named by J. W. L. Glaisher, “On a class of definite integrals,” Philosophical Magazine 42 (1871), 294-302. Via
scipy.special.erf().Examples
>>> import math >>> round(float(error_function(1 / math.sqrt(2))), 4) # the "68%" of the 68-95-99.7 rule 0.6827
- mathematicskit.special_functions.euler_product(s, prime_bound)[source]#
Euler’s product \(\prod_{p \le P} (1 - p^{-s})^{-1}\) over primes up to
prime_bound.Euler (1737) showed that for \(s > 1\) this product converges to \(\zeta(s)\) as \(P\to\infty\), by unique prime factorization. See L. Euler, “Variae observationes circa series infinitas,” Commentarii Academiae Scientiarum Petropolitanae 9 (1744), 160-188 (presented 1737).
- Parameters:
- Returns:
float
Examples
>>> round(euler_product(2.0, 10), 6) # primes 2, 3, 5, 7 1.595052
- mathematicskit.special_functions.fft_numpy(x)[source]#
The discrete Fourier transform via
numpy.fft.fft().The primary API for actually computing a DFT in this domain: a highly optimized mixed-radix FFT (handling any length, not just powers of 2). See NumPy’s FFT documentation.
Examples
>>> import numpy as np >>> np.allclose(fft_numpy(np.array([1.0, 0.0, -1.0, 0.0])), np.fft.fft([1.0, 0.0, -1.0, 0.0])) True
- mathematicskit.special_functions.fft_radix2(x)[source]#
The Cooley-Tukey radix-2 fast Fourier transform.
Recursively splits the DFT of length \(N\) (a power of 2) into two DFTs of length \(N/2\) (even- and odd-indexed samples), combined via the “butterfly” identity \(X_k = E_k + e^{-2\pi ik/N}O_k\), \(X_{k+N/2} = E_k - e^{-2\pi ik/N}O_k\) – giving \(O(N\log N)\) time overall. Requires
len(x)to be a power of 2. See Cooley & Tukey (1965), Math. Comp. 19, and Cormen et al., Introduction to Algorithms, 3rd ed., Ch. 30.2.- Parameters:
x (
ndarray) – Complex (or real) input signal;len(x)must be a power of 2.- Return type:
- Returns:
ndarray, complex
Examples
>>> import numpy as np >>> x = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]) >>> np.allclose(fft_radix2(x), np.fft.fft(x)) True
- mathematicskit.special_functions.fresnel_integrals(t)[source]#
The Fresnel integrals \(S(t) = \int_0^t \sin(\pi\tau^2/2)\,d\tau\) and \(C(t) = \int_0^t \cos(\pi\tau^2/2)\,d\tau\).
The parametric curve \((C(t), S(t))\) is the Euler (Cornu) spiral, which winds into the points \(\pm(\tfrac12, \tfrac12)\) as \(t\to\pm\infty\). Via
scipy.special.fresnel(). See A. Fresnel, “Mémoire sur la diffraction de la lumière,” Mémoires de l’Académie Royale des Sciences de l’Institut de France 5 (1826), 339-475.Examples
>>> r = fresnel_integrals(1.0) >>> round(float(r.c), 6), round(float(r.s), 6) (0.779893, 0.438259)
- mathematicskit.special_functions.gamma_function(x)[source]#
The gamma function \(\Gamma(x) = \int_0^\infty t^{x-1}e^{-t}\,dt\).
Extends the factorial to real (and complex) arguments: \(\Gamma(n) = (n-1)!\) for positive integers
n. Viascipy.special.gamma(). See Abramowitz & Stegun, Handbook of Mathematical Functions, Sec. 6.1.Examples
>>> round(float(gamma_function(5.0)), 6) # 4! = 24 24.0 >>> round(float(gamma_function(0.5)), 6) # sqrt(pi) 1.772454
- mathematicskit.special_functions.hermite_polynomial(n, x)[source]#
The (physicists’) Hermite polynomial \(H_n(x)\), orthogonal on \((-\infty,\infty)\) with weight \(e^{-x^2}\).
Via
numpy.polynomial.hermite.hermval(). Arise as the eigenfunctions of the quantum harmonic oscillator. See Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., Ch. 18.3.Examples
>>> float(hermite_polynomial(0, 2.0)) 1.0 >>> float(hermite_polynomial(1, 2.0)) # H_1(x) = 2x 4.0 >>> round(float(hermite_polynomial(2, 1.0)), 6) # H_2(x) = 4x^2 - 2 2.0
- mathematicskit.special_functions.hypergeometric_2f1(a, b, c, z)[source]#
Gauss’s hypergeometric function \({}_2F_1(a, b; c; z) = \sum_{n=0}^\infty \frac{(a)_n (b)_n}{(c)_n}\frac{z^n}{n!}\).
\((q)_n = q(q+1)\cdots(q+n-1)\) is the rising factorial. The series converges for \(|z| < 1\);
scipy.special.hyp2f1()uses analytic continuation elsewhere. Gauss’s summation theorem (1812) evaluates it at \(z = 1\) for \(\operatorname{Re}(c - a - b) > 0\):\[{}_2F_1(a, b; c; 1) = \frac{\Gamma(c)\Gamma(c-a-b)}{\Gamma(c-a)\Gamma(c-b)}.\]See C. F. Gauss, “Disquisitiones generales circa seriem infinitam,” Commentationes Societatis Regiae Scientiarum Gottingensis Recentiores 2 (1813).
- Parameters:
- Returns:
float or ndarray
Examples
>>> import math >>> # ln(1 + z) = z * 2F1(1, 1; 2; -z) >>> round(0.5 * float(hypergeometric_2f1(1.0, 1.0, 2.0, -0.5)), 12) == round(math.log(1.5), 12) True
- mathematicskit.special_functions.inner_product(f, g, weight, a, b)[source]#
The weighted inner product \(\langle f, g\rangle = \int_a^b f(x)g(x)w(x)\,dx\), via
scipy.integrate.quad().Used to numerically verify an orthogonal polynomial family’s defining property directly: \(\langle p_m, p_n\rangle_w = 0\) for \(m \neq n\). See Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., Ch. 18.1.
- Parameters:
- Return type:
- Returns:
float
Examples
>>> from mathematicskit.special_functions.systems.orthogonal_polynomials import legendre_polynomial >>> p2 = lambda x: legendre_polynomial(2, x) >>> p3 = lambda x: legendre_polynomial(3, x) >>> abs(inner_product(p2, p3, lambda x: 1.0, -1.0, 1.0)) < 1e-10 True
- mathematicskit.special_functions.jacobi_elliptic_functions(u, m)[source]#
The Jacobi elliptic functions \(\operatorname{sn}\), \(\operatorname{cn}\), \(\operatorname{dn}\).
Defined by inverting the incomplete elliptic integral of the first kind: if \(u = \int_0^\varphi d\theta/\sqrt{1 - m\sin^2\theta}\) then \(\operatorname{sn} u = \sin\varphi\), \(\operatorname{cn} u = \cos\varphi\), \(\operatorname{dn} u = \sqrt{1 - m\sin^2\varphi}\). They are doubly periodic; on the real axis, \(\operatorname{sn}\) and \(\operatorname{cn}\) have period \(4K(m)\). Via
scipy.special.ellipj(). See C. G. J. Jacobi, Fundamenta nova theoriae functionum ellipticarum (1829).- Parameters:
- Returns:
JacobiEllipticResult
Examples
>>> r = jacobi_elliptic_functions(0.7, 0.5) >>> round(float(r.sn**2 + r.cn**2), 12) 1.0
- mathematicskit.special_functions.laguerre_polynomial(n, x)[source]#
The (simple) Laguerre polynomial \(L_n(x)\), orthogonal on \([0,\infty)\) with weight \(e^{-x}\).
Via
numpy.polynomial.laguerre.lagval(). Arise as the radial part of the hydrogen atom’s wavefunctions. See Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., Ch. 18.5.- Parameters:
n (
int)x (float or array-like of float,
x >= 0)
- Returns:
float or ndarray
Examples
>>> float(laguerre_polynomial(0, 3.0)) 1.0 >>> round(float(laguerre_polynomial(1, 1.0)), 6) # L_1(x) = 1 - x 0.0
- mathematicskit.special_functions.lambert_w(z, branch=0)[source]#
The Lambert W function: the solution \(w\) of \(we^w = z\).
Via
scipy.special.lambertw(). For real \(z \ge -1/e\) the principal branch (branch=0) is real, and for \(-1/e \le z < 0\) the lower branchbranch=-1is a second real solution; in those cases the imaginary part (identically zero) is dropped and a real result is returned. See Corless et al. (1996).- Parameters:
- Returns:
float, complex, or ndarray
Examples
>>> import math >>> round(float(lambert_w(math.e)), 12) # 1 * e^1 = e 1.0 >>> # both real branches at z = -ln(2)/2: W_0 = -ln 2, W_{-1} = -ln 4 >>> round(float(lambert_w(-math.log(2) / 2)), 12) == round(-math.log(2), 12) True >>> round(float(lambert_w(-math.log(2) / 2, branch=-1)), 12) == round(-math.log(4), 12) True
- mathematicskit.special_functions.legendre_polynomial(n, x)[source]#
The Legendre polynomial \(P_n(x)\), orthogonal on \([-1,1]\) with weight 1.
Via
numpy.polynomial.legendre.legval()on the basis coefficients for degreen. Solutions of Legendre’s differential equation; arise e.g. in the multipole expansion of the electrostatic/gravitational potential. See Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., Ch. 18.2.- Parameters:
n (
int)x (float or array-like of float, in
[-1, 1])
- Returns:
float or ndarray
Examples
>>> import numpy as np >>> round(float(legendre_polynomial(2, 1.0)), 6) # P_n(1) = 1 for all n 1.0 >>> np.round(legendre_polynomial(2, np.array([-1.0, 0.0, 1.0])), 6) array([ 1. , -0.5, 1. ])
- mathematicskit.special_functions.log_gamma_function(x)[source]#
The (natural) log of the gamma function, \(\ln|\Gamma(x)|\).
Via
scipy.special.gammaln()– numerically stable for largex, where \(\Gamma(x)\) itself overflowsfloat64long before its logarithm does. See NIST Digital Library of Mathematical Functions, Sec. 5.4.Examples
>>> round(float(log_gamma_function(171.0)), 4) 706.5731
- mathematicskit.special_functions.mathieu_characteristic_a(m, q)[source]#
Characteristic value \(a_m(q)\) for the even Mathieu function \(\operatorname{ce}_m\).
At \(q = 0\) the equation is \(y'' + ay = 0\) and \(a_m(0) = m^2\). Via
scipy.special.mathieu_a().Examples
>>> float(mathieu_characteristic_a(3, 0.0)) 9.0
- mathematicskit.special_functions.mathieu_characteristic_b(m, q)[source]#
Characteristic value \(b_m(q)\) for the odd Mathieu function \(\operatorname{se}_m\).
\(b_m(0) = m^2\). Via
scipy.special.mathieu_b().Examples
>>> float(mathieu_characteristic_b(2, 0.0)) 4.0
- mathematicskit.special_functions.mathieu_even(m, q, x)[source]#
The even Mathieu function \(\operatorname{ce}_m(x, q)\), with
xin radians.Normalized so that \(\operatorname{ce}_m(x, 0) = \cos mx\) for \(m \ge 1\) (and \(1/\sqrt2\) for \(m = 0\)). Via
scipy.special.mathieu_cem().- Parameters:
- Returns:
float or ndarray
Examples
>>> import math >>> round(float(mathieu_even(2, 0.0, 0.3)), 12) == round(math.cos(0.6), 12) True
- mathematicskit.special_functions.mathieu_odd(m, q, x)[source]#
The odd Mathieu function \(\operatorname{se}_m(x, q)\), with
xin radians.Normalized so that \(\operatorname{se}_m(x, 0) = \sin mx\). Via
scipy.special.mathieu_sem().- Parameters:
- Returns:
float or ndarray
Examples
>>> import math >>> round(float(mathieu_odd(1, 0.0, 0.3)), 12) == round(math.sin(0.3), 12) True
- mathematicskit.special_functions.riemann_zeta(s)[source]#
The Riemann zeta function \(\zeta(s) = \sum_{n=1}^\infty n^{-s}\) (for \(s > 1\)), analytically continued.
Via
scipy.special.zeta(), real arguments only. Euler’s 1735 solution of the Basel problem gives \(\zeta(2) = \pi^2/6\); the continuation gives values such as \(\zeta(0) = -1/2\) and \(\zeta(-1) = -1/12\), and the “trivial zeros” \(\zeta(-2n) = 0\). See NIST Digital Library of Mathematical Functions, Ch. 25.- Parameters:
s (float or array-like of float) – Real argument, \(s \ne 1\).
- Returns:
float or ndarray
Examples
>>> import math >>> round(float(riemann_zeta(2.0)), 12) == round(math.pi**2 / 6, 12) True >>> round(float(riemann_zeta(-1.0)), 12) -0.083333333333
- mathematicskit.special_functions.stirling_factorial(n)[source]#
Stirling’s approximation \(n! \approx \sqrt{2\pi n}\,(n/e)^n\).
The leading term of Stirling’s series (J. Stirling, Methodus Differentialis, 1730, Prop. 28; A. de Moivre, Miscellanea Analytica, 1730). The relative error is about \(1/(12n)\), so the ratio to \(n!\) tends to 1 even though the absolute error grows without bound.
Examples
>>> import math >>> round(float(stirling_factorial(10.0)) / math.factorial(10), 4) 0.9917
- mathematicskit.special_functions.stirling_log_gamma(x, terms=3)[source]#
Stirling’s asymptotic series for \(\ln\Gamma(x)\), truncated after
termscorrections.\[\ln\Gamma(x) \sim \left(x - \tfrac12\right)\ln x - x + \tfrac12\ln(2\pi) + \sum_{k=1}^{K} \frac{B_{2k}}{2k(2k-1)\,x^{2k-1}},\]where \(B_{2k}\) are the Bernoulli numbers (from
scipy.special.bernoulli()). The series is asymptotic, not convergent: for fixed \(x\), adding terms helps only up to roughly \(K \approx \pi x\), after which the error grows. See NIST Digital Library of Mathematical Functions, Eq. 5.11.1.- Parameters:
- Returns:
float or ndarray
Examples
>>> abs(float(stirling_log_gamma(10.0, terms=3)) - float(log_gamma_function(10.0))) < 1e-9 True