Special functions

GSL Shell has moved

The project, its documentation and the release packages are now hosted on Github. Please find here the user manual.

The library includes routines for calculating the values of Airy functions, Bessel functions, Clausen functions, Coulomb wave functions, Coupling coefficients, the Dawson function, Debye functions, Dilogarithms, Elliptic integrals, Jacobi elliptic functions, Error functions, Exponential integrals, Fermi-Dirac functions, Gamma functions, Gegenbauer functions, Hypergeometric functions, Laguerre functions, Legendre functions and Spherical Harmonics, the Psi (Digamma) Function, Synchrotron functions, Transport functions, Trigonometric functions and Zeta functions. Any error returned by the special function (such as invalid input domains, etc...) will be signaled.

Overview

Airy Functions

sf.airyAi(x)
sf.airyBi(x)

The Airy functions Ai(x) and Bi(x) are defined by the integral representations,

Ai(x) = {1 \over \pi} \int_0^\infty \cos({1 \over 3} t^3 + xt) dt

Bi(x) = {1 \over \pi} \int_0^\infty \left( e^{-{1 \over 3} t^3} + \sin(1/3 \, t^3 + xt) \right) dt

For further information see Abramowitz & Stegun, Section 10.4.

sf.airyAi_scaled(x)
sf.airyBi_scaled(x)

These routines compute a scaled version of the Airy function S_A(x) Ai(x). For x>0 the scaling factor S_A(x) is

\exp(+(2/3) x^(3/2))

and is 1 for x<0.

sf.airyAi_deriv(x)
sf.airyAi_deriv_scaled(x)
sf.airyBi_deriv(x)
sf.airyBi_deriv_scaled(x)

These routines compute the Airy function derivatives and their scaled version with the same scaling factor as the airyAi_scaled and airyBi_scaled versions.

sf.airyAi_zero(n)
sf.airyAi_deriv_zero(n)
sf.airyBi_zero(n)
sf.airyBi_deriv_zero(n)

Return the n-th zero of the respective functions or their derivatives.

_images/examples-airy-functions-plot.png

Airy functions Ai and Bi.

Bessel Functions

sf.besselJ(n, x)

These routines compute the regular cylindrical Bessel function of n-th order, Jn(x)

_images/sf-besselJ-functions.png

Bessel functions J0(red), J1(green), J2(blue)

sf.besselJ_zero(n, s)

Return the s-th zero of the Bessel Jn function.

sf.besselY(n, x)

These routines compute the irregular cylindrical Bessel function of n-th order, Yn(x)

_images/sf-bessel-Y.png

Bessel functions Y0(red), Y1(green), Y2(blue)

sf.besselYnu(nu, x)

Compute the irregular cylindrical Bessel function of fractional order \nu, Y_\nu(x).

sf.besselI(n, x)
sf.besselI_scaled(n, x)

Regular modified cylindrical Bessel function of n-th order and their scaled version.

sf.besselInu(n, x)
sf.besselInu_scaled(n, x)

These routines compute the (scaled) regular modified Bessel function of fractional order \nu with I_{\nu}(x) for x>0, nu>0.

sf.besselK(n, x)
sf.besselK_scaled(n, x)

Irregular modified cylindrical Bessel function of order n and their scaled version.

sf.besselKnu(nu, x)
sf.bessellnKnu(nu, x)
sf.besselKnu_scaled(nu, x)

Compute the (scaled or logarithm) irregular modified Bessel function of fractional order \nu, K_\nu(x) for x>0, nu>0.

sf.besselj(l, x)

Compute the regular spherical Bessel function of l-th order.

sf.bessely(l, x)

Compute the irregular spherical Bessel function of l-th order.

sf.besseli_scaled(l, x)

Compute the scaled regular modified spherical Bessel function of l-th order.

sf.besselk(l, x)

Compute the irregular modified spherical Bessel function of l-th order.

Clausen Function

