diff options
Diffstat (limited to 'libm/float/hyp2f1f.c')
-rw-r--r-- | libm/float/hyp2f1f.c | 442 |
1 files changed, 0 insertions, 442 deletions
diff --git a/libm/float/hyp2f1f.c b/libm/float/hyp2f1f.c deleted file mode 100644 index 01fe54928..000000000 --- a/libm/float/hyp2f1f.c +++ /dev/null @@ -1,442 +0,0 @@ -/* hyp2f1f.c - * - * Gauss hypergeometric function F - * 2 1 - * - * - * SYNOPSIS: - * - * float a, b, c, x, y, hyp2f1f(); - * - * y = hyp2f1f( 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-6 of the nearest integer. - * - * ACCURACY: - * - * Relative error (-1 < x < 1): - * arithmetic domain # trials peak rms - * IEEE 0,3 30000 5.8e-4 4.3e-6 - */ - -/* hyp2f1 */ - - -/* -Cephes Math Library Release 2.2: November, 1992 -Copyright 1984, 1987, 1992 by Stephen L. Moshier -Direct inquiries to 30 Frost Street, Cambridge, MA 02140 -*/ - - -#include <math.h> - -#define EPS 1.0e-5 -#define EPS2 1.0e-5 -#define ETHRESH 1.0e-5 - -extern float MAXNUMF, MACHEPF; - -#define fabsf(x) ( (x) < 0 ? -(x) : (x) ) - -#ifdef ANSIC -float powf(float, float); -static float hys2f1f(float, float, float, float, float *); -static float hyt2f1f(float, float, float, float, float *); -float gammaf(float), logf(float), expf(float), psif(float); -float floorf(float); -#else -float powf(), gammaf(), logf(), expf(), psif(); -float floorf(); -static float hyt2f1f(), hys2f1f(); -#endif - -#define roundf(x) (floorf((x)+(float )0.5)) - - - - -float hyp2f1f( float aa, float bb, float cc, float xx ) -{ -float a, b, c, x; -float d, d1, d2, e; -float p, q, r, s, y, ax; -float ia, ib, ic, id, err; -int flag, i, aid; - -a = aa; -b = bb; -c = cc; -x = xx; -err = 0.0; -ax = fabsf(x); -s = 1.0 - x; -flag = 0; -ia = roundf(a); /* nearest integer to a */ -ib = roundf(b); - -if( a <= 0 ) - { - if( fabsf(a-ia) < EPS ) /* a is a negative integer */ - flag |= 1; - } - -if( b <= 0 ) - { - if( fabsf(b-ib) < EPS ) /* b is a negative integer */ - flag |= 2; - } - -if( ax < 1.0 ) - { - if( fabsf(b-c) < EPS ) /* b = c */ - { - y = powf( s, -a ); /* s to the -a power */ - goto hypdon; - } - if( fabsf(a-c) < EPS ) /* a = c */ - { - y = powf( s, -b ); /* s to the -b power */ - goto hypdon; - } - } - - - -if( c <= 0.0 ) - { - ic = roundf(c); /* nearest integer to c */ - if( fabsf(c-ic) < EPS ) /* c is a negative integer */ - { - /* check if termination before explosion */ - if( (flag & 1) && (ia > ic) ) - goto hypok; - if( (flag & 2) && (ib > ic) ) - goto hypok; - goto hypdiv; - } - } - -if( flag ) /* function is a polynomial */ - goto hypok; - -if( ax > 1.0 ) /* series diverges */ - goto hypdiv; - -p = c - a; -ia = roundf(p); -if( (ia <= 0.0) && (fabsf(p-ia) < EPS) ) /* negative int c - a */ - flag |= 4; - -r = c - b; -ib = roundf(r); /* nearest integer to r */ -if( (ib <= 0.0) && (fabsf(r-ib) < EPS) ) /* negative int c - b */ - flag |= 8; - -d = c - a - b; -id = roundf(d); /* nearest integer to d */ -q = fabsf(d-id); - -if( fabsf(ax-1.0) < EPS ) /* |x| == 1.0 */ - { - if( x > 0.0 ) - { - if( flag & 12 ) /* negative int c-a or c-b */ - { - if( d >= 0.0 ) - goto hypf; - else - goto hypdiv; - } - if( d <= 0.0 ) - goto hypdiv; - y = gammaf(c)*gammaf(d)/(gammaf(p)*gammaf(r)); - goto hypdon; - } - - if( d <= -1.0 ) - goto hypdiv; - } - -/* Conditionally make d > 0 by recurrence on c - * AMS55 #15.2.27 - */ -if( d < 0.0 ) - { -/* Try the power series first */ - y = hyt2f1f( a, b, c, x, &err ); - if( err < ETHRESH ) - goto hypdon; -/* Apply the recurrence if power series fails */ - err = 0.0; - aid = 2 - id; - e = c + aid; - d2 = hyp2f1f(a,b,e,x); - d1 = hyp2f1f(a,b,e+1.0,x); - q = a + b + 1.0; - for( i=0; i<aid; i++ ) - { - r = e - 1.0; - y = (e*(r-(2.0*e-q)*x)*d2 + (e-a)*(e-b)*x*d1)/(e*r*s); - e = r; - d1 = d2; - d2 = y; - } - goto hypdon; - } - - -if( flag & 12 ) - goto hypf; /* negative integer c-a or c-b */ - -hypok: -y = hyt2f1f( a, b, c, x, &err ); - -hypdon: -if( err > ETHRESH ) - { - mtherr( "hyp2f1", PLOSS ); -/* printf( "Estimated err = %.2e\n", err );*/ - } -return(y); - -/* The transformation for c-a or c-b negative integer - * AMS55 #15.3.3 - */ -hypf: -y = powf( s, d ) * hys2f1f( c-a, c-b, c, x, &err ); -goto hypdon; - -/* The alarm exit */ -hypdiv: -mtherr( "hyp2f1f", OVERFLOW ); -return( MAXNUMF ); -} - - - - -/* Apply transformations for |x| near 1 - * then call the power series - */ -static float hyt2f1f( float aa, float bb, float cc, float xx, float *loss ) -{ -float a, b, c, x; -float p, q, r, s, t, y, d, err, err1; -float ax, id, d1, d2, e, y1; -int i, aid; - -a = aa; -b = bb; -c = cc; -x = xx; -err = 0.0; -s = 1.0 - x; -if( x < -0.5 ) - { - if( b > a ) - y = powf( s, -a ) * hys2f1f( a, c-b, c, -x/s, &err ); - - else - y = powf( s, -b ) * hys2f1f( c-a, b, c, -x/s, &err ); - - goto done; - } - - - -d = c - a - b; -id = roundf(d); /* nearest integer to d */ - -if( x > 0.8 ) -{ - -if( fabsf(d-id) > EPS2 ) /* test for integer c-a-b */ - { -/* Try the power series first */ - y = hys2f1f( a, b, c, x, &err ); - if( err < ETHRESH ) - goto done; -/* If power series fails, then apply AMS55 #15.3.6 */ - q = hys2f1f( a, b, 1.0-d, s, &err ); - q *= gammaf(d) /(gammaf(c-a) * gammaf(c-b)); - r = powf(s,d) * hys2f1f( c-a, c-b, d+1.0, s, &err1 ); - r *= gammaf(-d)/(gammaf(a) * gammaf(b)); - y = q + r; - - q = fabsf(q); /* estimate cancellation error */ - r = fabsf(r); - if( q > r ) - r = q; - err += err1 + (MACHEPF*r)/y; - - y *= gammaf(c); - goto done; - } -else - { -/* Psi function expansion, AMS55 #15.3.10, #15.3.11, #15.3.12 */ - if( id >= 0.0 ) - { - e = d; - d1 = d; - d2 = 0.0; - aid = id; - } - else - { - e = -d; - d1 = 0.0; - d2 = d; - aid = -id; - } - - ax = logf(s); - - /* sum for t = 0 */ - y = psif(1.0) + psif(1.0+e) - psif(a+d1) - psif(b+d1) - ax; - y /= gammaf(e+1.0); - - p = (a+d1) * (b+d1) * s / gammaf(e+2.0); /* Poch for t=1 */ - t = 1.0; - do - { - r = psif(1.0+t) + psif(1.0+t+e) - psif(a+t+d1) - - psif(b+t+d1) - ax; - q = p * r; - y += q; - p *= s * (a+t+d1) / (t+1.0); - p *= (b+t+d1) / (t+1.0+e); - t += 1.0; - } - while( fabsf(q/y) > EPS ); - - - if( id == 0.0 ) - { - y *= gammaf(c)/(gammaf(a)*gammaf(b)); - goto psidon; - } - - y1 = 1.0; - - if( aid == 1 ) - goto nosum; - - t = 0.0; - p = 1.0; - for( i=1; i<aid; i++ ) - { - r = 1.0-e+t; - p *= s * (a+t+d2) * (b+t+d2) / r; - t += 1.0; - p /= t; - y1 += p; - } - - -nosum: - p = gammaf(c); - y1 *= gammaf(e) * p / (gammaf(a+d1) * gammaf(b+d1)); - y *= p / (gammaf(a+d2) * gammaf(b+d2)); - if( (aid & 1) != 0 ) - y = -y; - - q = powf( s, id ); /* s to the id power */ - if( id > 0.0 ) - y *= q; - else - y1 *= q; - - y += y1; -psidon: - goto done; - } -} - - -/* Use defining power series if no special cases */ -y = hys2f1f( a, b, c, x, &err ); - -done: -*loss = err; -return(y); -} - - - - - -/* Defining power series expansion of Gauss hypergeometric function */ - -static float hys2f1f( float aa, float bb, float cc, float xx, float *loss ) -{ -int i; -float a, b, c, x; -float f, g, h, k, m, s, u, umax; - - -a = aa; -b = bb; -c = cc; -x = xx; -i = 0; -umax = 0.0; -f = a; -g = b; -h = c; -k = 0.0; -s = 1.0; -u = 1.0; - -do - { - if( fabsf(h) < EPS ) - return( MAXNUMF ); - m = k + 1.0; - u = u * ((f+k) * (g+k) * x / ((h+k) * m)); - s += u; - k = fabsf(u); /* remember largest term summed */ - if( k > umax ) - umax = k; - k = m; - if( ++i > 10000 ) /* should never happen */ - { - *loss = 1.0; - return(s); - } - } -while( fabsf(u/s) > MACHEPF ); - -/* return estimated relative error */ -*loss = (MACHEPF*umax)/fabsf(s) + (MACHEPF*i); - -return(s); -} - - |