Skip to main content

Special Functions

The functions in this section take a number argument, the type that includes \pm\infty and NaN alongside the finite numbers. They are not restricted to finite input the way the number theory functions are:

  • An infinite argument is accepted. Where the function has a limit there, that limit is the value — \operatorname{erf}(\infty) = 1, \operatorname{erfc}(\infty) = 0. Otherwise the expression stays symbolic, as ["Gamma", "PositiveInfinity"] does.
  • A NaN argument gives NaN. This happens under plain evaluate(), not only under N(): NaN is not an exact value, so there is nothing to hold symbolically.

A function may still be infinite at a finite point — K(1) = \infty — in which case the result is PositiveInfinity.

Erf(z:number)

Evaluate to the error function of a complex number.

The error function is an odd function ( 

 \operatorname{erf} -z = -
\operatorname{erf} z
) that is used in statistics to calculate probabilities of normally distributed events.

The formula for the error function of a complex number is:

\operatorname{erf} z = \frac{2}{\sqrt{\pi}} \int_0^z e^{-t^2} dt

where z is a complex number.

Erfc(z:number)

Evaluate to the complementary error function of a complex number.

It is defined as \operatorname{erfc} z = 1 - \operatorname {erf} z.

ErfInv(x:number)

Evaluate to the inverse error function of a real number -1 < x < 1

Outside that interval the value is not defined: ["ErfInv", 2] and ["ErfInv", "PositiveInfinity"] both evaluate to NaN.

It is defined as 

\operatorname{erf} \left(\operatorname{erf} ^{-1}x\right)
= x
.

Factorial(n)

n!
$$$n!$$
["Factorial", 5]
// -> 120

For a non-integer n, n! is \Gamma(n+1). An exact non-integer argument stays symbolic under evaluate(), and .N() or a float argument gives the value to ce.precision digits:

["Factorial", ["Rational", 5, 2]]
// -> (5/2)!

["Factorial", 2.5]
// -> 3.32335097044784255118

Above machine precision, a factorial too large for an exact integer is a big decimal under .N():

["N", ["Factorial", ["Power", 10, 6]]]
// -> 8.26393168833124006238e+5565708

Factorial2(n)

The double factorial of n: 

 n!! = n \cdot (n-2) \cdot (n-4) \times
\cdots
, that is the product of all the positive integers up to n that have the same parity (odd or even) as n.

n!!
$$$n!!$$
["Factorial2", 5]
// -> 15

It can also be written in terms of the \Gamma function:

n!! = 2^{\frac{n}{2}+\frac{1}{4}(1-\cos(\pi n))}\pi^{\frac{1}{4}(\cos(\pi
n)-1)}\Gamma\left(\frac{n}{2}+1\right)

This is not the same as the factorial of the factorial of n (i.e. ((n!)!)).

Reference

Gamma(z)

\Gamma(n) = (n-1)!
$$$\Gamma(n) = (n-1)!$$

The Gamma Function is an extension of the factorial function, with its argument shifted by 1, to real and complex numbers.

\operatorname{\Gamma}\left(z\right) = \int\limits_{0}^{\infty} t^{z-1}
\mathrm{e}^{-t} \, \mathrm{d}t
["Gamma", 5]
// 24

Gamma(s, z)

\Gamma(s, z)
$$$\Gamma(s, z)$$

With two arguments, Gamma is the upper incomplete Gamma function. The lower limit of the integral is the second argument z instead of 0:

\operatorname{\Gamma}\left(s, z\right) = \int\limits_{z}^{\infty} t^{s-1}
\mathrm{e}^{-t} \, \mathrm{d}t

The order s and the lower limit z may be real or complex, including negative and fractional orders. The two-argument form is evaluated numerically with .N() and otherwise stays symbolic; \Gamma(s, 0) reduces to \Gamma(s). (This matches the Gamma[s, z] convention of Mathematica.)

["N", ["Gamma", 2, 1]]
// 0.7357588823428849

GammaLn(z)

\ln(\Gamma(z))
$$$\ln(\Gamma(z))$$

GammaLn(z) is the principal logarithm of \Gamma(z), as gammaln in MATLAB and SciPy. It is not the LogGamma of Mathematica: that function is the analytic continuation of \ln\Gamma, which the Compute Engine provides as LogGamma. The two agree for a positive real argument and differ by a multiple of 2\pi i elsewhere:

