/* acoshf.c
*
* Inverse hyperbolic cosine
*
*
*
* SYNOPSIS:
*
* float x, y, acoshf();
*
* y = acoshf( x );
*
*
*
* DESCRIPTION:
*
* Returns inverse hyperbolic cosine of argument.
*
* If 1 <= x < 1.5, a polynomial approximation
*
* sqrt(z) * P(z)
*
* where z = x-1, is used. Otherwise,
*
* acosh(x) = log( x + sqrt( (x-1)(x+1) ).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 1,3 100000 1.8e-7 3.9e-8
* IEEE 1,2000 100000 3.0e-8
*
*
* ERROR MESSAGES:
*
* message condition value returned
* acoshf domain |x| < 1 0.0
*
*/
/* airy.c
*
* Airy function
*
*
*
* SYNOPSIS:
*
* float x, ai, aip, bi, bip;
* int airyf();
*
* airyf( x, _&ai, _&aip, _&bi, _&bip );
*
*
*
* DESCRIPTION:
*
* Solution of the differential equation
*
* y"(x) = xy.
*
* The function returns the two independent solutions Ai, Bi
* and their first derivatives Ai'(x), Bi'(x).
*
* Evaluation is by power series summation for small x,
* by rational minimax approximations for large x.
*
*
*
* ACCURACY:
* Error criterion is absolute when function <= 1, relative
* when function > 1, except * denotes relative error criterion.
* For large negative x, the absolute error increases as x^1.5.
* For large positive x, the relative error increases as x^1.5.
*
* Arithmetic domain function # trials peak rms
* IEEE -10, 0 Ai 50000 7.0e-7 1.2e-7
* IEEE 0, 10 Ai 50000 9.9e-6* 6.8e-7*
* IEEE -10, 0 Ai' 50000 2.4e-6 3.5e-7
* IEEE 0, 10 Ai' 50000 8.7e-6* 6.2e-7*
* IEEE -10, 10 Bi 100000 2.2e-6 2.6e-7
* IEEE -10, 10 Bi' 50000 2.2e-6 3.5e-7
*
*/
/* asinf.c
*
* Inverse circular sine
*
*
*
* SYNOPSIS:
*
* float x, y, asinf();
*
* y = asinf( x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle between -pi/2 and +pi/2 whose sine is x.
*
* A polynomial of the form x + x**3 P(x**2)
* is used for |x| in the interval [0, 0.5]. If |x| > 0.5 it is
* transformed by the identity
*
* asin(x) = pi/2 - 2 asin( sqrt( (1-x)/2 ) ).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -1, 1 100000 2.5e-7 5.0e-8
*
*
* ERROR MESSAGES:
*
* message condition value returned
* asinf domain |x| > 1 0.0
*
*/
/* acosf()
*
* Inverse circular cosine
*
*
*
* SYNOPSIS:
*
* float x, y, acosf();
*
* y = acosf( x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle between -pi/2 and +pi/2 whose cosine
* is x.
*
* Analytically, acos(x) = pi/2 - asin(x). However if |x| is
* near 1, there is cancellation error in subtracting asin(x)
* from pi/2. Hence if x < -0.5,
*
* acos(x) = pi - 2.0 * asin( sqrt((1+x)/2) );
*
* or if x > +0.5,
*
* acos(x) = 2.0 * asin( sqrt((1-x)/2) ).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -1, 1 100000 1.4e-7 4.2e-8
*
*
* ERROR MESSAGES:
*
* message condition value returned
* acosf domain |x| > 1 0.0
*/
/* asinhf.c
*
* Inverse hyperbolic sine
*
*
*
* SYNOPSIS:
*
* float x, y, asinhf();
*
* y = asinhf( x );
*
*
*
* DESCRIPTION:
*
* Returns inverse hyperbolic sine of argument.
*
* If |x| < 0.5, the function is approximated by a rational
* form x + x**3 P(x)/Q(x). Otherwise,
*
* asinh(x) = log( x + sqrt(1 + x*x) ).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -3,3 100000 2.4e-7 4.1e-8
*
*/
/* atanf.c
*
* Inverse circular tangent
* (arctangent)
*
*
*
* SYNOPSIS:
*
* float x, y, atanf();
*
* y = atanf( x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle between -pi/2 and +pi/2 whose tangent
* is x.
*
* Range reduction is from four intervals into the interval
* from zero to tan( pi/8 ). A polynomial approximates
* the function in this basic interval.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10, 10 100000 1.9e-7 4.1e-8
*
*/
/* atan2f()
*
* Quadrant correct inverse circular tangent
*
*
*
* SYNOPSIS:
*
* float x, y, z, atan2f();
*
* z = atan2f( y, x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle whose tangent is y/x.
* Define compile time symbol ANSIC = 1 for ANSI standard,
* range -PI < z <= +PI, args (y,x); else ANSIC = 0 for range
* 0 to 2PI, args (x,y).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10, 10 100000 1.9e-7 4.1e-8
* See atan.c.
*
*/
/* atanhf.c
*
* Inverse hyperbolic tangent
*
*
*
* SYNOPSIS:
*
* float x, y, atanhf();
*
* y = atanhf( x );
*
*
*
* DESCRIPTION:
*
* Returns inverse hyperbolic tangent of argument in the range
* MINLOGF to MAXLOGF.
*
* If |x| < 0.5, a polynomial approximation is used.
* Otherwise,
* atanh(x) = 0.5 * log( (1+x)/(1-x) ).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -1,1 100000 1.4e-7 3.1e-8
*
*/
/* bdtrf.c
*
* Binomial distribution
*
*
*
* SYNOPSIS:
*
* int k, n;
* float p, y, bdtrf();
*
* y = bdtrf( k, n, p );
*
*
*
* DESCRIPTION:
*
* Returns the sum of the terms 0 through k of the Binomial
* probability density:
*
* k
* -- ( n ) j n-j
* > ( ) p (1-p)
* -- ( j )
* j=0
*
* The terms are not summed directly; instead the incomplete
* beta integral is employed, according to the formula
*
* y = bdtr( k, n, p ) = incbet( n-k, k+1, 1-p ).
*
* The arguments must be positive, with p ranging from 0 to 1.
*
*
*
* ACCURACY:
*
* Relative error (p varies from 0 to 1):
* arithmetic domain # trials peak rms
* IEEE 0,100 2000 6.9e-5 1.1e-5
*
* ERROR MESSAGES:
*
* message condition value returned
* bdtrf domain k < 0 0.0
* n < k
* x < 0, x > 1
*
*/
/* bdtrcf()
*
* Complemented binomial distribution
*
*
*
* SYNOPSIS:
*
* int k, n;
* float p, y, bdtrcf();
*
* y = bdtrcf( k, n, p );
*
*
*
* DESCRIPTION:
*
* Returns the sum of the terms k+1 through n of the Binomial
* probability density:
*
* n
* -- ( n ) j n-j
* > ( ) p (1-p)
* -- ( j )
* j=k+1
*
* The terms are not summed directly; instead the incomplete
* beta integral is employed, according to the formula
*
* y = bdtrc( k, n, p ) = incbet( k+1, n-k, p ).
*
* The arguments must be positive, with p ranging from 0 to 1.
*
*
*
* ACCURACY:
*
* Relative error (p varies from 0 to 1):
* arithmetic domain # trials peak rms
* IEEE 0,100 2000 6.0e-5 1.2e-5
*
* ERROR MESSAGES:
*
* message condition value returned
* bdtrcf domain x<0, x>1, n<k 0.0
*/
/* bdtrif()
*
* Inverse binomial distribution
*
*
*
* SYNOPSIS:
*
* int k, n;
* float p, y, bdtrif();
*
* p = bdtrf( k, n, y );
*
*
*
* DESCRIPTION:
*
* Finds the event probability p such that the sum of the
* terms 0 through k of the Binomial probability density
* is equal to the given cumulative probability y.
*
* This is accomplished using the inverse beta integral
* function and the relation
*
* 1 - p = incbi( n-k, k+1, y ).
*
*
*
*
* ACCURACY:
*
* Relative error (p varies from 0 to 1):
* arithmetic domain # trials peak rms
* IEEE 0,100 2000 3.5e-5 3.3e-6
*
* ERROR MESSAGES:
*
* message condition value returned
* bdtrif domain k < 0, n <= k 0.0
* x < 0, x > 1
*
*/
/* betaf.c
*
* Beta function
*
*
*
* SYNOPSIS:
*
* float a, b, y, betaf();
*
* y = betaf( a, b );
*
*
*
* DESCRIPTION:
*
* - -
* | (a) | (b)
* beta( a, b ) = -----------.
* -
* | (a+b)
*
* For large arguments the logarithm of the function is
* evaluated using lgam(), then exponentiated.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,30 10000 4.0e-5 6.0e-6
* IEEE -20,0 10000 4.9e-3 5.4e-5
*
* ERROR MESSAGES:
*
* message condition value returned
* betaf overflow log(beta) > MAXLOG 0.0
* a or b <0 integer 0.0
*
*/
/* cbrtf.c
*
* Cube root
*
*
*
* SYNOPSIS:
*
* float x, y, cbrtf();
*
* y = cbrtf( x );
*
*
*
* DESCRIPTION:
*
* Returns the cube root of the argument, which may be negative.
*
* Range reduction involves determining the power of 2 of
* the argument. A polynomial of degree 2 applied to the
* mantissa, and multiplication by the cube root of 1, 2, or 4
* approximates the root to within about 0.1%. Then Newton's
* iteration is used to converge to an accurate result.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,1e38 100000 7.6e-8 2.7e-8
*
*/
/* chbevlf.c
*
* Evaluate Chebyshev series
*
*
*
* SYNOPSIS:
*
* int N;
* float x, y, coef[N], chebevlf();
*
* y = chbevlf( x, coef, N );
*
*
*
* DESCRIPTION:
*
* Evaluates the series
*
* N-1
* - '
* y = > coef[i] T (x/2)
* - i
* i=0
*
* of Chebyshev polynomials Ti at argument x/2.
*
* Coefficients are stored in reverse order, i.e. the zero
* order term is last in the array. Note N is the number of
* coefficients, not the order.
*
* If coefficients are for the interval a to b, x must
* have been transformed to x -> 2(2x - b - a)/(b-a) before
* entering the routine. This maps x from (a, b) to (-1, 1),
* over which the Chebyshev polynomials are defined.
*
* If the coefficients are for the inverted interval, in
* which (a, b) is mapped to (1/b, 1/a), the transformation
* required is x -> 2(2ab/x - b - a)/(b-a). If b is infinity,
* this becomes x -> 4a/x - 1.
*
*
*
* SPEED:
*
* Taking advantage of the recurrence properties of the
* Chebyshev polynomials, the routine requires one more
* addition per loop than evaluating a nested polynomial of
* the same degree.
*
*/
/* chdtrf.c
*
* Chi-square distribution
*
*
*
* SYNOPSIS:
*
* float df, x, y, chdtrf();
*
* y = chdtrf( df, x );
*
*
*
* DESCRIPTION:
*
* Returns the area under the left hand tail (from 0 to x)
* of the Chi square probability density function with
* v degrees of freedom.
*
*
* inf.
* -
* 1 | | v/2-1 -t/2
* P( x | v ) = ----------- | t e dt
* v/2 - | |
* 2 | (v/2) -
* x
*
* where x is the Chi-square variable.
*
* The incomplete gamma integral is used, according to the
* formula
*
* y = chdtr( v, x ) = igam( v/2.0, x/2.0 ).
*
*
* The arguments must both be positive.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,100 5000 3.2e-5 5.0e-6
*
* ERROR MESSAGES:
*
* message condition value returned
* chdtrf domain x < 0 or v < 1 0.0
*/
/* chdtrcf()
*
* Complemented Chi-square distribution
*
*
*
* SYNOPSIS:
*
* float v, x, y, chdtrcf();
*
* y = chdtrcf( v, x );
*
*
*
* DESCRIPTION:
*
* Returns the area under the right hand tail (from x to
* infinity) of the Chi square probability density function
* with v degrees of freedom:
*
*
* inf.
* -
* 1 | | v/2-1 -t/2
* P( x | v ) = ----------- | t e dt
* v/2 - | |
* 2 | (v/2) -
* x
*
* where x is the Chi-square variable.
*
* The incomplete gamma integral is used, according to the
* formula
*
* y = chdtr( v, x ) = igamc( v/2.0, x/2.0 ).
*
*
* The arguments must both be positive.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,100 5000 2.7e-5 3.2e-6
*
* ERROR MESSAGES:
*
* message condition value returned
* chdtrc domain x < 0 or v < 1 0.0
*/
/* chdtrif()
*
* Inverse of complemented Chi-square distribution
*
*
*
* SYNOPSIS:
*
* float df, x, y, chdtrif();
*
* x = chdtrif( df, y );
*
*
*
*
* DESCRIPTION:
*
* Finds the Chi-square argument x such that the integral
* from x to infinity of the Chi-square density is equal
* to the given cumulative probability y.
*
* This is accomplished using the inverse gamma integral
* function and the relation
*
* x/2 = igami( df/2, y );
*
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,100 10000 2.2e-5 8.5e-7
*
* ERROR MESSAGES:
*
* message condition value returned
* chdtri domain y < 0 or y > 1 0.0
* v < 1
*
*/
/* clogf.c
*
* Complex natural logarithm
*
*
*
* SYNOPSIS:
*
* void clogf();
* cmplxf z, w;
*
* clogf( &z, &w );
*
*
*
* DESCRIPTION:
*
* Returns complex logarithm to the base e (2.718...) of
* the complex argument x.
*
* If z = x + iy, r = sqrt( x**2 + y**2 ),
* then
* w = log(r) + i arctan(y/x).
*
* The arctangent ranges from -PI to +PI.
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10,+10 30000 1.9e-6 6.2e-8
*
* Larger relative error can be observed for z near 1 +i0.
* In IEEE arithmetic the peak absolute error is 3.1e-7.
*
*/
/* cexpf()
*
* Complex exponential function
*
*
*
* SYNOPSIS:
*
* void cexpf();
* cmplxf z, w;
*
* cexpf( &z, &w );
*
*
*
* DESCRIPTION:
*
* Returns the exponential of the complex argument z
* into the complex result w.
*
* If
* z = x + iy,
* r = exp(x),
*
* then
*
* w = r cos y + i r sin y.
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10,+10 30000 1.4e-7 4.5e-8
*
*/
/* csinf()
*
* Complex circular sine
*
*
*
* SYNOPSIS:
*
* void csinf();
* cmplxf z, w;
*
* csinf( &z, &w );
*
*
*
* DESCRIPTION:
*
* If
* z = x + iy,
*
* then
*
* w = sin x cosh y + i cos x sinh y.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10,+10 30000 1.9e-7 5.5e-8
*
*/
/* ccosf()
*
* Complex circular cosine
*
*
*
* SYNOPSIS:
*
* void ccosf();
* cmplxf z, w;
*
* ccosf( &z, &w );
*
*
*
* DESCRIPTION:
*
* If
* z = x + iy,
*
* then
*
* w = cos x cosh y - i sin x sinh y.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10,+10 30000 1.8e-7 5.5e-8
*/
/* ctanf()
*
* Complex circular tangent
*
*
*
* SYNOPSIS:
*
* void ctanf();
* cmplxf z, w;
*
* ctanf( &z, &w );
*
*
*
* DESCRIPTION:
*
* If
* z = x + iy,
*
* then
*
* sin 2x + i sinh 2y
* w = --------------------.
* cos 2x + cosh 2y
*
* On the real axis the denominator is zero at odd multiples
* of PI/2. The denominator is evaluated by its Taylor
* series near these points.
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10,+10 30000 3.3e-7 5.1e-8
*/
/* ccotf()
*
* Complex circular cotangent
*
*
*
* SYNOPSIS:
*
* void ccotf();
* cmplxf z, w;
*
* ccotf( &z, &w );
*
*
*
* DESCRIPTION:
*
* If
* z = x + iy,
*
* then
*
* sin 2x - i sinh 2y
* w = --------------------.
* cosh 2y - cos 2x
*
* On the real axis, the denominator has zeros at even
* multiples of PI/2. Near these points it is evaluated
* by a Taylor series.
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10,+10 30000 3.6e-7 5.7e-8
* Also tested by ctan * ccot = 1 + i0.
*/
/* casinf()
*
* Complex circular arc sine
*
*
*
* SYNOPSIS:
*
* void casinf();
* cmplxf z, w;
*
* casinf( &z, &w );
*
*
*
* DESCRIPTION:
*
* Inverse complex sine:
*
* 2
* w = -i clog( iz + csqrt( 1 - z ) ).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10,+10 30000 1.1e-5 1.5e-6
* Larger relative error can be observed for z near zero.
*
*/
/* cacosf()
*
* Complex circular arc cosine
*
*
*
* SYNOPSIS:
*
* void cacosf();
* cmplxf z, w;
*
* cacosf( &z, &w );
*
*
*
* DESCRIPTION:
*
*
* w = arccos z = PI/2 - arcsin z.
*
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10,+10 30000 9.2e-6 1.2e-6
*
*/
/* catan()
*
* Complex circular arc tangent
*
*
*
* SYNOPSIS:
*
* void catan();
* cmplxf z, w;
*
* catan( &z, &w );
*
*
*
* DESCRIPTION:
*
* If
* z = x + iy,
*
* then
* 1 ( 2x )
* Re w = - arctan(-----------) + k PI
* 2 ( 2 2)
* (1 - x - y )
*
* ( 2 2)
* 1 (x + (y+1) )
* Im w = - log(------------)
* 4 ( 2 2)
* (x + (y-1) )
*
* Where k is an arbitrary integer.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10,+10 30000 2.3e-6 5.2e-8
*
*/
/* cmplxf.c
*
* Complex number arithmetic
*
*
*
* SYNOPSIS:
*
* typedef struct {
* float r; real part
* float i; imaginary part
* }cmplxf;
*
* cmplxf *a, *b, *c;
*
* caddf( a, b, c ); c = b + a
* csubf( a, b, c ); c = b - a
* cmulf( a, b, c ); c = b * a
* cdivf( a, b, c ); c = b / a
* cnegf( c ); c = -c
* cmovf( b, c ); c = b
*
*
*
* DESCRIPTION:
*
* Addition:
* c.r = b.r + a.r
* c.i = b.i + a.i
*
* Subtraction:
* c.r = b.r - a.r
* c.i = b.i - a.i
*
* Multiplication:
* c.r = b.r * a.r - b.i * a.i
* c.i = b.r * a.i + b.i * a.r
*
* Division:
* d = a.r * a.r + a.i * a.i
* c.r = (b.r * a.r + b.i * a.i)/d
* c.i = (b.i * a.r - b.r * a.i)/d
* ACCURACY:
*
* In DEC arithmetic, the test (1/z) * z = 1 had peak relative
* error 3.1e-17, rms 1.2e-17. The test (y/z) * (z/y) = 1 had
* peak relative error 8.3e-17, rms 2.1e-17.
*
* Tests in the rectangle {-10,+10}:
* Relative error:
* arithmetic function # trials peak rms
* IEEE cadd 30000 5.9e-8 2.6e-8
* IEEE csub 30000 6.0e-8 2.6e-8
* IEEE cmul 30000 1.1e-7 3.7e-8
* IEEE cdiv 30000 2.1e-7 5.7e-8
*/
/* cabsf()
*
* Complex absolute value
*
*
*
* SYNOPSIS:
*
* float cabsf();
* cmplxf z;
* float a;
*
* a = cabsf( &z );
*
*
*
* DESCRIPTION:
*
*
* If z = x + iy
*
* then
*
* a = sqrt( x**2 + y**2 ).
*
* Overflow and underflow are avoided by testing the magnitudes
* of x and y before squaring. If either is outside half of
* the floating point full scale range, both are rescaled.
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10,+10 30000 1.2e-7 3.4e-8
*/
/* csqrtf()
*
* Complex square root
*
*
*
* SYNOPSIS:
*
* void csqrtf();
* cmplxf z, w;
*
* csqrtf( &z, &w );
*
*
*
* DESCRIPTION:
*
*
* If z = x + iy, r = |z|, then
*
* 1/2
* Im w = [ (r - x)/2 ] ,
*
* Re w = y / 2 Im w.
*
*
* Note that -w is also a square root of z. The solution
* reported is always in the upper half plane.
*
* Because of the potential for cancellation error in r - x,
* the result is sharpened by doing a Heron iteration
* (see sqrt.c) in complex arithmetic.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10,+10 100000 1.8e-7 4.2e-8
*
*/
/* coshf.c
*
* Hyperbolic cosine
*
*
*
* SYNOPSIS:
*
* float x, y, coshf();
*
* y = coshf( x );
*
*
*
* DESCRIPTION:
*
* Returns hyperbolic cosine of argument in the range MINLOGF to
* MAXLOGF.
*
* cosh(x) = ( exp(x) + exp(-x) )/2.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE +-MAXLOGF 100000 1.2e-7 2.8e-8
*
*
* ERROR MESSAGES:
*
* message condition value returned
* coshf overflow |x| > MAXLOGF MAXNUMF
*
*
*/
/* dawsnf.c
*
* Dawson's Integral
*
*
*
* SYNOPSIS:
*
* float x, y, dawsnf();
*
* y = dawsnf( x );
*
*
*
* DESCRIPTION:
*
* Approximates the integral
*
* x
* -
* 2 | | 2
* dawsn(x) = exp( -x ) | exp( t ) dt
* | |
* -
* 0
*
* Three different rational approximations are employed, for
* the intervals 0 to 3.25; 3.25 to 6.25; and 6.25 up.
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,10 50000 4.4e-7 6.3e-8
*
*
*/
/* ellief.c
*
* Incomplete elliptic integral of the second kind
*
*
*
* SYNOPSIS:
*
* float phi, m, y, ellief();
*
* y = ellief( phi, m );
*
*
*
* DESCRIPTION:
*
* Approximates the integral
*
*
* phi
* -
* | |
* | 2
* E(phi\m) = | sqrt( 1 - m sin t ) dt
* |
* | |
* -
* 0
*
* of amplitude phi and modulus m, using the arithmetic -
* geometric mean algorithm.
*
*
*
* ACCURACY:
*
* Tested at random arguments with phi in [0, 2] and m in
* [0, 1].
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,2 10000 4.5e-7 7.4e-8
*
*
*/
/* ellikf.c
*
* Incomplete elliptic integral of the first kind
*
*
*
* SYNOPSIS:
*
* float phi, m, y, ellikf();
*
* y = ellikf( phi, m );
*
*
*
* DESCRIPTION:
*
* Approximates the integral
*
*
*
* phi
* -
* | |
* | dt
* F(phi\m) = | ------------------
* | 2
* | | sqrt( 1 - m sin t )
* -
* 0
*
* of amplitude phi and modulus m, using the arithmetic -
* geometric mean algorithm.
*
*
*
*
* ACCURACY:
*
* Tested at random points with phi in [0, 2] and m in
* [0, 1].
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,2 10000 2.9e-7 5.8e-8
*
*
*/
/* ellpef.c
*
* Complete elliptic integral of the second kind
*
*
*
* SYNOPSIS:
*
* float m1, y, ellpef();
*
* y = ellpef( m1 );
*
*
*
* DESCRIPTION:
*
* Approximates the integral
*
*
* pi/2
* -
* | | 2
* E(m) = | sqrt( 1 - m sin t ) dt
* | |
* -
* 0
*
* Where m = 1 - m1, using the approximation
*
* P(x) - x log x Q(x).
*
* Though there are no singularities, the argument m1 is used
* rather than m for compatibility with ellpk().
*
* E(1) = 1; E(0) = pi/2.
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0, 1 30000 1.1e-7 3.9e-8
*
*
* ERROR MESSAGES:
*
* message condition value returned
* ellpef domain x<0, x>1 0.0
*
*/
/* ellpjf.c
*
* Jacobian Elliptic Functions
*
*
*
* SYNOPSIS:
*
* float u, m, sn, cn, dn, phi;
* int ellpj();
*
* ellpj( u, m, _&sn, _&cn, _&dn, _&phi );
*
*
*
* DESCRIPTION:
*
*
* Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m),
* and dn(u|m) of parameter m between 0 and 1, and real
* argument u.
*
* These functions are periodic, with quarter-period on the
* real axis equal to the complete elliptic integral
* ellpk(1.0-m).
*
* Relation to incomplete elliptic integral:
* If u = ellik(phi,m), then sn(u|m) = sin(phi),
* and cn(u|m) = cos(phi). Phi is called the amplitude of u.
*
* Computation is by means of the arithmetic-geometric mean
* algorithm, except when m is within 1e-9 of 0 or 1. In the
* latter case with m close to 1, the approximation applies
* only for phi < pi/2.
*
* ACCURACY:
*
* Tested at random points with u between 0 and 10, m between
* 0 and 1.
*
* Absolute error (* = relative error):
* arithmetic function # trials peak rms
* IEEE sn 10000 1.7e-6 2.2e-7
* IEEE cn 10000 1.6e-6 2.2e-7
* IEEE dn 10000 1.4e-3 1.9e-5
* IEEE phi 10000 3.9e-7* 6.7e-8*
*
* Peak error observed in consistency check using addition
* theorem for sn(u+v) was 4e-16 (absolute). Also tested by
* the above relation to the incomplete elliptic integral.
* Accuracy deteriorates when u is large.
*
*/
/* ellpkf.c
*
* Complete elliptic integral of the first kind
*
*
*
* SYNOPSIS:
*
* float m1, y, ellpkf();
*
* y = ellpkf( m1 );
*
*
*
* DESCRIPTION:
*
* Approximates the integral
*
*
*
* pi/2
* -
* | |
* | dt
* K(m) = | ------------------
* | 2
* | | sqrt( 1 - m sin t )
* -
* 0
*
* where m = 1 - m1, using the approximation
*
* P(x) - log x Q(x).
*
* The argument m1 is used rather than m so that the logarithmic
* singularity at m = 1 will be shifted to the origin; this
* preserves maximum accuracy.
*
* K(0) = pi/2.
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,1 30000 1.3e-7 3.4e-8
*
* ERROR MESSAGES:
*
* message condition value returned
* ellpkf domain x<0, x>1 0.0
*
*/
/* exp10f.c
*
* Base 10 exponential function
* (Common antilogarithm)
*
*
*
* SYNOPSIS:
*
* float x, y, exp10f();
*
* y = exp10f( x );
*
*
*
* DESCRIPTION:
*
* Returns 10 raised to the x power.
*
* Range reduction is accomplished by expressing the argument
* as 10**x = 2**n 10**f, with |f| < 0.5 log10(2).
* A polynomial approximates 10**f.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -38,+38 100000 9.8e-8 2.8e-8
*
* ERROR MESSAGES:
*
* message condition value returned
* exp10 underflow x < -MAXL10 0.0
* exp10 overflow x > MAXL10 MAXNUM
*
* IEEE single arithmetic: MAXL10 = 38.230809449325611792.
*
*/
/* exp2f.c
*
* Base 2 exponential function
*
*
*
* SYNOPSIS:
*
* float x, y, exp2f();
*
* y = exp2f( x );
*
*
*
* DESCRIPTION:
*
* Returns 2 raised to the x power.
*
* Range reduction is accomplished by separating the argument
* into an integer k and fraction f such that
* x k f
* 2 = 2 2.
*
* A polynomial approximates 2**x in the basic range [-0.5, 0.5].
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -127,+127 100000 1.7e-7 2.8e-8
*
*
* See exp.c for comments on error amplification.
*
*
* ERROR MESSAGES:
*
* message condition value returned
* exp underflow x < -MAXL2 0.0
* exp overflow x > MAXL2 MAXNUMF
*
* For IEEE arithmetic, MAXL2 = 127.
*/
/* expf.c
*
* Exponential function
*
*
*
* SYNOPSIS:
*
* float x, y, expf();
*
* y = expf( x );
*
*
*
* DESCRIPTION:
*
* Returns e (2.71828...) raised to the x power.
*
* Range reduction is accomplished by separating the argument
* into an integer k and fraction f such that
*
* x k f
* e = 2 e.
*
* A polynomial is used to approximate exp(f)
* in the basic range [-0.5, 0.5].
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE +- MAXLOG 100000 1.7e-7 2.8e-8
*
*
* Error amplification in the exponential function can be
* a serious matter. The error propagation involves
* exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
* which shows that a 1 lsb error in representing X produces
* a relative error of X times 1 lsb in the function.
* While the routine gives an accurate result for arguments
* that are exactly represented by a double precision
* computer number, the result contains amplified roundoff
* error for large arguments not exactly represented.
*
*
* ERROR MESSAGES:
*
* message condition value returned
* expf underflow x < MINLOGF 0.0
* expf overflow x > MAXLOGF MAXNUMF
*
*/
/* expnf.c
*
* Exponential integral En
*
*
*
* SYNOPSIS:
*
* int n;
* float x, y, expnf();
*
* y = expnf( n, x );
*
*
*
* DESCRIPTION:
*
* Evaluates the exponential integral
*
* inf.
* -
* | | -xt
* | e
* E (x) = | ---- dt.
* n | n
* | | t
* -
* 1
*
*
* Both n and x must be nonnegative.
*
* The routine employs either a power series, a continued
* fraction, or an asymptotic formula depending on the
* relative values of n and x.
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0, 30 10000 5.6e-7 1.2e-7
*
*/
/* facf.c
*
* Factorial function
*
*
*
* SYNOPSIS:
*
* float y, facf();
* int i;
*
* y = facf( i );
*
*
*
* DESCRIPTION:
*
* Returns factorial of i = 1 * 2 * 3 * ... * i.
* fac(0) = 1.0.
*
* Due to machine arithmetic bounds the largest value of
* i accepted is 33 in single precision arithmetic.
* Greater values, or negative ones,
* produce an error message and return MAXNUM.
*
*
*
* ACCURACY:
*
* For i < 34 the values are simply tabulated, and have
* full machine accuracy.
*
*/
/* fdtrf.c
*
* F distribution
*
*
*
* SYNOPSIS:
*
* int df1, df2;
* float x, y, fdtrf();
*
* y = fdtrf( df1, df2, x );
*
*
*
* DESCRIPTION:
*
* Returns the area from zero to x under the F density
* function (also known as Snedcor's density or the
* variance ratio density). This is the density
* of x = (u1/df1)/(u2/df2), where u1 and u2 are random
* variables having Chi square distributions with df1
* and df2 degrees of freedom, respectively.
*
* The incomplete beta integral is used, according to the
* formula
*
* P(x) = incbet( df1/2, df2/2, (df1*x/(df2 + df1*x) ).
*
*
* The arguments a and b are greater than zero, and x
* x is nonnegative.
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,100 5000 2.2e-5 1.1e-6
*
* ERROR MESSAGES:
*
* message condition value returned
* fdtrf domain a<0, b<0, x<0 0.0
*
*/
/* fdtrcf()
*
* Complemented F distribution
*
*
*
* SYNOPSIS:
*
* int df1, df2;
* float x, y, fdtrcf();
*
* y = fdtrcf( df1, df2, x );
*
*
*
* DESCRIPTION:
*
* Returns the area from x to infinity under the F density
* function (also known as Snedcor's density or the
* variance ratio density).
*
*
* inf.
* -
* 1 | | a-1 b-1
* 1-P(x) = ------ | t (1-t) dt
* B(a,b) | |
* -
* x
*
* (See fdtr.c.)
*
* The incomplete beta integral is used, according to the
* formula
*
* P(x) = incbet( df2/2, df1/2, (df2/(df2 + df1*x) ).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,100 5000 7.3e-5 1.2e-5
*
* ERROR MESSAGES:
*
* message condition value returned
* fdtrcf domain a<0, b<0, x<0 0.0
*
*/
/* fdtrif()
*
* Inverse of complemented F distribution
*
*
*
* SYNOPSIS:
*
* float df1, df2, x, y, fdtrif();
*
* x = fdtrif( df1, df2, y );
*
*
*
*
* DESCRIPTION:
*
* Finds the F density argument x such that the integral
* from x to infinity of the F density is equal to the
* given probability y.
*
* This is accomplished using the inverse beta integral
* function and the relations
*
* z = incbi( df2/2, df1/2, y )
* x = df2 (1-z) / (df1 z).
*
* Note: the following relations hold for the inverse of
* the uncomplemented F distribution:
*
* z = incbi( df1/2, df2/2, y )
* x = df2 z / (df1 (1-z)).
*
*
*
* ACCURACY:
*
* arithmetic domain # trials peak rms
* Absolute error:
* IEEE 0,100 5000 4.0e-5 3.2e-6
* Relative error:
* IEEE 0,100 5000 1.2e-3 1.8e-5
*
* ERROR MESSAGES:
*
* message condition value returned
* fdtrif domain y <= 0 or y > 1 0.0
* v < 1
*
*/
/* ceilf()
* floorf()
* frexpf()
* ldexpf()
*
* Single precision floating point numeric utilities
*
*
*
* SYNOPSIS:
*
* float x, y;
* float ceilf(), floorf(), frexpf(), ldexpf();
* int expnt, n;
*
* y = floorf(x);
* y = ceilf(x);
* y = frexpf( x, &expnt );
* y = ldexpf( x, n );
*
*
*
* DESCRIPTION:
*
* All four routines return a single precision floating point
* result.
*
* sfloor() returns the largest integer less than or equal to x.
* It truncates toward minus infinity.
*
* sceil() returns the smallest integer greater than or equal
* to x. It truncates toward plus infinity.
*
* sfrexp() extracts the exponent from x. It returns an integer
* power of two to expnt and the significand between 0.5 and 1
* to y. Thus x = y * 2**expn.
*
* sldexp() multiplies x by 2**n.
*
* These functions are part of the standard C run time library
* for many but not all C compilers. The ones supplied are
* written in C for either DEC or IEEE arithmetic. They should
* be used only if your compiler library does not already have
* them.
*
* The IEEE versions assume that denormal numbers are implemented
* in the arithmetic. Some modifications will be required if
* the arithmetic has abrupt rather than gradual underflow.
*/
/* fresnlf.c
*
* Fresnel integral
*
*
*
* SYNOPSIS:
*
* float x, S, C;
* void fresnlf();
*
* fresnlf( x, _&S, _&C );
*
*
* DESCRIPTION:
*
* Evaluates the Fresnel integrals
*
* x
* -
* | |
* C(x) = | cos(pi/2 t**2) dt,
* | |
* -
* 0
*
* x
* -
* | |
* S(x) = | sin(pi/2 t**2) dt.
* | |
* -
* 0
*
*
* The integrals are evaluated by power series for small x.
* For x >= 1 auxiliary functions f(x) and g(x) are employed
* such that
*
* C(x) = 0.5 + f(x) sin( pi/2 x**2 ) - g(x) cos( pi/2 x**2 )
* S(x) = 0.5 - f(x) cos( pi/2 x**2 ) - g(x) sin( pi/2 x**2 )
*
*
*
* ACCURACY:
*
* Relative error.
*
* Arithmetic function domain # trials peak rms
* IEEE S(x) 0, 10 30000 1.1e-6 1.9e-7
* IEEE C(x) 0, 10 30000 1.1e-6 2.0e-7
*/
/* gammaf.c
*
* Gamma function
*
*
*
* SYNOPSIS:
*
* float x, y, gammaf();
* extern int sgngamf;
*
* y = gammaf( x );
*
*
*
* DESCRIPTION:
*
* Returns gamma function of the argument. The result is
* correctly signed, and the sign (+1 or -1) is also
* returned in a global (extern) variable named sgngamf.
* This same variable is also filled in by the logarithmic
* gamma function lgam().
*
* Arguments between 0 and 10 are reduced by recurrence and the
* function is approximated by a polynomial function covering
* the interval (2,3). Large arguments are handled by Stirling's
* formula. Negative arguments are made positive using
* a reflection formula.
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,-33 100,000 5.7e-7 1.0e-7
* IEEE -33,0 100,000 6.1e-7 1.2e-7
*
*
*/
/* lgamf()
*
* Natural logarithm of gamma function
*
*
*
* SYNOPSIS:
*
* float x, y, lgamf();
* extern int sgngamf;
*
* y = lgamf( x );
*
*
*
* DESCRIPTION:
*
* Returns the base e (2.718...) logarithm of the absolute
* value of the gamma function of the argument.
* The sign (+1 or -1) of the gamma function is returned in a
* global (extern) variable named sgngamf.
*
* For arguments greater than 6.5, the logarithm of the gamma
* function is approximated by the logarithmic version of
* Stirling's formula. Arguments between 0 and +6.5 are reduced by
* by recurrence to the interval [.75,1.25] or [1.5,2.5] of a rational
* approximation. The cosecant reflection formula is employed for
* arguments less than zero.
*
* Arguments greater than MAXLGM = 2.035093e36 return MAXNUM and an
* error message.
*
*
*
* ACCURACY:
*
*
*
* arithmetic domain # trials peak rms
* IEEE -100,+100 500,000 7.4e-7 6.8e-8
* The error criterion was relative when the function magnitude
* was greater than one but absolute when it was less than one.
* The routine has low relative error for positive arguments.
*
* The following test used the relative error criterion.
* IEEE -2, +3 100000 4.0e-7 5.6e-8
*
*/
/* gdtrf.c
*
* Gamma distribution function
*
*
*
* SYNOPSIS:
*
* float a, b, x, y, gdtrf();
*
* y = gdtrf( a, b, x );
*
*
*
* DESCRIPTION:
*
* Returns the integral from zero to x of the gamma probability
* density function:
*
*
* x
* b -
* a | | b-1 -at
* y = ----- | t e dt
* - | |
* | (b) -
* 0
*
* The incomplete gamma integral is used, according to the
* relation
*
* y = igam( b, ax ).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,100 5000 5.8e-5 3.0e-6
*
* ERROR MESSAGES:
*
* message condition value returned
* gdtrf domain x < 0 0.0
*
*/
/* gdtrcf.c
*
* Complemented gamma distribution function
*
*
*
* SYNOPSIS:
*
* float a, b, x, y, gdtrcf();
*
* y = gdtrcf( a, b, x );
*
*
*
* DESCRIPTION:
*
* Returns the integral from x to infinity of the gamma
* probability density function:
*
*
* inf.
* b -
* a | | b-1 -at
* y = ----- | t e dt
* - | |
* | (b) -
* x
*
* The incomplete gamma integral is used, according to the
* relation
*
* y = igamc( b, ax ).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,100 5000 9.1e-5 1.5e-5
*
* ERROR MESSAGES:
*
* message condition value returned
* gdtrcf domain x < 0 0.0
*
*/
/* hyp2f1f.c
*
* Gauss hypergeometric function F
* 2 1
*
*
* SYNOPSIS:
*
* float a, b, c, x, y, hyp2f1f();
*
* y = hyp2f1f( a, b, c, x );
*
*
* DESCRIPTION:
*
*
* hyp2f1( a, b, c, x ) = F ( a, b; c; x )
* 2 1
*
* inf.
* - a(a+1)...(a+k) b(b+1)...(b+k) k+1
* = 1 + > ----------------------------- x .
* - c(c+1)...(c+k) (k+1)!
* k = 0
*
* Cases addressed are
* Tests and escapes for negative integer a, b, or c
* Linear transformation if c - a or c - b negative integer
* Special case c = a or c = b
* Linear transformation for x near +1
* Transformation for x < -0.5
* Psi function expansion if x > 0.5 and c - a - b integer
* Conditionally, a recurrence on c to make c-a-b > 0
*
* |x| > 1 is rejected.
*
* The parameters a, b, c are considered to be integer
* valued if they are within 1.0e-6 of the nearest integer.
*
* ACCURACY:
*
* Relative error (-1 < x < 1):
* arithmetic domain # trials peak rms
* IEEE 0,3 30000 5.8e-4 4.3e-6
*/
/* hypergf.c
*
* Confluent hypergeometric function
*
*
*
* SYNOPSIS:
*
* float a, b, x, y, hypergf();
*
* y = hypergf( a, b, x );
*
*
*
* DESCRIPTION:
*
* Computes the confluent hypergeometric function
*
* 1 2
* a x a(a+1) x
* F ( a,b;x ) = 1 + ---- + --------- + ...
* 1 1 b 1! b(b+1) 2!
*
* Many higher transcendental functions are special cases of
* this power series.
*
* As is evident from the formula, b must not be a negative
* integer or zero unless a is an integer with 0 >= a > b.
*
* The routine attempts both a direct summation of the series
* and an asymptotic expansion. In each case error due to
* roundoff, cancellation, and nonconvergence is estimated.
* The result with smaller estimated error is returned.
*
*
*
* ACCURACY:
*
* Tested at random points (a, b, x), all three variables
* ranging from 0 to 30.
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,5 10000 6.6e-7 1.3e-7
* IEEE 0,30 30000 1.1e-5 6.5e-7
*
* Larger errors can be observed when b is near a negative
* integer or zero. Certain combinations of arguments yield
* serious cancellation error in the power series summation
* and also are not in the region of near convergence of the
* asymptotic series. An error message is printed if the
* self-estimated relative error is greater than 1.0e-3.
*
*/
/* i0f.c
*
* Modified Bessel function of order zero
*
*
*
* SYNOPSIS:
*
* float x, y, i0();
*
* y = i0f( x );
*
*
*
* DESCRIPTION:
*
* Returns modified Bessel function of order zero of the
* argument.
*
* The function is defined as i0(x) = j0( ix ).
*
* The range is partitioned into the two intervals [0,8] and
* (8, infinity). Chebyshev polynomial expansions are employed
* in each interval.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,30 100000 4.0e-7 7.9e-8
*
*/
/* i0ef.c
*
* Modified Bessel function of order zero,
* exponentially scaled
*
*
*
* SYNOPSIS:
*
* float x, y, i0ef();
*
* y = i0ef( x );
*
*
*
* DESCRIPTION:
*
* Returns exponentially scaled modified Bessel function
* of order zero of the argument.
*
* The function is defined as i0e(x) = exp(-|x|) j0( ix ).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,30 100000 3.7e-7 7.0e-8
* See i0f().
*
*/
/* i1f.c
*
* Modified Bessel function of order one
*
*
*
* SYNOPSIS:
*
* float x, y, i1f();
*
* y = i1f( x );
*
*
*
* DESCRIPTION:
*
* Returns modified Bessel function of order one of the
* argument.
*
* The function is defined as i1(x) = -i j1( ix ).
*
* The range is partitioned into the two intervals [0,8] and
* (8, infinity). Chebyshev polynomial expansions are employed
* in each interval.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0, 30 100000 1.5e-6 1.6e-7
*
*
*/
/* i1ef.c
*
* Modified Bessel function of order one,
* exponentially scaled
*
*
*
* SYNOPSIS:
*
* float x, y, i1ef();
*
* y = i1ef( x );
*
*
*
* DESCRIPTION:
*
* Returns exponentially scaled modified Bessel function
* of order one of the argument.
*
* The function is defined as i1(x) = -i exp(-|x|) j1( ix ).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0, 30 30000 1.5e-6 1.5e-7
* See i1().
*
*/
/* igamf.c
*
* Incomplete gamma integral
*
*
*
* SYNOPSIS:
*
* float a, x, y, igamf();
*
* y = igamf( a, x );
*
*
*
* DESCRIPTION:
*
* The function is defined by
*
* x
* -
* 1 | | -t a-1
* igam(a,x) = ----- | e t dt.
* - | |
* | (a) -
* 0
*
*
* In this implementation both arguments must be positive.
* The integral is evaluated by either a power series or
* continued fraction expansion, depending on the relative
* values of a and x.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,30 20000 7.8e-6 5.9e-7
*
*/
/* igamcf()
*
* Complemented incomplete gamma integral
*
*
*
* SYNOPSIS:
*
* float a, x, y, igamcf();
*
* y = igamcf( a, x );
*
*
*
* DESCRIPTION:
*
* The function is defined by
*
*
* igamc(a,x) = 1 - igam(a,x)
*
* inf.
* -
* 1 | | -t a-1
* = ----- | e t dt.
* - | |
* | (a) -
* x
*
*
* In this implementation both arguments must be positive.
* The integral is evaluated by either a power series or
* continued fraction expansion, depending on the relative
* values of a and x.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,30 30000 7.8e-6 5.9e-7
*
*/
/* igamif()
*
* Inverse of complemented imcomplete gamma integral
*
*
*
* SYNOPSIS:
*
* float a, x, y, igamif();
*
* x = igamif( a, y );
*
*
*
* DESCRIPTION:
*
* Given y, the function finds x such that
*
* igamc( a, x ) = y.
*
* Starting with the approximate value
*
* 3
* x = a t
*
* where
*
* t = 1 - d - ndtri(y) sqrt(d)
*
* and
*
* d = 1/9a,
*
* the routine performs up to 10 Newton iterations to find the
* root of igamc(a,x) - y = 0.
*
*
* ACCURACY:
*
* Tested for a ranging from 0 to 100 and x from 0 to 1.
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,100 5000 1.0e-5 1.5e-6
*
*/
/* incbetf.c
*
* Incomplete beta integral
*
*
* SYNOPSIS:
*
* float a, b, x, y, incbetf();
*
* y = incbetf( a, b, x );
*
*
* DESCRIPTION:
*
* Returns incomplete beta integral of the arguments, evaluated
* from zero to x. The function is defined as
*
* x
* - -
* | (a+b) | | a-1 b-1
* ----------- | t (1-t) dt.
* - - | |
* | (a) | (b) -
* 0
*
* The domain of definition is 0 <= x <= 1. In this
* implementation a and b are restricted to positive values.
* The integral from x to 1 may be obtained by the symmetry
* relation
*
* 1 - incbet( a, b, x ) = incbet( b, a, 1-x ).
*
* The integral is evaluated by a continued fraction expansion.
* If a < 1, the function calls itself recursively after a
* transformation to increase a to a+1.
*
* ACCURACY:
*
* Tested at random points (a,b,x) with a and b in the indicated
* interval and x between 0 and 1.
*
* arithmetic domain # trials peak rms
* Relative error:
* IEEE 0,30 10000 3.7e-5 5.1e-6
* IEEE 0,100 10000 1.7e-4 2.5e-5
* The useful domain for relative error is limited by underflow
* of the single precision exponential function.
* Absolute error:
* IEEE 0,30 100000 2.2e-5 9.6e-7
* IEEE 0,100 10000 6.5e-5 3.7e-6
*
* Larger errors may occur for extreme ratios of a and b.
*
* ERROR MESSAGES:
* message condition value returned
* incbetf domain x<0, x>1 0.0
*/
/* incbif()
*
* Inverse of imcomplete beta integral
*
*
*
* SYNOPSIS:
*
* float a, b, x, y, incbif();
*
* x = incbif( a, b, y );
*
*
*
* DESCRIPTION:
*
* Given y, the function finds x such that
*
* incbet( a, b, x ) = y.
*
* the routine performs up to 10 Newton iterations to find the
* root of incbet(a,b,x) - y = 0.
*
*
* ACCURACY:
*
* Relative error:
* x a,b
* arithmetic domain domain # trials peak rms
* IEEE 0,1 0,100 5000 2.8e-4 8.3e-6
*
* Overflow and larger errors may occur for one of a or b near zero
* and the other large.
*/
/* ivf.c
*
* Modified Bessel function of noninteger order
*
*
*
* SYNOPSIS:
*
* float v, x, y, ivf();
*
* y = ivf( v, x );
*
*
*
* DESCRIPTION:
*
* Returns modified Bessel function of order v of the
* argument. If x is negative, v must be integer valued.
*
* The function is defined as Iv(x) = Jv( ix ). It is
* here computed in terms of the confluent hypergeometric
* function, according to the formula
*
* v -x
* Iv(x) = (x/2) e hyperg( v+0.5, 2v+1, 2x ) / gamma(v+1)
*
* If v is a negative integer, then v is replaced by -v.
*
*
* ACCURACY:
*
* Tested at random points (v, x), with v between 0 and
* 30, x between 0 and 28.
* arithmetic domain # trials peak rms
* Relative error:
* IEEE 0,15 3000 4.7e-6 5.4e-7
* Absolute error (relative when function > 1)
* IEEE 0,30 5000 8.5e-6 1.3e-6
*
* Accuracy is diminished if v is near a negative integer.
* The useful domain for relative error is limited by overflow
* of the single precision exponential function.
*
* See also hyperg.c.
*
*/
/* j0f.c
*
* Bessel function of order zero
*
*
*
* SYNOPSIS:
*
* float x, y, j0f();
*
* y = j0f( x );
*
*
*
* DESCRIPTION:
*
* Returns Bessel function of order zero of the argument.
*
* The domain is divided into the intervals [0, 2] and
* (2, infinity). In the first interval the following polynomial
* approximation is used:
*
*
* 2 2 2
* (w - r ) (w - r ) (w - r ) P(w)
* 1 2 3
*
* 2
* where w = x and the three r's are zeros of the function.
*
* In the second interval, the modulus and phase are approximated
* by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x)
* and Phase(x) = x + 1/x R(1/x^2) - pi/4. The function is
*
* j0(x) = Modulus(x) cos( Phase(x) ).
*
*
*
* ACCURACY:
*
* Absolute error:
* arithmetic domain # trials peak rms
* IEEE 0, 2 100000 1.3e-7 3.6e-8
* IEEE 2, 32 100000 1.9e-7 5.4e-8
*
*/
/* y0f.c
*
* Bessel function of the second kind, order zero
*
*
*
* SYNOPSIS:
*
* float x, y, y0f();
*
* y = y0f( x );
*
*
*
* DESCRIPTION:
*
* Returns Bessel function of the second kind, of order
* zero, of the argument.
*
* The domain is divided into the intervals [0, 2] and
* (2, infinity). In the first interval a rational approximation
* R(x) is employed to compute
*
* 2 2 2
* y0(x) = (w - r ) (w - r ) (w - r ) R(x) + 2/pi ln(x) j0(x).
* 1 2 3
*
* Thus a call to j0() is required. The three zeros are removed
* from R(x) to improve its numerical stability.
*
* In the second interval, the modulus and phase are approximated
* by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x)
* and Phase(x) = x + 1/x S(1/x^2) - pi/4. Then the function is
*
* y0(x) = Modulus(x) sin( Phase(x) ).
*
*
*
*
* ACCURACY:
*
* Absolute error, when y0(x) < 1; else relative error:
*
* arithmetic domain # trials peak rms
* IEEE 0, 2 100000 2.4e-7 3.4e-8
* IEEE 2, 32 100000 1.8e-7 5.3e-8
*
*/
/* j1f.c
*
* Bessel function of order one
*
*
*
* SYNOPSIS:
*
* float x, y, j1f();
*
* y = j1f( x );
*
*
*
* DESCRIPTION:
*
* Returns Bessel function of order one of the argument.
*
* The domain is divided into the intervals [0, 2] and
* (2, infinity). In the first interval a polynomial approximation
* 2
* (w - r ) x P(w)
* 1
* 2
* is used, where w = x and r is the first zero of the function.
*
* In the second interval, the modulus and phase are approximated
* by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x)
* and Phase(x) = x + 1/x R(1/x^2) - 3pi/4. The function is
*
* j0(x) = Modulus(x) cos( Phase(x) ).
*
*
*
* ACCURACY:
*
* Absolute error:
* arithmetic domain # trials peak rms
* IEEE 0, 2 100000 1.2e-7 2.5e-8
* IEEE 2, 32 100000 2.0e-7 5.3e-8
*
*
*/
/* y1.c
*
* Bessel function of second kind of order one
*
*
*
* SYNOPSIS:
*
* double x, y, y1();
*
* y = y1( x );
*
*
*
* DESCRIPTION:
*
* Returns Bessel function of the second kind of order one
* of the argument.
*
* The domain is divided into the intervals [0, 2] and
* (2, infinity). In the first interval a rational approximation
* R(x) is employed to compute
*
* 2
* y0(x) = (w - r ) x R(x^2) + 2/pi (ln(x) j1(x) - 1/x) .
* 1
*
* Thus a call to j1() is required.
*
* In the second interval, the modulus and phase are approximated
* by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x)
* and Phase(x) = x + 1/x S(1/x^2) - 3pi/4. Then the function is
*
* y0(x) = Modulus(x) sin( Phase(x) ).
*
*
*
*
* ACCURACY:
*
* Absolute error:
* arithmetic domain # trials peak rms
* IEEE 0, 2 100000 2.2e-7 4.6e-8
* IEEE 2, 32 100000 1.9e-7 5.3e-8
*
* (error criterion relative when |y1| > 1).
*
*/
/* jnf.c
*
* Bessel function of integer order
*
*
*
* SYNOPSIS:
*
* int n;
* float x, y, jnf();
*
* y = jnf( n, x );
*
*
*
* DESCRIPTION:
*
* Returns Bessel function of order n, where n is a
* (possibly negative) integer.
*
* The ratio of jn(x) to j0(x) is computed by backward
* recurrence. First the ratio jn/jn-1 is found by a
* continued fraction expansion. Then the recurrence
* relating successive orders is applied until j0 or j1 is
* reached.
*
* If n = 0 or 1 the routine for j0 or j1 is called
* directly.
*
*
*
* ACCURACY:
*
* Absolute error:
* arithmetic range # trials peak rms
* IEEE 0, 15 30000 3.6e-7 3.6e-8
*
*
* Not suitable for large n or x. Use jvf() instead.
*
*/
/* jvf.c
*
* Bessel function of noninteger order
*
*
*
* SYNOPSIS:
*
* float v, x, y, jvf();
*
* y = jvf( v, x );
*
*
*
* DESCRIPTION:
*
* Returns Bessel function of order v of the argument,
* where v is real. Negative x is allowed if v is an integer.
*
* Several expansions are included: the ascending power
* series, the Hankel expansion, and two transitional
* expansions for large v. If v is not too large, it
* is reduced by recurrence to a region of best accuracy.
*
* The single precision routine accepts negative v, but with
* reduced accuracy.
*
*
*
* ACCURACY:
* Results for integer v are indicated by *.
* Error criterion is absolute, except relative when |jv()| > 1.
*
* arithmetic domain # trials peak rms
* v x
* IEEE 0,125 0,125 30000 2.0e-6 2.0e-7
* IEEE -17,0 0,125 30000 1.1e-5 4.0e-7
* IEEE -100,0 0,125 3000 1.5e-4 7.8e-6
*/
/* k0f.c
*
* Modified Bessel function, third kind, order zero
*
*
*
* SYNOPSIS:
*
* float x, y, k0f();
*
* y = k0f( x );
*
*
*
* DESCRIPTION:
*
* Returns modified Bessel function of the third kind
* of order zero of the argument.
*
* The range is partitioned into the two intervals [0,8] and
* (8, infinity). Chebyshev polynomial expansions are employed
* in each interval.
*
*
*
* ACCURACY:
*
* Tested at 2000 random points between 0 and 8. Peak absolute
* error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15.
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0, 30 30000 7.8e-7 8.5e-8
*
* ERROR MESSAGES:
*
* message condition value returned
* K0 domain x <= 0 MAXNUM
*
*/
/* k0ef()
*
* Modified Bessel function, third kind, order zero,
* exponentially scaled
*
*
*
* SYNOPSIS:
*
* float x, y, k0ef();
*
* y = k0ef( x );
*
*
*
* DESCRIPTION:
*
* Returns exponentially scaled modified Bessel function
* of the third kind of order zero of the argument.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0, 30 30000 8.1e-7 7.8e-8
* See k0().
*
*/
/* k1f.c
*
* Modified Bessel function, third kind, order one
*
*
*
* SYNOPSIS:
*
* float x, y, k1f();
*
* y = k1f( x );
*
*
*
* DESCRIPTION:
*
* Computes the modified Bessel function of the third kind
* of order one of the argument.
*
* The range is partitioned into the two intervals [0,2] and
* (2, infinity). Chebyshev polynomial expansions are employed
* in each interval.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0, 30 30000 4.6e-7 7.6e-8
*
* ERROR MESSAGES:
*
* message condition value returned
* k1 domain x <= 0 MAXNUM
*
*/
/* k1ef.c
*
* Modified Bessel function, third kind, order one,
* exponentially scaled
*
*
*
* SYNOPSIS:
*
* float x, y, k1ef();
*
* y = k1ef( x );
*
*
*
* DESCRIPTION:
*
* Returns exponentially scaled modified Bessel function
* of the third kind of order one of the argument:
*
* k1e(x) = exp(x) * k1(x).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0, 30 30000 4.9e-7 6.7e-8
* See k1().
*
*/
/* knf.c
*
* Modified Bessel function, third kind, integer order
*
*
*
* SYNOPSIS:
*
* float x, y, knf();
* int n;
*
* y = knf( n, x );
*
*
*
* DESCRIPTION:
*
* Returns modified Bessel function of the third kind
* of order n of the argument.
*
* The range is partitioned into the two intervals [0,9.55] and
* (9.55, infinity). An ascending power series is used in the
* low range, and an asymptotic expansion in the high range.
*
*
*
* ACCURACY:
*
* Absolute error, relative when function > 1:
* arithmetic domain # trials peak rms
* IEEE 0,30 10000 2.0e-4 3.8e-6
*
* Error is high only near the crossover point x = 9.55
* between the two expansions used.
*/
/* log10f.c
*
* Common logarithm
*
*
*
* SYNOPSIS:
*
* float x, y, log10f();
*
* y = log10f( x );
*
*
*
* DESCRIPTION:
*
* Returns logarithm to the base 10 of x.
*
* The argument is separated into its exponent and fractional
* parts. The logarithm of the fraction is approximated by
*
* log(1+x) = x - 0.5 x**2 + x**3 P(x).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 100000 1.3e-7 3.4e-8
* IEEE 0, MAXNUMF 100000 1.3e-7 2.6e-8
*
* In the tests over the interval [0, MAXNUM], the logarithms
* of the random arguments were uniformly distributed over
* [-MAXL10, MAXL10].
*
* ERROR MESSAGES:
*
* log10f singularity: x = 0; returns -MAXL10
* log10f domain: x < 0; returns -MAXL10
* MAXL10 = 38.230809449325611792
*/
/* log2f.c
*
* Base 2 logarithm
*
*
*
* SYNOPSIS:
*
* float x, y, log2f();
*
* y = log2f( x );
*
*
*
* DESCRIPTION:
*
* Returns the base 2 logarithm of x.
*
* The argument is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the base e
* logarithm of the fraction is approximated by
*
* log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
*
* Otherwise, setting z = 2(x-1)/x+1),
*
* log(x) = z + z**3 P(z)/Q(z).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE exp(+-88) 100000 1.1e-7 2.4e-8
* IEEE 0.5, 2.0 100000 1.1e-7 3.0e-8
*
* In the tests over the interval [exp(+-88)], the logarithms
* of the random arguments were uniformly distributed.
*
* ERROR MESSAGES:
*
* log singularity: x = 0; returns MINLOGF/log(2)
* log domain: x < 0; returns MINLOGF/log(2)
*/
/* logf.c
*
* Natural logarithm
*
*
*
* SYNOPSIS:
*
* float x, y, logf();
*
* y = logf( x );
*
*
*
* DESCRIPTION:
*
* Returns the base e (2.718...) logarithm of x.
*
* The argument is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the logarithm
* of the fraction is approximated by
*
* log(1+x) = x - 0.5 x**2 + x**3 P(x)
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 100000 7.6e-8 2.7e-8
* IEEE 1, MAXNUMF 100000 2.6e-8
*
* In the tests over the interval [1, MAXNUM], the logarithms
* of the random arguments were uniformly distributed over
* [0, MAXLOGF].
*
* ERROR MESSAGES:
*
* logf singularity: x = 0; returns MINLOG
* logf domain: x < 0; returns MINLOG
*/
/* mtherr.c
*
* Library common error handling routine
*
*
*
* SYNOPSIS:
*
* char *fctnam;
* int code;
* void mtherr();
*
* mtherr( fctnam, code );
*
*
*
* DESCRIPTION:
*
* This routine may be called to report one of the following
* err