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 print an error message when an error occurs. By default this printing is disabled. To enable such messages use errprint(1) 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 (from the cephesmodule). The output is the previous state. With errprint(0) no error messages are shown; the default is errprint(1). If no argument is given the current state of the flag is returned and no change occurs.
errstate with errstate(**state): –> operations in following block use given state.

Available functions

Airy functions

airy (x[, out1, out2, out3, ...) (Ai,Aip,Bi,Bip)=airy(z) calculates the Airy functions and their derivatives evaluated at real or complex number z. The Airy functions Ai and Bi are two independent solutions of y’‘(x)=xy. Aip and Bip are the first derivatives evaluated at x of Ai and Bi respectively.
airye (x[, out1, out2, out3, ...) (Aie,Aipe,Bie,Bipe)=airye(z) calculates the exponentially scaled Airy functions and their derivatives evaluated at real or complex number z. airye(z)[0:1] = airy(z)[0:1] * exp(2.0/3.0*z*sqrt(z)) airye(z)[2:3] = airy(z)[2:3] * exp(-abs((2.0/3.0*z*sqrt(z)).real))
ai_zeros (nt) Compute the zeros of Airy Functions Ai(x) and Ai’(x), a and a’ respectively, and the associated values of Ai(a’) and Ai’(a).
bi_zeros (nt) Compute the zeros of Airy Functions Bi(x) and Bi’(x), b and b’ respectively, and the associated values of Ai(b’) and Ai’(b).

Elliptic Functions and Integrals

ellipj (x1, x2[, out1, out2, ...) (sn,cn,dn,ph)=ellipj(u,m) calculates the Jacobian elliptic functions of parameter m between 0 and 1, and real u. The returned functions are often written sn(u|m), cn(u|m), and dn(u|m). The value of ph is such that if u = ellik(ph,m), then sn(u|m) = sin(ph) and cn(u|m) = cos(ph).
ellipk (x[, out]) y=ellipk(m) returns the complete integral of the first kind: integral(1/sqrt(1-m*sin(t)**2),t=0..pi/2)
ellipkinc (x1, x2[, out]) y=ellipkinc(phi,m) returns the incomplete elliptic integral of the first kind: integral(1/sqrt(1-m*sin(t)**2),t=0..phi)
ellipe (x[, out]) y=ellipe(m) returns the complete integral of the second kind: integral(sqrt(1-m*sin(t)**2),t=0..pi/2)
ellipeinc (x1, x2[, out]) y=ellipeinc(phi,m) returns the incomplete elliptic integral of the second kind: integral(sqrt(1-m*sin(t)**2),t=0..phi)

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 v at complex z: jve(v,z) = jv(v,z) * exp(-abs(z.imag))
yn (x1, x2[, out]) y=yn(n,x) returns the Bessel function of the second kind of integer order n at x.
yv (x1, x2[, out]) y=yv(v,z) returns the Bessel function of the second kind of real order v at complex z.
yve (x1, x2[, out]) y=yve(v,z) returns the exponentially scaled Bessel function of the second kind of real order v at complex z: yve(v,z) = yv(v,z) * exp(-abs(z.imag))
kn (x1, x2[, out]) y=kn(n,x) returns the modified Bessel function of the second kind (sometimes called the third kind) for integer order n at x.
kv (x1, x2[, out]) y=kv(v,z) returns the modified Bessel function of the second kind (sometimes called the third kind) for real order v at complex z.
kve (x1, x2[, out]) y=kve(v,z) returns the exponentially scaled, modified Bessel function of the second kind (sometimes called the third kind) for real order v at complex z: kve(v,z) = kv(v,z) * exp(z)
iv (x1, x2[, out]) y=iv(v,z) returns the modified Bessel function of real order v of z. If z is of real type and negative, v must be integer valued.
ive (x1, x2[, out]) y=ive(v,z) returns the exponentially scaled modified Bessel function of real order v and complex z: ive(v,z) = iv(v,z) * exp(-abs(z.real))
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 kind for real order v and complex argument z: hankel1e(v,z) = hankel1(v,z) * exp(-1j * z)
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 kind for real order v and complex argument z: hankel1e(v,z) = hankel1(v,z) * exp(1j * z)

The following is not an universal function:

lmbda (v, x) Compute sequence of lambda functions with arbitrary order v and their derivatives. Lv0(x)..Lv(x) are computed with v0=v-int(v).

Zeros of Bessel Functions

These are not universal functions:

jnjnp_zeros (nt) Compute nt (<=1200) zeros of the bessel functions Jn and Jn’ and arange them in order of their magnitudes.
jnyn_zeros (n, nt) Compute nt zeros of the Bessel functions Jn(x), Jn’(x), Yn(x), and Yn’(x), respectively. Returns 4 arrays of length nt.
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 of Y0’(z0) = -Y1(z0) at each zero.
y1_zeros (nt[, complex]) Returns nt (complex or real) zeros of Y1(z), z1, and the value of Y1’(z1) = Y0(z1) at each zero.
y1p_zeros (nt[, complex]) Returns nt (complex or real) zeros of Y1’(z), z1’, and the value of Y1(z1’) at each zero.

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 of order 0 at x. i0e(x) = exp(-|x|) * i0(x).
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 of order 0 at x. i1e(x) = exp(-|x|) * i1(x).
k0 (x[, out]) y=k0(x) returns the modified Bessel function of the second kind (sometimes called the third kind) of order 0 at x.
k0e (x[, out]) y=k0e(x) returns the exponentially scaled modified Bessel function of the second kind (sometimes called the third kind) of order 0 at x. k0e(x) = exp(x) * k0(x).
k1 (x[, out]) y=i1(x) returns the modified Bessel function of the second kind (sometimes called the third kind) of order 1 at x.
k1e (x[, out]) y=k1e(x) returns the exponentially scaled modified Bessel function of the second kind (sometimes called the third kind) of order 1 at x. k1e(x) = exp(x) * k1(x)

Integrals of Bessel Functions

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

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 all orders up to and including n.
sph_yn (n, z) Compute the spherical Bessel function yn(z) and its derivative for all orders up to and including n.
sph_jnyn (n, z) Compute the spherical Bessel functions, jn(z) and yn(z) and their derivatives for all orders up to and including n.
sph_in (n, z) Compute the spherical Bessel function in(z) and its derivative for all orders up to and including n.
sph_kn (n, z) Compute the spherical Bessel function kn(z) and its derivative for all orders up to and including n.
sph_inkn (n, z) Compute the spherical Bessel functions, in(z) and kn(z) and their derivatives for all orders up to and including n.

Ricatti-Bessel Functions

These are not universal functions:

riccati_jn (n, x) Compute the Ricatti-Bessel function of the first kind and its derivative for all orders up to and including n.
riccati_yn (n, x) Compute the Ricatti-Bessel function of the second kind and its derivative for all orders up to and including n.

Struve Functions

struve (x1, x2[, out]) y=struve(v,x) returns the Struve function Hv(x) of order v at x, x must be positive unless v is an integer.
modstruve (x1, x2[, out]) y=modstruve(v,x) returns the modified Struve function Lv(x) of order v at x, x must be positive unless v is an integer and it is recommended that |v|<=20.
itstruve0 (x[, out]) y=itstruve0(x) returns the integral of the Struve function of order 0 from 0 to x: integral(H0(t), t=0..x).
it2struve0 (x[, out]) y=it2struve0(x) returns the integral of the Struve function of order 0 divided by t from x to infinity: integral(H0(t)/t, t=x..inf).
itmodstruve0 (x[, out]) y=itmodstruve0(x) returns the integral of the modified Struve function of order 0 from 0 to x: integral(L0(t), t=0..x).

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 Binomial probability density: sum(nCj p**j (1-p)**(n-j),j=0..k)
bdtrc (x1, x2, x3[, out]) y=bdtrc(k,n,p) returns the sum of the terms k+1 through n of the Binomial probability density: sum(nCj p**j (1-p)**(n-j), j=k+1..n)
bdtri (x1, x2, x3[, out]) p=bdtri(k,n,y) finds the probability p such that the sum of the terms 0 through k of the Binomial probability density is equal to the given cumulative probability y.
btdtr (x1, x2, x3[, out]) y=btdtr(a,b,x) returns the area from zero to x under the beta density function: gamma(a+b)/(gamma(a)*gamma(b)))*integral(t**(a-1) (1-t)**(b-1), t=0..x). SEE ALSO betainc
btdtri (x1, x2, x3[, out]) x=btdtri(a,b,p) returns the pth quantile of the beta distribution. It is effectively the inverse of btdtr returning the value of x for which btdtr(a,b,x) = p. SEE ALSO betaincinv
fdtr (x1, x2, x3[, out]) y=fdtr(dfn,dfd,x) returns the area from zero to x under the F density function (also known as Snedcor’s density or the variance ratio density). This is the density of X = (unum/dfn)/(uden/dfd), where unum and uden are random variables having Chi square distributions with dfn and dfd degrees of freedom, respectively.
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 fdtr(dfn,dfd,x)=p.
gdtr (x1, x2, x3[, out]) y=gdtr(a,b,x) returns the integral from zero to x of the gamma probability density function: a**b / gamma(b) * integral(t**(b-1) exp(-at),t=0..x). The arguments a and b are used differently here than in other definitions.
gdtrc (x1, x2, x3[, out]) y=gdtrc(a,b,x) returns the integral from x to infinity of the gamma probability density function. SEE gdtr, gdtri
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 negative binomial distribution: sum((n+j-1)Cj p**n (1-p)**j,j=0..k). In a sequence of Bernoulli trials this is the probability that k or fewer failures precede the nth success.
nbdtrc (x1, x2, x3[, out]) y=nbdtrc(k,n,p) returns the sum of the terms k+1 to infinity of the negative binomial distribution.
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 distribution: sum(exp(-m) * m**j / j!, j=0..k) = gammaincc( k+1, m). Arguments must both be positive and k an integer.
pdtrc (x1, x2[, out]) y=pdtrc(k,m) returns the sum of the terms from k+1 to infinity of the Poisson distribution: sum(exp(-m) * m**j / j!, j=k+1..inf) = gammainc( k+1, m). Arguments must both be positive and k an integer.
pdtri (x1, x2[, out]) m=pdtri(k,y) returns the Poisson variable m such that the sum from 0 to k of the Poisson density is equal to the given probability y: calculated by gammaincinv( k+1, y). k must be a nonnegative integer and y between 0 and 1.
stdtr (x1, x2[, out]) p=stdtr(df,t) returns the integral from minus infinity to t of the Student t distribution with df > 0 degrees of freedom: gamma((df+1)/2)/(sqrt(df*pi)*gamma(df/2)) * integral((1+x**2/df)**(-df/2-1/2), x=-inf..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 square probability density function with v degrees of freedom: 1/(2**(v/2) * gamma(v/2)) * integral(t**(v/2-1) * exp(-t/2), t=0..x)
chdtrc (x1, x2[, out]) p=chdtrc(v,x) returns the area under the right hand tail (from x to infinity) of the Chi square probability density function with v degrees of freedom: 1/(2**(v/2) * gamma(v/2)) * integral(t**(v/2-1) * exp(-t/2), t=x..inf)
chdtri (x1, x2[, out]) x=chdtri(v,p) returns the argument x such that chdtrc(v,x) is equal to p.
ndtr (x[, out]) y=ndtr(x) returns the area under the standard Gaussian probability density function, integrated from minus infinity to x: 1/sqrt(2*pi) * integral(exp(-t**2 / 2),t=-inf..x)
ndtri (x[, out]) x=ndtri(y) returns the argument x for which the area udnder the Gaussian probability density function (integrated from minus infinity to x) is equal to y.
smirnov (x1, x2[, out]) y=smirnov(n,e) returns the exact Kolmogorov-Smirnov complementary cumulative distribution function (Dn+ or Dn-) for a one-sided test of equality between an empirical and a theoretical distribution. It is equal to the probability that the maximum difference between a theoretical distribution and an empirical one based on n samples is greater than e.
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 function of Kolmogorov’s limiting distribution (Kn* for large n) of a two-sided test for equality between an empirical and a theoretical distribution. It is equal to the (limit as n->infinity of the) probability that sqrt(n) * max absolute deviation > y.
kolmogi (x[, out]) y=kolmogi(p) returns y such that kolmogorov(y) = p
tklmbda (x1, x2[, out])

Error Function and Fresnel Integrals

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

These are not universal functions:

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 m and nonnegative degree v: |x|<=1.
sph_harm () Compute spherical harmonics.

These are not universal functions:

lpn (n, z) Compute sequence of Legendre functions of the first kind (polynomials), Pn(z) and derivatives for all degrees from 0 to n (inclusive).
lqn (n, z) Compute sequence of Legendre functions of the second kind, Qn(z) and derivatives for all degrees from 0 to n (inclusive).
lpmn (m, n, z) Associated Legendre functions of the first kind, Pmn(z) and its derivative, Pmn’(z) of order m and degree n. Returns two arrays of size (m+1,n+1) containing Pmn(z) and Pmn’(z) for all orders from 0..m and degrees from 0..n.
lqmn (m, n, z) Associated Legendre functions of the second kind, Qmn(z) and its derivative, Qmn’(z) of order m and degree n. Returns two arrays of size (m+1,n+1) containing Qmn(z) and Qmn’(z) for all orders from 0..m and degrees from 0..n.

Orthogonal polynomials

These functions all return a polynomial class which can then be evaluated: vals = chebyt(n)(x).

The class also has an attribute ‘weights’ which return 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.

Warning

Evaluating large-order polynomials using these functions can be numerically unstable.

The reason is that the functions below return polynomials as numpy.poly1d objects, which represent the polynomial in terms of their coefficients, and this can result to loss of precision when the polynomial terms are summed.

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

Hypergeometric Functions

hyp2f1 (x1, x2, x3, x4[, out]) y=hyp2f1(a,b,c,z) returns the gauss hypergeometric function ( 2F1(a,b;c;z) ).
hyp1f1 (x1, x2, x3[, out]) y=hyp1f1(a,b,x) returns the confluent hypergeometeric function ( 1F1(a,b;x) ) evaluated at the values a, b, and x.
hyperu (x1, x2, x3[, out]) y=hyperu(a,b,x) returns the confluent hypergeometric function of the second kind U(a,b,x).
hyp0f1 (v, z) Confluent hypergeometric limit function 0F1. Limit as q->infinity of 1F1(q;a;z/q)
hyp2f0 (x1, x2, x3, x4[, out1, ...) (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 can be either 1 or 2.
hyp1f2 (x1, x2, x3, x4[, out1, ...) (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, ...) (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 d and the derivative, Dv’(x) in dp.
pbvv (x1, x2[, out1, out2]) (v,vp)=pbvv(v,x) returns (v,vp) with the parabolic cylinder function Vv(x) in v and the derivative, Vv’(x) in vp.
pbwa (x1, x2[, out1, out2]) (w,wp)=pbwa(a,x) returns (w,wp) with the parabolic cylinder function W(a,x) in w and the derivative, W’(a,x) in wp. May not be accurate for large (>5) arguments in a and/or x.

These are not universal functions:

pbdv_seq (v, x) Compute sequence of parabolic cylinder functions Dv(x) and their derivatives for Dv0(x)..Dv(x) with v0=v-int(v).
pbvv_seq (v, x) Compute sequence of parabolic cylinder functions Dv(x) and their derivatives for Dv0(x)..Dv(x) with v0=v-int(v).
pbdn_seq (n, z) Compute sequence of parabolic cylinder functions Dn(z) and their derivatives for D0(z)..Dn(z).

Spheroidal Wave Functions

pro_ang1 (x1, x2, x3, x4[, out1, ...) (s,sp)=pro_ang1(m,n,c,x) computes the prolate sheroidal angular function of the first kind and its derivative (with respect to x) for mode paramters m>=0 and n>=m, spheroidal parameter c and |x|<1.0.
pro_rad1 (x1, x2, x3, x4[, out1, ...) (s,sp)=pro_rad1(m,n,c,x) computes the prolate sheroidal radial function of the first kind and its derivative (with respect to x) for mode paramters m>=0 and n>=m, spheroidal parameter c and |x|<1.0.
pro_rad2 (x1, x2, x3, x4[, out1, ...) (s,sp)=pro_rad2(m,n,c,x) computes the prolate sheroidal radial function of the second kind and its derivative (with respect to x) for mode paramters m>=0 and n>=m, spheroidal parameter c and |x|<1.0.
obl_ang1 (x1, x2, x3, x4[, out1, ...) (s,sp)=obl_ang1(m,n,c,x) computes the oblate sheroidal angular function of the first kind and its derivative (with respect to x) for mode paramters m>=0 and n>=m, spheroidal parameter c and |x|<1.0.
obl_rad1 (x1, x2, x3, x4[, out1, ...) (s,sp)=obl_rad1(m,n,c,x) computes the oblate sheroidal radial function of the first kind and its derivative (with respect to x) for mode paramters m>=0 and n>=m, spheroidal parameter c and |x|<1.0.
obl_rad2 (x1, x2, x3, x4[, out1, ...) (s,sp)=obl_rad2(m,n,c,x) computes the oblate sheroidal radial function of the second kind and its derivative (with respect to x) for mode paramters m>=0 and n>=m, spheroidal parameter c and |x|<1.0.
pro_cv (x1, x2, x3[, out]) cv=pro_cv(m,n,c) computes the characteristic value of prolate spheroidal wave functions of order m,n (n>=m) and spheroidal parameter c.
obl_cv (x1, x2, x3[, out]) cv=obl_cv(m,n,c) computes the characteristic value of oblate spheroidal wave functions of order m,n (n>=m) and spheroidal parameter c.
pro_cv_seq (m, n, c) Compute a sequence of characteristic values for the prolate spheroidal wave functions for mode m and n’=m..n and spheroidal parameter c.
obl_cv_seq (m, n, c) Compute a sequence of characteristic values for the oblate spheroidal wave functions for mode m and n’=m..n and spheroidal parameter c.

The following functions require pre-computed characteristic value:

pro_ang1_cv (x1, x2, x3, x4, x5[, ...) (s,sp)=pro_ang1_cv(m,n,c,cv,x) computes the prolate sheroidal angular function of the first kind and its derivative (with respect to x) for mode paramters m>=0 and n>=m, spheroidal parameter c and |x|<1.0. Requires pre-computed characteristic value.
pro_rad1_cv (x1, x2, x3, x4, x5[, ...) (s,sp)=pro_rad1_cv(m,n,c,cv,x) computes the prolate sheroidal radial function of the first kind and its derivative (with respect to x) for mode paramters m>=0 and n>=m, spheroidal parameter c and |x|<1.0. Requires pre-computed characteristic value.
pro_rad2_cv (x1, x2, x3, x4, x5[, ...) (s,sp)=pro_rad2_cv(m,n,c,cv,x) computes the prolate sheroidal radial function of the second kind and its derivative (with respect to x) for mode paramters m>=0 and n>=m, spheroidal parameter c and |x|<1.0. Requires pre-computed characteristic value.
obl_ang1_cv (x1, x2, x3, x4, x5[, ...) (s,sp)=obl_ang1_cv(m,n,c,cv,x) computes the oblate sheroidal angular function of the first kind and its derivative (with respect to x) for mode paramters m>=0 and n>=m, spheroidal parameter c and |x|<1.0. Requires pre-computed characteristic value.
obl_rad1_cv (x1, x2, x3, x4, x5[, ...) (s,sp)=obl_rad1_cv(m,n,c,cv,x) computes the oblate sheroidal radial function of the first kind and its derivative (with respect to x) for mode paramters m>=0 and n>=m, spheroidal parameter c and |x|<1.0. Requires pre-computed characteristic value.
obl_rad2_cv (x1, x2, x3, x4, x5[, ...) (s,sp)=obl_rad2_cv(m,n,c,cv,x) computes the oblate sheroidal radial function of the second kind and its derivative (with respect to x) for mode paramters m>=0 and n>=m, spheroidal parameter c and |x|<1.0. Requires pre-computed characteristic value.

Kelvin Functions

kelvin (x[, out1, out2, out3, ...) (Be, Ke, Bep, Kep)=kelvin(x) returns the tuple (Be, Ke, Bep, Kep) which containes complex numbers representing the real and imaginary Kelvin functions and their derivatives evaluated at x. For example, kelvin(x)[0].real = ber x and kelvin(x)[0].imag = bei x with similar relationships for ker and kei.
kelvin_zeros (nt) Compute nt zeros of all the kelvin functions returned in a length 8 tuple of arrays of length nt. The tuple containse the arrays of zeros of (ber, bei, ker, kei, ber’, bei’, ker’, kei’)
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 non-negative x and n: integral(exp(-x*t) / t**n, t=1..inf).
exp1 (x[, out]) y=exp1(z) returns the exponential integral (n=1) of complex argument z: integral(exp(-z*t)/t,t=1..inf).
expi (x[, out]) y=expi(x) returns an exponential integral of argument x defined as integral(exp(t)/t,t=-inf..x). See expn for a different exponential integral.
wofz (x[, out]) y=wofz(z) returns the value of the fadeeva function for complex argument z: exp(-z**2)*erfc(-i*z)
dawsn (x[, out]) y=dawsn(x) returns dawson’s integral: exp(-x**2) * integral(exp(t**2),t=0..x).
shichi (x[, out1, out2]) (shi,chi)=shichi(x) returns the hyperbolic sine and cosine integrals: integral(sinh(t)/t,t=0..x) and eul + ln x + integral((cosh(t)-1)/t,t=0..x) where eul is Euler’s Constant.
sici (x[, out1, out2]) (si,ci)=sici(x) returns in si the integral of the sinc function from 0 to x: integral(sin(t)/t,t=0..x). It returns in ci the cosine integral: eul + ln x + integral((cos(t) - 1)/t,t=0..x).
spence (x[, out]) y=spence(x) returns the dilogarithm integral: -integral(log t / (t-1),t=1..x)
zeta (x1, x2[, out]) y=zeta(x,q) returns the Riemann zeta function of two arguments: sum((k+q)**(-x),k=0..inf)
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 (s)econds in radians.
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 floating point result. If x ends in 0.5 exactly, the nearest even integer is chosen.