sf.clausen(x)

The Clausen function is defined by the following integral,

Cl_2(x) = - \int_0^x dt \log(2 \sin(t/2))

It is related to the dilogarithm by Cl_2(\theta) = \Im Li_2(\exp(i\theta)).

Coulomb Functions

sf.hydrogenicR_1(Z, r)

These routines compute the lowest-order normalized hydrogenic bound state radial wavefunction

R_1 := 2Z \sqrt{Z} \exp(-Z r)

sf.hydrogenicR(n, l, Z, r)

These routines compute the n-th normalized hydrogenic bound state radial wave function,

R_n := 2 (Z^{3/2}/n^2) \sqrt{(n-l-1)!/(n+l)!} \exp(-Z r/n) (2Zr/n)^lL^{2l+1}_{n-l-1}(2Zr/n)

where L^a_b(x) is the generalized Laguerre polynomial (see Laguerre Functions). The normalization is chosen such that the wave function \psi is given by \psi(n,l,r) = R_n Y_{lm}.

sf.coulomb_wave_FG(eta, x, L_F, k)

This function computes the Coulomb wave functions F_L(\eta,x), G_{L-k}(\eta,x) and their derivatives F'_L(\eta,x), G'_{L-k}(\eta,x) with respect to x. The parameters are restricted to L, L-k > -1/2, x > 0 and integer k.

Note that L itself is not restricted to being an integer. The results are returned as:

F, Fp, G, Gp, exp_F, exp_G = coulomb_wave_FG(eta, x, L_F,k)

with Fp and Gp being the derivates. If an overflow occurs, GSL_EOVRFLW is returned as an error and scaling exponents are stored in the return values exp_F, exp_G.

sf.coulomb_CL(L, eta)

This function computes the Coulomb wave function normalization constant C_L(\eta) for L > -1.

Coupling Coefficients

sf.coupling_3j(two_ja, two_jb, two_jc, two_ma, two_mb, two_mc)

These routines compute the Wigner 3-j coefficient,

\begin{pmatrix}
   j_a & j_b & j_c \\
   m_a & m_b & m_c
\end{pmatrix}

where the arguments are given in half-integer units, ja = two_ja/2, ma = two_ma/2, etc.

sf.coupling_6j(two_ja, two_jb, two_jc, two_ma, two_mb, two_mc)

These routines compute the Wigner 6-j coefficient,

\begin{pmatrix}
   j_a & j_b & j_c \\
   j_d & j_e & j_f
\end{pmatrix}

where the arguments are given in half-integer units, ja = two_ja/2, ma = two_ma/2, etc.

sf.coupling_9j(two_ja, two_jb, two_jc, two_jd, two_je, two_jf, two_jg, two_jh, two_ji)

These routines compute the Wigner 9-j coefficient,

\begin{pmatrix}
   j_a & j_b & j_c \\
   j_d & j_e & j_f \\
   j_g & j_h & j_i
\end{pmatrix}

where the arguments are given in half-integer units, ja = two_ja/2, ma = two_ma/2, etc.

Dawson Function

sf.dawson(x)

The Dawson integral is defined by

\exp(-x^2) \int_0^x dt \exp(t^2)

A table of Dawson’s integral can be found in Abramowitz & Stegun, Table 7.5.

Debye Functions

sf.debye(n, x)

The Debye functions D_n(x) are defined by the following integral,

D_n(x) = {n \over x^n} \int_0^x dt {t^n \over e^t - 1}

For further information see Abramowitz & Stegun, Section 27.1.

Dilogarithms

sf.dilog(x)

These routines compute the dilogarithm for a real argument. In Lewin’s notation this is Li2(x), the real part of the dilogarithm of a real x. It is defined by the integral representation

Li_2(x) = - \Re \int_0^x ds {\log(1-s) \over s}

Note that Im(Li_2(x)) = 0 for x <= 1, and -|pgr| log(x) for x > 1.