["N", ["GammaLn", -1.5]]
// -> 0.860047015376481014511

["N", ["LogGamma", -1.5]]
// -> 0.8600470153764792 - 6.283185307179586i

LogGamma(z)

The analytic continuation of \ln\Gamma(z), with its branch cut on (-\infty, 0]. This is the LogGamma of Mathematica and the loggamma of mpmath. It is exact where \Gamma is exact, and the poles at the non-positive integers are +\infty. A real argument is computed to ce.precision digits. A complex argument is computed in double precision.

["LogGamma", 5]
// -> 3ln(2) + ln(3)

["LogGamma", ["Rational", 1, 2]]
// -> 1/2 * ln(pi)

["N", ["LogGamma", ["Complex", -2.5, 1.5]]]
// -> -3.7175134511917927 - 7.713065525834192i

PolyGamma(-1, z) is LogGamma(z).

BarnesG(z)

The Barnes G-function, with G(z+1) = \Gamma(z)\,G(z) and G(1) = 1. At a positive integer n, BarnesG(n) is the superfactorial \prod_{k=0}^{n-2} k!, an exact integer. A real argument is computed to ce.precision digits. A complex argument is computed in double precision.

["BarnesG", 5]
// -> 12

["BarnesG", 10]
// -> 5056584744960000

["N", ["BarnesG", ["Rational", 1, 2]]]
// -> 0.603244281209446206191

LogBarnesG(z)

The logarithm of the Barnes G-function, continued with LogGamma as in Mathematica. LogBarnesG(0) is -\infty.

["N", ["LogBarnesG", -2.5]]
// -> -2.574748476853131 + 18.84955592153876i

Zeta(s)

\zeta(s)
$$$\zeta(s)$$

The Riemann zeta function, defined for complex numbers with real part greater than 1 as:

\zeta(s) = \sum_{n=1}^{\infty} \frac{1}{n^s}

and extended to other values by analytic continuation.

["Zeta", 2]
// ➔ π²/6

Zeta(s, a)

The two-argument form is the generalized zeta function, as Mathematica defines it:

\zeta(s, a) = \sum_{n=0}^{\infty} \frac{1}{|n + a|^s}

with the term where n + a = 0 omitted. For a > 0 it is equal to the Hurwitz zeta function HurwitzZeta(s, a), which sums (n + a)^{-s} with the principal power. For a \le 0 the two differ: Zeta(2.976, -2.664) is the real 29.9024, while HurwitzZeta(2.976, -2.664) is 22.6003 - 0.2754i.

DirichletEta(s)

\eta(s)
$$$\eta(s)$$

The Dirichlet eta function, the alternating form of \zeta:

\eta(s) = \sum_{n=1}^{\infty} \frac{(-1)^{n-1}}{n^s} = (1 - 2^{1-s})\,\zeta(s)

It is entire. The values at the integers are exact. A real s is computed to ce.precision digits, also next to s = 1. A complex s is computed in double precision.

["DirichletEta", 2]
// -> π²/12

["DirichletEta", 1]
// -> ln(2)

["N", ["DirichletEta", ["Rational", 1, 2]]]
// -> 0.604898643421630370247

DirichletBeta(s)

\beta(s)
$$$\beta(s)$$

The Dirichlet beta function:

\beta(s) = \sum_{n=0}^{\infty} \frac{(-1)^n}{(2n+1)^s}

It is entire. The values at the integers are exact. A real s is computed to ce.precision digits. A complex s is computed in double precision.

["DirichletBeta", 1]
// -> π/4

["DirichletBeta", 2]
// -> CatalanConstant

["DirichletBeta", 5]
// -> 5π⁵/1536

StieltjesGamma(n)

StieltjesGamma(n, a)

\gamma_n(a)
$$$\gamma_n(a)$$

The generalized Stieltjes constant \gamma_n(a), the coefficient of the Laurent series of the Hurwitz zeta function at its pole:

\zeta(s, a) = \frac{1}{s - 1} + \sum_{n=0}^{\infty} \frac{(-1)^n}{n!}\,\gamma_n(a)\,(s - 1)^n

