PC Professionell 2004 December

home *** CD-ROM | disk | FTP | other *** search

/ PC Professionell 2004 December / PCpro_2004_12.ISO / files / webserver / tsw / TSW_3.4.0.exe / Apache2 / perl / Cephes.pod < prev next >

Wrap

Text File | 2002-09-09 | 115.2 KB | 4,776 lines

=head1 NAME Math::Cephes - perl interface to the cephes math library =head1 SYNOPSIS use Math::Cephes qw(:all); =head1 DESCRIPTION This module provides an interface to over 150 functions of the cephes math library of Stephen Moshier. No functions are exported by default, but rather must be imported explicitly, as in use Math::Cephes qw(sin cos); There are a number of export tags defined which allow importing groups of functions: =over 4 =item use Math::Cephes qw(:constants); imports the variables $PI : 3.14159265358979323846 # pi $PIO2 : 1.57079632679489661923 # pi/2 $PIO4 : 0.785398163397448309616 # pi/4 $SQRT2 : 1.41421356237309504880 # sqrt(2) $SQRTH : 0.707106781186547524401 # sqrt(2)/2 $LOG2E : 1.4426950408889634073599 # 1/log(2) $SQ2OPI : 0.79788456080286535587989 # sqrt( 2/pi ) $LOGE2 : 0.693147180559945309417 # log(2) $LOGSQ2 : 0.346573590279972654709 # log(2)/2 $THPIO4 : 2.35619449019234492885 # 3*pi/4 $TWOOPI : 0.636619772367581343075535 # 2/pi As well, there are 4 machine-specific numbers available: $MACHEP : machine roundoff error $MAXLOG : maximum log on the machine $MINLOG : minimum log on the machine $MAXNUM : largest number represented =item use Math::Cephes qw(:trigs); imports acos: Inverse circular cosine asin: Inverse circular sine atan: Inverse circular tangent (arctangent) atan2: Quadrant correct inverse circular tangent cos: Circular cosine cosdg: Circular cosine of angle in degrees cot: Circular cotangent cotdg: Circular cotangent of argument in degrees hypot: hypotenuse associated with the sides of a right triangle radian: Degrees, minutes, seconds to radians sin: Circular sine sindg: Circular sine of angle in degrees tan: Circular tangent tandg: Circular tangent of argument in degrees cosm1: Relative error approximations for function arguments near unity =item use Math::Cephes qw(:hypers); imports acosh: Inverse hyperbolic cosine asinh: Inverse hyperbolic sine atanh: Inverse hyperbolic tangent cosh: Hyperbolic cosine sinh: Hyperbolic sine tanh: Hyperbolic tangent =item use Math::Cephes qw(:explog); imports exp: Exponential function expxx: exp(x*x) exp10: Base 10 exponential function (Common antilogarithm) exp2: Base 2 exponential function log: Natural logarithm log10: Common logarithm log2: Base 2 logarithm log1p,expm1: Relative error approximations for function arguments near unity. =item use Math::Cephes qw(:cmplx); imports new_cmplx: create a new complex number object cabs: Complex absolute value cacos: Complex circular arc cosine cacosh: Complex inverse hyperbolic cosine casin: Complex circular arc sine casinh: Complex inverse hyperbolic sine catan: Complex circular arc tangent catanh: Complex inverse hyperbolic tangent ccos: Complex circular cosine ccosh: Complex hyperbolic cosine ccot: Complex circular cotangent cexp: Complex exponential function clog: Complex natural logarithm cadd: add two complex numbers csub: subtract two complex numbers cmul: multiply two complex numbers cdiv: divide two complex numbers cmov: copy one complex number to another cneg: negate a complex number cpow: Complex power function csin: Complex circular sine csinh: Complex hyperbolic sine csqrt: Complex square root ctan: Complex circular tangent ctanh: Complex hyperbolic tangent =item use Math::Cephes qw(:utils); imports cbrt: Cube root ceil: ceil drand: Pseudorandom number generator fabs: Absolute value fac: Factorial function floor: floor frexp: frexp ldexp: multiplies x by 2**n. lrand: Pseudorandom number generator lsqrt: Integer square root pow: Power function powi: Real raised to integer power round: Round double to nearest or even integer valued double sqrt: Square root =item use Math::Cephes qw(:bessels); imports i0: Modified Bessel function of order zero i0e: Modified Bessel function of order zero, exponentially scaled i1: Modified Bessel function of order one i1e: Modified Bessel function of order one, exponentially scaled iv: Modified Bessel function of noninteger order j0: Bessel function of order zero j1: Bessel function of order one jn: Bessel function of integer order jv: Bessel function of noninteger order k0: Modified Bessel function, third kind, order zero k0e: Modified Bessel function, third kind, order zero, exponentially scaled k1: Modified Bessel function, third kind, order one k1e: Modified Bessel function, third kind, order one, exponentially scaled kn: Modified Bessel function, third kind, integer order y0: Bessel function of the second kind, order zero y1: Bessel function of second kind of order one yn: Bessel function of second kind of integer order yv: Bessel function Yv with noninteger v =item use Math::Cephes qw(:dists); imports bdtr: Binomial distribution bdtrc: Complemented binomial distribution bdtri: Inverse binomial distribution btdtr: Beta distribution chdtr: Chi-square distribution chdtrc: Complemented Chi-square distribution chdtri: Inverse of complemented Chi-square distribution fdtr: F distribution fdtrc: Complemented F distribution fdtri: Inverse of complemented F distribution gdtr: Gamma distribution function gdtrc: Complemented gamma distribution function nbdtr: Negative binomial distribution nbdtrc: Complemented negative binomial distribution nbdtri: Functional inverse of negative binomial distribution ndtr: Normal distribution function ndtri: Inverse of Normal distribution function pdtr: Poisson distribution pdtrc: Complemented poisson distribution pdtri: Inverse Poisson distribution stdtr: Student's t distribution stdtri: Functional inverse of Student's t distribution =item use Math::Cephes qw(:gammas); imports fac: Factorial function gamma: Gamma function igam: Incomplete gamma integral igamc: Complemented incomplete gamma integral igami: Inverse of complemented imcomplete gamma integral psi: Psi (digamma) function rgamma: Reciprocal gamma function =item use Math::Cephes qw(:betas); imports beta: Beta function incbet: Incomplete beta integral incbi: Inverse of imcomplete beta integral lbeta: Natural logarithm of |beta| =item use Math::Cephes qw(:elliptics); imports ellie: Incomplete elliptic integral of the second kind ellik: Incomplete elliptic integral of the first kind ellpe: Complete elliptic integral of the second kind ellpj: Jacobian Elliptic Functions ellpk: Complete elliptic integral of the first kind =item use Math::Cephes qw(:hypergeometrics); imports hyp2f0: Gauss hypergeometric function F hyp2f1: Gauss hypergeometric function F hyperg: Confluent hypergeometric function onef2: Hypergeometric function 1F2 threef0: Hypergeometric function 3F0 =item use Math::Cephes qw(:misc); imports airy: Airy function bernum: Bernoulli numbers dawsn: Dawson's Integral ei: Exponential integral erf: Error function erfc: Complementary error function expn: Exponential integral En fresnl: Fresnel integral plancki: Integral of Planck's black body radiation formula polylog: Polylogarithm function shichi: Hyperbolic sine and cosine integrals sici: Sine and cosine integrals simpson: Simpson's rule to find an integral spence: Dilogarithm struve: Struve function vecang: angle between two vectors zeta: Riemann zeta function of two arguments zetac: Riemann zeta function =item use Math::Cephes qw(:fract); imports new_fract: create a new fraction object radd: add two fractions rmul: multiply two fractions rsub: subtracttwo fractions rdiv: divide two fractions euclid: finds the greatest common divisor =back =head1 FUNCTIONS A description of the various functions available follows. =over 4 =item I<acosh>: Inverse hyperbolic cosine SYNOPSIS: # double x, y, acosh(); $y = acosh( $x ); DESCRIPTION: Returns inverse hyperbolic cosine of argument. If 1 <= x < 1.5, a rational approximation sqrt(z) * P(z)/Q(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 DEC 1,3 30000 4.2e-17 1.1e-17 IEEE 1,3 30000 4.6e-16 8.7e-17 ERROR MESSAGES: message condition value returned acosh domain |x| < 1 NAN =item I<airy>: Airy function SYNOPSIS: # double x, ai, aiprime, bi, biprime; # int airy(); ($flag, $ai, $aiprime, $bi, $biprime) = airy( $x ); 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 10000 1.6e-15 2.7e-16 IEEE 0, 10 Ai 10000 2.3e-14* 1.8e-15* IEEE -10, 0 Ai' 10000 4.6e-15 7.6e-16 IEEE 0, 10 Ai' 10000 1.8e-14* 1.5e-15* IEEE -10, 10 Bi 30000 4.2e-15 5.3e-16 IEEE -10, 10 Bi' 30000 4.9e-15 7.3e-16 DEC -10, 0 Ai 5000 1.7e-16 2.8e-17 DEC 0, 10 Ai 5000 2.1e-15* 1.7e-16* DEC -10, 0 Ai' 5000 4.7e-16 7.8e-17 DEC 0, 10 Ai' 12000 1.8e-15* 1.5e-16* DEC -10, 10 Bi 10000 5.5e-16 6.8e-17 DEC -10, 10 Bi' 7000 5.3e-16 8.7e-17 =item I<radian>: Degrees, minutes, seconds to radians SYNOPSIS: # double d, m, s, radian(); $r = radian( $d, $m, $s ); DESCRIPTION: Converts an angle of degrees, minutes, seconds to radians. =item I<hypot>: returns the hypotenuse associated with the sides of a right triangle SYNOPSIS: # double a, b, c, hypot(); $c = hypot( $a, $b ); DESCRIPTION: Calculates the hypotenuse associated with the sides of a right triangle, according to c = sqrt( a**2 + b**2) =item I<asin>: Inverse circular sine SYNOPSIS: # double x, y, asin(); $y = asin( $x ); DESCRIPTION: Returns radian angle between -pi/2 and +pi/2 whose sine is x. A rational function of the form x + x**3 P(x**2)/Q(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 DEC -1, 1 40000 2.6e-17 7.1e-18 IEEE -1, 1 10^6 1.9e-16 5.4e-17 ERROR MESSAGES: message condition value returned asin domain |x| > 1 NAN =item I<acos>: Inverse circular cosine SYNOPSIS: # double x, y, acos(); $y = acos( $x ); DESCRIPTION: Returns radian angle between 0 and pi 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 DEC -1, 1 50000 3.3e-17 8.2e-18 IEEE -1, 1 10^6 2.2e-16 6.5e-17 ERROR MESSAGES: message condition value returned asin domain |x| > 1 NAN =item I<asinh>: Inverse hyperbolic sine SYNOPSIS: # double x, y, asinh(); $y = asinh( $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 DEC -3,3 75000 4.6e-17 1.1e-17 IEEE -1,1 30000 3.7e-16 7.8e-17 IEEE 1,3 30000 2.5e-16 6.7e-17 =item I<atan>: Inverse circular tangent (arctangent) SYNOPSIS: # double x, y, atan(); $y = atan( $x ); DESCRIPTION: Returns radian angle between -pi/2 and +pi/2 whose tangent is x. Range reduction is from three intervals into the interval from zero to 0.66. The approximant uses a rational function of degree 4/5 of the form x + x**3 P(x)/Q(x). ACCURACY: Relative error: arithmetic domain # trials peak rms DEC -10, 10 50000 2.4e-17 8.3e-18 IEEE -10, 10 10^6 1.8e-16 5.0e-17 =item I<atan2>: Quadrant correct inverse circular tangent SYNOPSIS: # double x, y, z, atan2(); $z = atan2( $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 10^6 2.5e-16 6.9e-17 See atan.c. =item I<atanh>: Inverse hyperbolic tangent SYNOPSIS: # double x, y, atanh(); $y = atanh( $x ); DESCRIPTION: Returns inverse hyperbolic tangent of argument in the range MINLOG to MAXLOG. If |x| < 0.5, the rational form x + x**3 P(x)/Q(x) is employed. Otherwise, atanh(x) = 0.5 * log( (1+x)/(1-x) ). ACCURACY: Relative error: arithmetic domain # trials peak rms DEC -1,1 50000 2.4e-17 6.4e-18 IEEE -1,1 30000 1.9e-16 5.2e-17 =item I<bdtr>: Binomial distribution SYNOPSIS: # int k, n; # double p, y, bdtr(); $y = bdtr( $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: Tested at random points (a,b,p), with p between 0 and 1. a,b Relative error: arithmetic domain # trials peak rms For p between 0.001 and 1: IEEE 0,100 100000 4.3e-15 2.6e-16 See also incbet.c. ERROR MESSAGES: message condition value returned bdtr domain k < 0 0.0 n < k x < 0, x > 1 =item I<bdtrc>: Complemented binomial distribution SYNOPSIS: # int k, n; # double p, y, bdtrc(); $y = bdtrc( $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: Tested at random points (a,b,p). a,b Relative error: arithmetic domain # trials peak rms For p between 0.001 and 1: IEEE 0,100 100000 6.7e-15 8.2e-16 For p between 0 and .001: IEEE 0,100 100000 1.5e-13 2.7e-15 ERROR MESSAGES: message condition value returned bdtrc domain x<0, x>1, n<k 0.0 =item I<bdtri>: Inverse binomial distribution SYNOPSIS: # int k, n; # double p, y, bdtri(); $p = bdtr( $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: Tested at random points (a,b,p). a,b Relative error: arithmetic domain # trials peak rms For p between 0.001 and 1: IEEE 0,100 100000 2.3e-14 6.4e-16 IEEE 0,10000 100000 6.6e-12 1.2e-13 For p between 10^-6 and 0.001: IEEE 0,100 100000 2.0e-12 1.3e-14 IEEE 0,10000 100000 1.5e-12 3.2e-14 See also incbi.c. ERROR MESSAGES: message condition value returned bdtri domain k < 0, n <= k 0.0 x < 0, x > 1 =item I<beta>: Beta function SYNOPSIS: # double a, b, y, beta(); $y = beta( $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 DEC 0,30 1700 7.7e-15 1.5e-15 IEEE 0,30 30000 8.1e-14 1.1e-14 ERROR MESSAGES: message condition value returned beta overflow log(beta) > MAXLOG 0.0 a or b <0 integer 0.0 =item I<lbeta>: Natural logarithm of |beta| SYNOPSIS: # double a, b; # double lbeta( a, b ); $y = lbeta( $a, $b); =item I<btdtr>: Beta distribution SYNOPSIS: # double a, b, x, y, btdtr(); $y = btdtr( $a, $b, $x ); DESCRIPTION: Returns the area from zero to x under the beta density function: x - - | (a+b) | | a-1 b-1 P(x) = ---------- | t (1-t) dt - - | | | (a) | (b) - 0 This function is identical to the incomplete beta integral function incbet(a, b, x). The complemented function is 1 - P(1-x) = incbet( b, a, x ); ACCURACY: See incbet.c. =item I<cbrt>: Cube root SYNOPSIS: # double x, y, cbrt(); $y = cbrt( $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 three times to converge to an accurate result. ACCURACY: Relative error: arithmetic domain # trials peak rms DEC -10,10 200000 1.8e-17 6.2e-18 IEEE 0,1e308 30000 1.5e-16 5.0e-17 =item I<chdtr>: Chi-square distribution SYNOPSIS: # double v, x, y, chdtr(); $y = chdtr( $v, $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: See igam(). ERROR MESSAGES: message condition value returned chdtr domain x < 0 or v < 1 0.0 =item I<chdtrc>: Complemented Chi-square distribution SYNOPSIS: # double v, x, y, chdtrc(); $y = chdtrc( $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 = chdtrc( v, x ) = igamc( v/2.0, x/2.0 ). The arguments must both be positive. ACCURACY: See igamc(). ERROR MESSAGES: message condition value returned chdtrc domain x < 0 or v < 1 0.0 =item I<chdtri>: Inverse of complemented Chi-square distribution SYNOPSIS: # double df, x, y, chdtri(); $x = chdtri( $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: See igami.c. ERROR MESSAGES: message condition value returned chdtri domain y < 0 or y > 1 0.0 v < 1 =item I<clog>: Complex natural logarithm SYNOPSIS: # void clog(); # cmplx z, w; $z = new_cmplx(2, 3); # $z = 2 + 3 i $w = new_cmplx(); clog($z, $w ); print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $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 DEC -10,+10 7000 8.5e-17 1.9e-17 IEEE -10,+10 30000 5.0e-15 1.1e-16 Larger relative error can be observed for z near 1 +i0. In IEEE arithmetic the peak absolute error is 5.2e-16, rms absolute error 1.0e-16. =item I<cexp>: Complex exponential function SYNOPSIS: # void cexp(); # cmplx z, w; $z = new_cmplx(2, 3); # $z = 2 + 3 i $w = new_cmplx(); cexp($z, $w ); print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $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 DEC -10,+10 8700 3.7e-17 1.1e-17 IEEE -10,+10 30000 3.0e-16 8.7e-17 =item I<csin>: Complex circular sine SYNOPSIS: # void csin(); # cmplx z, w; $z = new_cmplx(2, 3); # $z = 2 + 3 i $w = new_cmplx(); csin($z, $w ); print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $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 DEC -10,+10 8400 5.3e-17 1.3e-17 IEEE -10,+10 30000 3.8e-16 1.0e-16 Also tested by csin(casin(z)) = z. =item I<ccos>: Complex circular cosine SYNOPSIS: # void ccos(); # cmplx z, w; $z = new_cmplx(2, 3); # $z = 2 + 3 i $w = new_cmplx(); ccos($z, $w ); print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $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 DEC -10,+10 8400 4.5e-17 1.3e-17 IEEE -10,+10 30000 3.8e-16 1.0e-16 =item I<ctan>: Complex circular tangent SYNOPSIS: # void ctan(); # cmplx z, w; $z = new_cmplx(2, 3); # $z = 2 + 3 i $w = new_cmplx(); ctan($z, $w ); print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $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 DEC -10,+10 5200 7.1e-17 1.6e-17 IEEE -10,+10 30000 7.2e-16 1.2e-16 Also tested by ctan * ccot = 1 and catan(ctan(z)) = z. =item I<ccot>: Complex circular cotangent SYNOPSIS: # void ccot(); # cmplx z, w; $z = new_cmplx(2, 3); # $z = 2 + 3 i $w = new_cmplx(); ccot($z, $w ); print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $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 DEC -10,+10 3000 6.5e-17 1.6e-17 IEEE -10,+10 30000 9.2e-16 1.2e-16 Also tested by ctan * ccot = 1 + i0. =item I<casin>: Complex circular arc sine SYNOPSIS: # void casin(); # cmplx z, w; $z = new_cmplx(2, 3); # $z = 2 + 3 i $w = new_cmplx(); casin($z, $w ); print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $w DESCRIPTION: Inverse complex sine: 2 w = -i clog( iz + csqrt( 1 - z ) ). ACCURACY: Relative error: arithmetic domain # trials peak rms DEC -10,+10 10100 2.1e-15 3.4e-16 IEEE -10,+10 30000 2.2e-14 2.7e-15 Larger relative error can be observed for z near zero. Also tested by csin(casin(z)) = z. =item I<cacos>: Complex circular arc cosine SYNOPSIS: # void cacos(); # cmplx z, w; $z = new_cmplx(2, 3); # $z = 2 + 3 i $w = new_cmplx(); cacos($z, $w ); print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $w DESCRIPTION: w = arccos z = PI/2 - arcsin z. ACCURACY: Relative error: arithmetic domain # trials peak rms DEC -10,+10 5200 1.6e-15 2.8e-16 IEEE -10,+10 30000 1.8e-14 2.2e-15 =item I<catan>: Complex circular arc tangent SYNOPSIS: # void catan(); # cmplx z, w; $z = new_cmplx(2, 3); # $z = 2 + 3 i $w = new_cmplx(); catan($z, $w ); print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $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 DEC -10,+10 5900 1.3e-16 7.8e-18 IEEE -10,+10 30000 2.3e-15 8.5e-17 The check catan( ctan(z) ) = z, with |x| and |y| < PI/2, had peak relative error 1.5e-16, rms relative error 2.9e-17. See also clog(). =item I<csinh>: Complex hyperbolic sine SYNOPSIS: # void csinh(); # cmplx z, w; $z = new_cmplx(2, 3); # $z = 2 + 3 i $w = new_cmplx(); csinh($z, $w ); print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $w DESCRIPTION: csinh z = (cexp(z) - cexp(-z))/2 = sinh x * cos y + i cosh x * sin y . ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE -10,+10 30000 3.1e-16 8.2e-17 =item I<casinh>: Complex inverse hyperbolic sine SYNOPSIS: # void casinh(); # cmplx z, w; $z = new_cmplx(2, 3); # $z = 2 + 3 i $w = new_cmplx(); casinh($z, $w ); print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $w print_new_cmplx($w); # prints $w as Re($w) + i Im($w) DESCRIPTION: casinh z = -i casin iz . ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE -10,+10 30000 1.8e-14 2.6e-15 =item I<ccosh>: Complex hyperbolic cosine SYNOPSIS: # void ccosh(); # cmplx z, w; $z = new_cmplx(2, 3); # $z = 2 + 3 i $w = new_cmplx(); ccosh($z, $w ); print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $w DESCRIPTION: ccosh(z) = cosh x cos y + i sinh x sin y . ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE -10,+10 30000 2.9e-16 8.1e-17 =item I<cacosh>: Complex inverse hyperbolic cosine SYNOPSIS: # void cacosh(); # cmplx z, w; $z = new_cmplx(2, 3); # $z = 2 + 3 i $w = new_cmplx(); cacosh($z, $w ); print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $w DESCRIPTION: acosh z = i acos z . ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE -10,+10 30000 1.6e-14 2.1e-15 =item I<ctanh>: Complex hyperbolic tangent SYNOPSIS: # void ctanh(); # cmplx z, w; $z = new_cmplx(2, 3); # $z = 2 + 3 i $w = new_cmplx(); ctanh($z, $w ); print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $w DESCRIPTION: tanh z = (sinh 2x + i sin 2y) / (cosh 2x + cos 2y) . ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE -10,+10 30000 1.7e-14 2.4e-16 =item I<catanh>: Complex inverse hyperbolic tangent SYNOPSIS: # void catanh(); # cmplx z, w; $z = new_cmplx(2, 3); # $z = 2 + 3 i $w = new_cmplx(); catanh($z, $w ); print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $w DESCRIPTION: Inverse tanh, equal to -i catan (iz); ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE -10,+10 30000 2.3e-16 6.2e-17 =item I<cpow>: Complex power function SYNOPSIS: # void cpow(); # cmplx a, z, w; $a = new_cmplx(5, 6); # $z = 5 + 6 i $z = new_cmplx(2, 3); # $z = 2 + 3 i $w = new_cmplx(); cpow($a, $z, $w ); print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $w DESCRIPTION: Raises complex A to the complex Zth power. Definition is per AMS55 # 4.2.8, analytically equivalent to cpow(a,z) = cexp(z clog(a)). ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE -10,+10 30000 9.4e-15 1.5e-15 =item I<cmplx>: Complex number arithmetic SYNOPSIS: # typedef struct { # double r; real part # double i; imaginary part # }cmplx; # cmplx *a, *b, *c; $a = new_cmplx(3, 5); # $a = 3 + 5 i $b = new_cmplx(2, 3); # $b = 2 + 3 i $c = new_cmplx(); cadd( $a, $b, $c ); # c = b + a csub( $a, $b, $c ); # c = b - a cmul( $a, $b, $c ); # c = b * a cdiv( $a, $b, $c ); # c = b / a cneg( $c ); # c = -c cmov( $b, $c ); # c = b print $c->{r}, ' ', $c->{i}; # prints real and imaginary parts of $c 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 DEC cadd 10000 1.4e-17 3.4e-18 IEEE cadd 100000 1.1e-16 2.7e-17 DEC csub 10000 1.4e-17 4.5e-18 IEEE csub 100000 1.1e-16 3.4e-17 DEC cmul 3000 2.3e-17 8.7e-18 IEEE cmul 100000 2.1e-16 6.9e-17 DEC cdiv 18000 4.9e-17 1.3e-17 IEEE cdiv 100000 3.7e-16 1.1e-16 =item I<cabs>: Complex absolute value SYNOPSIS: # double a, cabs(); # cmplx z; $z = new_cmplx(2, 3); # $z = 2 + 3 i $a = cabs( $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 DEC -30,+30 30000 3.2e-17 9.2e-18 IEEE -10,+10 100000 2.7e-16 6.9e-17 =item I<csqrt>: Complex square root SYNOPSIS: # void csqrt(); # cmplx z, w; $z = new_cmplx(2, 3); # $z = 2 + 3 i $w = new_cmplx(); csqrt($z, $w ); print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $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 root chosen 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 DEC -10,+10 25000 3.2e-17 9.6e-18 IEEE -10,+10 100000 3.2e-16 7.7e-17 2 Also tested by csqrt( z ) = z, and tested by arguments close to the real axis. =item I<machconst>: Globally declared constants SYNOPSIS: extern double nameofconstant; DESCRIPTION: This file contains a number of mathematical constants and also some needed size parameters of the computer arithmetic. The values are supplied as arrays of hexadecimal integers for IEEE arithmetic; arrays of octal constants for DEC arithmetic; and in a normal decimal scientific notation for other machines. The particular notation used is determined by a symbol (DEC, IBMPC, or UNK) defined in the include file mconf.h. The default size parameters are as follows. For DEC and UNK modes: MACHEP = 1.38777878078144567553E-17 2**-56 MAXLOG = 8.8029691931113054295988E1 log(2**127) MINLOG = -8.872283911167299960540E1 log(2**-128) MAXNUM = 1.701411834604692317316873e38 2**127 For IEEE arithmetic (IBMPC): MACHEP = 1.11022302462515654042E-16 2**-53 MAXLOG = 7.09782712893383996843E2 log(2**1024) MINLOG = -7.08396418532264106224E2 log(2**-1022) MAXNUM = 1.7976931348623158E308 2**1024 These lists are subject to change. =item I<cosh>: Hyperbolic cosine SYNOPSIS: # double x, y, cosh(); $y = cosh( $x ); DESCRIPTION: Returns hyperbolic cosine of argument in the range MINLOG to MAXLOG. cosh(x) = ( exp(x) + exp(-x) )/2. ACCURACY: Relative error: arithmetic domain # trials peak rms DEC +- 88 50000 4.0e-17 7.7e-18 IEEE +-MAXLOG 30000 2.6e-16 5.7e-17 ERROR MESSAGES: message condition value returned cosh overflow |x| > MAXLOG MAXNUM =item I<dawsn>: Dawson's Integral SYNOPSIS: # double x, y, dawsn(); $y = dawsn( $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 10000 6.9e-16 1.0e-16 DEC 0,10 6000 7.4e-17 1.4e-17 =item I<drand>: Pseudorandom number generator SYNOPSIS: # double y, drand(); ($flag, $y) = drand( ); DESCRIPTION: Yields a random number 1.0 <= y < 2.0. The three-generator congruential algorithm by Brian Wichmann and David Hill (BYTE magazine, March, 1987, pp 127-8) is used. The period, given by them, is 6953607871644. Versions invoked by the different arithmetic compile time options DEC, IBMPC, and MIEEE, produce approximately the same sequences, differing only in the least significant bits of the numbers. The UNK option implements the algorithm as recommended in the BYTE article. It may be used on all computers. However, the low order bits of a double precision number may not be adequately random, and may vary due to arithmetic implementation details on different computers. The other compile options generate an additional random integer that overwrites the low order bits of the double precision number. This reduces the period by a factor of two but tends to overcome the problems mentioned. =item I<ellie>: Incomplete elliptic integral of the second kind SYNOPSIS: # double phi, m, y, ellie(); $y = ellie( $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 [-10, 10] and m in [0, 1]. Relative error: arithmetic domain # trials peak rms DEC 0,2 2000 1.9e-16 3.4e-17 IEEE -10,10 150000 3.3e-15 1.4e-16 =item I<ellik>: Incomplete elliptic integral of the first kind SYNOPSIS: # double phi, m, y, ellik(); $y = ellik( $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 m in [0, 1] and phi as indicated. Relative error: arithmetic domain # trials peak rms IEEE -10,10 200000 7.4e-16 1.0e-16 =item I<ellpe>: Complete elliptic integral of the second kind SYNOPSIS: # double m1, y, ellpe(); $y = ellpe( $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 DEC 0, 1 13000 3.1e-17 9.4e-18 IEEE 0, 1 10000 2.1e-16 7.3e-17 ERROR MESSAGES: message condition value returned ellpe domain x<0, x>1 0.0 =item I<ellpj>: Jacobian Elliptic Functions SYNOPSIS: # double u, m, sn, cn, dn, phi; # int ellpj(); ($flag, $sn, $cn, $dn, $phi) = ellpj( $u, $m ); 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 DEC sn 1800 4.5e-16 8.7e-17 IEEE phi 10000 9.2e-16* 1.4e-16* IEEE sn 50000 4.1e-15 4.6e-16 IEEE cn 40000 3.6e-15 4.4e-16 IEEE dn 10000 1.3e-12 1.8e-14 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. =item I<ellpk>: Complete elliptic integral of the first kind SYNOPSIS: # double m1, y, ellpk(); $y = ellpk( $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 DEC 0,1 16000 3.5e-17 1.1e-17 IEEE 0,1 30000 2.5e-16 6.8e-17 ERROR MESSAGES: message condition value returned ellpk domain x<0, x>1 0.0 =item I<euclid>: Rational arithmetic routines SYNOPSIS: # typedef struct # { # double n; numerator # double d; denominator # }fract; $a = new_fract(3, 4); # a = 3 / 4 $b = new_fract(2, 3); # b = 2 / 3 $c = new_fract(); radd( $a, $b, $c ); # c = b + a rsub( $a, $b, $c ); # c = b - a rmul( $a, $b, $c ); # c = b * a rdiv( $a, $b, $c ); # c = b / a print $c->{n}, ' ', $c->{d}; # prints numerator and denominator of $c ($gcd, $m_reduced, $n_reduced) = euclid($m, $n); # returns the greatest common divisor of $m and $n, as well as # the result of reducing $m and $n by $gcd Arguments of the routines are pointers to the structures. The double precision numbers are assumed, without checking, to be integer valued. Overflow conditions are reported. =item I<exp>: Exponential function SYNOPSIS: # double x, y, exp(); $y = exp( $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 Pade' form 1 + 2x P(x**2)/( Q(x**2) - P(x**2) ) of degree 2/3 is used to approximate exp(f) in the basic interval [-0.5, 0.5]. ACCURACY: Relative error: arithmetic domain # trials peak rms DEC +- 88 50000 2.8e-17 7.0e-18 IEEE +- 708 40000 2.0e-16 5.6e-17 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 exp underflow x < MINLOG 0.0 exp overflow x > MAXLOG INFINITY =item I<expxx>: exp(x*x) # double x, y, expxx(); # int sign; $y = expxx( $x, $sign ); DESCRIPTION: Computes y = exp(x*x) while suppressing error amplification that would ordinarily arise from the inexactness of the exponential argument x*x. If sign < 0, exp(-x*x) is returned. If sign > 0, or omitted, exp(x*x) is returned. ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE -26.6, 26.6 10^7 3.9e-16 8.9e-17 =item I<exp10>: Base 10 exponential function (Common antilogarithm) SYNOPSIS: # double x, y, exp10(); $y = exp10( $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). The Pade' form 1 + 2x P(x**2)/( Q(x**2) - P(x**2) ) is used to approximate 10**f. ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE -307,+307 30000 2.2e-16 5.5e-17 Test result from an earlier version (2.1): DEC -38,+38 70000 3.1e-17 7.0e-18 ERROR MESSAGES: message condition value returned exp10 underflow x < -MAXL10 0.0 exp10 overflow x > MAXL10 MAXNUM DEC arithmetic: MAXL10 = 38.230809449325611792. IEEE arithmetic: MAXL10 = 308.2547155599167. =item I<exp2>: Base 2 exponential function SYNOPSIS: # double x, y, exp2(); $y = exp2( $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 Pade' form 1 + 2x P(x**2) / (Q(x**2) - x P(x**2) ) approximates 2**x in the basic range [-0.5, 0.5]. ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE -1022,+1024 30000 1.8e-16 5.4e-17 See exp.c for comments on error amplification. ERROR MESSAGES: message condition value returned exp underflow x < -MAXL2 0.0 exp overflow x > MAXL2 MAXNUM For DEC arithmetic, MAXL2 = 127. For IEEE arithmetic, MAXL2 = 1024. =item I<ei>: Exponential integral SYNOPSIS: #double x, y, ei(); $y = ei( $x ); DESCRIPTION: x - t | | e Ei(x) = -|- --- dt . | | t - -inf Not defined for x <= 0. See also expn.c. ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE 0,100 50000 8.6e-16 1.3e-16 =item I<expn>: Exponential integral En SYNOPSIS: # int n; # double x, y, expn(); $y = expn( $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 DEC 0, 30 5000 2.0e-16 4.6e-17 IEEE 0, 30 10000 1.7e-15 3.6e-16 =item I<fabs>: Absolute value SYNOPSIS: # double x, y; $y = fabs( $x ); DESCRIPTION: Returns the absolute value of the argument. =item I<fac>: Factorial function SYNOPSIS: # double y, fac(); # int i; $y = fac( $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 DEC arithmetic or 170 in IEEE 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. If i > 55, fac(i) = gamma(i+1); see gamma.c. Relative error: arithmetic domain peak IEEE 0, 170 1.4e-15 DEC 0, 33 1.4e-17 =item I<fdtr>: F distribution SYNOPSIS: # int df1, df2; # double x, y, fdtr(); $y = fdtr( $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 is nonnegative. ACCURACY: Tested at random points (a,b,x). x a,b Relative error: arithmetic domain domain # trials peak rms IEEE 0,1 0,100 100000 9.8e-15 1.7e-15 IEEE 1,5 0,100 100000 6.5e-15 3.5e-16 IEEE 0,1 1,10000 100000 2.2e-11 3.3e-12 IEEE 1,5 1,10000 100000 1.1e-11 1.7e-13 See also incbet.c. ERROR MESSAGES: message condition value returned fdtr domain a<0, b<0, x<0 0.0 =item I<fdtrc>: Complemented F distribution SYNOPSIS: # int df1, df2; # double x, y, fdtrc(); $y = fdtrc( $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 The incomplete beta integral is used, according to the formula P(x) = incbet( df2/2, df1/2, df2/(df2 + df1*x) ). ACCURACY: Tested at random points (a,b,x) in the indicated intervals. x a,b Relative error: arithmetic domain domain # trials peak rms IEEE 0,1 1,100 100000 3.7e-14 5.9e-16 IEEE 1,5 1,100 100000 8.0e-15 1.6e-15 IEEE 0,1 1,10000 100000 1.8e-11 3.5e-13 IEEE 1,5 1,10000 100000 2.0e-11 3.0e-12 See also incbet.c. ERROR MESSAGES: message condition value returned fdtrc domain a<0, b<0, x<0 0.0 =item I<fdtri>: Inverse of complemented F distribution SYNOPSIS: # int df1, df2; # double x, p, fdtri(); $x = fdtri( $df1, $df2, $p ); 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 p. This is accomplished using the inverse beta integral function and the relations z = incbi( df2/2, df1/2, p ) 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, p ) x = df2 z / (df1 (1-z)). ACCURACY: Tested at random points (a,b,p). a,b Relative error: arithmetic domain # trials peak rms For p between .001 and 1: IEEE 1,100 100000 8.3e-15 4.7e-16 IEEE 1,10000 100000 2.1e-11 1.4e-13 For p between 10^-6 and 10^-3: IEEE 1,100 50000 1.3e-12 8.4e-15 IEEE 1,10000 50000 3.0e-12 4.8e-14 See also fdtrc.c. ERROR MESSAGES: message condition value returned fdtri domain p <= 0 or p > 1 0.0 v < 1 =item I<ceil>: ceil ceil() returns the smallest integer greater than or equal to x. It truncates toward plus infinity. SYNOPSIS: # double x, y, ceil(); $y = ceil( $x ); =item I<floor>: floor floor() returns the largest integer less than or equal to x. It truncates toward minus infinity. SYNOPSIS: # double x, y, floor(); $y = floor( $x ); =item I<frexp>: frexp frexp() 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. SYNOPSIS: # double x, y, frexp(); # int expnt; ($y, $expnt) = frexp( $x ); =item I<ldexp>: multiplies x by 2**n. SYNOPSIS: # double x, y, ldexp(); # int n; $y = ldexp( $x, $n ); =item I<fresnl>: Fresnel integral SYNOPSIS: # double x, S, C; # void fresnl(); ($flag, $S, $C) = fresnl( $x ); 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 a power series for x < 1. 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 10000 2.0e-15 3.2e-16 IEEE C(x) 0, 10 10000 1.8e-15 3.3e-16 DEC S(x) 0, 10 6000 2.2e-16 3.9e-17 DEC C(x) 0, 10 5000 2.3e-16 3.9e-17 =item I<gamma>: Gamma function SYNOPSIS: # double x, y, gamma(); # extern int sgngam; $y = gamma( $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 sgngam. This variable is also filled in by the logarithmic gamma function lgam(). Arguments |x| <= 34 are reduced by recurrence and the function approximated by a rational function of degree 6/7 in the interval (2,3). Large arguments are handled by Stirling's formula. Large negative arguments are made positive using a reflection formula. ACCURACY: Relative error: arithmetic domain # trials peak rms DEC -34, 34 10000 1.3e-16 2.5e-17 IEEE -170,-33 20000 2.3e-15 3.3e-16 IEEE -33, 33 20000 9.4e-16 2.2e-16 IEEE 33, 171.6 20000 2.3e-15 3.2e-16 Error for arguments outside the test range will be larger owing to error amplification by the exponential function. =item I<lgam>: Natural logarithm of gamma function SYNOPSIS: # double x, y, lgam(); # extern int sgngam; $y = lgam( $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 sgngam. For arguments greater than 13, the logarithm of the gamma function is approximated by the logarithmic version of Stirling's formula using a polynomial approximation of degree 4. Arguments between -33 and +33 are reduced by recurrence to the interval [2,3] of a rational approximation. The cosecant reflection formula is employed for arguments less than -33. Arguments greater than MAXLGM return MAXNUM and an error message. MAXLGM = 2.035093e36 for DEC arithmetic or 2.556348e305 for IEEE arithmetic. ACCURACY: arithmetic domain # trials peak rms DEC 0, 3 7000 5.2e-17 1.3e-17 DEC 2.718, 2.035e36 5000 3.9e-17 9.9e-18 IEEE 0, 3 28000 5.4e-16 1.1e-16 IEEE 2.718, 2.556e305 40000 3.5e-16 8.3e-17 The error criterion was relative when the function magnitude was greater than one but absolute when it was less than one. The following test used the relative error criterion, though at certain points the relative error could be much higher than indicated. IEEE -200, -4 10000 4.8e-16 1.3e-16 =item I<gdtr>: Gamma distribution function SYNOPSIS: # double a, b, x, y, gdtr(); $y = gdtr( $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: See igam(). ERROR MESSAGES: message condition value returned gdtr domain x < 0 0.0 =item I<gdtrc>: Complemented gamma distribution function SYNOPSIS: # double a, b, x, y, gdtrc(); $y = gdtrc( $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: See igamc(). ERROR MESSAGES: message condition value returned gdtrc domain x < 0 0.0 =item I<hyp2f0>: Gauss hypergeometric function 2F0 SYNOPSIS: # double a, b, x, value, *err; # int type; /* determines what converging factor to use */ ($value, $err) = hyp2f0( $a, $b, $x, $type ) =item I<hyp2f1>: Gauss hypergeometric function 2F1 SYNOPSIS: # double a, b, c, x, y, hyp2f1(); $y = hyp2f1( $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-14 of the nearest integer (1.0e-13 for IEEE arithmetic). ACCURACY: Relative error (-1 < x < 1): arithmetic domain # trials peak rms IEEE -1,7 230000 1.2e-11 5.2e-14 Several special cases also tested with a, b, c in the range -7 to 7. ERROR MESSAGES: A "partial loss of precision" message is printed if the internally estimated relative error exceeds 1^-12. A "singularity" message is printed on overflow or in cases not addressed (such as x < -1). =item I<hyperg>: Confluent hypergeometric function SYNOPSIS: # double a, b, x, y, hyperg(); $y = hyperg( $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 DEC 0,30 2000 1.2e-15 1.3e-16 IEEE 0,30 30000 1.8e-14 1.1e-15 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-12. =item I<i0>: Modified Bessel function of order zero SYNOPSIS: # double x, y, i0(); $y = i0( $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 DEC 0,30 6000 8.2e-17 1.9e-17 IEEE 0,30 30000 5.8e-16 1.4e-16 =item I<i0e>: Modified Bessel function of order zero, exponentially scaled SYNOPSIS: # double x, y, i0e(); $y = i0e( $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 30000 5.4e-16 1.2e-16 See i0(). =item I<i1>: Modified Bessel function of order one SYNOPSIS: # double x, y, i1(); $y = i1( $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 DEC 0, 30 3400 1.2e-16 2.3e-17 IEEE 0, 30 30000 1.9e-15 2.1e-16 =item I<i1e>: Modified Bessel function of order one, exponentially scaled SYNOPSIS: # double x, y, i1e(); $y = i1e( $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 2.0e-15 2.0e-16 See i1(). =item I<igam>: Incomplete gamma integral SYNOPSIS: # double a, x, y, igam(); $y = igam( $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 200000 3.6e-14 2.9e-15 IEEE 0,100 300000 9.9e-14 1.5e-14 =item I<igamc>: Complemented incomplete gamma integral SYNOPSIS: # double a, x, y, igamc(); $y = igamc( $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: Tested at random a, x. a x Relative error: arithmetic domain domain # trials peak rms IEEE 0.5,100 0,100 200000 1.9e-14 1.7e-15 IEEE 0.01,0.5 0,100 200000 1.4e-13 1.6e-15 =item I<igami>: Inverse of complemented imcomplete gamma integral SYNOPSIS: # double a, x, p, igami(); $x = igami( $a, $p ); DESCRIPTION: Given p, the function finds x such that igamc( a, x ) = p. Starting with the approximate value 3 x = a t where t = 1 - d - ndtri(p) sqrt(d) and d = 1/9a, the routine performs up to 10 Newton iterations to find the root of igamc(a,x) - p = 0. ACCURACY: Tested at random a, p in the intervals indicated. a p Relative error: arithmetic domain domain # trials peak rms IEEE 0.5,100 0,0.5 100000 1.0e-14 1.7e-15 IEEE 0.01,0.5 0,0.5 100000 9.0e-14 3.4e-15 IEEE 0.5,10000 0,0.5 20000 2.3e-13 3.8e-14 =item I<incbet>: Incomplete beta integral SYNOPSIS: # double a, b, x, y, incbet(); $y = incbet( $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 or, when b*x is small, by a power series. ACCURACY: Tested at uniformly distributed random points (a,b,x) with a and b in "domain" and x between 0 and 1. Relative error arithmetic domain # trials peak rms IEEE 0,5 10000 6.9e-15 4.5e-16 IEEE 0,85 250000 2.2e-13 1.7e-14 IEEE 0,1000 30000 5.3e-12 6.3e-13 IEEE 0,10000 250000 9.3e-11 7.1e-12 IEEE 0,100000 10000 8.7e-10 4.8e-11 Outputs smaller than the IEEE gradual underflow threshold were excluded from these statistics. ERROR MESSAGES: message condition value returned incbet domain x<0, x>1 0.0 incbet underflow 0.0 =item I<incbi>: Inverse of imcomplete beta integral SYNOPSIS: # double a, b, x, y, incbi(); $x = incbi( $a, $b, $y ); DESCRIPTION: Given y, the function finds x such that incbet( a, b, x ) = y . The routine performs interval halving or 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 .5,10000 50000 5.8e-12 1.3e-13 IEEE 0,1 .25,100 100000 1.8e-13 3.9e-15 IEEE 0,1 0,5 50000 1.1e-12 5.5e-15 VAX 0,1 .5,100 25000 3.5e-14 1.1e-15 With a and b constrained to half-integer or integer values: IEEE 0,1 .5,10000 50000 5.8e-12 1.1e-13 IEEE 0,1 .5,100 100000 1.7e-14 7.9e-16 With a = .5, b constrained to half-integer or integer values: IEEE 0,1 .5,10000 10000 8.3e-11 1.0e-11 =item I<iv>: Modified Bessel function of noninteger order SYNOPSIS: # double v, x, y, iv(); $y = iv( $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. Relative error: arithmetic domain # trials peak rms DEC 0,30 2000 3.1e-15 5.4e-16 IEEE 0,30 10000 1.7e-14 2.7e-15 Accuracy is diminished if v is near a negative integer. See also hyperg.c. =item I<j0>: Bessel function of order zero SYNOPSIS: # double x, y, j0(); $y = j0( $x ); DESCRIPTION: Returns Bessel function of order zero of the argument. The domain is divided into the intervals [0, 5] and (5, infinity). In the first interval the following rational approximation is used: 2 2 (w - r ) (w - r ) P (w) / Q (w) 1 2 3 8 2 where w = x and the two r's are zeros of the function. In the second interval, the Hankel asymptotic expansion is employed with two rational functions of degree 6/6 and 7/7. ACCURACY: Absolute error: arithmetic domain # trials peak rms DEC 0, 30 10000 4.4e-17 6.3e-18 IEEE 0, 30 60000 4.2e-16 1.1e-16 =item I<y0>: Bessel function of the second kind, order zero SYNOPSIS: # double x, y, y0(); $y = y0( $x ); DESCRIPTION: Returns Bessel function of the second kind, of order zero, of the argument. The domain is divided into the intervals [0, 5] and (5, infinity). In the first interval a rational approximation R(x) is employed to compute y0(x) = R(x) + 2 * log(x) * j0(x) / PI. Thus a call to j0() is required. In the second interval, the Hankel asymptotic expansion is employed with two rational functions of degree 6/6 and 7/7. ACCURACY: Absolute error, when y0(x) < 1; else relative error: arithmetic domain # trials peak rms DEC 0, 30 9400 7.0e-17 7.9e-18 IEEE 0, 30 30000 1.3e-15 1.6e-16 =item I<j1>: Bessel function of order one SYNOPSIS: # double x, y, j1(); $y = j1( $x ); DESCRIPTION: Returns Bessel function of order one of the argument. The domain is divided into the intervals [0, 8] and (8, infinity). In the first interval a 24 term Chebyshev expansion is used. In the second, the asymptotic trigonometric representation is employed using two rational functions of degree 5/5. ACCURACY: Absolute error: arithmetic domain # trials peak rms DEC 0, 30 10000 4.0e-17 1.1e-17 IEEE 0, 30 30000 2.6e-16 1.1e-16 =item I<y1>: 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, 8] and (8, infinity). In the first interval a 25 term Chebyshev expansion is used, and a call to j1() is required. In the second, the asymptotic trigonometric representation is employed using two rational functions of degree 5/5. ACCURACY: Absolute error: arithmetic domain # trials peak rms DEC 0, 30 10000 8.6e-17 1.3e-17 IEEE 0, 30 30000 1.0e-15 1.3e-16 (error criterion relative when |y1| > 1). =item I<jn>: Bessel function of integer order SYNOPSIS: # int n; # double x, y, jn(); $y = jn( $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 DEC 0, 30 5500 6.9e-17 9.3e-18 IEEE 0, 30 5000 4.4e-16 7.9e-17 Not suitable for large n or x. Use jv() instead. =item I<jv>: Bessel function of noninteger order SYNOPSIS: # double v, x, y, jv(); $y = jv( $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 transitional expansions give 12D accuracy for v > 500. ACCURACY: Results for integer v are indicated by *, where x and v both vary from -125 to +125. Otherwise, x ranges from 0 to 125, v ranges as indicated by "domain." Error criterion is absolute, except relative when |jv()| > 1. arithmetic v domain x domain # trials peak rms IEEE 0,125 0,125 100000 4.6e-15 2.2e-16 IEEE -125,0 0,125 40000 5.4e-11 3.7e-13 IEEE 0,500 0,500 20000 4.4e-15 4.0e-16 Integer v: IEEE -125,125 -125,125 50000 3.5e-15* 1.9e-16* =item I<k0>: Modified Bessel function, third kind, order zero SYNOPSIS: # double x, y, k0(); $y = k0( $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 DEC 0, 30 3100 1.3e-16 2.1e-17 IEEE 0, 30 30000 1.2e-15 1.6e-16 ERROR MESSAGES: message condition value returned K0 domain x <= 0 MAXNUM =item I<k0e>: Modified Bessel function, third kind, order zero, exponentially scaled SYNOPSIS: # double x, y, k0e(); $y = k0e( $x ); DESCRIPTION: Returns exponentially scaled modified Bessel function of the third kind of order zero of the argument. k0e(x) = exp(x) * k0(x). ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE 0, 30 30000 1.4e-15 1.4e-16 See k0(). =item I<k1>: Modified Bessel function, third kind, order one SYNOPSIS: # double x, y, k1(); $y = k1( $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 DEC 0, 30 3300 8.9e-17 2.2e-17 IEEE 0, 30 30000 1.2e-15 1.6e-16 ERROR MESSAGES: message condition value returned k1 domain x <= 0 MAXNUM =item I<k1e>: Modified Bessel function, third kind, order one, exponentially scaled SYNOPSIS: # double x, y, k1e(); $y = k1e( $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 7.8e-16 1.2e-16 See k1(). =item I<kn>: Modified Bessel function, third kind, integer order SYNOPSIS: # double x, y, kn(); # int n; $y = kn( $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: Relative error: arithmetic domain # trials peak rms DEC 0,30 3000 1.3e-9 5.8e-11 IEEE 0,30 90000 1.8e-8 3.0e-10 Error is high only near the crossover point x = 9.55 between the two expansions used. =item I<log>: Natural logarithm SYNOPSIS: # double x, y, log(); $y = log( $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)/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 0.5, 2.0 150000 1.44e-16 5.06e-17 IEEE +-MAXNUM 30000 1.20e-16 4.78e-17 DEC 0, 10 170000 1.8e-17 6.3e-18 In the tests over the interval [+-MAXNUM], the logarithms of the random arguments were uniformly distributed over [0, MAXLOG]. ERROR MESSAGES: log singularity: x = 0; returns -INFINITY log domain: x < 0; returns NAN =item I<log10>: Common logarithm SYNOPSIS: # double x, y, log10(); $y = log10( $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)/Q(x). ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE 0.5, 2.0 30000 1.5e-16 5.0e-17 IEEE 0, MAXNUM 30000 1.4e-16 4.8e-17 DEC 1, MAXNUM 50000 2.5e-17 6.0e-18 In the tests over the interval [1, MAXNUM], the logarithms of the random arguments were uniformly distributed over [0, MAXLOG]. ERROR MESSAGES: log10 singularity: x = 0; returns -INFINITY log10 domain: x < 0; returns NAN =item I<log2>: Base 2 logarithm SYNOPSIS: # double x, y, log2(); $y = log2( $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 0.5, 2.0 30000 2.0e-16 5.5e-17 IEEE exp(+-700) 40000 1.3e-16 4.6e-17 In the tests over the interval [exp(+-700)], the logarithms of the random arguments were uniformly distributed. ERROR MESSAGES: log2 singularity: x = 0; returns -INFINITY log2 domain: x < 0; returns NAN =item I<lrand>: Pseudorandom number generator SYNOPSIS: long y, lrand(); $y = lrand( ); DESCRIPTION: Yields a long integer random number. The three-generator congruential algorithm by Brian Wichmann and David Hill (BYTE magazine, March, 1987, pp 127-8) is used. The period, given by them, is 6953607871644. =item I<lsqrt>: Integer square root SYNOPSIS: long x, y; long lsqrt(); $y = lsqrt( $x ); DESCRIPTION: Returns a long integer square root of the long integer argument. The computation is by binary long division. The largest possible result is lsqrt(2,147,483,647) = 46341. If x < 0, the square root of |x| is returned, and an error message is printed. ACCURACY: An extra, roundoff, bit is computed; hence the result is the nearest integer to the actual square root. NOTE: only DEC arithmetic is currently supported. =item I<mtherr>: Library common error handling routine SYNOPSIS: char *fctnam; # int code; # int mtherr(); mtherr( $fctnam, $code ); DESCRIPTION: This routine may be called to report one of the following error conditions (in the include file mconf.h). Mnemonic Value Significance DOMAIN 1 argument domain error SING 2 function singularity OVERFLOW 3 overflow range error UNDERFLOW 4 underflow range error TLOSS 5 total loss of precision PLOSS 6 partial loss of precision EDOM 33 Unix domain error code ERANGE 34 Unix range error code The default version of the file prints the function name, passed to it by the pointer fctnam, followed by the error condition. The display is directed to the standard output device. The routine then returns to the calling program. Users may wish to modify the program to abort by calling exit() under severe error conditions such as domain errors. Since all error conditions pass control to this function, the display may be easily changed, eliminated, or directed to an error logging device. SEE ALSO: mconf.h =item I<nbdtr>: Negative binomial distribution SYNOPSIS: # int k, n; # double p, y, nbdtr(); $y = nbdtr( $k, $n, $p ); DESCRIPTION: Returns the sum of the terms 0 through k of the negative binomial distribution: k -- ( n+j-1 ) n j > ( ) p (1-p) -- ( j ) j=0 In a sequence of Bernoulli trials, this is the probability that k or fewer failures precede the nth success. The terms are not computed individually; instead the incomplete beta integral is employed, according to the formula y = nbdtr( k, n, p ) = incbet( n, k+1, p ). The arguments must be positive, with p ranging from 0 to 1. ACCURACY: Tested at random points (a,b,p), with p between 0 and 1. a,b Relative error: arithmetic domain # trials peak rms IEEE 0,100 100000 1.7e-13 8.8e-15 See also incbet.c. =item I<nbdtrc>: Complemented negative binomial distribution SYNOPSIS: # int k, n; # double p, y, nbdtrc(); $y = nbdtrc( $k, $n, $p ); DESCRIPTION: Returns the sum of the terms k+1 to infinity of the negative binomial distribution: inf -- ( n+j-1 ) n j > ( ) p (1-p) -- ( j ) j=k+1 The terms are not computed individually; instead the incomplete beta integral is employed, according to the formula y = nbdtrc( k, n, p ) = incbet( k+1, n, 1-p ). The arguments must be positive, with p ranging from 0 to 1. ACCURACY: Tested at random points (a,b,p), with p between 0 and 1. a,b Relative error: arithmetic domain # trials peak rms IEEE 0,100 100000 1.7e-13 8.8e-15 See also incbet.c. =item I<nbdtrc>: Complemented negative binomial distribution SYNOPSIS: # int k, n; # double p, y, nbdtrc(); $y = nbdtrc( $k, $n, $p ); DESCRIPTION: Returns the sum of the terms k+1 to infinity of the negative binomial distribution: inf -- ( n+j-1 ) n j > ( ) p (1-p) -- ( j ) j=k+1 The terms are not computed individually; instead the incomplete beta integral is employed, according to the formula y = nbdtrc( k, n, p ) = incbet( k+1, n, 1-p ). The arguments must be positive, with p ranging from 0 to 1. ACCURACY: See incbet.c. =item I<nbdtri>: Functional inverse of negative binomial distribution SYNOPSIS: # int k, n; # double p, y, nbdtri(); $p = nbdtri( $k, $n, $y ); DESCRIPTION: Finds the argument p such that nbdtr(k,n,p) is equal to y. ACCURACY: Tested at random points (a,b,y), with y between 0 and 1. a,b Relative error: arithmetic domain # trials peak rms IEEE 0,100 100000 1.5e-14 8.5e-16 See also incbi.c. =item I<ndtr>: Normal distribution function SYNOPSIS: # double x, y, ndtr(); $y = ndtr( $x ); DESCRIPTION: Returns the area under the Gaussian probability density function, integrated from minus infinity to x: x - 1 | | 2 ndtr(x) = --------- | exp( - t /2 ) dt sqrt(2pi) | | - -inf. = ( 1 + erf(z) ) / 2 where z = x/sqrt(2). Computation is via the functions erf and erfc. ACCURACY: Relative error: arithmetic domain # trials peak rms DEC -13,0 8000 2.1e-15 4.8e-16 IEEE -13,0 30000 3.4e-14 6.7e-15 ERROR MESSAGES: message condition value returned erfc underflow x > 37.519379347 0.0 =item I<erf>: Error function SYNOPSIS: # double x, y, erf(); $y = erf( $x ); DESCRIPTION: The integral is x - 2 | | 2 erf(x) = -------- | exp( - t ) dt. sqrt(pi) | | - 0 The magnitude of x is limited to 9.231948545 for DEC arithmetic; 1 or -1 is returned outside this range. For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise erf(x) = 1 - erfc(x). ACCURACY: Relative error: arithmetic domain # trials peak rms DEC 0,1 14000 4.7e-17 1.5e-17 IEEE 0,1 30000 3.7e-16 1.0e-16 =item I<erfc>: Complementary error function SYNOPSIS: # double x, y, erfc(); $y = erfc( $x ); DESCRIPTION: 1 - erf(x) = inf. - 2 | | 2 erfc(x) = -------- | exp( - t ) dt sqrt(pi) | | - x For small x, erfc(x) = 1 - erf(x); otherwise rational approximations are computed. ACCURACY: Relative error: arithmetic domain # trials peak rms DEC 0, 9.2319 12000 5.1e-16 1.2e-16 IEEE 0,26.6417 30000 5.7e-14 1.5e-14 ERROR MESSAGES: message condition value returned erfc underflow x > 9.231948545 (DEC) 0.0 =item I<ndtri>: Inverse of Normal distribution function SYNOPSIS: # double x, y, ndtri(); $x = ndtri( $y ); DESCRIPTION: Returns the argument, x, for which the area under the Gaussian probability density function (integrated from minus infinity to x) is equal to y. For small arguments 0 < y < exp(-2), the program computes z = sqrt( -2.0 * log(y) ); then the approximation is x = z - log(z)/z - (1/z) P(1/z) / Q(1/z). There are two rational functions P/Q, one for 0 < y < exp(-32) and the other for y up to exp(-2). For larger arguments, w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)). ACCURACY: Relative error: arithmetic domain # trials peak rms DEC 0.125, 1 5500 9.5e-17 2.1e-17 DEC 6e-39, 0.135 3500 5.7e-17 1.3e-17 IEEE 0.125, 1 20000 7.2e-16 1.3e-16 IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17 ERROR MESSAGES: message condition value returned ndtri domain x <= 0 -MAXNUM ndtri domain x >= 1 MAXNUM =item I<pdtr>: Poisson distribution SYNOPSIS: # int k; # double m, y, pdtr(); $y = pdtr( $k, $m ); DESCRIPTION: Returns the sum of the first k terms of the Poisson distribution: k j -- -m m > e -- -- j! j=0 The terms are not summed directly; instead the incomplete gamma integral is employed, according to the relation y = pdtr( k, m ) = igamc( k+1, m ). The arguments must both be positive. ACCURACY: See igamc(). =item I<pdtrc>: Complemented poisson distribution SYNOPSIS: # int k; # double m, y, pdtrc(); $y = pdtrc( $k, $m ); DESCRIPTION: Returns the sum of the terms k+1 to infinity of the Poisson distribution: inf. j -- -m m > e -- -- j! j=k+1 The terms are not summed directly; instead the incomplete gamma integral is employed, according to the formula y = pdtrc( k, m ) = igam( k+1, m ). The arguments must both be positive. ACCURACY: See igam.c. =item I<pdtri>: Inverse Poisson distribution SYNOPSIS: # int k; # double m, y, pdtr(); $m = pdtri( $k, $y ); DESCRIPTION: Finds the Poisson variable x such that the integral from 0 to x of the Poisson density is equal to the given probability y. This is accomplished using the inverse gamma integral function and the relation m = igami( k+1, y ). ACCURACY: See igami.c. ERROR MESSAGES: message condition value returned pdtri domain y < 0 or y >= 1 0.0 k < 0 =item I<pow>: Power function SYNOPSIS: # double x, y, z, pow(); $z = pow( $x, $y ); DESCRIPTION: Computes x raised to the yth power. Analytically, x**y = exp( y log(x) ). Following Cody and Waite, this program uses a lookup table of 2**-i/16 and pseudo extended precision arithmetic to obtain an extra three bits of accuracy in both the logarithm and the exponential. ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE -26,26 30000 4.2e-16 7.7e-17 DEC -26,26 60000 4.8e-17 9.1e-18 1/26 < x < 26, with log(x) uniformly distributed. -26 < y < 26, y uniformly distributed. IEEE 0,8700 30000 1.5e-14 2.1e-15 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed. ERROR MESSAGES: message condition value returned pow overflow x**y > MAXNUM INFINITY pow underflow x**y < 1/MAXNUM 0.0 pow domain x<0 and y noninteger 0.0 =item I<powi>: Real raised to integer power SYNOPSIS: # double x, y, powi(); # int n; $y = powi( $x, $n ); DESCRIPTION: Returns argument x raised to the nth power. The routine efficiently decomposes n as a sum of powers of two. The desired power is a product of two-to-the-kth powers of x. Thus to compute the 32767 power of x requires 28 multiplications instead of 32767 multiplications. ACCURACY: Relative error: arithmetic x domain n domain # trials peak rms DEC .04,26 -26,26 100000 2.7e-16 4.3e-17 IEEE .04,26 -26,26 50000 2.0e-15 3.8e-16 IEEE 1,2 -1022,1023 50000 8.6e-14 1.6e-14 Returns MAXNUM on overflow, zero on underflow. =item I<psi>: Psi (digamma) function SYNOPSIS: # double x, y, psi(); $y = psi( $x ); DESCRIPTION: d - psi(x) = -- ln | (x) dx is the logarithmic derivative of the gamma function. For integer x, n-1 - psi(n) = -EUL + > 1/k. - k=1 This formula is used for 0 < n <= 10. If x is negative, it is transformed to a positive argument by the reflection formula psi(1-x) = psi(x) + pi cot(pi x). For general positive x, the argument is made greater than 10 using the recurrence psi(x+1) = psi(x) + 1/x. Then the following asymptotic expansion is applied: inf. B - 2k psi(x) = log(x) - 1/2x - > ------- - 2k k=1 2k x where the B2k are Bernoulli numbers. ACCURACY: Relative error (except absolute when |psi| < 1): arithmetic domain # trials peak rms DEC 0,30 2500 1.7e-16 2.0e-17 IEEE 0,30 30000 1.3e-15 1.4e-16 IEEE -30,0 40000 1.5e-15 2.2e-16 ERROR MESSAGES: message condition value returned psi singularity x integer <=0 MAXNUM =item I<rgamma>: Reciprocal gamma function SYNOPSIS: # double x, y, rgamma(); $y = rgamma( $x ); DESCRIPTION: Returns one divided by the gamma function of the argument. The function is approximated by a Chebyshev expansion in the interval [0,1]. Range reduction is by recurrence for arguments between -34.034 and +34.84425627277176174. 1/MAXNUM is returned for positive arguments outside this range. For arguments less than -34.034 the cosecant reflection formula is applied; lograrithms are employed to avoid unnecessary overflow. The reciprocal gamma function has no singularities, but overflow and underflow may occur for large arguments. These conditions return either MAXNUM or 1/MAXNUM with appropriate sign. ACCURACY: Relative error: arithmetic domain # trials peak rms DEC -30,+30 4000 1.2e-16 1.8e-17 IEEE -30,+30 30000 1.1e-15 2.0e-16 For arguments less than -34.034 the peak error is on the order of 5e-15 (DEC), excepting overflow or underflow. =item I<round>: Round double to nearest or even integer valued double SYNOPSIS: # double x, y, round(); $y = round( $x ); DESCRIPTION: Returns the nearest integer to x as a double precision floating point result. If x ends in 0.5 exactly, the nearest even integer is chosen. ACCURACY: If x is greater than 1/(2*MACHEP), its closest machine representation is already an integer, so rounding does not change it. =item I<shichi>: Hyperbolic sine and cosine integrals SYNOPSIS: # double x, Chi, Shi, shichi(); ($flag, $Shi, $Chi) = shichi( $x ); DESCRIPTION: Approximates the integrals x - | | cosh t - 1 Chi(x) = eul + ln x + | ----------- dt, | | t - 0 x - | | sinh t Shi(x) = | ------ dt | | t - 0 where eul = 0.57721566490153286061 is Euler's constant. The integrals are evaluated by power series for x < 8 and by Chebyshev expansions for x between 8 and 88. For large x, both functions approach exp(x)/2x. Arguments greater than 88 in magnitude return MAXNUM. ACCURACY: Test interval 0 to 88. Relative error: arithmetic function # trials peak rms DEC Shi 3000 9.1e-17 IEEE Shi 30000 6.9e-16 1.6e-16 Absolute error, except relative when |Chi| > 1: DEC Chi 2500 9.3e-17 IEEE Chi 30000 8.4e-16 1.4e-16 =item I<sici>: Sine and cosine integrals SYNOPSIS: # double x, Ci, Si, sici(); ($flag, $Si, $Ci) = sici( $x ); DESCRIPTION: Evaluates the integrals x - | cos t - 1 Ci(x) = eul + ln x + | --------- dt, | t - 0 x - | sin t Si(x) = | ----- dt | t - 0 where eul = 0.57721566490153286061 is Euler's constant. The integrals are approximated by rational functions. For x > 8 auxiliary functions f(x) and g(x) are employed such that Ci(x) = f(x) sin(x) - g(x) cos(x) Si(x) = pi/2 - f(x) cos(x) - g(x) sin(x) ACCURACY: Test interval = [0,50]. Absolute error, except relative when > 1: arithmetic function # trials peak rms IEEE Si 30000 4.4e-16 7.3e-17 IEEE Ci 30000 6.9e-16 5.1e-17 DEC Si 5000 4.4e-17 9.0e-18 DEC Ci 5300 7.9e-17 5.2e-18 =item I<sin>: Circular sine SYNOPSIS: # double x, y, sin(); $y = sin( $x ); DESCRIPTION: Range reduction is into intervals of pi/4. The reduction error is nearly eliminated by contriving an extended precision modular arithmetic. Two polynomial approximating functions are employed. Between 0 and pi/4 the sine is approximated by x + x**3 P(x**2). Between pi/4 and pi/2 the cosine is represented as 1 - x**2 Q(x**2). ACCURACY: Relative error: arithmetic domain # trials peak rms DEC 0, 10 150000 3.0e-17 7.8e-18 IEEE -1.07e9,+1.07e9 130000 2.1e-16 5.4e-17 ERROR MESSAGES: message condition value returned sin total loss x > 1.073741824e9 0.0 Partial loss of accuracy begins to occur at x = 2**30 = 1.074e9. The loss is not gradual, but jumps suddenly to about 1 part in 10e7. Results may be meaningless for x > 2**49 = 5.6e14. The routine as implemented flags a TLOSS error for x > 2**30 and returns 0.0. =item I<cos>: Circular cosine SYNOPSIS: # double x, y, cos(); $y = cos( $x ); DESCRIPTION: Range reduction is into intervals of pi/4. The reduction error is nearly eliminated by contriving an extended precision modular arithmetic. Two polynomial approximating functions are employed. Between 0 and pi/4 the cosine is approximated by 1 - x**2 Q(x**2). Between pi/4 and pi/2 the sine is represented as x + x**3 P(x**2). ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE -1.07e9,+1.07e9 130000 2.1e-16 5.4e-17 DEC 0,+1.07e9 17000 3.0e-17 7.2e-18 =item I<sindg>: Circular sine of angle in degrees SYNOPSIS: # double x, y, sindg(); $y = sindg( $x ); DESCRIPTION: Range reduction is into intervals of 45 degrees. Two polynomial approximating functions are employed. Between 0 and pi/4 the sine is approximated by x + x**3 P(x**2). Between pi/4 and pi/2 the cosine is represented as 1 - x**2 P(x**2). ACCURACY: Relative error: arithmetic domain # trials peak rms DEC +-1000 3100 3.3e-17 9.0e-18 IEEE +-1000 30000 2.3e-16 5.6e-17 ERROR MESSAGES: message condition value returned sindg total loss x > 8.0e14 (DEC) 0.0 x > 1.0e14 (IEEE) =item I<cosdg>: Circular cosine of angle in degrees SYNOPSIS: # double x, y, cosdg(); $y = cosdg( $x ); DESCRIPTION: Range reduction is into intervals of 45 degrees. Two polynomial approximating functions are employed. Between 0 and pi/4 the cosine is approximated by 1 - x**2 P(x**2). Between pi/4 and pi/2 the sine is represented as x + x**3 P(x**2). ACCURACY: Relative error: arithmetic domain # trials peak rms DEC +-1000 3400 3.5e-17 9.1e-18 IEEE +-1000 30000 2.1e-16 5.7e-17 See also sin(). =item I<sinh>: Hyperbolic sine SYNOPSIS: # double x, y, sinh(); $y = sinh( $x ); DESCRIPTION: Returns hyperbolic sine of argument in the range MINLOG to MAXLOG. The range is partitioned into two segments. If |x| <= 1, a rational function of the form x + x**3 P(x)/Q(x) is employed. Otherwise the calculation is sinh(x) = ( exp(x) - exp(-x) )/2. ACCURACY: Relative error: arithmetic domain # trials peak rms DEC +- 88 50000 4.0e-17 7.7e-18 IEEE +-MAXLOG 30000 2.6e-16 5.7e-17 =item I<spence>: Dilogarithm SYNOPSIS: # double x, y, spence(); $y = spence( $x ); DESCRIPTION: Computes the integral x - | | log t spence(x) = - | ----- dt | | t - 1 - 1 for x >= 0. A rational approximation gives the integral in the interval (0.5, 1.5). Transformation formulas for 1/x and 1-x are employed outside the basic expansion range. ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE 0,4 30000 3.9e-15 5.4e-16 DEC 0,4 3000 2.5e-16 4.5e-17 =item I<sqrt>: Square root SYNOPSIS: # double x, y, sqrt(); $y = sqrt( $x ); DESCRIPTION: Returns the square root of x. Range reduction involves isolating the power of two of the argument and using a polynomial approximation to obtain a rough value for the square root. Then Heron's iteration is used three times to converge to an accurate value. ACCURACY: Relative error: arithmetic domain # trials peak rms DEC 0, 10 60000 2.1e-17 7.9e-18 IEEE 0,1.7e308 30000 1.7e-16 6.3e-17 ERROR MESSAGES: message condition value returned sqrt domain x < 0 0.0 =item I<stdtr>: Student's t distribution SYNOPSIS: # double t, stdtr(); short k; $y = stdtr( $k, $t ); DESCRIPTION: Computes the integral from minus infinity to t of the Student t distribution with integer k > 0 degrees of freedom: t - | | - | 2 -(k+1)/2 | ( (k+1)/2 ) | ( x ) ---------------------- | ( 1 + --- ) dx - | ( k ) sqrt( k pi ) | ( k/2 ) | | | - -inf. Relation to incomplete beta integral: 1 - stdtr(k,t) = 0.5 * incbet( k/2, 1/2, z ) where z = k/(k + t**2). For t < -2, this is the method of computation. For higher t, a direct method is derived from integration by parts. Since the function is symmetric about t=0, the area under the right tail of the density is found by calling the function with -t instead of t. ACCURACY: Tested at random 1 <= k <= 25. The "domain" refers to t. Relative error: arithmetic domain # trials peak rms IEEE -100,-2 50000 5.9e-15 1.4e-15 IEEE -2,100 500000 2.7e-15 4.9e-17 =item I<stdtri>: Functional inverse of Student's t distribution SYNOPSIS: # double p, t, stdtri(); # int k; $t = stdtri( $k, $p ); DESCRIPTION: Given probability p, finds the argument t such that stdtr(k,t) is equal to p. ACCURACY: Tested at random 1 <= k <= 100. The "domain" refers to p: Relative error: arithmetic domain # trials peak rms IEEE .001,.999 25000 5.7e-15 8.0e-16 IEEE 10^-6,.001 25000 2.0e-12 2.9e-14 =item I<struve>: Struve function SYNOPSIS: # double v, x, y, struve(); $y = struve( $v, $x ); DESCRIPTION: Computes the Struve function Hv(x) of order v, argument x. Negative x is rejected unless v is an integer. ACCURACY: Not accurately characterized, but spot checked against tables. =item I<plancki>: Integral of Planck's black body radiation formula SYNOPSIS: # double lambda, T, y, plancki() $y = plancki( $lambda, $T ); DESCRIPTION: Evaluates the definite integral, from wavelength 0 to lambda, of Planck's radiation formula -5 c1 lambda E = ------------------ c2/(lambda T) e - 1 Physical constants c1 = 3.7417749e-16 and c2 = 0.01438769 are built in to the function program. They are scaled to provide a result in watts per square meter. Argument T represents temperature in degrees Kelvin; lambda is wavelength in meters. The integral is expressed in closed form, in terms of polylogarithms (see polylog.c). The total area under the curve is (-1/8) (42 zeta(4) - 12 pi^2 zeta(2) + pi^4 ) c1 (T/c2)^4 = (pi^4 / 15) c1 (T/c2)^4 = 5.6705032e-8 T^4 where sigma = 5.6705032e-8 W m^2 K^-4 is the Stefan-Boltzmann constant. ACCURACY: The left tail of the function experiences some relative error amplification in computing the dominant term exp(-c2/(lambda T)). For the right-hand tail see planckc, below. Relative error. The domain refers to lambda T / c2. arithmetic domain # trials peak rms IEEE 0.1, 10 50000 7.1e-15 5.4e-16 =item I<polylog>: polylogarithm function SYNOPSIS: # double x, y, polylog(); # int n; $y = polylog( $n, $x ); The polylogarithm of order n is defined by the series inf k - x Li (x) = > --- . n - n k=1 k For x = 1, inf - 1 Li (1) = > --- = Riemann zeta function (n) . n - n k=1 k When n = 2, the function is the dilogarithm, related to Spence's integral: x 1-x - - | | -ln(1-t) | | ln t Li (x) = | -------- dt = | ------ dt = spence(1-x) . 2 | | t | | 1 - t - - 0 1 ACCURACY: Relative error: arithmetic domain n # trials peak rms IEEE 0, 1 2 50000 6.2e-16 8.0e-17 IEEE 0, 1 3 100000 2.5e-16 6.6e-17 IEEE 0, 1 4 30000 1.7e-16 4.9e-17 IEEE 0, 1 5 30000 5.1e-16 7.8e-17 =item I<bernum>: Bernoulli numbers SYNOPSIS: ($num, $den) = bernum( $n); ($num_array, $den_array) = bernum(); DESCRIPTION: This calculates the Bernoulli numbers, up to 30th order. If called with an integer argument, the numerator and denominator of that Bernoulli number is returned; if called with no argument, two array references representing the numerator and denominators of the first 30 Bernoulli numbers are returned. =item I<simpson>: Simpson's rule to find an integral SYNOPSIS: $result = simpson(\&fun, $a, $b, $abs_err, $rel_err, $nmax); sub fun { my $x = shift; return cos($x)*exp($x); } DESCRIPTION: This evaluates the area under the graph of a function, represented in a subroutine, from $a to $b, using an 8-point Newton-Cotes formula. The routine divides up the interval into equal segments, evaluates the integral, then compares that to the result with double the number of segments. If the two results agree, to within an absolute error $abs_err or a relative error $rel_err, the result is returned; otherwise, the number of segments is doubled again, and the results compared. This continues until the desired accuracy is attained, or until the maximum number of iterations $nmax is reached. =item I<vecang>: angle between two vectors SYNOPSIS: # double p[3], q[3], vecang(); $y = vecang( $p, $q ); DESCRIPTION: For two vectors p, q, the angle A between them is given by p.q / (|p| |q|) = cos A . where "." represents inner product, "|x|" the length of vector x. If the angle is small, an expression in sin A is preferred. Set r = q - p. Then p.q = p.p + p.r , |p|^2 = p.p , |q|^2 = p.p + 2 p.r + r.r , p.p^2 + 2 p.p p.r + p.r^2 cos^2 A = ---------------------------- p.p (p.p + 2 p.r + r.r) p.p + 2 p.r + p.r^2 / p.p = --------------------------- , p.p + 2 p.r + r.r sin^2 A = 1 - cos^2 A r.r - p.r^2 / p.p = -------------------- p.p + 2 p.r + r.r = (r.r - p.r^2 / p.p) / q.q . ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE -1, 1 10^6 1.7e-16 4.2e-17 =item I<onef2>: Hypergeometric function 1F2 SYNOPSIS: # double a, b, c, x, value; # double *err; ($value, $err) = onef2( $a, $b, $c, $x) ACCURACY: Not accurately characterized, but spot checked against tables. =item I<threef0>: Hypergeometric function 3F0 SYNOPSIS: # double a, b, c, x, value; # double *err; ($value, $err) = threef0( $a, $b, $c, $x ) ACCURACY: Not accurately characterized, but spot checked against tables. =item I<yv>: Bessel function Yv with noninteger v SYNOPSIS: # double v, x; # double yv( v, x ); $y = yv( $v, $x ); ACCURACY: Not accurately characterized, but spot checked against tables. =item I<tan>: Circular tangent SYNOPSIS: # double x, y, tan(); $y = tan( $x ); DESCRIPTION: Returns the circular tangent of the radian argument x. Range reduction is modulo pi/4. A rational function x + x**3 P(x**2)/Q(x**2) is employed in the basic interval [0, pi/4]. ACCURACY: Relative error: arithmetic domain # trials peak rms DEC +-1.07e9 44000 4.1e-17 1.0e-17 IEEE +-1.07e9 30000 2.9e-16 8.1e-17 ERROR MESSAGES: message condition value returned tan total loss x > 1.073741824e9 0.0 =item I<cot>: Circular cotangent SYNOPSIS: # double x, y, cot(); $y = cot( $x ); DESCRIPTION: Returns the circular cotangent of the radian argument x. Range reduction is modulo pi/4. A rational function x + x**3 P(x**2)/Q(x**2) is employed in the basic interval [0, pi/4]. ACCURACY: Relative error: arithmetic domain # trials peak rms IEEE +-1.07e9 30000 2.9e-16 8.2e-17 ERROR MESSAGES: message condition value returned cot total loss x > 1.073741824e9 0.0 cot singularity x = 0 INFINITY =item I<tandg>: Circular tangent of argument in degrees SYNOPSIS: # double x, y, tandg(); $y = tandg( $x ); DESCRIPTION: Returns the circular tangent of the argument x in degrees. Range reduction is modulo pi/4. A rational function x + x**3 P(x**2)/Q(x**2) is employed in the basic interval [0, pi/4]. ACCURACY: Relative error: arithmetic domain # trials peak rms DEC 0,10 8000 3.4e-17 1.2e-17 IEEE 0,10 30000 3.2e-16 8.4e-17 ERROR MESSAGES: message condition value returned tandg total loss x > 8.0e14 (DEC) 0.0 x > 1.0e14 (IEEE) tandg singularity x = 180 k + 90 MAXNUM =item I<cotdg>: Circular cotangent of argument in degrees SYNOPSIS: # double x, y, cotdg(); $y = cotdg( $x ); DESCRIPTION: Returns the circular cotangent of the argument x in degrees. Range reduction is modulo pi/4. A rational function x + x**3 P(x**2)/Q(x**2) is employed in the basic interval [0, pi/4]. ERROR MESSAGES: message condition value returned cotdg total loss x > 8.0e14 (DEC) 0.0 x > 1.0e14 (IEEE) cotdg singularity x = 180 k MAXNUM =item I<tanh>: Hyperbolic tangent SYNOPSIS: # double x, y, tanh(); $y = tanh( $x ); DESCRIPTION: Returns hyperbolic tangent of argument in the range MINLOG to MAXLOG. A rational function is used for |x| < 0.625. The form x + x**3 P(x)/Q(x) of Cody _& Waite is employed. Otherwise, tanh(x) = sinh(x)/cosh(x) = 1 - 2/(exp(2x) + 1). ACCURACY: Relative error: arithmetic domain # trials peak rms DEC -2,2 50000 3.3e-17 6.4e-18 IEEE -2,2 30000 2.5e-16 5.8e-17 =item I<unity>: Relative error approximations for function arguments near unity. SYNOPSIS: # log1p(x) = log(1+x) $y = log1p( $x ); # expm1(x) = exp(x) - 1 $y = expm1( $x ); # cosm1(x) = cos(x) - 1 $y = cosm1( $x ); =item I<yn>: Bessel function of second kind of integer order SYNOPSIS: # double x, y, yn(); # int n; $y = yn( $n, $x ); DESCRIPTION: Returns Bessel function of order n, where n is a (possibly negative) integer. The function is evaluated by forward recurrence on n, starting with values computed by the routines y0() and y1(). If n = 0 or 1 the routine for y0 or y1 is called directly. ACCURACY: Absolute error, except relative when y > 1: arithmetic domain # trials peak rms DEC 0, 30 2200 2.9e-16 5.3e-17 IEEE 0, 30 30000 3.4e-15 4.3e-16 ERROR MESSAGES: message condition value returned yn singularity x = 0 MAXNUM yn overflow MAXNUM Spot checked against tables for x, n between 0 and 100. =item I<zeta>: Riemann zeta function of two arguments SYNOPSIS: # double x, q, y, zeta(); $y = zeta( $x, $q ); DESCRIPTION: inf. - -x zeta(x,q) = > (k+q) - k=0 where x > 1 and q is not a negative integer or zero. The Euler-Maclaurin summation formula is used to obtain the expansion n - -x zeta(x,q) = > (k+q) - k=1 1-x inf. B x(x+1)...(x+2j) (n+q) 1 - 2j + --------- - ------- + > -------------------- x-1 x - x+2j+1 2(n+q) j=1 (2j)! (n+q) where the B2j are Bernoulli numbers. Note that (see zetac.c) zeta(x,1) = zetac(x) + 1. ACCURACY: REFERENCE: Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals, Series, and Products, p. 1073; Academic Press, 1980. =item I<zetac>: Riemann zeta function SYNOPSIS: # double x, y, zetac(); $y = zetac( $x ); DESCRIPTION: inf. - -x zetac(x) = > k , x > 1, - k=2 is related to the Riemann zeta function by Riemann zeta(x) = zetac(x) + 1. Extension of the function definition for x < 1 is implemented. Zero is returned for x > log2(MAXNUM). An overflow error may occur for large negative x, due to the gamma function in the reflection formula. ACCURACY: Tabulated values have full machine accuracy. Relative error: arithmetic domain # trials peak rms IEEE 1,50 10000 9.8e-16 1.3e-16 DEC 1,50 2000 1.1e-16 1.9e-17 =back =head1 TODO =over 4 =item * Make the configuration of mconf.h automatic. =item * Include the rest of the routines in the cephes library, such as polynomial and matrix manipulation functions; this will involve writing typemaps for arrays. =back =head1 BUGS Please report any to Randy Kobes <randy@theoryx5.uwinnipeg.ca> =head1 SEE ALSO For interfaces to programs which can do symbolic manipulation, see L<PDL>, L<Math::Pari>, and L<Math::ematica>. For a command line interface to the routines of I<Math::Cephes>, see the included C<pmath> script. For a different interface to the fraction and complex number routines, see L<Math::Cephes::Fraction> and L<Math::Cephes::Complex>. For an interface to some polynomial routines, see L<Math::Cephes::Polynomial>, and for some matrix routines, see L<Math::Cephes::Matrix>. =head1 COPYRIGHT The C code for the Cephes Math Library is Copyright 1984, 1987, 1989, 2002 by Stephen L. Moshier, and is available at http://www.netlib.org/cephes/. Direct inquiries to 30 Frost Street, Cambridge, MA 02140. The file arrays.c included here to handle passing arrays into and out of C routines comes from the PGPLOT module of Karl Glazebrook <kgb@zzoepp.aao.gov.au>. The perl interface is copyright 2000, 2002 by Randy Kobes. This library is free software; you can redistribute it and/or modify it under the same terms as Perl itself. =cut