Special functions (scipy.special)

Nearly all of the functions below are universal functions and follow broadcasting and automatic array-looping rules. Exceptions are noted.

Error handling

Errors are handled by returning nans, or other appropriate values. Some of the special function routines will emit warnings when an error occurs. By default this is disabled. To enable such messages use errprint(1), and to disable such messages use errprint(0).

Example:

>>> print scipy.special.bdtr(-1,10,0.3)
>>> scipy.special.errprint(1)
>>> print scipy.special.bdtr(-1,10,0.3)
errprint errprint({flag}) sets the error printing flag for special functions

Available functions

Airy functions

airy(x[, out1, out2, out3, out4]) (Ai,Aip,Bi,Bip)=airy(z) calculates the Airy functions and their derivatives
airye(x[, out1, out2, out3, out4]) (Aie,Aipe,Bie,Bipe)=airye(z) calculates the exponentially scaled Airy functions and
ai_zeros(nt) Compute the zeros of Airy Functions Ai(x) and Ai’(x), a and a’
bi_zeros(nt) Compute the zeros of Airy Functions Bi(x) and Bi’(x), b and b’

Elliptic Functions and Integrals

ellipj(x1, x2[, out1, out2, out3, out4]) (sn,cn,dn,ph)=ellipj(u,m) calculates the Jacobian elliptic functions of
ellipk(...) This function is rather imprecise around m==1.
ellipkm1(x[, out]) y=ellipkm1(1 - m) returns the complete integral of the first kind:
ellipkinc(x1, x2[, out]) y=ellipkinc(phi,m) returns the incomplete elliptic integral of the first
ellipe(x[, out]) y=ellipe(m) returns the complete integral of the second kind:
ellipeinc(x1, x2[, out]) y=ellipeinc(phi,m) returns the incomplete elliptic integral of the

Bessel Functions

jn(x1, x2[, out]) y=jv(v,z) returns the Bessel function of real order v at complex z.
jv(x1, x2[, out]) y=jv(v,z) returns the Bessel function of real order v at complex z.
jve(x1, x2[, out]) y=jve(v,z) returns the exponentially scaled Bessel function of real order
yn(x1, x2[, out]) y=yn(n,x) returns the Bessel function of the second kind of integer
yv(x1, x2[, out]) y=yv(v,z) returns the Bessel function of the second kind of real
yve(x1, x2[, out]) y=yve(v,z) returns the exponentially scaled Bessel function of the second
kn(x1, x2[, out]) y=kn(n,x) returns the modified Bessel function of the second kind (sometimes called the third kind) for
kv(x1, x2[, out]) y=kv(v,z) returns the modified Bessel function of the second kind (sometimes called the third kind) for
kve(x1, x2[, out]) y=kve(v,z) returns the exponentially scaled, modified Bessel function
iv(x1, x2[, out]) y=iv(v,z) returns the modified Bessel function of real order v of
ive(x1, x2[, out]) y=ive(v,z) returns the exponentially scaled modified Bessel function of
hankel1(x1, x2[, out]) y=hankel1(v,z) returns the Hankel function of the first kind for real order v and complex argument z.
hankel1e(x1, x2[, out]) y=hankel1e(v,z) returns the exponentially scaled Hankel function of the first
hankel2(x1, x2[, out]) y=hankel2(v,z) returns the Hankel function of the second kind for real order v and complex argument z.
hankel2e(x1, x2[, out]) y=hankel2e(v,z) returns the exponentially scaled Hankel function of the second

The following is not an universal function:

lmbda(v, x) Compute sequence of lambda functions with arbitrary order v and their derivatives.

Zeros of Bessel Functions

These are not universal functions:

jnjnp_zeros(nt) Compute nt (<=1200) zeros of the bessel functions Jn and Jn’
jnyn_zeros(n, nt) Compute nt zeros of the Bessel functions Jn(x), Jn’(x), Yn(x), and
jn_zeros(n, nt) Compute nt zeros of the Bessel function Jn(x).
jnp_zeros(n, nt) Compute nt zeros of the Bessel function Jn’(x).
yn_zeros(n, nt) Compute nt zeros of the Bessel function Yn(x).
ynp_zeros(n, nt) Compute nt zeros of the Bessel function Yn’(x).
y0_zeros(nt[, complex]) Returns nt (complex or real) zeros of Y0(z), z0, and the value
y1_zeros(nt[, complex]) Returns nt (complex or real) zeros of Y1(z), z1, and the value
y1p_zeros(nt[, complex]) Returns nt (complex or real) zeros of Y1’(z), z1’, and the value