StieltjesGamma(n) is \gamma_n(1), the Stieltjes constant of the Riemann zeta function. StieltjesGamma(0) is EulerGamma, and StieltjesGamma(0, a) is -PolyGamma(0, a). A non-positive integer a is a pole. A real a > 0 is computed to ce.precision digits. A complex a or a negative non-integer a is computed in double precision. An order above 30 stays unevaluated.

["StieltjesGamma", 0]
// -> EulerGamma

["N", ["StieltjesGamma", 1]]
// -> -0.0728158454836767248606

["StieltjesGamma", 0, "a"]
// -> -PolyGamma(0, a)

ClausenCl(n, θ)

\operatorname{Cl}_n(\theta)
$$$\operatorname{Cl}_n(\theta)$$

The Clausen function of integer order n \ge 1 and real \theta:

\operatorname{Cl}_n(\theta) = \sum_{k=1}^{\infty} \frac{\sin(k\theta)}{k^n}
\text{ for even } n, \qquad
\operatorname{Cl}_n(\theta) = \sum_{k=1}^{\infty} \frac{\cos(k\theta)}{k^n}
\text{ for odd } n

These are the clsin and clcos functions of mpmath. The values at the special points are exact. A real \theta is computed to ce.precision digits for an order up to 40. A non-integer or non-positive order, and a symbolic \theta, stay unevaluated.

["ClausenCl", 2, ["Divide", "Pi", 2]]
// -> CatalanConstant

["ClausenCl", 3, 0]
// -> Zeta(3)

["N", ["ClausenCl", 2, 1]]
// -> 1.01395913236076850429

DirichletCharacter(k, j, n)

The Dirichlet character \chi_{k,j}(n) modulo k, in the indexing of the Wolfram Language: j = 1 is the principal character, and j runs up to \varphi(k). The value is 0 when \gcd(n, k) > 1. The function threads over a list.

["DirichletCharacter", 5, 2, 2]
// -> i

["DirichletCharacter", 7, 3, 3]
// -> e^(2/3i * pi)

["DirichletCharacter", 5, 2, 5]
// -> 0

DirichletL(k, j, s)

The Dirichlet L-function L(s, \chi_{k,j}) of the character DirichletCharacter(k, j, ·). The value at a non-positive integer s is exact. For a real character, the value at s = 1 is exact (from the class number formula). The principal character has a pole at s = 1. A real s of a real character is computed to ce.precision digits. A character with complex values is computed in double precision. A modulus above 1000 stays unevaluated.

["DirichletL", 4, 2, 1]
// -> π/4

["DirichletL", 3, 2, 1]
// -> sqrt(3)/9 * pi

["DirichletL", 5, 1, 1]
// -> ~oo

["N", ["DirichletL", 3, 2, 1.5]]
// -> 0.703968244868733261668

Beta(a, b)

\Beta(a, b)
$$$\Beta(a, b)$$

The Euler beta function, defined as:

\operatorname{B}(a, b) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}

It can also be expressed as an integral:

\operatorname{B}(a, b) = \int_0^1 t^{a-1}(1-t)^{b-1} \, dt
["Beta", 2, 3]
// ➔ 1/12

GammaRegularized(a, z)

The upper regularized incomplete gamma function, defined as:

Q(a, z) = \frac{\Gamma(a, z)}{\Gamma(a)} = \frac{1}{\Gamma(a)}\int_z^\infty t^{a-1}e^{-t} \, dt

Exact arguments stay symbolic, and special values fold: GammaRegularized(a, 0) is 1 and GammaRegularized(1, z) is e^{-z}. Use N() for a numeric value, at machine or arbitrary precision.

The cumulative distribution function of a Poisson distribution evaluates to this function: 

\operatorname{CDF}(\operatorname{Poisson}(\lambda), k) =
Q(k+1, \lambda)
.

["N", ["GammaRegularized", 3, 5]]
// ➔ 0.12465201948308115

BetaRegularized(x, a, b)

The regularized incomplete beta function, defined as:

I_x(a, b) = \frac{1}{\operatorname{B}(a, b)}\int_0^x t^{a-1}(1-t)^{b-1} \, dt

Exact arguments stay symbolic, and the endpoints fold: BetaRegularized(0, a, b) is 0 and BetaRegularized(1, a, b) is 1. Use N() for a numeric value, at machine or arbitrary precision.

