summaryrefslogtreecommitdiff
path: root/libm/ldouble/incbetl.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/ldouble/incbetl.c')
-rw-r--r--libm/ldouble/incbetl.c406
1 files changed, 0 insertions, 406 deletions
diff --git a/libm/ldouble/incbetl.c b/libm/ldouble/incbetl.c
deleted file mode 100644
index fc85ead4c..000000000
--- a/libm/ldouble/incbetl.c
+++ /dev/null
@@ -1,406 +0,0 @@
-/* incbetl.c
- *
- * Incomplete beta integral
- *
- *
- * SYNOPSIS:
- *
- * long double a, b, x, y, incbetl();
- *
- * y = incbetl( a, b, x );
- *
- *
- * DESCRIPTION:
- *
- * Returns incomplete beta integral of the arguments, evaluated
- * from zero to x. The function is defined as
- *
- * x
- * - -
- * | (a+b) | | a-1 b-1
- * ----------- | t (1-t) dt.
- * - - | |
- * | (a) | (b) -
- * 0
- *
- * The domain of definition is 0 <= x <= 1. In this
- * implementation a and b are restricted to positive values.
- * The integral from x to 1 may be obtained by the symmetry
- * relation
- *
- * 1 - incbet( a, b, x ) = incbet( b, a, 1-x ).
- *
- * The integral is evaluated by a continued fraction expansion
- * or, when b*x is small, by a power series.
- *
- * ACCURACY:
- *
- * Tested at random points (a,b,x) with x between 0 and 1.
- * arithmetic domain # trials peak rms
- * IEEE 0,5 20000 4.5e-18 2.4e-19
- * IEEE 0,100 100000 3.9e-17 1.0e-17
- * Half-integer a, b:
- * IEEE .5,10000 100000 3.9e-14 4.4e-15
- * Outputs smaller than the IEEE gradual underflow threshold
- * were excluded from these statistics.
- *
- * ERROR MESSAGES:
- *
- * message condition value returned
- * incbetl domain x<0, x>1 0.0
- */
-
-
-/*
-Cephes Math Library, Release 2.3: January, 1995
-Copyright 1984, 1995 by Stephen L. Moshier
-*/
-
-#include <math.h>
-
-#define MAXGAML 1755.455L
-static long double big = 9.223372036854775808e18L;
-static long double biginv = 1.084202172485504434007e-19L;
-extern long double MACHEPL, MINLOGL, MAXLOGL;
-
-#ifdef ANSIPROT
-extern long double gammal ( long double );
-extern long double lgaml ( long double );
-extern long double expl ( long double );
-extern long double logl ( long double );
-extern long double fabsl ( long double );
-extern long double powl ( long double, long double );
-static long double incbcfl( long double, long double, long double );
-static long double incbdl( long double, long double, long double );
-static long double pseriesl( long double, long double, long double );
-#else
-long double gammal(), lgaml(), expl(), logl(), fabsl(), powl();
-static long double incbcfl(), incbdl(), pseriesl();
-#endif
-
-long double incbetl( aa, bb, xx )
-long double aa, bb, xx;
-{
-long double a, b, t, x, w, xc, y;
-int flag;
-
-if( aa <= 0.0L || bb <= 0.0L )
- goto domerr;
-
-if( (xx <= 0.0L) || ( xx >= 1.0L) )
- {
- if( xx == 0.0L )
- return( 0.0L );
- if( xx == 1.0L )
- return( 1.0L );
-domerr:
- mtherr( "incbetl", DOMAIN );
- return( 0.0L );
- }
-
-flag = 0;
-if( (bb * xx) <= 1.0L && xx <= 0.95L)
- {
- t = pseriesl(aa, bb, xx);
- goto done;
- }
-
-w = 1.0L - xx;
-
-/* Reverse a and b if x is greater than the mean. */
-if( xx > (aa/(aa+bb)) )
- {
- flag = 1;
- a = bb;
- b = aa;
- xc = xx;
- x = w;
- }
-else
- {
- a = aa;
- b = bb;
- xc = w;
- x = xx;
- }
-
-if( flag == 1 && (b * x) <= 1.0L && x <= 0.95L)
- {
- t = pseriesl(a, b, x);
- goto done;
- }
-
-/* Choose expansion for optimal convergence */
-y = x * (a+b-2.0L) - (a-1.0L);
-if( y < 0.0L )
- w = incbcfl( a, b, x );
-else
- w = incbdl( a, b, x ) / xc;
-
-/* Multiply w by the factor
- a b _ _ _
- x (1-x) | (a+b) / ( a | (a) | (b) ) . */
-
-y = a * logl(x);
-t = b * logl(xc);
-if( (a+b) < MAXGAML && fabsl(y) < MAXLOGL && fabsl(t) < MAXLOGL )
- {
- t = powl(xc,b);
- t *= powl(x,a);
- t /= a;
- t *= w;
- t *= gammal(a+b) / (gammal(a) * gammal(b));
- goto done;
- }
-else
- {
- /* Resort to logarithms. */
- y += t + lgaml(a+b) - lgaml(a) - lgaml(b);
- y += logl(w/a);
- if( y < MINLOGL )
- t = 0.0L;
- else
- t = expl(y);
- }
-
-done:
-
-if( flag == 1 )
- {
- if( t <= MACHEPL )
- t = 1.0L - MACHEPL;
- else
- t = 1.0L - t;
- }
-return( t );
-}
-
-/* Continued fraction expansion #1
- * for incomplete beta integral
- */
-
-static long double incbcfl( a, b, x )
-long double a, b, x;
-{
-long double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
-long double k1, k2, k3, k4, k5, k6, k7, k8;
-long double r, t, ans, thresh;
-int n;
-
-k1 = a;
-k2 = a + b;
-k3 = a;
-k4 = a + 1.0L;
-k5 = 1.0L;
-k6 = b - 1.0L;
-k7 = k4;
-k8 = a + 2.0L;
-
-pkm2 = 0.0L;
-qkm2 = 1.0L;
-pkm1 = 1.0L;
-qkm1 = 1.0L;
-ans = 1.0L;
-r = 1.0L;
-n = 0;
-thresh = 3.0L * MACHEPL;
-do
- {
-
- xk = -( x * k1 * k2 )/( k3 * k4 );
- pk = pkm1 + pkm2 * xk;
- qk = qkm1 + qkm2 * xk;
- pkm2 = pkm1;
- pkm1 = pk;
- qkm2 = qkm1;
- qkm1 = qk;
-
- xk = ( x * k5 * k6 )/( k7 * k8 );
- pk = pkm1 + pkm2 * xk;
- qk = qkm1 + qkm2 * xk;
- pkm2 = pkm1;
- pkm1 = pk;
- qkm2 = qkm1;
- qkm1 = qk;
-
- if( qk != 0.0L )
- r = pk/qk;
- if( r != 0.0L )
- {
- t = fabsl( (ans - r)/r );
- ans = r;
- }
- else
- t = 1.0L;
-
- if( t < thresh )
- goto cdone;
-
- k1 += 1.0L;
- k2 += 1.0L;
- k3 += 2.0L;
- k4 += 2.0L;
- k5 += 1.0L;
- k6 -= 1.0L;
- k7 += 2.0L;
- k8 += 2.0L;
-
- if( (fabsl(qk) + fabsl(pk)) > big )
- {
- pkm2 *= biginv;
- pkm1 *= biginv;
- qkm2 *= biginv;
- qkm1 *= biginv;
- }
- if( (fabsl(qk) < biginv) || (fabsl(pk) < biginv) )
- {
- pkm2 *= big;
- pkm1 *= big;
- qkm2 *= big;
- qkm1 *= big;
- }
- }
-while( ++n < 400 );
-mtherr( "incbetl", PLOSS );
-
-cdone:
-return(ans);
-}
-
-
-/* Continued fraction expansion #2
- * for incomplete beta integral
- */
-
-static long double incbdl( a, b, x )
-long double a, b, x;
-{
-long double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
-long double k1, k2, k3, k4, k5, k6, k7, k8;
-long double r, t, ans, z, thresh;
-int n;
-
-k1 = a;
-k2 = b - 1.0L;
-k3 = a;
-k4 = a + 1.0L;
-k5 = 1.0L;
-k6 = a + b;
-k7 = a + 1.0L;
-k8 = a + 2.0L;
-
-pkm2 = 0.0L;
-qkm2 = 1.0L;
-pkm1 = 1.0L;
-qkm1 = 1.0L;
-z = x / (1.0L-x);
-ans = 1.0L;
-r = 1.0L;
-n = 0;
-thresh = 3.0L * MACHEPL;
-do
- {
-
- xk = -( z * k1 * k2 )/( k3 * k4 );
- pk = pkm1 + pkm2 * xk;
- qk = qkm1 + qkm2 * xk;
- pkm2 = pkm1;
- pkm1 = pk;
- qkm2 = qkm1;
- qkm1 = qk;
-
- xk = ( z * k5 * k6 )/( k7 * k8 );
- pk = pkm1 + pkm2 * xk;
- qk = qkm1 + qkm2 * xk;
- pkm2 = pkm1;
- pkm1 = pk;
- qkm2 = qkm1;
- qkm1 = qk;
-
- if( qk != 0.0L )
- r = pk/qk;
- if( r != 0.0L )
- {
- t = fabsl( (ans - r)/r );
- ans = r;
- }
- else
- t = 1.0L;
-
- if( t < thresh )
- goto cdone;
-
- k1 += 1.0L;
- k2 -= 1.0L;
- k3 += 2.0L;
- k4 += 2.0L;
- k5 += 1.0L;
- k6 += 1.0L;
- k7 += 2.0L;
- k8 += 2.0L;
-
- if( (fabsl(qk) + fabsl(pk)) > big )
- {
- pkm2 *= biginv;
- pkm1 *= biginv;
- qkm2 *= biginv;
- qkm1 *= biginv;
- }
- if( (fabsl(qk) < biginv) || (fabsl(pk) < biginv) )
- {
- pkm2 *= big;
- pkm1 *= big;
- qkm2 *= big;
- qkm1 *= big;
- }
- }
-while( ++n < 400 );
-mtherr( "incbetl", PLOSS );
-
-cdone:
-return(ans);
-}
-
-/* Power series for incomplete gamma integral.
- Use when b*x is small. */
-
-static long double pseriesl( a, b, x )
-long double a, b, x;
-{
-long double s, t, u, v, n, t1, z, ai;
-
-ai = 1.0L / a;
-u = (1.0L - b) * x;
-v = u / (a + 1.0L);
-t1 = v;
-t = u;
-n = 2.0L;
-s = 0.0L;
-z = MACHEPL * ai;
-while( fabsl(v) > z )
- {
- u = (n - b) * x / n;
- t *= u;
- v = t / (a + n);
- s += v;
- n += 1.0L;
- }
-s += t1;
-s += ai;
-
-u = a * logl(x);
-if( (a+b) < MAXGAML && fabsl(u) < MAXLOGL )
- {
- t = gammal(a+b)/(gammal(a)*gammal(b));
- s = s * t * powl(x,a);
- }
-else
- {
- t = lgaml(a+b) - lgaml(a) - lgaml(b) + u + logl(s);
- if( t < MINLOGL )
- s = 0.0L;
- else
- s = expl(t);
- }
-return(s);
-}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-06 21:04:58 +0000