Note that Abramowitz & Stegun refer to the Spence integral S(x)= Li_2(1-x) as the dilogarithm rather than Li_2(x).

sf.cdilog(z)

Compute the dilogarithm for a complex argument.

Elliptic Integrals

The Legendre forms of elliptic integrals F(\phi,k), E(\phi,k) and \Pi(\phi,k,n) are defined by,

F(\phi,k) = \int_0^\phi dt 1/\sqrt((1 - k^2 \sin^2(t)))

E(\phi,k) = \int_0^\phi dt   \sqrt((1 - k^2 \sin^2(t)))

Pi(\phi,k,n) = \int_0^\phi dt 1/((1 + n \sin^2(t))\sqrt(1 - k^2 \sin^2(t)))

The complete Legendre forms are denoted by K(k) = F(\pi/2, k) and E(k) = E(\pi/2, k). The notation used here is based on Carlson, Numerische Mathematik 33 (1979) 1 and differs slightly from that used by Abramowitz & Stegun, where the functions are given in terms of the parameter m = k^2 and n is replaced by -n.

The Carlson symmetric forms of elliptical integrals RC(x,y), RD(x,y,z), RF(x,y,z) and RJ(x,y,z,p) are defined by,

RC(x,y) = 1/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1)

RD(x,y,z) = 3/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-3/2)

RF(x,y,z) = 1/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-1/2)

RJ(x,y,z,p) = 3/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-1/2) (t+p)^(-1)

sf.ellint_D(phi, k, n)

These functions compute the incomplete elliptic integral D(\phi,k) which is defined through the Carlson form RD(x,y,z) by the following relation,

D(\phi,k,n) = (1/3)(\sin(\phi))^3 RD (1-\sin^2(\phi), 1-k^2 \sin^2(\phi), 1).

The argument n is not used and will be removed in a future release.

sf.ellint_E(phi, lk)

These routines compute the incomplete elliptic integral E(\phi,k) to the accuracy specified by the mode variable mode. Note that Abramowitz & Stegun define this function in terms of the parameter m = k^2.

sf.ellint_F(phi, k)

These routines compute the incomplete elliptic integral F(\phi,k) to the accuracy specified by the mode variable mode. Note that Abramowitz & Stegun define this function in terms of the parameter m = k^2.

sf.ellint_P(phi, k, n)

These routines compute the incomplete elliptic integral \Pi(\phi,k,n) to the accuracy specified by the mode variable mode. Note that Abramowitz & Stegun define this function in terms of the parameters m = k^2 and \sin^2(\alpha) = k^2, with the change of sign n \to -n.

sf.ellint_RC(x, y)

These routines compute the incomplete elliptic integral RC(x,y) to the accuracy specified by the mode variable mode.

sf.ellint_RD(x, y, z)

These routines compute the incomplete elliptic integral RD(x,y,z) to the accuracy specified by the mode variable mode.

sf.ellint_RF(x, y, z)

These routines compute the incomplete elliptic integral RF(x,y,z) to the accuracy specified by the mode variable mode.

sf.ellint_RJ(x, y, z, p)

These routines compute the incomplete elliptic integral RJ(x,y,z,p) to the accuracy specified by the mode variable mode.

Elliptic Function

sf.elljac(u, m)

This function computes the Jacobian elliptic functions sn(u|m), cn(u|m), dn(u|m) by descending Landen transformations. It returns sn, cn and dn as:

sn, cn, dn = elljac(u, m)

Error Functions

sf.erf(x)

The error function erf(x)

{2 \over \sqrt{\pi}} \int_0^x dt e^{-t^2}

sf.erfc(x)
sf.log_erfc(x)

These routines compute the (logarithmic) complementary error function erfc(x) = 1 - erf(x) = (2/\sqrt(\pi)) \int_x^\infty \exp(-t^2).

sf.erf_Q(1)

