1 files changed, 386 insertions, 0 deletions
diff --git a/libm/double/hyperg.c b/libm/double/hyperg.c
new file mode 100644
index 000000000..36a3f9781
--- /dev/null
+++ b/libm/double/hyperg.c
@@ -0,0 +1,386 @@
+/*							hyperg.c
+ *
+ *	Confluent hypergeometric function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double a, b, x, y, hyperg();
+ *
+ * y = hyperg( 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
+ *    DEC       0,30         2000       1.2e-15     1.3e-16
+ qtst1:
+ 21800   max =  1.4200E-14   rms =  1.0841E-15  ave = -5.3640E-17 
+ ltstd:
+ 25500   max = 1.2759e-14   rms = 3.7155e-16  ave = 1.5384e-18 
+ *    IEEE      0,30        30000       1.8e-14     1.1e-15
+ *
+ * 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-12.
+ *
+ */
+
+/*							hyperg.c */
+
+
+/*
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
+*/
+
+#include <math.h>
+
+#ifdef ANSIPROT
+extern double exp ( double );
+extern double log ( double );
+extern double gamma ( double );
+extern double lgam ( double );
+extern double fabs ( double );
+double hyp2f0 ( double, double, double, int, double * );
+static double hy1f1p(double, double, double, double *);
+static double hy1f1a(double, double, double, double *);
+double hyperg (double, double, double);
+#else
+double exp(), log(), gamma(), lgam(), fabs(), hyp2f0();
+static double hy1f1p();
+static double hy1f1a();
+double hyperg();
+#endif
+extern double MAXNUM, MACHEP;
+
+double hyperg( a, b, x)
+double a, b, x;
+{
+double asum, psum, acanc, pcanc, temp;
+
+/* See if a Kummer transformation will help */
+temp = b - a;
+if( fabs(temp) < 0.001 * fabs(a) )
+	return( exp(x) * hyperg( temp, b, -x )  );
+
+
+psum = hy1f1p( a, b, x, &pcanc );
+if( pcanc < 1.0e-15 )
+	goto done;
+
+
+/* try asymptotic series */
+
+asum = hy1f1a( a, b, x, &acanc );
+
+
+/* Pick the result with less estimated error */
+
+if( acanc < pcanc )
+	{
+	pcanc = acanc;
+	psum = asum;
+	}
+
+done:
+if( pcanc > 1.0e-12 )
+	mtherr( "hyperg", PLOSS );
+
+return( psum );
+}
+
+
+
+
+/* Power series summation for confluent hypergeometric function		*/
+
+
+static double hy1f1p( a, b, x, err )
+double a, b, x;
+double *err;
+{
+double n, a0, sum, t, u, temp;
+double an, bn, maxt, pcanc;
+
+
+/* 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 > MACHEP )
+	{
+	if( bn == 0 )			/* check bn first since if both	*/
+		{
+		mtherr( "hyperg", SING );
+		return( MAXNUM );	/* 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 = fabs(u);
+	if( (temp > 1.0 ) && (maxt > (MAXNUM/temp)) )
+		{
+		pcanc = 1.0;	/* estimate 100% error */
+		goto blowup;
+		}
+
+	a0 *= u;
+	sum += a0;
+	t = fabs(a0);
+	if( t > maxt )
+		maxt = t;
+/*
+	if( (maxt/fabs(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 /= fabs(sum);
+maxt *= MACHEP; 	/* this way avoids multiply overflow */
+pcanc = fabs( MACHEP * 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 double hy1f1a( a, b, x, err )
+double a, b, x;
+double *err;
+{
+double h1, h2, t, u, temp, acanc, asum, err1, err2;
+
+if( x == 0 )
+	{
+	acanc = 1.0;
+	asum = MAXNUM;
+	goto adone;
+	}
+temp = log( fabs(x) );
+t = x + temp * (a-b);
+u = -temp * a;
+
+if( b > 0 )
+	{
+	temp = lgam(b);
+	t += temp;
+	u += temp;
+	}
+
+h1 = hyp2f0( a, a-b+1, -1.0/x, 1, &err1 );
+
+temp = exp(u) / gamma(b-a);
+h1 *= temp;
+err1 *= temp;
+
+h2 = hyp2f0( b-a, 1.0-a, 1.0/x, 2, &err2 );
+
+if( a < 0 )
+	temp = exp(t) / gamma(a);
+else
+	temp = exp( t - lgam(a) );
+
+h2 *= temp;
+err2 *= temp;
+
+if( x < 0.0 )
+	asum = h1;
+else
+	asum = h2;
+
+acanc = fabs(err1) + fabs(err2);
+
+
+if( b < 0 )
+	{
+	temp = gamma(b);
+	asum *= temp;
+	acanc *= fabs(temp);
+	}
+
+
+if( asum != 0.0 )
+	acanc /= fabs(asum);
+
+acanc *= 30.0;	/* fudge factor, since error of asymptotic formula
+		 * often seems this much larger than advertised */
+
+adone:
+
+
+*err = acanc;
+return( asum );
+}
+
+/*							hyp2f0()	*/
+
+double hyp2f0( a, b, x, type, err )
+double a, b, x;
+int type;	/* determines what converging factor to use */
+double *err;
+{
+double a0, alast, t, tlast, maxt;
+double n, an, bn, u, sum, temp;
+
+an = a;
+bn = b;
+a0 = 1.0e0;
+alast = 1.0e0;
+sum = 0.0;
+n = 1.0e0;
+t = 1.0e0;
+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 = fabs(u);
+	if( (temp > 1.0 ) && (maxt > (MAXNUM/temp)) )
+		goto error;
+
+	a0 *= u;
+	t = fabs(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.0e0;
+	bn += 1.0e0;
+	n += 1.0e0;
+	if( t > maxt )
+		maxt = t;
+	}
+while( t > MACHEP );
+
+
+pdone:	/* series converged! */
+
+/* estimate error due to roundoff and cancellation */
+*err = fabs(  MACHEP * (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 = MACHEP * (n + maxt)  +  fabs ( a0 );
+
+
+done:
+sum += alast;
+return( sum );
+
+/* series blew up: */
+error:
+*err = MAXNUM;
+mtherr( "hyperg", TLOSS );
+return( sum );
+}