The cumulative distribution function of a binomial distribution evaluates to this function: 

\operatorname{CDF}(\operatorname{Binomial}(n, p), k) =
I_{1-p}(n-k, k+1)
.

["N", ["BetaRegularized", 0.5, 2, 3]]
// ➔ 0.6875

LambertW(x)

LambertW(x, k)

\operatorname{W}(x)
$$$\operatorname{W}(x)$$
\operatorname{W}_{k}(x)
$$$\operatorname{W}_{k}(x)$$

The Lambert W function, also called the product logarithm. It is the inverse function of f(w) = w e^w.

For a given value x, W_k(x) is the value w on branch k such that w e^w = x. The branch k is an integer. When it is not given, it is 0, the principal branch W_0.

W_k(x) is real only on the two real branches:

  • W_0(x) for x \ge -1/e
  • W_{-1}(x) for -1/e \le x < 0

On all other branches, and for all other arguments, the value is complex.

Argument order. The branch is the second argument: ["LambertW", x, k]. This is the order of mpmath (lambertw(z, k)), SciPy, SymPy, Julia and Fungrim. Mathematica (ProductLog[k, z], or LambertW[k, z] in Wolfram|Alpha), Maple, Sage and MATLAB put the branch first. When you copy an expression from one of these systems, swap the two arguments.

A branch that is not an integer is a type error. For example, ["LambertW", -1, -0.1] (the Wolfram order, copied without a swap) gives an incompatible-type error. But when both arguments are integers, a swap is not an error: ["LambertW", 1, 2] is W_2(1) \approx -2.4016 + 10.7763i, and ["LambertW", 2, 1] is W_1(2) \approx -0.8343 + 4.5303i.

To make the order clear, name the branch. The parameters are named z and branch, and named arguments can be in any order:

["N", ["LambertW", -0.1, ["NamedArgument", "'branch'", -1]]]
// ➔ -3.57715206395729721841

["N", ["LambertW", ["NamedArgument", "'branch'", -1], ["NamedArgument", "'z'", -0.1]]]
// ➔ -3.57715206395729721841

In Epsil, this is lambertW(-0.1, branch: -1).

The derivative of the Lambert W function is:

\frac{d}{dx} W_k(x) = \frac{W_k(x)}{x(1 + W_k(x))}
["LambertW", 1]
// ➔ Ω ≈ 0.5671 (the Omega constant)

["N", ["LambertW", -0.1, -1]]
// ➔ -3.57715206395729721841

Fresnel Integrals​

The Fresnel integrals arise in the description of near-field diffraction and in the geometry of the Cornu spiral (Euler spiral).

FresnelS(x)

\operatorname{FresnelS}(x)
$$$\operatorname{FresnelS}(x)$$

The Fresnel S integral:

S(x) = \int_0^x \sin\!\left(\frac{\pi t^2}{2}\right) dt

It is an odd function (S(-x) = -S(x)) with asymptotic value S(\infty) = \tfrac{1}{2}.

["FresnelS", 1]
// ➔ 0.4383

FresnelC(x)

\operatorname{FresnelC}(x)
$$$\operatorname{FresnelC}(x)$$

The Fresnel C integral:

C(x) = \int_0^x \cos\!\left(\frac{\pi t^2}{2}\right) dt

It is an odd function (C(-x) = -C(x)) with asymptotic value C(\infty) = \tfrac{1}{2}.

["FresnelC", 1]
// ➔ 0.7799

Bessel Functions​

Bessel functions are solutions to Bessel's differential equation:

x^2 \frac{d^2 y}{dx^2} + x \frac{dy}{dx} + (x^2 - n^2)y = 0

They arise in problems with cylindrical or spherical symmetry.

BesselJ, BesselY, BesselI and BesselK of an integer order evaluate numerically at a real or a complex argument under .N(), and under evaluate() when the argument is a float. On the negative real axis, their branch cut, BesselY and BesselK take the limit from above, as Ln(-2) does. A non-integer order stays symbolic. Compiled code computes these functions for a real argument only.

["N", ["BesselJ", 3, ["Complex", 0, 0.7]]]
// ➔ -0.007367373607628006i

["N", ["BesselY", 0, -2]]
// ➔ 0.5103756726497453 + 0.44778155828247124i

