summaryrefslogtreecommitdiff
path: root/libm/double/hyp2f1.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/double/hyp2f1.c')
-rw-r--r--libm/double/hyp2f1.c460
1 files changed, 0 insertions, 460 deletions
diff --git a/libm/double/hyp2f1.c b/libm/double/hyp2f1.c
deleted file mode 100644
index f2e93106c..000000000
--- a/libm/double/hyp2f1.c
+++ /dev/null
@@ -1,460 +0,0 @@
-/* hyp2f1.c
- *
- * Gauss hypergeometric function F
- * 2 1
- *
- *
- * SYNOPSIS:
- *
- * double a, b, c, x, y, hyp2f1();
- *
- * y = hyp2f1( 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-14 of the nearest integer
- * (1.0e-13 for IEEE arithmetic).
- *
- * ACCURACY:
- *
- *
- * Relative error (-1 < x < 1):
- * arithmetic domain # trials peak rms
- * IEEE -1,7 230000 1.2e-11 5.2e-14
- *
- * Several special cases also tested with a, b, c in
- * the range -7 to 7.
- *
- * ERROR MESSAGES:
- *
- * A "partial loss of precision" message is printed if
- * the internally estimated relative error exceeds 1^-12.
- * A "singularity" message is printed on overflow or
- * in cases not addressed (such as x < -1).
- */
-
-/* hyp2f1 */
-
-
-/*
-Cephes Math Library Release 2.8: June, 2000
-Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier
-*/
-
-
-#include <math.h>
-
-#ifdef DEC
-#define EPS 1.0e-14
-#define EPS2 1.0e-11
-#endif
-
-#ifdef IBMPC
-#define EPS 1.0e-13
-#define EPS2 1.0e-10
-#endif
-
-#ifdef MIEEE
-#define EPS 1.0e-13
-#define EPS2 1.0e-10
-#endif
-
-#ifdef UNK
-#define EPS 1.0e-13
-#define EPS2 1.0e-10
-#endif
-
-#define ETHRESH 1.0e-12
-
-#ifdef ANSIPROT
-extern double fabs ( double );
-extern double pow ( double, double );
-extern double round ( double );
-extern double gamma ( double );
-extern double log ( double );
-extern double exp ( double );
-extern double psi ( double );
-static double hyt2f1(double, double, double, double, double *);
-static double hys2f1(double, double, double, double, double *);
-double hyp2f1(double, double, double, double);
-#else
-double fabs(), pow(), round(), gamma(), log(), exp(), psi();
-static double hyt2f1();
-static double hys2f1();
-double hyp2f1();
-#endif
-extern double MAXNUM, MACHEP;
-
-double hyp2f1( a, b, c, x )
-double a, b, c, x;
-{
-double d, d1, d2, e;
-double p, q, r, s, y, ax;
-double ia, ib, ic, id, err;
-int flag, i, aid;
-
-err = 0.0;
-ax = fabs(x);
-s = 1.0 - x;
-flag = 0;
-ia = round(a); /* nearest integer to a */
-ib = round(b);
-
-if( a <= 0 )
- {
- if( fabs(a-ia) < EPS ) /* a is a negative integer */
- flag |= 1;
- }
-
-if( b <= 0 )
- {
- if( fabs(b-ib) < EPS ) /* b is a negative integer */
- flag |= 2;
- }
-
-if( ax < 1.0 )
- {
- if( fabs(b-c) < EPS ) /* b = c */
- {
- y = pow( s, -a ); /* s to the -a power */
- goto hypdon;
- }
- if( fabs(a-c) < EPS ) /* a = c */
- {
- y = pow( s, -b ); /* s to the -b power */
- goto hypdon;
- }
- }
-
-
-
-if( c <= 0.0 )
- {
- ic = round(c); /* nearest integer to c */
- if( fabs(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 = round(p); /* nearest integer to c-a */
-if( (ia <= 0.0) && (fabs(p-ia) < EPS) ) /* negative int c - a */
- flag |= 4;
-
-r = c - b;
-ib = round(r); /* nearest integer to c-b */
-if( (ib <= 0.0) && (fabs(r-ib) < EPS) ) /* negative int c - b */
- flag |= 8;
-
-d = c - a - b;
-id = round(d); /* nearest integer to d */
-q = fabs(d-id);
-
-/* Thanks to Christian Burger <BURGER@DMRHRZ11.HRZ.Uni-Marburg.DE>
- * for reporting a bug here. */
-if( fabs(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 = gamma(c)*gamma(d)/(gamma(p)*gamma(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 = hyt2f1( 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 = hyp2f1(a,b,e,x);
- d1 = hyp2f1(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 = hyt2f1( 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 = pow( s, d ) * hys2f1( c-a, c-b, c, x, &err );
-goto hypdon;
-
-/* The alarm exit */
-hypdiv:
-mtherr( "hyp2f1", OVERFLOW );
-return( MAXNUM );
-}
-
-
-
-
-
-
-/* Apply transformations for |x| near 1
- * then call the power series
- */
-static double hyt2f1( a, b, c, x, loss )
-double a, b, c, x;
-double *loss;
-{
-double p, q, r, s, t, y, d, err, err1;
-double ax, id, d1, d2, e, y1;
-int i, aid;
-
-err = 0.0;
-s = 1.0 - x;
-if( x < -0.5 )
- {
- if( b > a )
- y = pow( s, -a ) * hys2f1( a, c-b, c, -x/s, &err );
-
- else
- y = pow( s, -b ) * hys2f1( c-a, b, c, -x/s, &err );
-
- goto done;
- }
-
-d = c - a - b;
-id = round(d); /* nearest integer to d */
-
-if( x > 0.9 )
-{
-if( fabs(d-id) > EPS ) /* test for integer c-a-b */
- {
-/* Try the power series first */
- y = hys2f1( a, b, c, x, &err );
- if( err < ETHRESH )
- goto done;
-/* If power series fails, then apply AMS55 #15.3.6 */
- q = hys2f1( a, b, 1.0-d, s, &err );
- q *= gamma(d) /(gamma(c-a) * gamma(c-b));
- r = pow(s,d) * hys2f1( c-a, c-b, d+1.0, s, &err1 );
- r *= gamma(-d)/(gamma(a) * gamma(b));
- y = q + r;
-
- q = fabs(q); /* estimate cancellation error */
- r = fabs(r);
- if( q > r )
- r = q;
- err += err1 + (MACHEP*r)/y;
-
- y *= gamma(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 = log(s);
-
- /* sum for t = 0 */
- y = psi(1.0) + psi(1.0+e) - psi(a+d1) - psi(b+d1) - ax;
- y /= gamma(e+1.0);
-
- p = (a+d1) * (b+d1) * s / gamma(e+2.0); /* Poch for t=1 */
- t = 1.0;
- do
- {
- r = psi(1.0+t) + psi(1.0+t+e) - psi(a+t+d1)
- - psi(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( fabs(q/y) > EPS );
-
-
- if( id == 0.0 )
- {
- y *= gamma(c)/(gamma(a)*gamma(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 = gamma(c);
- y1 *= gamma(e) * p / (gamma(a+d1) * gamma(b+d1));
-
- y *= p / (gamma(a+d2) * gamma(b+d2));
- if( (aid & 1) != 0 )
- y = -y;
-
- q = pow( 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 = hys2f1( a, b, c, x, &err );
-
-done:
-*loss = err;
-return(y);
-}
-
-
-
-
-
-/* Defining power series expansion of Gauss hypergeometric function */
-
-static double hys2f1( a, b, c, x, loss )
-double a, b, c, x;
-double *loss; /* estimates loss of significance */
-{
-double f, g, h, k, m, s, u, umax;
-int i;
-
-i = 0;
-umax = 0.0;
-f = a;
-g = b;
-h = c;
-s = 1.0;
-u = 1.0;
-k = 0.0;
-do
- {
- if( fabs(h) < EPS )
- {
- *loss = 1.0;
- return( MAXNUM );
- }
- m = k + 1.0;
- u = u * ((f+k) * (g+k) * x / ((h+k) * m));
- s += u;
- k = fabs(u); /* remember largest term summed */
- if( k > umax )
- umax = k;
- k = m;
- if( ++i > 10000 ) /* should never happen */
- {
- *loss = 1.0;
- return(s);
- }
- }
-while( fabs(u/s) > MACHEP );
-
-/* return estimated relative error */
-*loss = (MACHEP*umax)/fabs(s) + (MACHEP*i);
-
-return(s);
-}
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-20 06:06:18 +0000