Faster versions of common Bessel Functions

j0(x[, out]) y=j0(x) returns the Bessel function of order 0 at x.
j1(x[, out]) y=j1(x) returns the Bessel function of order 1 at x.
y0(x[, out]) y=y0(x) returns the Bessel function of the second kind of order 0 at x.
y1(x[, out]) y=y1(x) returns the Bessel function of the second kind of order 1 at x.
i0(x[, out]) y=i0(x) returns the modified Bessel function of order 0 at x.
i0e(x[, out]) y=i0e(x) returns the exponentially scaled modified Bessel function
i1(x[, out]) y=i1(x) returns the modified Bessel function of order 1 at x.
i1e(x[, out]) y=i1e(x) returns the exponentially scaled modified Bessel function
k0(x[, out]) y=k0(x) returns the modified Bessel function of the second kind (sometimes called the third kind) of
k0e(x[, out]) y=k0e(x) returns the exponentially scaled modified Bessel function
k1(x[, out]) y=i1(x) returns the modified Bessel function of the second kind (sometimes called the third kind) of
k1e(x[, out]) y=k1e(x) returns the exponentially scaled modified Bessel function

Integrals of Bessel Functions

itj0y0(x[, out1, out2]) (ij0,iy0)=itj0y0(x) returns simple integrals from 0 to x of the zeroth order
it2j0y0(x[, out1, out2]) (ij0,iy0)=it2j0y0(x) returns the integrals int((1-j0(t))/t,t=0..x) and
iti0k0(x[, out1, out2]) (ii0,ik0)=iti0k0(x) returns simple integrals from 0 to x of the zeroth order
it2i0k0(x[, out1, out2]) (ii0,ik0)=it2i0k0(x) returns the integrals int((i0(t)-1)/t,t=0..x) and
besselpoly(x1, x2, x3[, out]) y=besselpoly(a,lam,nu) returns the value of the integral:

Derivatives of Bessel Functions

jvp(v, z[, n]) Return the nth derivative of Jv(z) with respect to z.
yvp(v, z[, n]) Return the nth derivative of Yv(z) with respect to z.
kvp(v, z[, n]) Return the nth derivative of Kv(z) with respect to z.
ivp(v, z[, n]) Return the nth derivative of Iv(z) with respect to z.
h1vp(v, z[, n]) Return the nth derivative of H1v(z) with respect to z.
h2vp(v, z[, n]) Return the nth derivative of H2v(z) with respect to z.

Spherical Bessel Functions

These are not universal functions:

sph_jn(n, z) Compute the spherical Bessel function jn(z) and its derivative for
sph_yn(n, z) Compute the spherical Bessel function yn(z) and its derivative for
sph_jnyn(n, z) Compute the spherical Bessel functions, jn(z) and yn(z) and their
sph_in(n, z) Compute the spherical Bessel function in(z) and its derivative for
sph_kn(n, z) Compute the spherical Bessel function kn(z) and its derivative for
sph_inkn(n, z) Compute the spherical Bessel functions, in(z) and kn(z) and their

Riccati-Bessel Functions

These are not universal functions:

riccati_jn(n, x) Compute the Ricatti-Bessel function of the first kind and its
riccati_yn(n, x) Compute the Ricatti-Bessel function of the second kind and its

Struve Functions

struve(x1, x2[, out]) y=struve(v,x) returns the Struve function Hv(x) of order v at x, x
modstruve(x1, x2[, out]) y=modstruve(v,x) returns the modified Struve function Lv(x) of order
itstruve0(x[, out]) y=itstruve0(x) returns the integral of the Struve function of order 0
it2struve0(x[, out]) y=it2struve0(x) returns the integral of the Struve function of order 0
itmodstruve0(x[, out]) y=itmodstruve0(x) returns the integral of the modified Struve function

Raw Statistical Functions

See also

scipy.stats: Friendly versions of these functions.

