From 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 22 Nov 2001 14:04:29 +0000 Subject: Totally rework the math library, this time based on the MacOs X math library (which is itself based on the math lib from FreeBSD). -Erik --- libm/float/hypergf.c | 384 --------------------------------------------------- 1 file changed, 384 deletions(-) delete mode 100644 libm/float/hypergf.c (limited to 'libm/float/hypergf.c') diff --git a/libm/float/hypergf.c b/libm/float/hypergf.c deleted file mode 100644 index 60d0eb4c5..000000000 --- a/libm/float/hypergf.c +++ /dev/null @@ -1,384 +0,0 @@ -/* 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