This package provides natural-syntax C++ wrapping for the large spectrum of special functions provided in the GNU Scientific Library (gsl). As of July 26, 2002, these wrappers correspond to release 1.2 of the gsl. The gsl itself exists as a separate external product available via the usual ups/upd mechanism. As such, the user must issue the command
before starting up either a compilation or a link. This is to insure that the environment variable $GSL_DIR is properly defined.
Tables specifying the entire list of function signatures, as well as the corresponding underlying routine in gsl, are provided. This page also explains the details of how these functions are used in C++.
The gsl Reference Manual documents these functions. There are other C sources in gsl. These are not provided with the package, but are described in the Reference Manual. ZOOM has not wrapped these for C++, and makes no claim as to whether the non-special-function portions will compile and correctly execute on ZOOM platforms under a C++ compiler. Only the Special Functions section is supported.
Structure of Classes in the Package
Syntax and Use of Most Functions
List of Basic Signatures for Each Function
List of Possible Exceptions for Each Class
Class and Struct Descriptions Generated by the doxygen Utility.
The structs are:
A second form, with
the same name but ending in
Relatively few of these functions logically ought to return something other than
a real. For example, the package defines
(In the case of routines returning arrays of results, the name of the Err
version has the word
Errors due to invalid input values or other difficulties
are dealt with via the
ZOOM Exceptions
mechanism, which if C++ exceptions are enabled can lead to throwing exceptions.
These exceptions are always derived from
Spherical Bessel Functions - regular, irregular, regular modified and irregular
modified
Bessel functions of non-integral order
Structure of Classes in the Package
The functions are grouped into several groups, each of which is implemented
as a struct containing a variety of specific functions. For instance,
there is an Bessel struct, accessed by including
SpecialFunctions/Bessel.h, which has functions like
Bessel::J0(double x). These functions are all static and the
structs have no data members.
Syntax and Use of Most Functions
Most of these functions logically ought to return a real number.
For such a function, two forms are provided. The basic form indeed
returns a double, which is the natural syntax. Most users will use
this simplest form.
#include "SpecialFunctions/Bessel.h"
using Bessel::J0;
double x = J0(5000.);
Err, returns a Dpair,
which is merely a structure holding two doubles: the usual
return value, and the maximum possible computational error.
#include "SpecialFunctions/Bessel.h"
using Bessel::J0Err;
Dpair xe = J0Err(5000.);
double x = xe.first;
double err = xe.second;
Dvector as a std::vector
of doubles. For clarity, we list here all the data types (other than the usual
primitive types, float, double and so on) used in various places in the package.
In the brief signatures below,
we will explicitly note when a return type is not the usual double.
typedef std::pair<double,double> Dpair;
typedef std::complex<double> Dcomplex;
typedef std::vector<double> Dvector;
typedef std::pair<Dcomplex,Dcomplex> Cpair;
Err tacked on before the word
Array.)
ZMxSpecFun
which in turn is derived from the basic ZMexception
class. Possible exceptions are listed (with brief explanations) in
SFExceptions.h.
List of Basic Signatures for Each Function
Airy
Bessel
Ordinary cylindrical Bessel functions - regular, irregular , regular modified
and irregular modified
Scale factor is exp{-| x |}
Zeroes of the Bessel functions
Arrays of Bessel functions for a range of order and fixed argument
ClaussenIntegral
Defined as - Integral(0,x) log(2sin(t/2)) dt
Coulomb
double* F, double* Fp, double* G, double* Gp, double* Fscale, double* Gscale )
Where the boolean return value indicates that the results are to be scaled
by exp{Fscale} and exp{Gscale} or not. Be careful with this. The result
is expressed in this way because it would otherwise overflow. Consequently,
if you simply form the product exp{Fscale} * Fp, for example, you will
trigger the very overflow it has worked so hard to avoid.
DawsonIntegral
Defined as exp(-x^2) * Integral(0,x) exp (t^2) dt
DebyeIntegral
Dn(x)
is defined as (n/x^n) * Integral(0,x) t^n / (exp(t)-1) dt
DiLogarithm
Defined as - Re [ Integral(0,x) log(1-t)/t dt ]
EllipticIntegrals
Defined as - Integral(0,pi/2) [sqrt(1-k^2*sin(t)^2)] dt
Defined as - Integral(0,pi/2) { dt /[sqrt(1-k^2*sin(t)^2)] }
Defined as - Integral(0,phi) [sqrt(1-k^2*sin(t)^2)] dt
Defined as - Integral(0,phi) { dt /[sqrt(1-k^2*sin(t)^2)] }
Defined as - Integral(0,phi) [sqrt(1-k^2*sin(t)^2)]/[1+n*sin(t)^2] dt
Definition not entirely clear!
where
RC(x,y) = 1/2 integral(0,infinity) { dt / [ sqrt(t+x) * (t+y) ] }
RD(x,y,z) = 3/2 integral(0,infinity) { dt / [ sqrt(t+x) * sqrt(t+y) * (t+z)^3/2 ] }
RF(x,y,z) = 1/2 integral(0,infinity) { dt / [ sqrt(t+x) * sqrt(t+y) * sqrt(t+z) ] }
RJ(x,y,z.p) = 3/2 integral(0,infinity) { dt / [ sqrt(t+x) * sqrt(t+y) * sqrt(t+z) * (t+p) ] }
ErrorFunction
Defined as 2/sqrt(pi) integral( x, infinity ) { exp[ -t^2 ] }
Defined as 2/sqrt(pi) integral( 0, x ) { exp[ -t^2 ] }
Defined as 1/sqrt(2*pi) ) { exp[ -t^2/2 ] }
Defined as integral( x, infinity ) ProbZ[ t ]
ExpIntegrals
E1(x) is Re ( Integral(1,infinity) exp(-xt)/t dt )
E2(x) is Re ( Integral(1,infinity) exp(-xt)/t^2 dt )
Ei(x) is the principal value of Integral(-x,infinity) exp(-t)/t dt )
Li(x) is the principal value of Integral(0,x) dt/log(t)
Si(x) is Integral(0,x) sin(t)/t dt; Shi(x) is the same for sinh;
Atanint(x) for atan.
Ci(x) is Integral(0,x) cos(t)/t dt; Chi(x) is the same for cosh
Ei3(x) is Integral(0,x) [ dt exp{ -t^3 } ]
Exponential
ExpM1(x) is exp { x } -1
RelExp(x) and RelExp2(x) are defined as (exp(x)-1)/x and (exp(x)-1-x)/x^2
These compute the function and its associated uncertainty for a given
uncertainty in the argument(s).
FermiDirac
FDj(x) is [1/Gamma(j+1)] Integral(0,infinity) t^j / (exp(t-x)+1)
dt
Gamma
Defined as Gamma(x) divided by the Stirling approximation to Gamma(x),
which is sqrt(2 PI/x) (x/e)^x
x^n/n!
Gamma(a+x)/Gamma(a)
[Gamma(a+x)/Gamma(a) - 1] / x
(1/Gamma(a)) Integral(x,infinity) t^(a-1) exp(-t) dt
(1/Gamma(a)) Integral(0,x,) avt^(a-1) exp(-t) dt
Beta_x(a,b)/Beta(a,b)
The first part of the result is ln | Gamma(x) |; the second part of the result
is plus or minus one indicating the sign of Gamma(x)
(for negative non-integers, Gamma(x) can be negative).
Gegenbauer
Hyperbolic
Note. There are no inverse hyperbolic functions for complex argument.
HyperGeometric
IntegerPowers
Pow6 ( double x ) -- Pow7 ( double x )
-- Pow8 ( double x ) -- Pow9 ( double x )
JacobiElliptics
The three values sn(u,m), cn(u,m) and dn(m) are returned via non-const
pointer arguments.
Users calling JacobiErr() should examine how the pairs of values and errors
are returned.
Laguerre
Legendre
Logarithm
Defined as ln(1+x)
Defined as ln(1+x)-x
PolyGamma
The poly-gamma function is defined as
psi( m, x ) = ( d / dx )[ ( m + 1 )
ln ( Gamma( x ) ] where psi( 0, x ) is called the di-gamma function.
Synchrotron
Synchrotron1 is defined as x * Integral(x,infinity)[ K5/3 (t) dt ]
Synchrotron2 is defined as x * K2/3 (x)
where Kn(x) is defined earlier.
Transport
where Tn = Integral(0,x)[ dt (t^n * e^t)/(e^t-1)^2 ]
Trigonometry
where Sinc(x) = sin(pi*x)/(pi*x)
Convert (r,theta) to (x,y)
Convert (x,y) to (r,theta)
Add or subtract n * 2pi so as to bring theta into the range [-pi,pi]
Add or subtract n * 2pi so as to bring theta into the range [0,2*pi]
Wigner
Note that 3J symbols and Clebsch Gordan coefficients have different
normalizations and different phase conventions.
int two_jg, int two_jh, int two_ji )
Zeta
Riemann zeta(s) = Sum[1/(k^s)] {k = 1,infinity}
Hurwicz zeta(s,q) = Sum[1/(k+q)^s] {k = 0,infinity}
eta(s) = [1-2^(1-s)]zeta(s)
| Class | Possible Exception |
| Airy | GSL_EUNDRFLW, GSL_EOVRFLW |
| Bessel | GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW, GSL_EINVAL |
| Chebyshev | GSL_ENOMEM, GSL_EFAILED |
| ClausenIntegral | None |
| Coulomb | GSL_EDOM, GSL_EMAXITER, GSL_EUNDRFLW, GSL_EOVRFLW, GSL_ERUNAWAY, GSL_ELOSS, GSL_ESANITY |
| DawsonIntegral | GSL_EUNDRFLW |
| DebyeIntegral | GSL_EDOM, GSL_EUNDRFLW |
| DiLogarithm | GSL_EMAXITER |
| EllipticIntegrals | GSL_EDOM |
| ErrorFunction | GSL_EUNDRFLW |
| ExpIntegrals | GSL_EDOM, GSL_EUNDRFLW, GSL_EOVRFLW |
| Exponential | GSL_EOVRFLW, GSL_EUNDRFLW |
| FermiDirac | GSL_EUNDRFLW, GSL_EOVRFLW |
| Gamma | GSL_EDOM, GSL_EROUND, GSL_ELOSS, GSL_EOVRFLW, GSL_EUNDRFLW |
| Gegenbauer | GSL_EDOM |
| HyperGeometric | GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW |
| IntegerPowers | None |
| JacobiElliptics | GSL_EDOM |
| Laguerre | GSL_EDOM, GSL_EOVRFLW |
| Legendre | GSL_EDOM, GSL_EOVRFLW |
| Logarithm | GSL_EDOM |
| PolyGamma | GSL_EDOM, GSL_ELOSS |
| Synchrotron | GSL_EDOM, GSL_EUNDRFLW |
| Transport | GSL_EDOM, GSL_EUNDRFLW |
| Trigonometry | GSL_EDOM, GSL_ELOSS GSL_EOVRFLW |
| Wigner | GSL_EDOM, GSL_EOVRFLW |
| Zeta | GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW |
 
The meanings and assigned severities of these exceptions are:
 
| Exception | Explanation | Severity |
| GSL_EDOM | Input value outside valid domain | ERROR |
| GSL_EFAILED | Generic (unclassified) failure | SEVERE |
| GSL_EINVAL | Invalid argument supplied | ERROR |
| GSL_ELOSS | Failure due to excessive precision loss | WARNING |
| GSL_EMAXITER | Maximum number of iterations exceeded | ERROR |
| GSL_ENOMEM | malloc failed | SEVERE |
| GSL_EOVRFLW | Overflow | WARNING |
| GSL_EROUND | Failure due to excessive roundoff | WARNING |
| GSL_ERUNAWAY | Iterative process out of control | ERROR |
| GSL_ESANITY | Sanity check failed | SEVERE |
| GSL_EUNDRFLW | Underflow | WARNING |