bdtr(x1, x2, x3[, out]) y=bdtr(k,n,p) returns the sum of the terms 0 through k of the
bdtrc(x1, x2, x3[, out]) y=bdtrc(k,n,p) returns the sum of the terms k+1 through n of the
bdtri(x1, x2, x3[, out]) p=bdtri(k,n,y) finds the probability p such that the sum of the
btdtr(x1, x2, x3[, out]) y=btdtr(a,b,x) returns the area from zero to x under the beta
btdtri(x1, x2, x3[, out]) x=btdtri(a,b,p) returns the pth quantile of the beta distribution. It is
fdtr(x1, x2, x3[, out]) y=fdtr(dfn,dfd,x) returns the area from zero to x under the F density
fdtrc(x1, x2, x3[, out]) y=fdtrc(dfn,dfd,x) returns the complemented F distribution function.
fdtri(x1, x2, x3[, out]) x=fdtri(dfn,dfd,p) finds the F density argument x such that
gdtr(x1, x2, x3[, out]) y=gdtr(a,b,x) returns the integral from zero to x of the gamma
gdtrc(x1, x2, x3[, out]) y=gdtrc(a,b,x) returns the integral from x to infinity of the gamma
gdtria(x1, x2, x3[, out])
gdtrib(x1, x2, x3[, out])
gdtrix(x1, x2, x3[, out])
nbdtr(x1, x2, x3[, out]) y=nbdtr(k,n,p) returns the sum of the terms 0 through k of the
nbdtrc(x1, x2, x3[, out]) y=nbdtrc(k,n,p) returns the sum of the terms k+1 to infinity of the
nbdtri(x1, x2, x3[, out]) p=nbdtri(k,n,y) finds the argument p such that nbdtr(k,n,p)=y.
pdtr(x1, x2[, out]) y=pdtr(k,m) returns the sum of the first k terms of the Poisson
pdtrc(x1, x2[, out]) y=pdtrc(k,m) returns the sum of the terms from k+1 to infinity of the
pdtri(x1, x2[, out]) m=pdtri(k,y) returns the Poisson variable m such that the sum
stdtr(x1, x2[, out]) p=stdtr(df,t) returns the integral from minus infinity to t of the Student t
stdtridf(x1, x2[, out]) t=stdtridf(p,t) returns the argument df such that stdtr(df,t) is equal to p.
stdtrit(x1, x2[, out]) t=stdtrit(df,p) returns the argument t such that stdtr(df,t) is equal to p.
chdtr(x1, x2[, out]) p=chdtr(v,x) Returns the area under the left hand tail (from 0 to x) of the Chi
chdtrc(x1, x2[, out]) p=chdtrc(v,x) returns the area under the right hand tail (from x to
chdtri(x1, x2[, out]) x=chdtri(v,p) returns the argument x such that chdtrc(v,x) is equal
ndtr(x[, out]) y=ndtr(x) returns the area under the standard Gaussian probability
ndtri(x[, out]) x=ndtri(y) returns the argument x for which the area udnder the
smirnov(x1, x2[, out]) y=smirnov(n,e) returns the exact Kolmogorov-Smirnov complementary
smirnovi(x1, x2[, out]) e=smirnovi(n,y) returns e such that smirnov(n,e) = y.
kolmogorov(x[, out]) p=kolmogorov(y) returns the complementary cumulative distribution
kolmogi(x[, out]) y=kolmogi(p) returns y such that kolmogorov(y) = p
tklmbda(x1, x2[, out])
logit(x[, out]) NULL
expit(x[, out]) NULL

Error Function and Fresnel Integrals

erf(x[, out]) y=erf(z) returns the error function of complex argument defined as
erfc(x[, out]) y=erfc(x) returns 1 - erf(x).
erfinv(y)
erfcinv(y)
fresnel(x[, out1, out2]) (ssa,cca)=fresnel(z) returns the fresnel sin and cos integrals: integral(sin(pi/2
fresnel_zeros(nt) Compute nt complex zeros of the sine and cosine fresnel integrals
modfresnelp(x[, out1, out2]) (fp,kp)=modfresnelp(x) returns the modified fresnel integrals F_+(x) and K_+(x)
modfresnelm(x[, out1, out2]) (fm,km)=modfresnelp(x) returns the modified fresnel integrals F_-(x) amd K_-(x)

These are not universal functions:

erf_zeros(nt) Compute nt complex zeros of the error function erf(z).
fresnelc_zeros(nt) Compute nt complex zeros of the cosine fresnel integral C(z).
fresnels_zeros(nt) Compute nt complex zeros of the sine fresnel integral S(z).

Legendre Functions

lpmv(x1, x2, x3[, out]) y=lpmv(m,v,x) returns the associated legendre function of integer order
sph_harm Compute spherical harmonics.

These are not universal functions:

lpn(n, z) Compute sequence of Legendre functions of the first kind (polynomials),
lqn(n, z) Compute sequence of Legendre functions of the second kind,
lpmn(m, n, z) Associated Legendre functions of the first kind, Pmn(z) and its
lqmn(m, n, z) Associated Legendre functions of the second kind, Qmn(z) and its

Orthogonal polynomials

The following functions evaluate values of orthogonal polynomials:

eval_legendre Evaluate Legendre polynomial at a point.
eval_chebyt Evaluate Chebyshev T polynomial at a point.
eval_chebyu Evaluate Chebyshev U polynomial at a point.
eval_chebyc Evaluate Chebyshev C polynomial at a point.
eval_chebys Evaluate Chebyshev S polynomial at a point.
eval_jacobi Evaluate Jacobi polynomial at a point.
eval_laguerre Evaluate Laguerre polynomial at a point.
eval_genlaguerre Evaluate generalized Laguerre polynomial at a point.
eval_hermite Evaluate Hermite polynomial at a point.
eval_hermitenorm Evaluate normalized Hermite polynomial at a point.
eval_gegenbauer Evaluate Gegenbauer polynomial at a point.
eval_sh_legendre Evaluate shifted Legendre polynomial at a point.
eval_sh_chebyt Evaluate shifted Chebyshev T polynomial at a point.
eval_sh_chebyu Evaluate shifted Chebyshev U polynomial at a point.
eval_sh_jacobi Evaluate shifted Jacobi polynomial at a point.

The functions below, in turn, return orthopoly1d objects, which functions similarly as numpy.poly1d. The orthopoly1d class also has an attribute weights which returns the roots, weights, and total weights for the appropriate form of Gaussian quadrature. These are returned in an n x 3 array with roots in the first column, weights in the second column, and total weights in the final column.

legendre(n[, monic]) Returns the nth order Legendre polynomial, P_n(x), orthogonal over
chebyt(n[, monic]) Return nth order Chebyshev polynomial of first kind, Tn(x). Orthogonal
chebyu(n[, monic]) Return nth order Chebyshev polynomial of second kind, Un(x). Orthogonal
chebyc(n[, monic]) Return nth order Chebyshev polynomial of first kind, Cn(x). Orthogonal
chebys(n[, monic]) Return nth order Chebyshev polynomial of second kind, Sn(x). Orthogonal
jacobi(n, alpha, beta[, monic]) Returns the nth order Jacobi polynomial, P^(alpha,beta)_n(x)
laguerre(n[, monic]) Return the nth order Laguerre polynoimal, L_n(x), orthogonal over
genlaguerre(n, alpha[, monic]) Returns the nth order generalized (associated) Laguerre polynomial,
hermite(n[, monic]) Return the nth order Hermite polynomial, H_n(x), orthogonal over
hermitenorm(n[, monic]) Return the nth order normalized Hermite polynomial, He_n(x), orthogonal
gegenbauer(n, alpha[, monic]) Return the nth order Gegenbauer (ultraspherical) polynomial,
sh_legendre(n[, monic]) Returns the nth order shifted Legendre polynomial, P^*_n(x), orthogonal
sh_chebyt(n[, monic]) Return nth order shifted Chebyshev polynomial of first kind, Tn(x).
sh_chebyu(n[, monic]) Return nth order shifted Chebyshev polynomial of second kind, Un(x).
sh_jacobi(n, p, q[, monic]) Returns the nth order Jacobi polynomial, G_n(p,q,x)

Warning

Large-order polynomials obtained from these functions are numerically unstable.

orthopoly1d objects are converted to poly1d, when doing arithmetic. numpy.poly1d works in power basis and cannot represent high-order polynomials accurately, which can cause significant inaccuracy.

Hypergeometric Functions

hyp2f1(x1, x2, x3, x4[, out]) y=hyp2f1(a,b,c,z) returns the gauss hypergeometric function
hyp1f1(x1, x2, x3[, out]) y=hyp1f1(a,b,x) returns the confluent hypergeometeric function
hyperu(x1, x2, x3[, out]) y=hyperu(a,b,x) returns the confluent hypergeometric function of the
hyp0f1(v, z) Confluent hypergeometric limit function 0F1.
hyp2f0(x1, x2, x3, x4[, out1, out2]) (y,err)=hyp2f0(a,b,x,type) returns (y,err) with the hypergeometric function 2F0 in y and an error estimate in err. The input type determines a convergence factor and
hyp1f2(x1, x2, x3, x4[, out1, out2]) (y,err)=hyp1f2(a,b,c,x) returns (y,err) with the hypergeometric function 1F2 in y and an error estimate in err.
hyp3f0(x1, x2, x3, x4[, out1, out2]) (y,err)=hyp3f0(a,b,c,x) returns (y,err) with the hypergeometric function 3F0 in y and an error estimate in err.

Parabolic Cylinder Functions

pbdv(x1, x2[, out1, out2]) (d,dp)=pbdv(v,x) returns (d,dp) with the parabolic cylinder function Dv(x) in
pbvv(x1, x2[, out1, out2]) (v,vp)=pbvv(v,x) returns (v,vp) with the parabolic cylinder function Vv(x) in
pbwa(x1, x2[, out1, out2]) (w,wp)=pbwa(a,x) returns (w,wp) with the parabolic cylinder function W(a,x) in

These are not universal functions:

pbdv_seq(v, x) Compute sequence of parabolic cylinder functions Dv(x) and
pbvv_seq(v, x) Compute sequence of parabolic cylinder functions Dv(x) and
pbdn_seq(n, z) Compute sequence of parabolic cylinder functions Dn(z) and

Spheroidal Wave Functions

pro_ang1(x1, x2, x3, x4[, out1, out2]) (s,sp)=pro_ang1(m,n,c,x) computes the prolate sheroidal angular function
pro_rad1(x1, x2, x3, x4[, out1, out2]) (s,sp)=pro_rad1(m,n,c,x) computes the prolate sheroidal radial function
pro_rad2(x1, x2, x3, x4[, out1, out2]) (s,sp)=pro_rad2(m,n,c,x) computes the prolate sheroidal radial function
obl_ang1(x1, x2, x3, x4[, out1, out2]) (s,sp)=obl_ang1(m,n,c,x) computes the oblate sheroidal angular function
obl_rad1(x1, x2, x3, x4[, out1, out2]) (s,sp)=obl_rad1(m,n,c,x) computes the oblate sheroidal radial function
obl_rad2(x1, x2, x3, x4[, out1, out2]) (s,sp)=obl_rad2(m,n,c,x) computes the oblate sheroidal radial function
pro_cv(x1, x2, x3[, out]) cv=pro_cv(m,n,c) computes the characteristic value of prolate spheroidal
obl_cv(x1, x2, x3[, out]) cv=obl_cv(m,n,c) computes the characteristic value of oblate spheroidal
pro_cv_seq(m, n, c) Compute a sequence of characteristic values for the prolate
obl_cv_seq(m, n, c) Compute a sequence of characteristic values for the oblate

The following functions require pre-computed characteristic value:

pro_ang1_cv(x1, x2, x3, x4, x5[, out1, out2]) (s,sp)=pro_ang1_cv(m,n,c,cv,x) computes the prolate sheroidal angular function
pro_rad1_cv(x1, x2, x3, x4, x5[, out1, out2]) (s,sp)=pro_rad1_cv(m,n,c,cv,x) computes the prolate sheroidal radial function
pro_rad2_cv(x1, x2, x3, x4, x5[, out1, out2]) (s,sp)=pro_rad2_cv(m,n,c,cv,x) computes the prolate sheroidal radial function
obl_ang1_cv(x1, x2, x3, x4, x5[, out1, out2]) (s,sp)=obl_ang1_cv(m,n,c,cv,x) computes the oblate sheroidal angular function
obl_rad1_cv(x1, x2, x3, x4, x5[, out1, out2]) (s,sp)=obl_rad1_cv(m,n,c,cv,x) computes the oblate sheroidal radial function
obl_rad2_cv(x1, x2, x3, x4, x5[, out1, out2]) (s,sp)=obl_rad2_cv(m,n,c,cv,x) computes the oblate sheroidal radial function

Kelvin Functions

kelvin(x[, out1, out2, out3, out4]) (Be, Ke, Bep, Kep)=kelvin(x) returns the tuple (Be, Ke, Bep, Kep) which containes
kelvin_zeros(nt) Compute nt zeros of all the kelvin functions returned in a length 8 tuple of arrays of length nt.
ber(x[, out]) y=ber(x) returns the Kelvin function ber x
bei(x[, out]) y=bei(x) returns the Kelvin function bei x
berp(x[, out]) y=berp(x) returns the derivative of the Kelvin function ber x
beip(x[, out]) y=beip(x) returns the derivative of the Kelvin function bei x
ker(x[, out]) y=ker(x) returns the Kelvin function ker x
kei(x[, out]) y=kei(x) returns the Kelvin function ker x
kerp(x[, out]) y=kerp(x) returns the derivative of the Kelvin function ker x
keip(x[, out]) y=keip(x) returns the derivative of the Kelvin function kei x

These are not universal functions:

ber_zeros(nt) Compute nt zeros of the kelvin function ber x
bei_zeros(nt) Compute nt zeros of the kelvin function bei x
berp_zeros(nt) Compute nt zeros of the kelvin function ber’ x
beip_zeros(nt) Compute nt zeros of the kelvin function bei’ x
ker_zeros(nt) Compute nt zeros of the kelvin function ker x
kei_zeros(nt) Compute nt zeros of the kelvin function kei x
kerp_zeros(nt) Compute nt zeros of the kelvin function ker’ x
keip_zeros(nt) Compute nt zeros of the kelvin function kei’ x

Other Special Functions

expn(x1, x2[, out]) y=expn(n,x) returns the exponential integral for integer n and
exp1(x[, out]) y=exp1(z) returns the exponential integral (n=1) of complex argument
expi(x[, out]) y=expi(x) returns an exponential integral of argument x defined as
wofz(x[, out]) y=wofz(z) returns the value of the fadeeva function for complex argument
dawsn(x[, out]) y=dawsn(x) returns dawson’s integral: exp(-x**2) *
shichi(x[, out1, out2]) (shi,chi)=shichi(x) returns the hyperbolic sine and cosine integrals:
sici(x[, out1, out2]) (si,ci)=sici(x) returns in si the integral of the sinc function from 0 to x:
spence(x[, out]) y=spence(x) returns the dilogarithm integral: -integral(log t /
lambertw(z[, k, tol]) Lambert W function.
zeta(x1, x2[, out]) y=zeta(x,q) returns the Riemann zeta function of two arguments:
zetac(x[, out]) y=zetac(x) returns 1.0 - the Riemann zeta function: sum(k**(-x), k=2..inf)

Convenience Functions

cbrt(x[, out]) y=cbrt(x) returns the real cube root of x.
exp10(x[, out]) y=exp10(x) returns 10 raised to the x power.
exp2(x[, out]) y=exp2(x) returns 2 raised to the x power.
radian(x1, x2, x3[, out]) y=radian(d,m,s) returns the angle given in (d)egrees, (m)inutes, and
cosdg(x[, out]) y=cosdg(x) calculates the cosine of the angle x given in degrees.
sindg(x[, out]) y=sindg(x) calculates the sine of the angle x given in degrees.
tandg(x[, out]) y=tandg(x) calculates the tangent of the angle x given in degrees.
cotdg(x[, out]) y=cotdg(x) calculates the cotangent of the angle x given in degrees.
log1p(x[, out]) y=log1p(x) calculates log(1+x) for use when x is near zero.
expm1(x[, out]) y=expm1(x) calculates exp(x) - 1 for use when x is near zero.
cosm1(x[, out]) y=calculates cos(x) - 1 for use when x is near zero.
round(x[, out]) y=Returns the nearest integer to x as a double precision