These routines compute the upper tail of the Gaussian probability function Q(x) = (1/\sqrt{2\pi}) \int_x^\infty dt \exp(-t^2/2

sf.erf_Z(1)

These routines compute the Gaussian probability density function Z(x) = (1/\sqrt{2\pi}) \exp(-x^2/2).

sf.hazard(x)

These routines compute the hazard function for the normal distribution.

Exponential Functions

sf.exp(x)

These routines provide an exponential function \exp(x).

sf.exp_err(x, dx)

This function exponentiates x with an associated absolute error dx.

sf.exp_mult(x, y)

These routines exponentiate x and multiply by the factor y to return the product y \exp(x).

sf.exp_mult_err(x, dx, y, dy)

This routine computes the product y \exp(x) for the quantities x, y with associated absolute errors dx, dy.

sf.expm1(x)

These routines compute the quantity \exp(x)-1 using an algorithm that is accurate for small x.

sf.exprel(x)

These routines compute the quantity (\exp(x)-1)/x using an algorithm that is accurate for small x. For small x the algorithm is based on the expansion (\exp(x)-1)/x = 1 + x/2 + x^2/(2*3) + x^3/(2*3*4) + \dots.

sf.exprel_2(x)

These routines compute the quantity 2(\exp(x)-1-x)/x^2 using an algorithm that is accurate for small x. For small x the algorithm is based on the expansion 2(\exp(x)-1-x)/x^2 = 1 + x/3 + x^2/(3*4) + x^3/(3*4*5) + \dots.

sf.exprel_n(n, 1)

These routines compute the N-relative exponential given by,

exprel_N(x) = N!/x^N (\exp(x) - \sum_{k=0}^{N-1} x^k/k!)
            = 1 + x/(N+1) + x^2/((N+1)(N+2)) + ...
            = 1F1 (1,1+N,x)

Exponential Integrals

sf.expint_E(n, x)

These routines compute the exponential integral E_n(x) of order n,

E_n(x) := \Re \int_1^\infty dt \exp(-xt)/t^n.

sf.expint_Ei(x)

These routines compute the exponential integral Ei(x),

Ei(x) := - PV(\int_{-x}^\infty dt \exp(-t)/t)

where PV denotes the principal value of the integral.

sf.Shi(x)

These routines compute the integral Shi(x) = \int_0^x dt \sinh(t)/t.

sf.Chi(x)

These routines compute the integral Chi(x) := \Re[ \gamma_E + \log(x) + \int_0^x dt (\cosh(t)-1)/t] , where \gamma_E is the Euler constant (available as the macro M_EULER).

sf.expint_3(x)

These routines compute the third-order exponential integral Ei_3(x) = \int_0^xdt \exp(-t^3) for x >= 0.

sf.Si(x)

These routines compute the Sine integral Si(x) = \int_0^x dt \sin(t)/t.

sf.Ci(x)

These routines compute the Cosine integral Ci(x) = -\int_x^\infty dt \cos(t)/t for x > 0.

sf.atanint(x)

These routines compute the Arctangent integral, which is defined as AtanInt(x) = \int_0^x dt \arctan(t)/t.

Fermi Dirac Function

sf.fermi_dirac(n, x)

The complete Fermi-Dirac integral Fn(x) is given by,

F_n(x) = {1 \over \Gamma(n+1)} \int_0^\infty dt {t^n \over e^{t-x} + 1}

Note that the Fermi-Dirac integral is sometimes defined without the normalization factor in other texts.

sf.fermi_dirac_inc(x, b)

These routines compute the incomplete Fermi-Dirac integral with an index of zero, F_0(x,b) = \ln(1 + e^{b-x}) - (b-x).

The incomplete Fermi-Dirac integral F_j(x,b) is given by,

F_j(x,b)   := (1/\Gamma(j+1)) \int_b^\infty dt (t^j / (\exp(t-x) + 1))

Gamma and Beta Functions

sf.Shi(x)

Compute the integral \textrm{Shi}(x) = \int_0^x dt \sinh(t)/t.

sf.Chi(x)

Compute the integral

\textrm{Chi}(x) = \Re \left( \gamma_E +
\log(x) + \int_0^x dt \dfrac{\cosh(t)-1}{t} \right)

where γE is the Euler constant.

sf.Si(x)

Compute the Sine integral \textrm{Si}(x) = \int_0^x dt \dfrac{\sin(t)}{t}.

sf.Ci(x)

Compute the Cosine integral \textrm{Ci}(x) = -\int_x^\infty dt \dfrac{\cos(t)}{t} for x > 0.

sf.atanint(x)

compute the Arctangent integral, which is defined as \textrm{AtanInt}(x) = \int_0^x dt \dfrac{\arctan(t)}{t}.

sf.fact(n)

Compute the factorial n!. The factorial is related to the Gamma function by n! = |Ggr| (n+1).

sf.doublefact(n)

Compute the double factorial n!! = n(n-2)(n-4) \dots.

sf.lnfact(n)

These routines compute the logarithm of the factorial of N, log(n!). The algorithm is faster than computing ln(Γ (n+1)).

sf.lndoublefact(n)

Compute the logarithm of the double factorial of N, log(n!!).

sf.choose(n, k)

Compute the combinatorial factor

\binom{n}{k} = \dfrac{n!}{k! (n-k)!}

sf.lnchoose(n, k)

Compute the logarithm of “n choose m”. This is equivalent to the sum \log(n!) - \log(m!) - \log((n-m)!).

sf.gamma(x)

Compute the Gamma function Γ (x), subject to x not being a negative integer or zero. The function is computed using the real Lanczos method.

sf.lngamma(x)

Compute the logarithm of the Gamma function, \log(\Gamma(x)), subject to x not being a negative integer or zero. For x<0 the real part of \log(\Gamma(x)) is returned, which is equivalent to \log(|\Gamma(x)|). The function is computed using the real Lanczos method.

sf.gammastar(x)

These routines compute the regulated Gamma Function \Gamma^*(x) for x > 0. The regulated gamma function is given by,

\Gamma^*(x) = \Gamma(x)/(\sqrt{2\pi} x^{(x-1/2)} \exp(-x))

            = (1 + (1/12x) + ...)  for x \to \infty

and is a useful suggestion of Temme.

sf.gammainv(x)

These routines compute the reciprocal of the gamma function, 1/\Gamma(x) using the real Lanczos method.

sf.lngammac(z)

This routine computes \log(\Gamma(z)) for complex z and z not a negative integer or zero, using the complex Lanczos method. The returned parameters are lnr = \log|\Gamma(z)| and arg = \arg(\Gamma(z)) in (-\pi,\pi]. Note that the phase part (arg) is not well-determined when |z| is very large, due to inevitable roundoff in restricting to (-pi,pi]. This will result in a GSL_ELOSS error when it occurs. The absolute value part (lnr), however, never suffers from loss of precision.

The functions returns:

lnr, arg = lngammac(z)
sf.beta(a, b)

These routines compute the Beta Function, B(a,b) = \Gamma(a)\Gamma(b)/\Gamma(a+b) subject to a and b not being negative integers.

sf.lnbeta(a, b)

These routines compute the logarithm of the Beta Function, \log(B(a,b)) subject to a and b not being negative integers.

sf.beta_inc(a, b, x)

These routines compute the normalized incomplete Beta function I_x(a,b)=B_x(a,b)/B(a,b) where B_x(a,b) = \int_0^x t^{a-1} (1-t)^{b-1} dt for 0 <= x <= 1. For a > 0, b > 0 the value is computed using a continued fraction expansion. For all other values it is computed using the relation I_x(a,b,x) = (1/a) x^a 2F1(a,1-b,a+1,x)/B(a,b).

sf.taylorcoeff(n, x)

These routines compute the Taylor coefficient x^n / n! for x >= 0, n >= 0.

sf.poch(a, x)

These routines compute the Pochhammer symbol (a)_x = \Gamma(a + x)/\Gamma(a). The Pochhammer symbol is also known as the Apell symbol and sometimes written as (a,x). When a and a+x are negative integers or zero, the limiting value of the ratio is returned.

sf.lnpoch(a, x)

These routines compute the logarithm of the Pochhammer symbol, \log((a)_x) = \log(\Gamma(a + x)/\Gamma(a)).

sf.pochrel(a, x)

These routines compute the relative Pochhammer symbol ((a)_x - 1)/x where (a)_x = \Gamma(a + x)/\Gamma(a).

sf.gamma_inc(a, x)

These functions compute the unnormalized incomplete Gamma Function \Gamma(a,x) = \int_x^\infty dt t^{a-1} \exp(-t) for a real and x >= 0.

sf.gamma_inc_Q(a, x)

These routines compute the normalized incomplete Gamma Function Q(a,x) = 1/\Gamma(a) \int_x^\infty dt t^{a-1} \exp(-t) for a > 0, x >= 0.

sf.gamma_inc_P(a, x)

These routines compute the complementary normalized incomplete Gamma Function P(a,x) = 1 - Q(a,x) = 1/\Gamma(a) \int_0^x dt t^{a-1} \exp(-t) for a > 0, x >= 0. Note that Abramowitz & Stegun call P(a,x) the incomplete gamma function (section 6.5).

Gegenbauer Functions

sf.gegenpoly(n, lambda, x)

These functions evaluate the Gegenbauer polynomial C^{(\lambda)}_n(x) for a specific value of n, lambda, x subject to \lambda > -1/2, n >= 0.

Hypergeometric functions

sf.hyperg0F1(a, b)

These routines compute the hypergeometric function 0F1(c,x).

sf.hyperg1F1(m, n, x)

These routines compute the confluent hypergeometric function 1F1(m,n,x) = M(m,n,x). The parameters m and n can be integer or real numbers.

sf.hypergU(m, n, x)

These routines compute the confluent hypergeometric function U(m,n,x). The parameters m and n can be integer or real numbers.

sf.hyperg2F1(a, b, c, x)

These routines compute the Gauss hypergeometric function 2F1(a,b,c,x) = F(a,b,c,x) for |x| < 1.

If the arguments (a,b,c,x) are too close to a singularity then the function can return the error code GSL_EMAXITER when the series approximation converges too slowly. This occurs in the region of x=1, c - a - b = m for integer m.

sf.hyperg2F1conj(a, c, x)

These routines compute the Gauss hypergeometric function 2F1(a, a^*, c, x) where a is complex parameter, c and x are real parameters with |x| < 1.

Laguerre Functions

sf.laguerre(n, a, x)

The generalized Laguerre polynomials are defined in terms of confluent hypergeometric functions as

L^a_n(x) = {(a+1)_n \over n!} {}_1 F_1(-n,a+1,x)

and are sometimes referred to as the associated Laguerre polynomials. They are related to the plain Laguerre polynomials L_n(x) by

L^0_n(x) = L_n(x)

and

L^k_n(x) = (-1)^k {d^k \over dx^k} L_{n+k}(x)

For more information see Abramowitz & Stegun, Chapter 22.

Lambert W Functions

Lambert’s W functions, W(x), are defined to be solutions of the equation W(x) \exp(W(x)) = x. This function has multiple branches for x < 0; however, it has only two real-valued branches. We define W_0(x) to be the principal branch, where W > -1 for x < 0, and W_{-1}(x) to be the other real branch, where W < -1 for x < 0.

sf.lambertW0(x)

These compute the principal branch of the Lambert W function, W_0(x).

sf.lambertWm1(x)

These compute the secondary real-valued branch of the Lambert W function, W_{-1}(x).

Legendre Functions and Spherical Harmonics

sf.legendreP(n, x)
_images/examples-legendre-polynomials.png

Legendre polynomials

These functions evaluate the Legendre polynomial P_l(x) for a specific value of l, x subject to l >= 0, |x| <= 1

sf.legendreQ(n, x)

These routines compute the Legendre function Q_l(x) for x > -1, x != 1 and l >= 0.

sf.legendrePlm(l, m, x)

These routines compute the associated Legendre polynomial P_l^m(x) for m >= 0, l >= m, |x| <= 1.

sf.legendresphPlm(l, m, x)

These routines compute the normalized associated Legendre polynomial \sqrt{(2l+1)/(4\pi)} \sqrt{(l-m)!/(l+m)!} P_l^m(x) suitable for use in spherical harmonics. The parameters must satisfy m >= 0, l >= m, |x| <= 1. Theses routines avoid the overflows that occur for the standard normalization of P_l^m(x).

sf.conicalP(n, lambda, x)

These routines compute the regular Spherical Conical Function P^{n}_{-1/2 + i \lambda}(x) for x > -1 where n in {-12, 0, 12, 1}

sf.conicalPsphreg(l, lambda, x)

These routines compute the Regular Spherical Conical Function P^{-1/2-l}_{-1/2 + i \lambda}(x) for x > -1, l >= -1.

sf.conicalPcylreg(m, lambda, x)

These routines compute the Regular Cylindrical Conical Function P^{-m}_{-1/2 + i \lambda}(x) for x > -1, m >= -1.

sf.legendre_H3d(l, lambda, eta)

These routines compute the l-th radial eigenfunction of the Laplacian on the 3-dimensional hyperbolic space \eta >= 0, l >= 0. In the flat limit this takes the form L^{H3d}_l(\lambda,\eta) = j_l(\lambda\eta).

Psi (Digamma) Functions

The polygamma functions of order n are defined by

\psi^{(n)}(x) = (d/dx)^n \psi(x) = (d/dx)^{n+1} \log(\Gamma(x))

where \psi(x) = \Gamma'(x)/\Gamma(x) is known as the digamma function.

sf.psi(x)

These routines compute the digamma function \psi(x) for general x, x ne 0.

sf.psi_1(x)

These routines compute the Trigamma function \psi'(x) for general x.

sf.1piy(y)

These routines compute the real part of the digamma function on the line 1+i y, \Re[\psi(1 + i y)].

sf.psi_n(n, x)

These routines compute the polygamma function \psi^{(n)}(x) for n >= 0, x > 0.

Synchrotron Functions

sf.synchrotron1(n, x)

These routines compute the first synchrotron function x \int_x^\infty dt K_{5/3}(t) for x >= 0.

sf.synchrotron2(n, x)

These routines compute the second synchrotron function x K_{2/3}(x) for x >= 0.

Transport Functions

The transport functions J(n,x) are defined by the integral representations J(n,x) := \int_0^x dt t^n e^t /(e^t - 1)^2.

sf.transport(n, x)

These routines compute the transport function J(n,x) with n \in {2,3,4,5}

Zeta Functions

sf.zeta(s)

The Riemann zeta function is defined by the infinite sum \zeta(s) = \sum_{k=1}^\infty k^{-s}. These routines compute the Riemann zeta function \zeta(s) for arbitrary s, s \ne 1.

sf.zetam1(s)

For large positive argument, the Riemann zeta function approaches one. In this region the fractional part is interesting, and therefore we need a function to evaluate it explicitly. These routines compute \zeta(s) - 1 for arbitrary s, s \ne 1.

sf.eta(s)

The eta function is defined by \eta(s) = (1-2^{1-s}) \zeta(s). These routines compute the eta function eta(s) for arbitrary s.

sf.hzeta(s, q)

These routines compute the Hurwitz zeta function \zeta(s,q) for s > 1, q > 0.