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, 460 insertions, 0 deletions
diff --git a/libm/double/hyp2f1.c b/libm/double/hyp2f1.c
new file mode 100644
index 000000000..f2e93106c
--- /dev/null
+++ b/libm/double/hyp2f1.c
@@ -0,0 +1,460 @@
+/* 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-06-01 05:09:13 +0000