["N", ["BesselJ", ["Rational", 1, 2], 1]]
// ➔ ["BesselJ", 0.5, 1]  (non-integer order: stays symbolic)

BesselJ(n, x)

J_n(x)
$$$J_n(x)$$

The Bessel function of the first kind of order n.

The derivative with respect to x is:

\frac{d}{dx} J_n(x) = \frac{1}{2}(J_{n-1}(x) - J_{n+1}(x))
["BesselJ", 0, 1]
// ➔ J₀(1) ≈ 0.7652

BesselY(n, x)

Y_n(x)
$$$Y_n(x)$$

The Bessel function of the second kind of order n, also called the Neumann function.

The derivative with respect to x is:

\frac{d}{dx} Y_n(x) = \frac{1}{2}(Y_{n-1}(x) - Y_{n+1}(x))
["BesselY", 0, 1]
// ➔ Y₀(1) ≈ 0.0883

BesselI(n, x)

I_n(x)
$$$I_n(x)$$

The modified Bessel function of the first kind of order n.

The derivative with respect to x is:

\frac{d}{dx} I_n(x) = \frac{1}{2}(I_{n-1}(x) + I_{n+1}(x))
["BesselI", 0, 1]
// ➔ I₀(1) ≈ 1.2661

BesselK(n, x)

K_n(x)
$$$K_n(x)$$

The modified Bessel function of the second kind of order n, also called the MacDonald function.

The derivative with respect to x is:

\frac{d}{dx} K_n(x) = -\frac{1}{2}(K_{n-1}(x) + K_{n+1}(x))
["BesselK", 0, 1]
// ➔ K₀(1) ≈ 0.4210

Airy Functions​

Airy functions are solutions to the Airy differential equation:

\frac{d^2 y}{dx^2} - xy = 0

They arise in physics, particularly in quantum mechanics and optics.

AiryAi, AiryBi, AiryAiPrime and AiryBiPrime evaluate numerically at a real or a complex argument under .N(), and under evaluate() when the argument is a float: ["N", ["AiryAi", ["Complex", 1, 1]]] is 0.060458308371838146 - 0.15188956587718141i. Compiled code computes them for a real argument only.

AiryAi(x)

\operatorname{Ai}(x)
$$$\operatorname{Ai}(x)$$

The Airy function of the first kind.

It is the solution to the Airy equation that decays exponentially for positive x and oscillates for negative x.

["AiryAi", 0]
// ➔ 1/(3^(2/3) Γ(2/3)) ≈ 0.3550

AiryBi(x)

\operatorname{Bi}(x)
$$$\operatorname{Bi}(x)$$

The Airy function of the second kind.

It is the solution to the Airy equation that grows exponentially for positive x and oscillates for negative x.

["AiryBi", 0]
// ➔ 1/(3^(1/6) Γ(2/3)) ≈ 0.6149

AiryAiPrime(x)

\operatorname{Ai}'(x)
$$$\operatorname{Ai}'(x)$$

The derivative of the Airy function of the first kind.

It satisfies \operatorname{Ai}''(x) = x\operatorname{Ai}(x), so differentiating an Airy expression stays within the Airy family.

["AiryAiPrime", 0]
// ➔ −1/(3^(1/3) Γ(1/3)) ≈ −0.2588

AiryBiPrime(x)

\operatorname{Bi}'(x)
$$$\operatorname{Bi}'(x)$$

The derivative of the Airy function of the second kind.

It satisfies \operatorname{Bi}''(x) = x\operatorname{Bi}(x).

["AiryBiPrime", 0]
// ➔ 3^(1/6)/Γ(1/3) ≈ 0.4483

Elliptic Integrals​

The complete elliptic integrals arise in computing the arc length of an ellipse, the period of a pendulum, and throughout number theory.

Convention: these functions use the parameter m = k^2, where k is the modulus. This matches Mathematica, mpmath and Fungrim. To evaluate in terms of the modulus k, pass k^2.

EllipticK(m)

K(m)
$$$K(m)$$

The complete elliptic integral of the first kind:

K(m) = \int_0^{\pi/2} \frac{d\theta}{\sqrt{1 - m \sin^2\theta}}

It is computed via the arithmetic-geometric mean: 

