diff options
Diffstat (limited to 'libm/ldouble/README.txt')
-rw-r--r-- | libm/ldouble/README.txt | 3502 |
1 files changed, 0 insertions, 3502 deletions
diff --git a/libm/ldouble/README.txt b/libm/ldouble/README.txt deleted file mode 100644 index 30fcaad36..000000000 --- a/libm/ldouble/README.txt +++ /dev/null @@ -1,3502 +0,0 @@ -/* acoshl.c - * - * Inverse hyperbolic cosine, long double precision - * - * - * - * SYNOPSIS: - * - * long double x, y, acoshl(); - * - * y = acoshl( x ); - * - * - * - * DESCRIPTION: - * - * Returns inverse hyperbolic cosine of argument. - * - * If 1 <= x < 1.5, a rational approximation - * - * sqrt(2z) * 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 - * IEEE 1,3 30000 2.0e-19 3.9e-20 - * - * - * ERROR MESSAGES: - * - * message condition value returned - * acoshl domain |x| < 1 0.0 - * - */ - -/* asinhl.c - * - * Inverse hyperbolic sine, long double precision - * - * - * - * SYNOPSIS: - * - * long double x, y, asinhl(); - * - * y = asinhl( x ); - * - * - * - * DESCRIPTION: - * - * Returns inverse hyperbolic sine of argument. - * - * If |x| < 0.5, the function is approximated by a rational - * form x + x**3 P(x)/Q(x). Otherwise, - * - * asinh(x) = log( x + sqrt(1 + x*x) ). - * - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * IEEE -3,3 30000 1.7e-19 3.5e-20 - * - */ - -/* asinl.c - * - * Inverse circular sine, long double precision - * - * - * - * SYNOPSIS: - * - * double x, y, asinl(); - * - * y = asinl( 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 - * IEEE -1, 1 30000 2.7e-19 4.8e-20 - * - * - * ERROR MESSAGES: - * - * message condition value returned - * asin domain |x| > 1 0.0 - * - */ -/* acosl() - * - * Inverse circular cosine, long double precision - * - * - * - * SYNOPSIS: - * - * double x, y, acosl(); - * - * y = acosl( x ); - * - * - * - * DESCRIPTION: - * - * Returns radian angle between -pi/2 and +pi/2 whose cosine - * is x. - * - * Analytically, acos(x) = pi/2 - asin(x). However if |x| is - * near 1, there is cancellation error in subtracting asin(x) - * from pi/2. Hence if x < -0.5, - * - * acos(x) = pi - 2.0 * asin( sqrt((1+x)/2) ); - * - * or if x > +0.5, - * - * acos(x) = 2.0 * asin( sqrt((1-x)/2) ). - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * IEEE -1, 1 30000 1.4e-19 3.5e-20 - * - * - * ERROR MESSAGES: - * - * message condition value returned - * asin domain |x| > 1 0.0 - */ - -/* atanhl.c - * - * Inverse hyperbolic tangent, long double precision - * - * - * - * SYNOPSIS: - * - * long double x, y, atanhl(); - * - * y = atanhl( x ); - * - * - * - * DESCRIPTION: - * - * Returns inverse hyperbolic tangent of argument in the range - * MINLOGL to MAXLOGL. - * - * 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 - * IEEE -1,1 30000 1.1e-19 3.3e-20 - * - */ - -/* atanl.c - * - * Inverse circular tangent, long double precision - * (arctangent) - * - * - * - * SYNOPSIS: - * - * long double x, y, atanl(); - * - * y = atanl( x ); - * - * - * - * DESCRIPTION: - * - * Returns radian angle between -pi/2 and +pi/2 whose tangent - * is x. - * - * Range reduction is from four intervals into the interval - * from zero to tan( pi/8 ). The approximant uses a rational - * function of degree 3/4 of the form x + x**3 P(x)/Q(x). - * - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * IEEE -10, 10 150000 1.3e-19 3.0e-20 - * - */ -/* atan2l() - * - * Quadrant correct inverse circular tangent, - * long double precision - * - * - * - * SYNOPSIS: - * - * long double x, y, z, atan2l(); - * - * z = atan2l( 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 60000 1.7e-19 3.2e-20 - * See atan.c. - * - */ - -/* bdtrl.c - * - * Binomial distribution - * - * - * - * SYNOPSIS: - * - * int k, n; - * long double p, y, bdtrl(); - * - * y = bdtrl( 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 (k,n,p) with a and b between 0 - * and 10000 and p between 0 and 1. - * Relative error: - * arithmetic domain # trials peak rms - * IEEE 0,10000 3000 1.6e-14 2.2e-15 - * - * ERROR MESSAGES: - * - * message condition value returned - * bdtrl domain k < 0 0.0 - * n < k - * x < 0, x > 1 - * - */ -/* bdtrcl() - * - * Complemented binomial distribution - * - * - * - * SYNOPSIS: - * - * int k, n; - * long double p, y, bdtrcl(); - * - * y = bdtrcl( 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: - * - * See incbet.c. - * - * ERROR MESSAGES: - * - * message condition value returned - * bdtrcl domain x<0, x>1, n<k 0.0 - */ -/* bdtril() - * - * Inverse binomial distribution - * - * - * - * SYNOPSIS: - * - * int k, n; - * long double p, y, bdtril(); - * - * p = bdtril( 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: - * - * See incbi.c. - * Tested at random k, n between 1 and 10000. The "domain" refers to p: - * Relative error: - * arithmetic domain # trials peak rms - * IEEE 0,1 3500 2.0e-15 8.2e-17 - * - * ERROR MESSAGES: - * - * message condition value returned - * bdtril domain k < 0, n <= k 0.0 - * x < 0, x > 1 - */ - - -/* btdtrl.c - * - * Beta distribution - * - * - * - * SYNOPSIS: - * - * long double a, b, x, y, btdtrl(); - * - * y = btdtrl( 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 - * - * - * The mean value of this distribution is a/(a+b). The variance - * is ab/[(a+b)^2 (a+b+1)]. - * - * This function is identical to the incomplete beta integral - * function, incbetl(a, b, x). - * - * The complemented function is - * - * 1 - P(1-x) = incbetl( b, a, x ); - * - * - * ACCURACY: - * - * See incbetl.c. - * - */ - -/* cbrtl.c - * - * Cube root, long double precision - * - * - * - * SYNOPSIS: - * - * long double x, y, cbrtl(); - * - * y = cbrtl( 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 - * IEEE .125,8 80000 7.0e-20 2.2e-20 - * IEEE exp(+-707) 100000 7.0e-20 2.4e-20 - * - */ - -/* chdtrl.c - * - * Chi-square distribution - * - * - * - * SYNOPSIS: - * - * long double df, x, y, chdtrl(); - * - * y = chdtrl( df, x ); - * - * - * - * DESCRIPTION: - * - * Returns the area under the left hand tail (from 0 to x) - * of the Chi square probability density function with - * v degrees of freedom. - * - * - * inf. - * - - * 1 | | v/2-1 -t/2 - * P( x | v ) = ----------- | t e dt - * v/2 - | | - * 2 | (v/2) - - * x - * - * where x is the Chi-square variable. - * - * The incomplete gamma integral is used, according to the - * formula - * - * y = chdtr( v, x ) = igam( v/2.0, x/2.0 ). - * - * - * The arguments must both be positive. - * - * - * - * ACCURACY: - * - * See igam(). - * - * ERROR MESSAGES: - * - * message condition value returned - * chdtr domain x < 0 or v < 1 0.0 - */ -/* chdtrcl() - * - * Complemented Chi-square distribution - * - * - * - * SYNOPSIS: - * - * long double v, x, y, chdtrcl(); - * - * y = chdtrcl( v, x ); - * - * - * - * DESCRIPTION: - * - * Returns the area under the right hand tail (from x to - * infinity) of the Chi square probability density function - * with v degrees of freedom: - * - * - * inf. - * - - * 1 | | v/2-1 -t/2 - * P( x | v ) = ----------- | t e dt - * v/2 - | | - * 2 | (v/2) - - * x - * - * where x is the Chi-square variable. - * - * The incomplete gamma integral is used, according to the - * formula - * - * y = chdtr( v, x ) = igamc( v/2.0, x/2.0 ). - * - * - * The arguments must both be positive. - * - * - * - * ACCURACY: - * - * See igamc(). - * - * ERROR MESSAGES: - * - * message condition value returned - * chdtrc domain x < 0 or v < 1 0.0 - */ -/* chdtril() - * - * Inverse of complemented Chi-square distribution - * - * - * - * SYNOPSIS: - * - * long double df, x, y, chdtril(); - * - * x = chdtril( 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 - * - */ - -/* clogl.c - * - * Complex natural logarithm - * - * - * - * SYNOPSIS: - * - * void clogl(); - * cmplxl z, w; - * - * clogl( &z, &w ); - * - * - * - * DESCRIPTION: - * - * Returns complex logarithm to the base e (2.718...) of - * the complex argument x. - * - * If z = x + iy, r = sqrt( x**2 + y**2 ), - * then - * w = log(r) + i arctan(y/x). - * - * The arctangent ranges from -PI to +PI. - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * 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. - */ - -/* cexpl() - * - * Complex exponential function - * - * - * - * SYNOPSIS: - * - * void cexpl(); - * cmplxl z, w; - * - * cexpl( &z, &w ); - * - * - * - * DESCRIPTION: - * - * Returns the exponential of the complex argument z - * into the complex result w. - * - * If - * z = x + iy, - * r = exp(x), - * - * then - * - * w = r cos y + i r sin y. - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * DEC -10,+10 8700 3.7e-17 1.1e-17 - * IEEE -10,+10 30000 3.0e-16 8.7e-17 - * - */ -/* csinl() - * - * Complex circular sine - * - * - * - * SYNOPSIS: - * - * void csinl(); - * cmplxl z, w; - * - * csinl( &z, &w ); - * - * - * - * DESCRIPTION: - * - * If - * z = x + iy, - * - * then - * - * w = sin x cosh y + i cos x sinh y. - * - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * 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. - * - */ -/* ccosl() - * - * Complex circular cosine - * - * - * - * SYNOPSIS: - * - * void ccosl(); - * cmplxl z, w; - * - * ccosl( &z, &w ); - * - * - * - * DESCRIPTION: - * - * If - * z = x + iy, - * - * then - * - * w = cos x cosh y - i sin x sinh y. - * - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * DEC -10,+10 8400 4.5e-17 1.3e-17 - * IEEE -10,+10 30000 3.8e-16 1.0e-16 - */ -/* ctanl() - * - * Complex circular tangent - * - * - * - * SYNOPSIS: - * - * void ctanl(); - * cmplxl z, w; - * - * ctanl( &z, &w ); - * - * - * - * DESCRIPTION: - * - * If - * z = x + iy, - * - * then - * - * sin 2x + i sinh 2y - * w = --------------------. - * cos 2x + cosh 2y - * - * On the real axis the denominator is zero at odd multiples - * of PI/2. The denominator is evaluated by its Taylor - * series near these points. - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * 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. - */ -/* ccotl() - * - * Complex circular cotangent - * - * - * - * SYNOPSIS: - * - * void ccotl(); - * cmplxl z, w; - * - * ccotl( &z, &w ); - * - * - * - * DESCRIPTION: - * - * If - * z = x + iy, - * - * then - * - * sin 2x - i sinh 2y - * w = --------------------. - * cosh 2y - cos 2x - * - * On the real axis, the denominator has zeros at even - * multiples of PI/2. Near these points it is evaluated - * by a Taylor series. - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * 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. - */ - -/* casinl() - * - * Complex circular arc sine - * - * - * - * SYNOPSIS: - * - * void casinl(); - * cmplxl z, w; - * - * casinl( &z, &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. - */ -/* cacosl() - * - * Complex circular arc cosine - * - * - * - * SYNOPSIS: - * - * void cacosl(); - * cmplxl z, w; - * - * cacosl( &z, &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 - */ - -/* catanl() - * - * Complex circular arc tangent - * - * - * - * SYNOPSIS: - * - * void catanl(); - * cmplxl z, w; - * - * catanl( &z, &w ); - * - * - * - * DESCRIPTION: - * - * If - * z = x + iy, - * - * then - * 1 ( 2x ) - * Re w = - arctan(-----------) + k PI - * 2 ( 2 2) - * (1 - x - y ) - * - * ( 2 2) - * 1 (x + (y+1) ) - * Im w = - log(------------) - * 4 ( 2 2) - * (x + (y-1) ) - * - * Where k is an arbitrary integer. - * - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * 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(). - */ - -/* cmplxl.c - * - * Complex number arithmetic - * - * - * - * SYNOPSIS: - * - * typedef struct { - * long double r; real part - * long double i; imaginary part - * }cmplxl; - * - * cmplxl *a, *b, *c; - * - * caddl( a, b, c ); c = b + a - * csubl( a, b, c ); c = b - a - * cmull( a, b, c ); c = b * a - * cdivl( a, b, c ); c = b / a - * cnegl( c ); c = -c - * cmovl( b, c ); c = b - * - * - * - * DESCRIPTION: - * - * Addition: - * c.r = b.r + a.r - * c.i = b.i + a.i - * - * Subtraction: - * c.r = b.r - a.r - * c.i = b.i - a.i - * - * Multiplication: - * c.r = b.r * a.r - b.i * a.i - * c.i = b.r * a.i + b.i * a.r - * - * Division: - * d = a.r * a.r + a.i * a.i - * c.r = (b.r * a.r + b.i * a.i)/d - * c.i = (b.i * a.r - b.r * a.i)/d - * ACCURACY: - * - * In DEC arithmetic, the test (1/z) * z = 1 had peak relative - * error 3.1e-17, rms 1.2e-17. The test (y/z) * (z/y) = 1 had - * peak relative error 8.3e-17, rms 2.1e-17. - * - * Tests in the rectangle {-10,+10}: - * Relative error: - * arithmetic function # trials peak rms - * 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 - */ - -/* cabsl() - * - * Complex absolute value - * - * - * - * SYNOPSIS: - * - * long double cabsl(); - * cmplxl z; - * long double a; - * - * 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 - */ -/* csqrtl() - * - * Complex square root - * - * - * - * SYNOPSIS: - * - * void csqrtl(); - * cmplxl z, w; - * - * csqrtl( &z, &w ); - * - * - * - * DESCRIPTION: - * - * - * If z = x + iy, r = |z|, then - * - * 1/2 - * Im w = [ (r - x)/2 ] , - * - * Re w = y / 2 Im w. - * - * - * Note that -w is also a square root of z. The 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. - */ - -/* coshl.c - * - * Hyperbolic cosine, long double precision - * - * - * - * SYNOPSIS: - * - * long double x, y, coshl(); - * - * y = coshl( x ); - * - * - * - * DESCRIPTION: - * - * Returns hyperbolic cosine of argument in the range MINLOGL to - * MAXLOGL. - * - * cosh(x) = ( exp(x) + exp(-x) )/2. - * - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * IEEE +-10000 30000 1.1e-19 2.8e-20 - * - * - * ERROR MESSAGES: - * - * message condition value returned - * cosh overflow |x| > MAXLOGL MAXNUML - * - * - */ - -/* elliel.c - * - * Incomplete elliptic integral of the second kind - * - * - * - * SYNOPSIS: - * - * long double phi, m, y, elliel(); - * - * y = elliel( 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 - * IEEE -10,10 50000 2.7e-18 2.3e-19 - * - * - */ - -/* ellikl.c - * - * Incomplete elliptic integral of the first kind - * - * - * - * SYNOPSIS: - * - * long double phi, m, y, ellikl(); - * - * y = ellikl( 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 30000 3.6e-18 4.1e-19 - * - * - */ - -/* ellpel.c - * - * Complete elliptic integral of the second kind - * - * - * - * SYNOPSIS: - * - * long double m1, y, ellpel(); - * - * y = ellpel( m1 ); - * - * - * - * DESCRIPTION: - * - * Approximates the integral - * - * - * pi/2 - * - - * | | 2 - * E(m) = | sqrt( 1 - m sin t ) dt - * | | - * - - * 0 - * - * Where m = 1 - m1, using the approximation - * - * P(x) - x log x Q(x). - * - * Though there are no singularities, the argument m1 is used - * rather than m for compatibility with ellpk(). - * - * E(1) = 1; E(0) = pi/2. - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * IEEE 0, 1 10000 1.1e-19 3.5e-20 - * - * - * ERROR MESSAGES: - * - * message condition value returned - * ellpel domain x<0, x>1 0.0 - * - */ - -/* ellpjl.c - * - * Jacobian Elliptic Functions - * - * - * - * SYNOPSIS: - * - * long double u, m, sn, cn, dn, phi; - * int ellpjl(); - * - * ellpjl( u, m, _&sn, _&cn, _&dn, _&phi ); - * - * - * - * DESCRIPTION: - * - * - * Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m), - * and dn(u|m) of parameter m between 0 and 1, and real - * argument u. - * - * These functions are periodic, with quarter-period on the - * real axis equal to the complete elliptic integral - * ellpk(1.0-m). - * - * Relation to incomplete elliptic integral: - * If u = ellik(phi,m), then sn(u|m) = sin(phi), - * and cn(u|m) = cos(phi). Phi is called the amplitude of u. - * - * Computation is by means of the arithmetic-geometric mean - * algorithm, except when m is within 1e-12 of 0 or 1. In the - * latter case with m close to 1, the approximation applies - * only for phi < pi/2. - * - * ACCURACY: - * - * Tested at random points with u between 0 and 10, m between - * 0 and 1. - * - * Absolute error (* = relative error): - * arithmetic function # trials peak rms - * IEEE sn 10000 1.7e-18 2.3e-19 - * IEEE cn 20000 1.6e-18 2.2e-19 - * IEEE dn 10000 4.7e-15 2.7e-17 - * IEEE phi 10000 4.0e-19* 6.6e-20* - * - * Accuracy deteriorates when u is large. - * - */ - -/* ellpkl.c - * - * Complete elliptic integral of the first kind - * - * - * - * SYNOPSIS: - * - * long double m1, y, ellpkl(); - * - * y = ellpkl( m1 ); - * - * - * - * DESCRIPTION: - * - * Approximates the integral - * - * - * - * pi/2 - * - - * | | - * | dt - * K(m) = | ------------------ - * | 2 - * | | sqrt( 1 - m sin t ) - * - - * 0 - * - * where m = 1 - m1, using the approximation - * - * P(x) - log x Q(x). - * - * The argument m1 is used rather than m so that the logarithmic - * singularity at m = 1 will be shifted to the origin; this - * preserves maximum accuracy. - * - * K(0) = pi/2. - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * IEEE 0,1 10000 1.1e-19 3.3e-20 - * - * ERROR MESSAGES: - * - * message condition value returned - * ellpkl domain x<0, x>1 0.0 - * - */ - -/* exp10l.c - * - * Base 10 exponential function, long double precision - * (Common antilogarithm) - * - * - * - * SYNOPSIS: - * - * long double x, y, exp10l() - * - * y = exp10l( 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 +-4900 30000 1.0e-19 2.7e-20 - * - * ERROR MESSAGES: - * - * message condition value returned |