This chapter describes the GSL special function library. 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. Each routine also computes an estimate of the numerical
error in the calculated value of the function.
The functions are declared in individual header files, such as
`gsl_sf_airy.h', `gsl_sf_bessel.h', etc. The complete set of
header files can be included using the file `gsl_sf.h'.
The special functions are available in two calling conventions, a
natural form which returns the numerical value of the function and
an error-handling form which returns an error code. The two types
of function provide alternative ways of accessing the same underlying
code.
The natural form returns only the value of the function and can be
used directly in mathematical expressions.. For example, the following
function call will compute the value of the Bessel function
J_0(x),
double y = gsl_sf_bessel_J0 (x);
There is no way to access an error code or to estimate the error using
this method. To allow access to this information the alternative
error-handling form stores the value and error in a modifiable argument,
gsl_sf_result result;
int status = gsl_sf_bessel_J0_e (x, &result);
The error-handling functions have the suffix _e. The returned
status value indicates error conditions such as overflow, underflow or
loss of precision. If there are no errors the error-handling functions
return GSL_SUCCESS.
The error handling form of the special functions always calculate an
error estimate along with the value of the result. Therefore,
structures are provided for amalgamating a value and error estimate.
These structures are declared in the header file `gsl_sf_result.h'.
The gsl_sf_result struct contains value and error fields.
The field val contains the value and the field err contains
an estimate of the absolute error in the value.
In some cases, an overflow or underflow can be detected and handled by a
function. In this case, it may be possible to return a scaling exponent
as well as an error/value pair in order to save the result from
exceeding the dynamic range of the built-in types. The
gsl_sf_result_e10 struct contains value and error fields as well
as an exponent field such that the actual result is obtained as
result * 10^(e10).
The goal of the library is to achieve double precision accuracy wherever
possible. However the cost of evaluating some special functions to
double precision can be significant, particularly where very high order
terms are required. In these cases a mode argument allows the
accuracy of the function to be reduced in order to improve performance.
The following precision levels are available for the mode argument,
GSL_PREC_DOUBLE
Double-precision, a relative accuracy of approximately
2 * 10^-16.
GSL_PREC_SINGLE
Single-precision, a relative accuracy of approximately
10^-7.
GSL_PREC_APPROX
Approximate values, a relative accuracy of approximately
5 * 10^-4.
The approximate mode provides the fastest evaluation at the lowest
accuracy.
Function: int gsl_sf_airy_Ai_scaled_e(double x, gsl_mode_t mode, gsl_sf_result * result)
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.
The routines described in this section compute the Cylindrical Bessel
functions J_n(x), Y_n(x), Modified cylindrical Bessel
functions I_n(x), K_n(x), Spherical Bessel functions
j_l(x), y_l(x), and Modified Spherical Bessel functions
i_l(x), k_l(x). For more information see Abramowitz & Stegun,
Chapters 9 and 10. The Bessel functions are defined in the header file
`gsl_sf_bessel.h'.
Function: int gsl_sf_bessel_J0_e(double x, gsl_sf_result * result)
These routines compute the regular cylindrical Bessel function of zeroth
order, J_0(x).
Function: double gsl_sf_bessel_J1(double x)
Function: int gsl_sf_bessel_J1_e(double x, gsl_sf_result * result)
These routines compute the regular cylindrical Bessel function of first
order, J_1(x).
Function: double gsl_sf_bessel_Jn(int n, double x)
Function: int gsl_sf_bessel_Jn_e(int n, double x, gsl_sf_result * result)
These routines compute the regular cylindrical Bessel function of
order n, J_n(x).
Function: int gsl_sf_bessel_Jn_array(int nmin, int nmax, double x, double result_array[])
This routine computes the values of the regular cylindrical Bessel
functions J_n(x) for n from nmin to nmax
inclusive, storing the results in the array result_array. The
values are computed using recurrence relations, for efficiency, and
therefore may differ slightly from the exact values.
Function: int gsl_sf_bessel_Y0_e(double x, gsl_sf_result * result)
These routines compute the irregular cylindrical Bessel function of zeroth
order, Y_0(x), for x>0.
Function: double gsl_sf_bessel_Y1(double x)
Function: int gsl_sf_bessel_Y1_e(double x, gsl_sf_result * result)
These routines compute the irregular cylindrical Bessel function of first
order, Y_1(x), for x>0.
Function: double gsl_sf_bessel_Yn(int n,double x)
Function: int gsl_sf_bessel_Yn_e(int n,double x, gsl_sf_result * result)
These routines compute the irregular cylindrical Bessel function of
order n, Y_n(x), for x>0.
Function: int gsl_sf_bessel_Yn_array(int nmin, int nmax, double x, double result_array[])
This routine computes the values of the irregular cylindrical Bessel
functions Y_n(x) for n from nmin to nmax
inclusive, storing the results in the array result_array. The
domain of the function is x>0. The values are computed using
recurrence relations, for efficiency, and therefore may differ slightly
from the exact values.
Function: int gsl_sf_bessel_I0_e(double x, gsl_sf_result * result)
These routines compute the regular modified cylindrical Bessel function
of zeroth order, I_0(x).
Function: double gsl_sf_bessel_I1(double x)
Function: int gsl_sf_bessel_I1_e(double x, gsl_sf_result * result)
These routines compute the regular modified cylindrical Bessel function
of first order, I_1(x).
Function: double gsl_sf_bessel_In(int n, double x)
Function: int gsl_sf_bessel_In_e(int n, double x, gsl_sf_result * result)
These routines compute the regular modified cylindrical Bessel function
of order n, I_n(x).
Function: int gsl_sf_bessel_In_array(int nmin, int nmax, double x, double result_array[])
This routine computes the values of the regular modified cylindrical
Bessel functions I_n(x) for n from nmin to
nmax inclusive, storing the results in the array
result_array. The start of the range nmin must be positive
or zero. The values are computed using recurrence relations, for
efficiency, and therefore may differ slightly from the exact values.
Function: int gsl_sf_bessel_I1_scaled_e(double x, gsl_sf_result * result)
These routines compute the scaled regular modified cylindrical Bessel
function of first order \exp(-|x|) I_1(x).
Function: double gsl_sf_bessel_In_scaled(int n, double x)
Function: int gsl_sf_bessel_In_scaled_e(int n, double x, gsl_sf_result * result)
These routines compute the scaled regular modified cylindrical Bessel
function of order n, \exp(-|x|) I_n(x)
Function: int gsl_sf_bessel_In_scaled_array(int nmin, int nmax, double x, double result_array[])
This routine computes the values of the scaled regular cylindrical
Bessel functions \exp(-|x|) I_n(x) for n from
nmin to nmax inclusive, storing the results in the array
result_array. The start of the range nmin must be positive
or zero. The values are computed using recurrence relations, for
efficiency, and therefore may differ slightly from the exact values.
Function: int gsl_sf_bessel_K0_e(double x, gsl_sf_result * result)
These routines compute the irregular modified cylindrical Bessel
function of zeroth order, K_0(x), for x > 0.
Function: double gsl_sf_bessel_K1(double x)
Function: int gsl_sf_bessel_K1_e(double x, gsl_sf_result * result)
These routines compute the irregular modified cylindrical Bessel
function of first order, K_1(x), for x > 0.
Function: double gsl_sf_bessel_Kn(int n, double x)
Function: int gsl_sf_bessel_Kn_e(int n, double x, gsl_sf_result * result)
These routines compute the irregular modified cylindrical Bessel
function of order n, K_n(x), for x > 0.
Function: int gsl_sf_bessel_Kn_array(int nmin, int nmax, double x, double result_array[])
This routine computes the values of the irregular modified cylindrical
Bessel functions K_n(x) for n from nmin to
nmax inclusive, storing the results in the array
result_array. The start of the range nmin must be positive
or zero. The domain of the function is x>0. The values are
computed using recurrence relations, for efficiency, and therefore
may differ slightly from the exact values.
Function: int gsl_sf_bessel_K1_scaled_e(double x, gsl_sf_result * result)
These routines compute the scaled irregular modified cylindrical Bessel
function of first order \exp(x) K_1(x) for x>0.
Function: double gsl_sf_bessel_Kn_scaled(int n, double x)
Function: int gsl_sf_bessel_Kn_scaled_e(int n, double x, gsl_sf_result * result)
These routines compute the scaled irregular modified cylindrical Bessel
function of order n, \exp(x) K_n(x), for x>0.
Function: int gsl_sf_bessel_Kn_scaled_array(int nmin, int nmax, double x, double result_array[])
This routine computes the values of the scaled irregular cylindrical
Bessel functions \exp(x) K_n(x) for n from nmin to
nmax inclusive, storing the results in the array
result_array. The start of the range nmin must be positive
or zero. The domain of the function is x>0. The values are
computed using recurrence relations, for efficiency, and therefore
may differ slightly from the exact values.
Function: int gsl_sf_bessel_jl_e(int l, double x, gsl_sf_result * result)
These routines compute the regular spherical Bessel function of
order l, j_l(x), for
l >= 0 and
x >= 0.
Function: int gsl_sf_bessel_jl_array(int lmax, double x, double result_array[])
This routine computes the values of the regular spherical Bessel
functions j_l(x) for l from 0 to lmax
inclusive for
lmax >= 0 and
x >= 0, storing the results in the array result_array.
The values are computed using recurrence relations, for
efficiency, and therefore may differ slightly from the exact values.
Function: int gsl_sf_bessel_jl_steed_array(int lmax, double x, double * jl_x_array)
This routine uses Steed's method to compute the values of the regular
spherical Bessel functions j_l(x) for l from 0 to
lmax inclusive for
lmax >= 0 and
x >= 0, storing the results in the array
result_array.
The Steed/Barnett algorithm is described in Comp. Phys. Comm. 21,
297 (1981). Steed's method is more stable than the
recurrence used in the other functions but is also slower.
Function: int gsl_sf_bessel_yl_e(int l, double x, gsl_sf_result * result)
These routines compute the irregular spherical Bessel function of
order l, y_l(x), for
l >= 0.
Function: int gsl_sf_bessel_yl_array(int lmax, double x, double result_array[])
This routine computes the values of the irregular spherical Bessel
functions y_l(x) for l from 0 to lmax
inclusive for
lmax >= 0, storing the results in the array result_array.
The values are computed using recurrence relations, for
efficiency, and therefore may differ slightly from the exact values.
The regular modified spherical Bessel functions i_l(x)
are related to the modified Bessel functions of fractional order,
i_l(x) = \sqrt{\pi/(2x)} I_{l+1/2}(x)
Function: int gsl_sf_bessel_il_scaled_e(int l, double x, gsl_sf_result * result)
These routines compute the scaled regular modified spherical Bessel
function of order l, \exp(-|x|) i_l(x)
Function: int gsl_sf_bessel_il_scaled_array(int lmax, double x, double result_array[])
This routine computes the values of the scaled regular modified
cylindrical Bessel functions \exp(-|x|) i_l(x) for l from
0 to lmax inclusive for
lmax >= 0, storing the results in
the array result_array.
The values are computed using recurrence relations, for
efficiency, and therefore may differ slightly from the exact values.
The irregular modified spherical Bessel functions k_l(x)
are related to the irregular modified Bessel functions of fractional order,
k_l(x) = \sqrt{\pi/(2x)} K_{l+1/2}(x).
Function: int gsl_sf_bessel_kl_scaled_e(int l, double x, gsl_sf_result * result)
These routines compute the scaled irregular modified spherical Bessel
function of order l, \exp(x) k_l(x), for x>0.
Function: int gsl_sf_bessel_kl_scaled_array(int lmax, double x, double result_array[])
This routine computes the values of the scaled irregular modified
spherical Bessel functions \exp(x) k_l(x) for l from
0 to lmax inclusive for
lmax >= 0 and x>0, storing the results in
the array result_array.
The values are computed using recurrence relations, for
efficiency, and therefore may differ slightly from the exact values.
Function: int gsl_sf_bessel_Jnu_e(double nu, double x, gsl_sf_result * result)
These routines compute the regular cylindrical Bessel function of
fractional order nu, J_\nu(x).
Function: int gsl_sf_bessel_sequence_Jnu_e(double nu, gsl_mode_t mode, size_t size, double v[])
This function computes the regular cylindrical Bessel function of
fractional order \nu, J_\nu(x), evaluated at a series of
x values. The array v of length size contains the
x values. They are assumed to be strictly ordered and positive.
The array is over-written with the values of J_\nu(x_i).
The Coulomb wave functions F_L(\eta,x), G_L(\eta,x) are
described in Abramowitz & Stegun, Chapter 14. Because there can be a
large dynamic range of values for these functions, overflows are handled
gracefully. If an overflow occurs, GSL_EOVRFLW is signalled and
exponent(s) are returned through the modifiable parameters exp_F,
exp_G. The full solution can be reconstructed from the following
relations,
Function: int gsl_sf_coulomb_wave_FG_e(double eta, double x, double L_F, int k, gsl_sf_result * F, gsl_sf_result * Fp, gsl_sf_result * G, gsl_sf_result * Gp, double * exp_F, double * exp_G)
This function computes the coulomb wave functions F_L(\eta,x),
G_{L-k}(\eta,x) and their derivatives with respect to x,
F'_L(\eta,x)
G'_{L-k}(\eta,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 stored in
the parameters F, G for the function values and Fp,
Gp for the derivative values. If an overflow occurs,
GSL_EOVRFLW is returned and scaling exponents are stored in
the modifiable parameters exp_F, exp_G.
Function: int gsl_sf_coulomb_wave_F_array(double L_min, int kmax, double eta, double x, double fc_array[], double * F_exponent)
This function computes the function F_L(eta,x) for L = Lmin
\dots Lmin + kmax storing the results in fc_array. In the case of
overflow the exponent is stored in F_exponent.
Function: int gsl_sf_coulomb_wave_FG_array(double L_min, int kmax, double eta, double x, double fc_array[], double gc_array[], double * F_exponent, double * G_exponent)
This function computes the functions F_L(\eta,x),
G_L(\eta,x) for L = Lmin \dots Lmin + kmax storing the
results in fc_array and gc_array. In the case of overflow the
exponents are stored in F_exponent and G_exponent.
Function: int gsl_sf_coulomb_wave_FGp_array(double L_min, int kmax, double eta, double x, double fc_array[], double fcp_array[], double gc_array[], double gcp_array[], double * F_exponent, double * G_exponent)
This function computes the functions F_L(\eta,x),
G_L(\eta,x) and their derivatives F'_L(\eta,x),
G'_L(\eta,x) for L = Lmin \dots Lmin + kmax storing the
results in fc_array, gc_array, fcp_array and gcp_array.
In the case of overflow the exponents are stored in F_exponent
and G_exponent.
Function: int gsl_sf_coulomb_wave_sphF_array(double L_min, int kmax, double eta, double x, double fc_array[], double F_exponent[])
This function computes the Coulomb wave function divided by the argument
F_L(\eta, x)/x for L = Lmin \dots Lmin + kmax, storing the
results in fc_array. In the case of overflow the exponent is
stored in F_exponent. This function reduces to spherical Bessel
functions in the limit \eta \to 0.
The Wigner 3-j, 6-j and 9-j symbols give the coupling coefficients for
combined angular momentum vectors. Since the arguments of the standard
coupling coefficient functions are integer or half-integer, the
arguments of the following functions are, by convention, integers equal
to twice the actual spin value. For information on the 3-j coefficients
see Abramowitz & Stegun, Section 27.9. The functions described in this
section are declared in the header file `gsl_sf_coupling.h'.
Function: double gsl_sf_coupling_9j(int two_ja, int two_jb, int two_jc, int two_jd, int two_je, int two_jf, int two_jg, int two_jh, int two_ji)
Function: int gsl_sf_coupling_9j_e(int two_ja, int two_jb, int two_jc, int two_jd, int two_je, int two_jf, int two_jg, int two_jh, int two_ji, gsl_sf_result * result)
These routines compute the Wigner 9-j coefficient,
(ja jb jc
jd je jf
jg jh ji)
where the arguments are given in half-integer units, ja =
two_ja/2, ma = two_ma/2, etc.
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. The Dawson functions are declared in the header file
`gsl_sf_dawson.h'.
Function: double gsl_sf_dawson(double x)
Function: int gsl_sf_dawson_e(double x, gsl_sf_result * result)
These routines compute the value of Dawson's integral for x.
The Debye functions are defined by the integral D_n(x) = n/x^n
\int_0^x dt (t^n/(e^t - 1)). For further information see Abramowitz &
Stegun, Section 27.1. The Debye functions are declared in the header
file `gsl_sf_debye.h'.
Function: double gsl_sf_debye_1(double x)
Function: int gsl_sf_debye_1_e(double x, gsl_sf_result * result)
These routines compute the first-order Debye function
D_1(x) = (1/x) \int_0^x dt (t/(e^t - 1)).
Function: double gsl_sf_debye_2(double x)
Function: int gsl_sf_debye_2_e(double x, gsl_sf_result * result)
These routines compute the second-order Debye function
D_2(x) = (2/x^2) \int_0^x dt (t^2/(e^t - 1)).
Function: double gsl_sf_debye_3(double x)
Function: int gsl_sf_debye_3_e(double x, gsl_sf_result * result)
These routines compute the third-order Debye function
D_3(x) = (3/x^3) \int_0^x dt (t^3/(e^t - 1)).
Function: double gsl_sf_debye_4(double x)
Function: int gsl_sf_debye_4_e(double x, gsl_sf_result * result)
These routines compute the fourth-order Debye function
D_4(x) = (4/x^4) \int_0^x dt (t^4/(e^t - 1)).
Function: int gsl_sf_dilog_e(double x, gsl_sf_result * result)
These routines compute the dilogarithm for a real argument. In Lewin's
notation this is Li_2(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) / s.
Note that \Im(Li_2(x)) = 0 for
x <= 1, and -\pi\log(x) for x > 1.
This function computes the full complex-valued dilogarithm for the
complex argument z = r \exp(i \theta). The real and imaginary
parts of the result are returned in result_re, result_im.
The following functions allow for the propagation of errors when
combining quantities by multiplication. The functions are declared in
the header file `gsl_sf_elementary.h'.
Function: int gsl_sf_multiply_e(double x, double y, gsl_sf_result * result)
This function multiplies x and y storing the product and its
associated error in result.
Function: int gsl_sf_multiply_err_e(double x, double dx, double y, double dy, gsl_sf_result * result)
This function multiplies x and y with associated absolute
errors dx and dy. The product
xy +/- xy \sqrt((dx/x)^2 +(dy/y)^2)
is stored in result.
The complete Legendre forms are denoted by K(k) = F(\pi/2, k) and
E(k) = E(\pi/2, k). Further information on the Legendre forms of
elliptic integrals can be found in Abramowitz & Stegun, Chapter 17. The
notation used here is based on Carlson, Numerische Mathematik 33
(1979) 1 and differs slightly from that used by Abramowitz & Stegun.
Function: int gsl_sf_exp_e(double x, gsl_sf_result * result)
These routines provide an exponential function \exp(x) using GSL
semantics and error checking.
Function: int gsl_sf_exp_e10_e(double x, gsl_sf_result_e10 * result)
This function computes the exponential \exp(x) using the
gsl_sf_result_e10 type to return a result with extended range.
This function may be useful if the value of \exp(x) would
overflow the numeric range of double.
Function: int gsl_sf_expm1_e(double x, gsl_sf_result * result)
These routines compute the quantity \exp(x)-1 using an algorithm
that is accurate for small x.
Function: double gsl_sf_exprel(double x)
Function: int gsl_sf_exprel_e(double x, gsl_sf_result * result)
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.
Function: double gsl_sf_exprel_2(double x)
Function: int gsl_sf_exprel_2_e(double x, gsl_sf_result * result)
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.
Function: double gsl_sf_exprel_n(int n, double x)
Function: int gsl_sf_exprel_n_e(int n, double x, gsl_sf_result * result)
These routines compute the N-relative exponential, which is the
n-th generalization of the functions gsl_sf_exprel and
gsl_sf_exprel2. The N-relative exponential is given by,
Function: int gsl_sf_exp_err_e(double x, double dx, gsl_sf_result * result)
This function exponentiates x with an associated absolute error
dx.
Function: int gsl_sf_exp_err_e10_e(double x, double dx, gsl_sf_result_e10 * result)
This functions exponentiate a quantity x with an associated absolute
error dx using the gsl_sf_result_e10 type to return a result with
extended range.
Function: int gsl_sf_exp_mult_err_e(double x, double dx, double y, double dy, gsl_sf_result * result)
This routine computes the product y \exp(x) for the quantities
x, y with associated absolute errors dx, dy.
Function: int gsl_sf_exp_mult_err_e10_e(double x, double dx, double y, double dy, gsl_sf_result_e10 * result)
This routine computes the product y \exp(x) for the quantities
x, y with associated absolute errors dx, dy using the
gsl_sf_result_e10 type to return a result with extended range.
Information on the exponential integrals can be found in Abramowitz &
Stegun, Chapter 5. These functions are declared in the header file
`gsl_sf_expint.h'.
Function: int gsl_sf_Shi_e(double x, gsl_sf_result * result)
These routines compute the integral Shi(x) = \int_0^x dt \sinh(t)/t.
Function: double gsl_sf_Chi(double x)
Function: int gsl_sf_Chi_e(double x, gsl_sf_result * result)
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).
The Gamma function is defined by the following integral,
\Gamma(x) = \int_0^t dt t^{x-1} \exp(-t)
Further information on the Gamma function can be found in Abramowitz &
Stegun, Chapter 6. The functions described in this section are declared
in the header file `gsl_sf_gamma.h'.
Function: double gsl_sf_gamma(double x)
Function: int gsl_sf_gamma_e(double x, gsl_sf_result * result)
These routines compute the Gamma function \Gamma(x), subject to x
not being a negative integer. The function is computed using the real
Lanczos method. The maximum value of x such that \Gamma(x) is not
considered an overflow is given by the macro GSL_SF_GAMMA_XMAX
and is 171.0.
Function: double gsl_sf_lngamma(double x)
Function: int gsl_sf_lngamma_e(double x, gsl_sf_result * result)
These routines compute the logarithm of the Gamma function,
\log(\Gamma(x)), subject to x not a being negative
integer. 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.
Function: int gsl_sf_lngamma_sgn_e(double x, gsl_sf_result * result_lg, double * sgn)
This routine computes the sign of the gamma function and the logarithm
its magnitude, subject to x not being a negative integer. The
function is computed using the real Lanczos method. The value of the
gamma function can be reconstructed using the relation \Gamma(x) =
sgn * \exp(resultlg).
Function: double gsl_sf_gammastar(double x)
Function: int gsl_sf_gammastar_e(double x, gsl_sf_result * result)
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.
Function: double gsl_sf_gammainv(double x)
Function: int gsl_sf_gammainv_e(double x, gsl_sf_result * result)
These routines compute the reciprocal of the gamma function,
1/\Gamma(x) using the real Lanczos method.
This routine computes \log(\Gamma(z)) for complex z=z_r+i
z_i and z not a negative integer, 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.
Function: double gsl_sf_taylorcoeff(int n, double x)
Function: int gsl_sf_taylorcoeff_e(int n, double x, gsl_sf_result * result)
These routines compute the Taylor coefficient x^n / n! for
x >= 0,
n >= 0.
Function: double gsl_sf_fact(unsigned int n)
Function: int gsl_sf_fact_e(unsigned int n, gsl_sf_result * result)
These routines compute the factorial n!. The factorial is
related to the Gamma function by n! = \Gamma(n+1).
Function: double gsl_sf_doublefact(unsigned int n)
Function: int gsl_sf_doublefact_e(unsigned int n, gsl_sf_result * result)
These routines compute the double factorial n!! = n(n-2)(n-4) \dots.
Function: double gsl_sf_lnfact(unsigned int n)
Function: int gsl_sf_lnfact_e(unsigned int n, gsl_sf_result * result)
These routines compute the logarithm of the factorial of n,
\log(n!). The algorithm is faster than computing
\ln(\Gamma(n+1)) via gsl_sf_lngamma for n < 170,
but defers for larger n.
Function: double gsl_sf_lndoublefact(unsigned int n)
Function: int gsl_sf_lndoublefact_e(unsigned int n, gsl_sf_result * result)
These routines compute the logarithm of the double factorial of n,
\log(n!!).
Function: double gsl_sf_choose(unsigned int n, unsigned int m)
Function: int gsl_sf_choose_e(unsigned int n, unsigned int m, gsl_sf_result * result)
These routines compute the combinatorial factor n choose m
= n!/(m!(n-m)!)
Function: double gsl_sf_lnchoose(unsigned int n, unsigned int m)
Function: int gsl_sf_lnchoose_e(unsigned int n, unsigned int m, gsl_sf_result * result)
These routines compute the logarithm of n choose m. This is
equivalent to the sum \log(n!) - \log(m!) - \log((n-m)!).
Function: double gsl_sf_poch(double a, double x)
Function: int gsl_sf_poch_e(double a, double x, gsl_sf_result * result)
These routines compute the Pochhammer symbol (a)_x := \Gamma(a +
x)/\Gamma(x), subject to a and a+x not being negative
integers. The Pochhammer symbol is also known as the Apell symbol.
Function: double gsl_sf_lnpoch(double a, double x)
Function: int gsl_sf_lnpoch_e(double a, double x, gsl_sf_result * result)
These routines compute the logarithm of the Pochhammer symbol,
\log((a)_x) = \log(\Gamma(a + x)/\Gamma(a)) for a > 0,
a+x > 0.
Function: int gsl_sf_lnpoch_sgn_e(double a, double x, gsl_sf_result * result, double * sgn)
These routines compute the sign of the Pochhammer symbol and the
logarithm of its magnitude. The computed parameters are result =
\log(|(a)_x|) and sgn = sgn((a)_x) where (a)_x :=
\Gamma(a + x)/\Gamma(a), subject to a, a+x not being
negative integers.
Function: double gsl_sf_pochrel(double a, double x)
Function: int gsl_sf_pochrel_e(double a, double x, gsl_sf_result * result)
These routines compute the relative Pochhammer symbol ((a,x) -
1)/x where (a,x) = (a)_x := \Gamma(a + x)/\Gamma(a).
Function: double gsl_sf_gamma_inc_Q(double a, double x)
Function: int gsl_sf_gamma_inc_Q_e(double a, double x, gsl_sf_result * result)
These routines compute the normalized incomplete Gamma Function
P(a,x) = 1/\Gamma(a) \int_x\infty dt t^{a-1} \exp(-t)
for a > 0,
x >= 0.
Function: double gsl_sf_gamma_inc_P(double a, double x)
Function: int gsl_sf_gamma_inc_P_e(double a, double x, gsl_sf_result * result)
These routines compute the complementary normalized incomplete Gamma Function
P(a,x) = 1/\Gamma(a) \int_0^x dt t^{a-1} \exp(-t)
for a > 0,
x >= 0.
Function: double gsl_sf_beta(double a, double b)
Function: int gsl_sf_beta_e(double a, double b, gsl_sf_result * result)
These routines compute the Beta Function, B(a,b) =
\Gamma(a)\Gamma(b)/\Gamma(a+b) for a > 0, b > 0.
Function: double gsl_sf_lnbeta(double a, double b)
Function: int gsl_sf_lnbeta_e(double a, double b, gsl_sf_result * result)
These routines compute the logarithm of the Beta Function, \log(B(a,b))
for a > 0, b > 0.
Function: double gsl_sf_beta_inc(double a, double b, double x)
Function: int gsl_sf_beta_inc_e(double a, double b, double x, gsl_sf_result * result)
These routines compute the normalize incomplete Beta function
B_x(a,b)/B(a,b) for a > 0, b > 0, and
0 <= x <= 1.
The Gegenbauer polynomials are defined in Abramowitz & Stegun, Chapter
22, where they are known as Ultraspherical polynomials. The functions
described in this section are declared in the header file
`gsl_sf_gegenbauer.h'.
Function: int gsl_sf_hyperg_0F1_e(double c, double x, gsl_sf_result * result)
These routines compute the hypergeometric function
0F1(c,x).
Function: double gsl_sf_hyperg_1F1_int(int m, int n, double x)
Function: int gsl_sf_hyperg_1F1_int_e(int m, int n, double x, gsl_sf_result * result)
These routines compute the confluent hypergeometric function
1F1(m,n,x) = M(m,n,x) for integer parameters m, n.
Function: double gsl_sf_hyperg_1F1(double a, double b, double x)
Function: int gsl_sf_hyperg_1F1_e(double a, double b, double x, gsl_sf_result * result)
These routines compute the confluent hypergeometric function
1F1(a,b,x) = M(a,b,x) for general parameters a, b.
Function: double gsl_sf_hyperg_U_int(int m, int n, double x)
Function: int gsl_sf_hyperg_U_int_e(int m, int n, double x, gsl_sf_result * result)
These routines compute the confluent hypergeometric function
U(m,n,x) for integer parameters m, n.
Function: int gsl_sf_hyperg_U_int_e10_e(int m, int n, double x, gsl_sf_result_e10 * result)
This routine computes the confluent hypergeometric function
U(m,n,x) for integer parameters m, n using the
gsl_sf_result_e10 type to return a result with extended range.
Function: double gsl_sf_hyperg_U(double a, double b, double x)
Function: int gsl_sf_hyperg_U_e(double a, double b, double x)
These routines compute the confluent hypergeometric function U(a,b,x).
Function: int gsl_sf_hyperg_U_e10_e(double a, double b, double x, gsl_sf_result_e10 * result)
This routine computes the confluent hypergeometric function
U(a,b,x) using the gsl_sf_result_e10 type to return a
result with extended range.
Function: double gsl_sf_hyperg_2F1(double a, double b, double c, double x)
Function: int gsl_sf_hyperg_2F1_e(double a, double b, double c, double x, gsl_sf_result * result)
These routines compute the Gauss hypergeometric function
2F1(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.
These routines compute the renormalized Gauss hypergeometric function
2F1(a_R + i a_I, a_R - i a_I, c, x) / \Gamma(c) for |x| < 1.
Function: double gsl_sf_hyperg_2F0(double a, double b, double x)
Function: int gsl_sf_hyperg_2F0_e(double a, double b, double x, gsl_sf_result * result)
These routines compute the hypergeometric function
2F0(a,b,x). The series representation
is a divergent hypergeometric series. However, for x < 0 we
have
2F0(a,b,x) = (-1/x)^a U(a,1+a-b,-1/x)
The Laguerre polynomials are defined in terms of confluent
hypergeometric functions as
L^a_n(x) = ((a+1)_n / n!) 1F1(-n,a+1,x). These functions are
declared in the header file `gsl_sf_laguerre.h'.
Function: double gsl_sf_laguerre_1(double a, double x)
Function: double gsl_sf_laguerre_2(double a, double x)
Function: double gsl_sf_laguerre_3(double a, double x)
Function: int gsl_sf_laguerre_1_e(double a, double x, gsl_sf_result * result)
Function: int gsl_sf_laguerre_2_e(double a, double x, gsl_sf_result * result)
Function: int gsl_sf_laguerre_3_e(double a, double x, gsl_sf_result * result)
These routines evaluate the generalized Laguerre polynomials
L^a_1(x), L^a_2(x), L^a_3(x) using explicit
representations.
Function: double gsl_sf_laguerre_n(const int n, const double a, const double x)
Function: int gsl_sf_laguerre_n_e(int n, double a, double x, gsl_sf_result * result)
Thse routines evaluate the generalized Laguerre polynomials
L^a_n(x) for a > -1,
n >= 0.
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. The Lambert functions are
declared in the header file `gsl_sf_lambert.h'.
Function: double gsl_sf_lambert_W0(double x)
Function: int gsl_sf_lambert_W0_e(double x, gsl_sf_result * result)
These compute the principal branch of the Lambert W function, W_0(x).
Function: double gsl_sf_lambert_Wm1(double x)
Function: int gsl_sf_lambert_Wm1_e(double x, gsl_sf_result * result)
These compute the secondary real-valued branch of the Lambert W function,
W_{-1}(x).
The Legendre Functions and Legendre Polynomials are described in
Abramowitz & Stegun, Chapter 8. These functions are declared in
the header file `gsl_sf_legendre.h'.
The following functions compute the associated Legendre Polynomials
P_l^m(x). Note that this function grows combinatorially with
l and can overflow for l larger than about 150. There is
no trouble for small m, but overflow occurs when m and
l are both large. Rather than allow overflows, these functions
refuse to calculate P_l^m(x) and return GSL_EOVRFLW when
they can sense that l and m are too big.
If you want to calculate a spherical harmonic, then do not use
these functions. Instead use gsl_sf_legendre_sphPlm() below,
which uses a similar recursion, but with the normalized functions.
Function: double gsl_sf_legendre_Plm(int l, int m, double x)
Function: int gsl_sf_legendre_Plm_e(int l, int m, double x, gsl_sf_result * result)
These routines compute the associated Legendre polynomial
P_l^m(x) for
m >= 0,
l >= m,
|x| <= 1.
Function: int gsl_sf_legendre_Plm_array(int lmax, int m, double x, double result_array[])
This function computes an array of Legendre polynomials
P_l^m(x) for
m >= 0,
l = |m|, ..., lmax,
|x| <= 1.
Function: double gsl_sf_legendre_sphPlm(int l, int m, double x)
Function: int gsl_sf_legendre_sphPlm_e(int l, int m, double x, gsl_sf_result * result)
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).
Function: int gsl_sf_legendre_sphPlm_array(int lmax, int m, double x, double result_array[])
This function computes an array of normalized associated Legendre functions
$\sqrt{(2l+1)/(4\pi)} \sqrt{(l-m)!/(l+m)!} P_l^m(x)$
for
m >= 0,
l = |m|, ..., lmax,
|x| <= 1.0
Function: int gsl_sf_legendre_array_size(const int lmax, const int m)
This functions returns the size of result_array[] needed for the array
versions of P_l^m(x), lmax - m + 1.
The following spherical functions are specializations of Legendre
functions which give the regular eigenfunctions of the Laplacian on a
3-dimensional hyperbolic space H3d. Of particular interest is
the flat limit, \lambda \to \infty, \eta \to 0,
\lambda\eta fixed.
Function: int gsl_sf_legendre_H3d_0_e(double lambda, double eta, gsl_sf_result * result)
These routines compute the zeroth radial eigenfunction of the Laplacian on the
3-dimensional hyperbolic space,
L^{H3d}_0(\lambda,\eta) := \sin(\lambda\eta)/(\lambda\sinh(\eta))
for
\eta >= 0.
In the flat limit this takes the form
L^{H3d}_0(\lambda,\eta) = j_0(\lambda\eta)
Function: int gsl_sf_legendre_H3d_1_e(double lambda, double eta, gsl_sf_result * result)
These routines compute the first radial eigenfunction of the Laplacian on
the 3-dimensional hyperbolic space,
L^{H3d}_1(\lambda,\eta) := 1/\sqrt{\lambda^2 + 1} \sin(\lambda \eta)/(\lambda \sinh(\eta)) (\coth(\eta) - \lambda \cot(\lambda\eta))
for
\eta >= 0.
In the flat limit this takes the form
L^{H3d}_1(\lambda,\eta) = j_1(\lambda\eta).
Function: int gsl_sf_legendre_H3d_e(int l, double lambda, double eta, gsl_sf_result * result)
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).
Function: int gsl_sf_legendre_H3d_array(int lmax, double lambda, double eta, double result_array[])
This function computes an array of radial eigenfunctions
L^{H3d}_l(\lambda, \eta)
for
0 <= l <= lmax.
Information on the properties of the Logarithm function can be found in
Abramowitz & Stegun, Chapter 4. The functions described in this section
are declared in the header file `gsl_sf_log.h'.
Function: double gsl_sf_log(double x)
Function: int gsl_sf_log_e(double x, gsl_sf_result * result)
These routines compute the logarithm of x, \log(x), for
x > 0.
Function: double gsl_sf_log_abs(double x)
Function: int gsl_sf_log_abs_e(double x, gsl_sf_result * result)
These routines compute the logarithm of the magnitude of x,
\log(|x|), for x \ne 0.
This routine computes the complex logarithm of z = z_r + i
z_i. The results are returned as lnr, theta such that
\exp(lnr + i \theta) = z_r + i z_i, where \theta lies in
the range [-\pi,\pi].
Function: double gsl_sf_log_1plusx(double x)
Function: int gsl_sf_log_1plusx_e(double x, gsl_sf_result * result)
These routines compute \log(1 + x) for x > -1 using an
algorithm that is accurate for small x.
Function: double gsl_sf_log_1plusx_mx(double x)
Function: int gsl_sf_log_1plusx_mx_e(double x, gsl_sf_result * result)
These routines compute \log(1 + x) - x for x > -1 using an
algorithm that is accurate for small x.
The following functions are equivalent to the function gsl_pow_int
(see section Small integer powers) with an error estimate. These functions are
declared in the header file `gsl_sf_pow_int.h'.
Function: double gsl_sf_pow_int(double x, int n)
Function: int gsl_sf_pow_int_e(double x, int n, gsl_sf_result * result)
These routines compute the power x^n for integer n. The
power is computed using the minimum number of multiplications. For
example, x^8 is computed as ((x^2)^2)^2, requiring only 3
multiplications. For reasons of efficiency, these functions do not
check for overflow or underflow conditions.
The polygamma functions of order m defined by
\psi^{(m)}(x) = (d/dx)^m \psi(x) = (d/dx)^{m+1} \log(\Gamma(x)),
where \psi(x) = \Gamma'(x)/\Gamma(x) is known as the digamma function.
These functions are declared in the header file `gsl_sf_psi.h'.
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.
They are declared in the header file `gsl_sf_transport.h'.
Function: double gsl_sf_transport_2(double x)
Function: int gsl_sf_transport_2_e(double x, gsl_sf_result * result)
These routines compute the transport function J(2,x).
Function: double gsl_sf_transport_3(double x)
Function: int gsl_sf_transport_3_e(double x, gsl_sf_result * result)
These routines compute the transport function J(3,x).
Function: double gsl_sf_transport_4(double x)
Function: int gsl_sf_transport_4_e(double x, gsl_sf_result * result)
These routines compute the transport function J(4,x).
Function: double gsl_sf_transport_5(double x)
Function: int gsl_sf_transport_5_e(double x, gsl_sf_result * result)
These routines compute the transport function J(5,x).
The library includes its own trigonometric functions in order to provide
consistency across platforms and reliable error estimates. These
functions are declared in the header file `gsl_sf_trig.h'.
This function converts the polar coordinates (r,theta) to
rectilinear coordinates (x,y), x = r\cos(\theta),
y = r\sin(\theta).
Function: int gsl_sf_rect_to_polar(double x, double y, gsl_sf_result * r, gsl_sf_result * theta)
This function converts the rectilinear coordinates (x,y) to
polar coordinates (r,theta), such that x =
r\cos(\theta), y = r\sin(\theta). The argument theta
lies in the range [-\pi, \pi].
The Riemann zeta function is defined in Abramowitz & Stegun, Section
23.2. The functions described in this section are declared in the
header file `gsl_sf_zeta.h'.
The next program computes the same quantity using the natural form of
the function. In this case the error term result.err and return
status are not accessible.
#include <stdio.h>
#include <gsl/gsl_sf_bessel.h>
int
main (void)
{
double x = 5.0;
double expected = -0.17759677131433830434739701;
double y = gsl_sf_bessel_J0 (x);
printf("J0(5.0) = %.18f\n", y);
printf("exact = %.18f\n", expected);
return 0;
}