From 1077fa4d772832f77a677ce7fb7c2d513b959e3f Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 10 May 2001 00:40:28 +0000 Subject: uClibc now has a math library. muahahahaha! -Erik --- libm/double/hyperg.c | 386 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 386 insertions(+) create mode 100644 libm/double/hyperg.c (limited to 'libm/double/hyperg.c') diff --git a/libm/double/hyperg.c b/libm/double/hyperg.c new file mode 100644 index 000000000..36a3f9781 --- /dev/null +++ b/libm/double/hyperg.c @@ -0,0 +1,386 @@ +/* hyperg.c + * + * 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 + qtst1: + 21800 max = 1.4200E-14 rms = 1.0841E-15 ave = -5.3640E-17 + ltstd: + 25500 max = 1.2759e-14 rms = 3.7155e-16 ave = 1.5384e-18 + * 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. + * + */ + +/* hyperg.c */ + + +/* +Cephes Math Library Release 2.8: June, 2000 +Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier +*/ + +#include + +#ifdef ANSIPROT +extern double exp ( double ); +extern double log ( double ); +extern double gamma ( double ); +extern double lgam ( double ); +extern double fabs ( double ); +double hyp2f0 ( double, double, double, int, double * ); +static double hy1f1p(double, double, double, double *); +static double hy1f1a(double, double, double, double *); +double hyperg (double, double, double); +#else +double exp(), log(), gamma(), lgam(), fabs(), hyp2f0(); +static double hy1f1p(); +static double hy1f1a(); +double hyperg(); +#endif +extern double MAXNUM, MACHEP; + +double hyperg( a, b, x) +double a, b, x; +{ +double asum, psum, acanc, pcanc, temp; + +/* See if a Kummer transformation will help */ +temp = b - a; +if( fabs(temp) < 0.001 * fabs(a) ) + return( exp(x) * hyperg( temp, b, -x ) ); + + +psum = hy1f1p( a, b, x, &pcanc ); +if( pcanc < 1.0e-15 ) + goto done; + + +/* try asymptotic series */ + +asum = hy1f1a( a, b, x, &acanc ); + + +/* Pick the result with less estimated error */ + +if( acanc < pcanc ) + { + pcanc = acanc; + psum = asum; + } + +done: +if( pcanc > 1.0e-12 ) + mtherr( "hyperg", PLOSS ); + +return( psum ); +} + + + + +/* Power series summation for confluent hypergeometric function */ + + +static double hy1f1p( a, b, x, err ) +double a, b, x; +double *err; +{ +double n, a0, sum, t, u, temp; +double an, bn, maxt, pcanc; + + +/* set up for power series summation */ +an = a; +bn = b; +a0 = 1.0; +sum = 1.0; +n = 1.0; +t = 1.0; +maxt = 0.0; + + +while( t > MACHEP ) + { + if( bn == 0 ) /* check bn first since if both */ + { + mtherr( "hyperg", SING ); + return( MAXNUM ); /* an and bn are zero it is */ + } + if( an == 0 ) /* a singularity */ + return( sum ); + if( n > 200 ) + goto pdone; + u = x * ( an / (bn * n) ); + + /* check for blowup */ + temp = fabs(u); + if( (temp > 1.0 ) && (maxt > (MAXNUM/temp)) ) + { + pcanc = 1.0; /* estimate 100% error */ + goto blowup; + } + + a0 *= u; + sum += a0; + t = fabs(a0); + if( t > maxt ) + maxt = t; +/* + if( (maxt/fabs(sum)) > 1.0e17 ) + { + pcanc = 1.0; + goto blowup; + } +*/ + an += 1.0; + bn += 1.0; + n += 1.0; + } + +pdone: + +/* estimate error due to roundoff and cancellation */ +if( sum != 0.0 ) + maxt /= fabs(sum); +maxt *= MACHEP; /* this way avoids multiply overflow */ +pcanc = fabs( MACHEP * n + maxt ); + +blowup: + +*err = pcanc; + +return( sum ); +} + + +/* hy1f1a() */ +/* asymptotic formula for hypergeometric function: + * + * ( -a + * -- ( |z| + * | (b) ( -------- 2f0( a, 1+a-b, -1/x ) + * ( -- + * ( | (b-a) + * + * + * x a-b ) + * e |x| ) + * + -------- 2f0( b-a, 1-a, 1/x ) ) + * -- ) + * | (a) ) + */ + +static double hy1f1a( a, b, x, err ) +double a, b, x; +double *err; +{ +double h1, h2, t, u, temp, acanc, asum, err1, err2; + +if( x == 0 ) + { + acanc = 1.0; + asum = MAXNUM; + goto adone; + } +temp = log( fabs(x) ); +t = x + temp * (a-b); +u = -temp * a; + +if( b > 0 ) + { + temp = lgam(b); + t += temp; + u += temp; + } + +h1 = hyp2f0( a, a-b+1, -1.0/x, 1, &err1 ); + +temp = exp(u) / gamma(b-a); +h1 *= temp; +err1 *= temp; + +h2 = hyp2f0( b-a, 1.0-a, 1.0/x, 2, &err2 ); + +if( a < 0 ) + temp = exp(t) / gamma(a); +else + temp = exp( t - lgam(a) ); + +h2 *= temp; +err2 *= temp; + +if( x < 0.0 ) + asum = h1; +else + asum = h2; + +acanc = fabs(err1) + fabs(err2); + + +if( b < 0 ) + { + temp = gamma(b); + asum *= temp; + acanc *= fabs(temp); + } + + +if( asum != 0.0 ) + acanc /= fabs(asum); + +acanc *= 30.0; /* fudge factor, since error of asymptotic formula + * often seems this much larger than advertised */ + +adone: + + +*err = acanc; +return( asum ); +} + +/* hyp2f0() */ + +double hyp2f0( a, b, x, type, err ) +double a, b, x; +int type; /* determines what converging factor to use */ +double *err; +{ +double a0, alast, t, tlast, maxt; +double n, an, bn, u, sum, temp; + +an = a; +bn = b; +a0 = 1.0e0; +alast = 1.0e0; +sum = 0.0; +n = 1.0e0; +t = 1.0e0; +tlast = 1.0e9; +maxt = 0.0; + +do + { + if( an == 0 ) + goto pdone; + if( bn == 0 ) + goto pdone; + + u = an * (bn * x / n); + + /* check for blowup */ + temp = fabs(u); + if( (temp > 1.0 ) && (maxt > (MAXNUM/temp)) ) + goto error; + + a0 *= u; + t = fabs(a0); + + /* terminating condition for asymptotic series */ + if( t > tlast ) + goto ndone; + + tlast = t; + sum += alast; /* the sum is one term behind */ + alast = a0; + + if( n > 200 ) + goto ndone; + + an += 1.0e0; + bn += 1.0e0; + n += 1.0e0; + if( t > maxt ) + maxt = t; + } +while( t > MACHEP ); + + +pdone: /* series converged! */ + +/* estimate error due to roundoff and cancellation */ +*err = fabs( MACHEP * (n + maxt) ); + +alast = a0; +goto done; + +ndone: /* series did not converge */ + +/* The following "Converging factors" are supposed to improve accuracy, + * but do not actually seem to accomplish very much. */ + +n -= 1.0; +x = 1.0/x; + +switch( type ) /* "type" given as subroutine argument */ +{ +case 1: + alast *= ( 0.5 + (0.125 + 0.25*b - 0.5*a + 0.25*x - 0.25*n)/x ); + break; + +case 2: + alast *= 2.0/3.0 - b + 2.0*a + x - n; + break; + +default: + ; +} + +/* estimate error due to roundoff, cancellation, and nonconvergence */ +*err = MACHEP * (n + maxt) + fabs ( a0 ); + + +done: +sum += alast; +return( sum ); + +/* series blew up: */ +error: +*err = MAXNUM; +mtherr( "hyperg", TLOSS ); +return( sum ); +} -- cgit v1.2.3