K(m) =
\dfrac{\pi}{2\operatorname{agm}(1, \sqrt{1-m})}
.

K(1) = \infty. For m > 1 the value is complex.

["EllipticK", 0.5]
// ➔ 1.85407467730137

EllipticE(m)

E(m)
$$$E(m)$$

The complete elliptic integral of the second kind:

E(m) = \int_0^{\pi/2} \sqrt{1 - m \sin^2\theta} \, d\theta

The perimeter of an ellipse with semi-major axis a and eccentricity e is 4aE(e^2).

E(1) = 1. For m > 1 the value is complex.

["EllipticE", 0.5]
// ➔ 1.35064388104768

AGM(a, b)

AGM(z)

\operatorname{agm}(a, b)
$$$\operatorname{agm}(a, b)$$

The arithmetic-geometric mean of two numbers: the common limit of the sequences 

a_{n+1} = \frac{a_n +
b_n}{2}
and b_{n+1} = \sqrt{a_n b_n}.

With a single argument, \operatorname{agm}(z) is shorthand for \operatorname{agm}(1, z).

["AGM", 1, 2]
// ➔ 1.45679103104691

Hypergeometric Functions​

Hypergeometric2F1(a, b, c, z)

{}_2F_1(a, b; c; z)
$$${}_2F_1(a, b; c; z)$$

The Gauss hypergeometric function, defined for |z| < 1 by the series:

{}_2F_1(a, b; c; z) = \sum_{n=0}^{\infty} \frac{(a)_n (b)_n}{(c)_n}
\frac{z^n}{n!}

where (q)_n is the Pochhammer symbol (rising factorial), and extended elsewhere by analytic continuation.

Many elementary and special functions are particular cases. For example \ln(1+z) = z \cdot {}_2F_1(1, 1; 2; -z) and K(m) = \frac{\pi}{2} \, {}_2F_1\big(\frac12, \frac12; 1; m\big).

If a or b is a non-positive integer the series terminates and the function is a polynomial in z, evaluable for any z. Otherwise the function evaluates numerically for z \le 1 (real) and for complex z within the unit disk and the Pfaff-transform region.

["Hypergeometric2F1", 1, 1, 2, 0.5]
// ➔ 1.38629436111989  (= 2 ln 2)

Hypergeometric1F1(a, b, z)

{}_1F_1(a; b; z)
$$${}_1F_1(a; b; z)$$

The Kummer confluent hypergeometric function, also written M(a, b, z):

{}_1F_1(a; b; z) = \sum_{n=0}^{\infty} \frac{(a)_n}{(b)_n} \frac{z^n}{n!}

It is an entire function of z and evaluates numerically for any real or complex z.

["Hypergeometric1F1", 1, 2, 2]
// ➔ 3.19452804946533  (= (e² − 1)/2)

Theta and Modular Functions​

JacobiTheta(j, z, τ)

\theta_j(z, \tau)
$$$\theta_j(z, \tau)$$

The Jacobi theta functions \theta_j(z, \tau) for j \in \{1, 2, 3, 4\}, with nome q = e^{i\pi\tau} and \operatorname{Im}(\tau) > 0:

\theta_3(z, \tau) = 1 + 2\sum_{n=1}^{\infty} q^{n^2} \cos(2n\pi z)

and similarly for \theta_1, \theta_2, \theta_4.

Convention: the trigonometric argument is a multiple of \pi z (the functions have period 1 in z), matching Fungrim and mpmath — not the classical convention with period \pi.

An optional fourth argument indicates the order of differentiation with respect to z; only order 0 currently evaluates numerically.

["JacobiTheta", 3, 0, "ImaginaryUnit"]
// ➔ 1.08643481121331  (= π^(1/4)/Γ(3/4))

DedekindEta(τ)

\eta(\tau)
$$$\eta(\tau)$$

The Dedekind eta function, defined on the upper half-plane (\operatorname{Im}(\tau) > 0) by:

\eta(\tau) = e^{i\pi\tau/12} \prod_{k=1}^{\infty} \left(1 - e^{2\pi i k
\tau}\right)

It is a modular form of weight \frac12, central to the theory of modular functions and integer partitions.

["DedekindEta", "ImaginaryUnit"]
// ➔ 0.768225422326057  (= Γ(1/4)/(2π^(3/4)))