summaryrefslogtreecommitdiff
path: root/libm/float/incbif.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/float/incbif.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD). -Erik
Diffstat (limited to 'libm/float/incbif.c')
-rw-r--r--libm/float/incbif.c197
1 files changed, 0 insertions, 197 deletions
diff --git a/libm/float/incbif.c b/libm/float/incbif.c
deleted file mode 100644
index 4d8c0652e..000000000
--- a/libm/float/incbif.c
+++ /dev/null
@@ -1,197 +0,0 @@
-/* incbif()
- *
- * Inverse of imcomplete beta integral
- *
- *
- *
- * SYNOPSIS:
- *
- * float a, b, x, y, incbif();
- *
- * x = incbif( a, b, y );
- *
- *
- *
- * DESCRIPTION:
- *
- * Given y, the function finds x such that
- *
- * incbet( a, b, x ) = y.
- *
- * the routine performs up to 10 Newton iterations to find the
- * root of incbet(a,b,x) - y = 0.
- *
- *
- * ACCURACY:
- *
- * Relative error:
- * x a,b
- * arithmetic domain domain # trials peak rms
- * IEEE 0,1 0,100 5000 2.8e-4 8.3e-6
- *
- * Overflow and larger errors may occur for one of a or b near zero
- * and the other large.
- */
-
-
-/*
-Cephes Math Library Release 2.2: July, 1992
-Copyright 1984, 1987, 1992 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
-*/
-
-#include <math.h>
-
-extern float MACHEPF, MINLOGF;
-
-#define fabsf(x) ( (x) < 0 ? -(x) : (x) )
-
-#ifdef ANSIC
-float incbetf(float, float, float);
-float ndtrif(float), expf(float), logf(float), sqrtf(float), lgamf(float);
-#else
-float incbetf();
-float ndtrif(), expf(), logf(), sqrtf(), lgamf();
-#endif
-
-float incbif( float aaa, float bbb, float yyy0 )
-{
-float aa, bb, yy0, a, b, y0;
-float d, y, x, x0, x1, lgm, yp, di;
-int i, rflg;
-
-
-aa = aaa;
-bb = bbb;
-yy0 = yyy0;
-if( yy0 <= 0 )
- return(0.0);
-if( yy0 >= 1.0 )
- return(1.0);
-
-/* approximation to inverse function */
-
-yp = -ndtrif(yy0);
-
-if( yy0 > 0.5 )
- {
- rflg = 1;
- a = bb;
- b = aa;
- y0 = 1.0 - yy0;
- yp = -yp;
- }
-else
- {
- rflg = 0;
- a = aa;
- b = bb;
- y0 = yy0;
- }
-
-
-if( (aa <= 1.0) || (bb <= 1.0) )
- {
- y = 0.5 * yp * yp;
- }
-else
- {
- lgm = (yp * yp - 3.0)* 0.16666666666666667;
- x0 = 2.0/( 1.0/(2.0*a-1.0) + 1.0/(2.0*b-1.0) );
- y = yp * sqrtf( x0 + lgm ) / x0
- - ( 1.0/(2.0*b-1.0) - 1.0/(2.0*a-1.0) )
- * (lgm + 0.833333333333333333 - 2.0/(3.0*x0));
- y = 2.0 * y;
- if( y < MINLOGF )
- {
- x0 = 1.0;
- goto under;
- }
- }
-
-x = a/( a + b * expf(y) );
-y = incbetf( a, b, x );
-yp = (y - y0)/y0;
-if( fabsf(yp) < 0.1 )
- goto newt;
-
-/* Resort to interval halving if not close enough */
-x0 = 0.0;
-x1 = 1.0;
-di = 0.5;
-
-for( i=0; i<20; i++ )
- {
- if( i != 0 )
- {
- x = di * x1 + (1.0-di) * x0;
- y = incbetf( a, b, x );
- yp = (y - y0)/y0;
- if( fabsf(yp) < 1.0e-3 )
- goto newt;
- }
-
- if( y < y0 )
- {
- x0 = x;
- di = 0.5;
- }
- else
- {
- x1 = x;
- di *= di;
- if( di == 0.0 )
- di = 0.5;
- }
- }
-
-if( x0 == 0.0 )
- {
-under:
- mtherr( "incbif", UNDERFLOW );
- goto done;
- }
-
-newt:
-
-x0 = x;
-lgm = lgamf(a+b) - lgamf(a) - lgamf(b);
-
-for( i=0; i<10; i++ )
- {
-/* compute the function at this point */
- if( i != 0 )
- y = incbetf(a,b,x0);
-/* compute the derivative of the function at this point */
- d = (a - 1.0) * logf(x0) + (b - 1.0) * logf(1.0-x0) + lgm;
- if( d < MINLOGF )
- {
- x0 = 0.0;
- goto under;
- }
- d = expf(d);
-/* compute the step to the next approximation of x */
- d = (y - y0)/d;
- x = x0;
- x0 = x0 - d;
- if( x0 <= 0.0 )
- {
- x0 = 0.0;
- goto under;
- }
- if( x0 >= 1.0 )
- {
- x0 = 1.0;
- goto under;
- }
- if( i < 2 )
- continue;
- if( fabsf(d/x0) < 256.0 * MACHEPF )
- goto done;
- }
-
-done:
-if( rflg )
- x0 = 1.0 - x0;
-return( x0 );
-}