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/float/hypergf.c | 384 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 384 insertions(+) create mode 100644 libm/float/hypergf.c (limited to 'libm/float/hypergf.c') diff --git a/libm/float/hypergf.c b/libm/float/hypergf.c new file mode 100644 index 000000000..60d0eb4c5 --- /dev/null +++ b/libm/float/hypergf.c @@ -0,0 +1,384 @@ +/* hypergf.c + * + * Confluent hypergeometric function + * + * + * + * SYNOPSIS: + * + * float a, b, x, y, hypergf(); + * + * y = hypergf( a, b, x ); + * + * + * + * DESCRIPTION: + * + * Computes the confluent hypergeometric function + * + * 1 2 + * a x a(a+1) x + * F ( a,b;x ) = 1 + ---- + --------- + ... + * 1 1 b 1! b(b+1) 2! + * + * Many higher transcendental functions are special cases of + * this power series. + * + * As is evident from the formula, b must not be a negative + * integer or zero unless a is an integer with 0 >= a > b. + * + * The routine attempts both a direct summation of the series + * and an asymptotic expansion. In each case error due to + * roundoff, cancellation, and nonconvergence is estimated. + * The result with smaller estimated error is returned. + * + * + * + * ACCURACY: + * + * Tested at random points (a, b, x), all three variables + * ranging from 0 to 30. + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0,5 10000 6.6e-7 1.3e-7 + * IEEE 0,30 30000 1.1e-5 6.5e-7 + * + * Larger errors can be observed when b is near a negative + * integer or zero. Certain combinations of arguments yield + * serious cancellation error in the power series summation + * and also are not in the region of near convergence of the + * asymptotic series. An error message is printed if the + * self-estimated relative error is greater than 1.0e-3. + * + */ + +/* hyperg.c */ + + +/* +Cephes Math Library Release 2.1: November, 1988 +Copyright 1984, 1987, 1988 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +#include + +extern float MAXNUMF, MACHEPF; + +#define fabsf(x) ( (x) < 0 ? -(x) : (x) ) + +#ifdef ANSIC +float expf(float); +float hyp2f0f(float, float, float, int, float *); +static float hy1f1af(float, float, float, float *); +static float hy1f1pf(float, float, float, float *); +float logf(float), gammaf(float), lgamf(float); +#else +float expf(), hyp2f0f(); +float logf(), gammaf(), lgamf(); +static float hy1f1pf(), hy1f1af(); +#endif + +float hypergf( float aa, float bb, float xx ) +{ +float a, b, x, asum, psum, acanc, pcanc, temp; + + +a = aa; +b = bb; +x = xx; +/* See if a Kummer transformation will help */ +temp = b - a; +if( fabsf(temp) < 0.001 * fabsf(a) ) + return( expf(x) * hypergf( temp, b, -x ) ); + +psum = hy1f1pf( a, b, x, &pcanc ); +if( pcanc < 1.0e-6 ) + goto done; + + +/* try asymptotic series */ + +asum = hy1f1af( a, b, x, &acanc ); + + +/* Pick the result with less estimated error */ + +if( acanc < pcanc ) + { + pcanc = acanc; + psum = asum; + } + +done: +if( pcanc > 1.0e-3 ) + mtherr( "hyperg", PLOSS ); + +return( psum ); +} + + + + +/* Power series summation for confluent hypergeometric function */ + + +static float hy1f1pf( float aa, float bb, float xx, float *err ) +{ +float a, b, x, n, a0, sum, t, u, temp; +float an, bn, maxt, pcanc; + +a = aa; +b = bb; +x = xx; +/* 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 > MACHEPF ) + { + if( bn == 0 ) /* check bn first since if both */ + { + mtherr( "hypergf", SING ); + return( MAXNUMF ); /* 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 = fabsf(u); + if( (temp > 1.0 ) && (maxt > (MAXNUMF/temp)) ) + { + pcanc = 1.0; /* estimate 100% error */ + goto blowup; + } + + a0 *= u; + sum += a0; + t = fabsf(a0); + if( t > maxt ) + maxt = t; +/* + if( (maxt/fabsf(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 /= fabsf(sum); +maxt *= MACHEPF; /* this way avoids multiply overflow */ +pcanc = fabsf( MACHEPF * 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 float hy1f1af( float aa, float bb, float xx, float *err ) +{ +float a, b, x, h1, h2, t, u, temp, acanc, asum, err1, err2; + +a = aa; +b = bb; +x = xx; +if( x == 0 ) + { + acanc = 1.0; + asum = MAXNUMF; + goto adone; + } +temp = logf( fabsf(x) ); +t = x + temp * (a-b); +u = -temp * a; + +if( b > 0 ) + { + temp = lgamf(b); + t += temp; + u += temp; + } + +h1 = hyp2f0f( a, a-b+1, -1.0/x, 1, &err1 ); + +temp = expf(u) / gammaf(b-a); +h1 *= temp; +err1 *= temp; + +h2 = hyp2f0f( b-a, 1.0-a, 1.0/x, 2, &err2 ); + +if( a < 0 ) + temp = expf(t) / gammaf(a); +else + temp = expf( t - lgamf(a) ); + +h2 *= temp; +err2 *= temp; + +if( x < 0.0 ) + asum = h1; +else + asum = h2; + +acanc = fabsf(err1) + fabsf(err2); + + +if( b < 0 ) + { + temp = gammaf(b); + asum *= temp; + acanc *= fabsf(temp); + } + + +if( asum != 0.0 ) + acanc /= fabsf(asum); + +acanc *= 30.0; /* fudge factor, since error of asymptotic formula + * often seems this much larger than advertised */ + +adone: + + +*err = acanc; +return( asum ); +} + +/* hyp2f0() */ + +float hyp2f0f(float aa, float bb, float xx, int type, float *err) +{ +float a, b, x, a0, alast, t, tlast, maxt; +float n, an, bn, u, sum, temp; + +a = aa; +b = bb; +x = xx; +an = a; +bn = b; +a0 = 1.0; +alast = 1.0; +sum = 0.0; +n = 1.0; +t = 1.0; +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 = fabsf(u); + if( (temp > 1.0 ) && (maxt > (MAXNUMF/temp)) ) + goto error; + + a0 *= u; + t = fabsf(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.0; + bn += 1.0; + n += 1.0; + if( t > maxt ) + maxt = t; + } +while( t > MACHEPF ); + + +pdone: /* series converged! */ + +/* estimate error due to roundoff and cancellation */ +*err = fabsf( MACHEPF * (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 = MACHEPF * (n + maxt) + fabsf( a0 ); + + +done: +sum += alast; +return( sum ); + +/* series blew up: */ +error: +*err = MAXNUMF; +mtherr( "hypergf", TLOSS ); +return( sum ); +} -- cgit v1.2.3