libm/float/igamf.c

   1 /*                                                      igamf.c
   2  *
   3  *      Incomplete gamma integral
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float a, x, y, igamf();
  10  *
  11  * y = igamf( a, x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * The function is defined by
  18  *
  19  *                           x
  20  *                            -
  21  *                   1       | |  -t  a-1
  22  *  igam(a,x)  =   -----     |   e   t   dt.
  23  *                  -      | |
  24  *                 | (a)    -
  25  *                           0
  26  *
  27  *
  28  * In this implementation both arguments must be positive.
  29  * The integral is evaluated by either a power series or
  30  * continued fraction expansion, depending on the relative
  31  * values of a and x.
  32  *
  33  *
  34  *
  35  * ACCURACY:
  36  *
  37  *                      Relative error:
  38  * arithmetic   domain     # trials      peak         rms
  39  *    IEEE      0,30        20000       7.8e-6      5.9e-7
  40  *
  41  */
  42 \f/*                                                     igamcf()
  43  *
  44  *      Complemented incomplete gamma integral
  45  *
  46  *
  47  *
  48  * SYNOPSIS:
  49  *
  50  * float a, x, y, igamcf();
  51  *
  52  * y = igamcf( a, x );
  53  *
  54  *
  55  *
  56  * DESCRIPTION:
  57  *
  58  * The function is defined by
  59  *
  60  *
  61  *  igamc(a,x)   =   1 - igam(a,x)
  62  *
  63  *                            inf.
  64  *                              -
  65  *                     1       | |  -t  a-1
  66  *               =   -----     |   e   t   dt.
  67  *                    -      | |
  68  *                   | (a)    -
  69  *                             x
  70  *
  71  *
  72  * In this implementation both arguments must be positive.
  73  * The integral is evaluated by either a power series or
  74  * continued fraction expansion, depending on the relative
  75  * values of a and x.
  76  *
  77  *
  78  *
  79  * ACCURACY:
  80  *
  81  *                      Relative error:
  82  * arithmetic   domain     # trials      peak         rms
  83  *    IEEE      0,30        30000       7.8e-6      5.9e-7
  84  *
  85  */
  86 \f
  87 /*
  88 Cephes Math Library Release 2.2: June, 1992
  89 Copyright 1985, 1987, 1992 by Stephen L. Moshier
  90 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  91 */
  92
  93 #include <math.h>
  94
  95 /* BIG = 1/MACHEPF */
  96 #define BIG   16777216.
  97
  98 extern float MACHEPF, MAXLOGF;
  99
 100 #ifdef ANSIC
 101 float lgamf(float), expf(float), logf(float), igamf(float, float);
 102 #else
 103 float lgamf(), expf(), logf(), igamf();
 104 #endif
 105
 106 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
 107
 108
 109
 110 float igamcf( float aa, float xx )
 111 {
 112 float a, x, ans, c, yc, ax, y, z;
 113 float pk, pkm1, pkm2, qk, qkm1, qkm2;
 114 float r, t;
 115 static float big = BIG;
 116
 117 a = aa;
 118 x = xx;
 119 if( (x <= 0) || ( a <= 0) )
 120         return( 1.0 );
 121
 122 if( (x < 1.0) || (x < a) )
 123         return( 1.0 - igamf(a,x) );
 124
 125 ax = a * logf(x) - x - lgamf(a);
 126 if( ax < -MAXLOGF )
 127         {
 128         mtherr( "igamcf", UNDERFLOW );
 129         return( 0.0 );
 130         }
 131 ax = expf(ax);
 132
 133 /* continued fraction */
 134 y = 1.0 - a;
 135 z = x + y + 1.0;
 136 c = 0.0;
 137 pkm2 = 1.0;
 138 qkm2 = x;
 139 pkm1 = x + 1.0;
 140 qkm1 = z * x;
 141 ans = pkm1/qkm1;
 142
 143 do
 144         {
 145         c += 1.0;
 146         y += 1.0;
 147         z += 2.0;
 148         yc = y * c;
 149         pk = pkm1 * z  -  pkm2 * yc;
 150         qk = qkm1 * z  -  qkm2 * yc;
 151         if( qk != 0 )
 152                 {
 153                 r = pk/qk;
 154                 t = fabsf( (ans - r)/r );
 155                 ans = r;
 156                 }
 157         else
 158                 t = 1.0;
 159         pkm2 = pkm1;
 160         pkm1 = pk;
 161         qkm2 = qkm1;
 162         qkm1 = qk;
 163         if( fabsf(pk) > big )
 164                 {
 165                 pkm2 *= MACHEPF;
 166                 pkm1 *= MACHEPF;
 167                 qkm2 *= MACHEPF;
 168                 qkm1 *= MACHEPF;
 169                 }
 170         }
 171 while( t > MACHEPF );
 172
 173 return( ans * ax );
 174 }
 175
 176
 177
 178 /* left tail of incomplete gamma function:
 179  *
 180  *          inf.      k
 181  *   a  -x   -       x
 182  *  x  e     >   ----------
 183  *           -     -
 184  *          k=0   | (a+k+1)
 185  *
 186  */
 187
 188 float igamf( float aa, float xx )
 189 {
 190 float a, x, ans, ax, c, r;
 191
 192 a = aa;
 193 x = xx;
 194 if( (x <= 0) || ( a <= 0) )
 195         return( 0.0 );
 196
 197 if( (x > 1.0) && (x > a ) )
 198         return( 1.0 - igamcf(a,x) );
 199
 200 /* Compute  x**a * exp(-x) / gamma(a)  */
 201 ax = a * logf(x) - x - lgamf(a);
 202 if( ax < -MAXLOGF )
 203         {
 204         mtherr( "igamf", UNDERFLOW );
 205         return( 0.0 );
 206         }
 207 ax = expf(ax);
 208
 209 /* power series */
 210 r = a;
 211 c = 1.0;
 212 ans = 1.0;
 213
 214 do
 215         {
 216         r += 1.0;
 217         c *= x/r;
 218         ans += c;
 219         }
 220 while( c/ans > MACHEPF );
 221
 222 return( ans * ax/a );
 223 }