summaryrefslogtreecommitdiff
path: root/libm/double/hyperg.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/double/hyperg.c')
-rw-r--r--libm/double/hyperg.c386
1 files changed, 0 insertions, 386 deletions
diff --git a/libm/double/hyperg.c b/libm/double/hyperg.c
deleted file mode 100644
index 36a3f9781..000000000
--- a/libm/double/hyperg.c
+++ /dev/null
@@ -1,386 +0,0 @@
-/* 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 );
-}
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-19 06:54:29 +0000