Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : IGAMI.C

 
Output of file : IGAMI.C contained in archive : CEPHES22.ZIP
/* igami()
*
* Inverse of complemented imcomplete gamma integral
*
*
*
* SYNOPSIS:
*
* double a, x, y, igami();
*
* x = igami( a, y );
*
*
*
* DESCRIPTION:
*
* Given y, the function finds x such that
*
* igamc( a, x ) = y.
*
* Starting with the approximate value
*
* 3
* x = a t
*
* where
*
* t = 1 - d - ndtri(y) sqrt(d)
*
* and
*
* d = 1/9a,
*
* the routine performs up to 10 Newton iterations to find the
* root of igamc(a,x) - y = 0.
*
*
* ACCURACY:
*
* Tested for a ranging from 0.5 to 30 and x from 0 to 0.5.
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC 0,0.5 3400 8.8e-16 1.3e-16
* IEEE 0,0.5 10000 1.1e-14 1.0e-15
*
*/

/*
Cephes Math Library Release 2.0: April, 1987
Copyright 1984, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

#include "mconf.h"

extern double MACHEP, MAXNUM, MAXLOG, MINLOG;

double igami( a, y0 )
double a, y0;
{
double d, y, x0, lgm;
int i;
double igamc();
double ndtri(), exp(), fabs(), log(), sqrt(), lgam();


/* approximation to inverse function */
d = 1.0/(9.0*a);
y = ( 1.0 - d - ndtri(y0) * sqrt(d) );
x0 = a * y * y * y;

lgm = lgam(a);

for( i=0; i<10; i++ )
{
if( x0 <= 0.0 )
{
mtherr( "igami", UNDERFLOW );
return(0.0);
}
y = igamc(a,x0);
/* compute the derivative of the function at this point */
d = (a - 1.0) * log(x0) - x0 - lgm;
if( d < -MAXLOG )
{
mtherr( "igami", UNDERFLOW );
goto done;
}
d = -exp(d);
/* compute the step to the next approximation of x */
if( d == 0.0 )
goto done;
d = (y - y0)/d;
x0 = x0 - d;
if( i < 3 )
continue;
if( fabs(d/x0) < 2.0 * MACHEP )
goto done;
}

done:
return( x0 );
}


  3 Responses to “Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : IGAMI.C

  1. Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!

  2. This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.

  3. But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/