lib/libm/gammaf.c

   1 /*                                                      gammaf.c
   2  *
   3  *      Gamma function
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, y, gammaf();
  10  * extern int sgngamf;
  11  *
  12  * y = gammaf( x );
  13  *
  14  *
  15  *
  16  * DESCRIPTION:
  17  *
  18  * Returns gamma function of the argument.  The result is
  19  * correctly signed, and the sign (+1 or -1) is also
  20  * returned in a global (extern) variable named sgngamf.
  21  * This same variable is also filled in by the logarithmic
  22  * gamma function lgam().
  23  *
  24  * Arguments between 0 and 10 are reduced by recurrence and the
  25  * function is approximated by a polynomial function covering
  26  * the interval (2,3).  Large arguments are handled by Stirling's
  27  * formula. Negative arguments are made positive using
  28  * a reflection formula.
  29  *
  30  *
  31  * ACCURACY:
  32  *
  33  *                      Relative error:
  34  * arithmetic   domain     # trials      peak         rms
  35  *    IEEE       0,-33      100,000     5.7e-7      1.0e-7
  36  *    IEEE       -33,0      100,000     6.1e-7      1.2e-7
  37  *
  38  *
  39  */\f
  40 /*                                                      lgamf()
  41  *
  42  *      Natural logarithm of gamma function
  43  *
  44  *
  45  *
  46  * SYNOPSIS:
  47  *
  48  * float x, y, lgamf();
  49  * extern int sgngamf;
  50  *
  51  * y = lgamf( x );
  52  *
  53  *
  54  *
  55  * DESCRIPTION:
  56  *
  57  * Returns the base e (2.718...) logarithm of the absolute
  58  * value of the gamma function of the argument.
  59  * The sign (+1 or -1) of the gamma function is returned in a
  60  * global (extern) variable named sgngamf.
  61  *
  62  * For arguments greater than 6.5, the logarithm of the gamma
  63  * function is approximated by the logarithmic version of
  64  * Stirling's formula.  Arguments between 0 and +6.5 are reduced by
  65  * by recurrence to the interval [.75,1.25] or [1.5,2.5] of a rational
  66  * approximation.  The cosecant reflection formula is employed for
  67  * arguments less than zero.
  68  *
  69  * Arguments greater than MAXLGM = 2.035093e36 return MAXNUM and an
  70  * error message.
  71  *
  72  *
  73  *
  74  * ACCURACY:
  75  *
  76  *
  77  *
  78  * arithmetic      domain        # trials     peak         rms
  79  *    IEEE        -100,+100       500,000    7.4e-7       6.8e-8
  80  * The error criterion was relative when the function magnitude
  81  * was greater than one but absolute when it was less than one.
  82  * The routine has low relative error for positive arguments.
  83  *
  84  * The following test used the relative error criterion.
  85  *    IEEE    -2, +3              100000     4.0e-7      5.6e-8
  86  *
  87  */
  88 \f
  89 /*                                                      gamma.c */
  90 /*      gamma function  */
  91
  92 /*
  93 Cephes Math Library Release 2.7:  July, 1998
  94 Copyright 1984, 1987, 1989, 1992, 1998 by Stephen L. Moshier
  95 */
  96
  97
  98 #include "mconf.h"
  99
 100 /* define MAXGAM 34.84425627277176174 */
 101
 102 /* Stirling's formula for the gamma function
 103  * gamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) ( 1 + 1/x P(1/x) )
 104  * .028 < 1/x < .1
 105  * relative error < 1.9e-11
 106  */
 107 static float STIR[] = {
 108 -2.705194986674176E-003,
 109  3.473255786154910E-003,
 110  8.333331788340907E-002,
 111 };
 112 static float MAXSTIR = 26.77;
 113 static float SQTPIF = 2.50662827463100050242; /* sqrt( 2 pi ) */
 114
 115 int sgngamf = 0;
 116 extern int sgngamf;
 117 extern float MAXLOGF, MAXNUMF, PIF;
 118
 119 #ifdef ANSIC
 120 float expf(float);
 121 float logf(float);
 122 float powf( float, float );
 123 float sinf(float);
 124 float gammaf(float);
 125 float floorf(float);
 126 static float stirf(float);
 127 float polevlf( float, float *, int );
 128 float p1evlf( float, float *, int );
 129 #else
 130 float expf(), logf(), powf(), sinf(), floorf();
 131 float polevlf(), p1evlf();
 132 static float stirf();
 133 #endif
 134
 135 /* Gamma function computed by Stirling's formula,
 136  * sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x))
 137  * The polynomial STIR is valid for 33 <= x <= 172.
 138  */
 139 #ifdef ANSIC
 140 static float stirf( float xx )
 141 #else
 142 static float stirf(xx)
 143 double xx;
 144 #endif
 145 {
 146 float x, y, w, v;
 147
 148 x = xx;
 149 w = 1.0/x;
 150 w = 1.0 + w * polevlf( w, STIR, 2 );
 151 y = expf( -x );
 152 if( x > MAXSTIR )
 153         { /* Avoid overflow in pow() */
 154         v = powf( x, 0.5 * x - 0.25 );
 155         y *= v;
 156         y *= v;
 157         }
 158 else
 159         {
 160         y = powf( x, x - 0.5 ) * y;
 161         }
 162 y = SQTPIF * y * w;
 163 return( y );
 164 }
 165
 166
 167 /* gamma(x+2), 0 < x < 1 */
 168 static float P[] = {
 169  1.536830450601906E-003,
 170  5.397581592950993E-003,
 171  4.130370201859976E-003,
 172  7.232307985516519E-002,
 173  8.203960091619193E-002,
 174  4.117857447645796E-001,
 175  4.227867745131584E-001,
 176  9.999999822945073E-001,
 177 };
 178
 179 #ifdef ANSIC
 180 float gammaf( float xx )
 181 #else
 182 float gammaf(xx)
 183 double xx;
 184 #endif
 185 {
 186 float p, q, x, z, nz;
 187 int i, direction, negative;
 188
 189 x = xx;
 190 sgngamf = 1;
 191 negative = 0;
 192 nz = 0.0;
 193 if( x < 0.0 )
 194         {
 195         negative = 1;
 196         q = -x;
 197         p = floorf(q);
 198         if( p == q )
 199                 goto goverf;
 200         i = p;
 201         if( (i & 1) == 0 )
 202                 sgngamf = -1;
 203         nz = q - p;
 204         if( nz > 0.5 )
 205                 {
 206                 p += 1.0;
 207                 nz = q - p;
 208                 }
 209         nz = q * sinf( PIF * nz );
 210         if( nz == 0.0 )
 211                 {
 212 goverf:
 213                 mtherr( "gamma", OVERFLOW );
 214                 return( sgngamf * MAXNUMF);
 215                 }
 216         if( nz < 0 )
 217                 nz = -nz;
 218         x = q;
 219         }
 220 if( x >= 10.0 )
 221         {
 222         z = stirf(x);
 223         }
 224 if( x < 2.0 )
 225         direction = 1;
 226 else
 227         direction = 0;
 228 z = 1.0;
 229 while( x >= 3.0 )
 230         {
 231         x -= 1.0;
 232         z *= x;
 233         }
 234 /*
 235 while( x < 0.0 )
 236         {
 237         if( x > -1.E-4 )
 238                 goto small;
 239         z *=x;
 240         x += 1.0;
 241         }
 242 */
 243 while( x < 2.0 )
 244         {
 245         if( x < 1.e-4 )
 246                 goto small;
 247         z *=x;
 248         x += 1.0;
 249         }
 250
 251 if( direction )
 252         z = 1.0/z;
 253
 254 if( x == 2.0 )
 255         return(z);
 256
 257 x -= 2.0;
 258 p = z * polevlf( x, P, 7 );
 259
 260 gdone:
 261
 262 if( negative )
 263         {
 264         p = sgngamf * PIF/(nz * p );
 265         }
 266 return(p);
 267
 268 small:
 269 if( x == 0.0 )
 270         {
 271         mtherr( "gamma", SING );
 272         return( MAXNUMF );
 273         }
 274 else
 275         {
 276         p = z / ((1.0 + 0.5772156649015329 * x) * x);
 277         goto gdone;
 278         }
 279 }
 280
 281
 282
 283
 284 /* log gamma(x+2), -.5 < x < .5 */
 285 static float B[] = {
 286  6.055172732649237E-004,
 287 -1.311620815545743E-003,
 288  2.863437556468661E-003,
 289 -7.366775108654962E-003,
 290  2.058355474821512E-002,
 291 -6.735323259371034E-002,
 292  3.224669577325661E-001,
 293  4.227843421859038E-001
 294 };
 295
 296 /* log gamma(x+1), -.25 < x < .25 */
 297 static float C[] = {
 298  1.369488127325832E-001,
 299 -1.590086327657347E-001,
 300  1.692415923504637E-001,
 301 -2.067882815621965E-001,
 302  2.705806208275915E-001,
 303 -4.006931650563372E-001,
 304  8.224670749082976E-001,
 305 -5.772156501719101E-001
 306 };
 307
 308 /* log( sqrt( 2*pi ) ) */
 309 static float LS2PI  =  0.91893853320467274178;
 310 #define MAXLGM 2.035093e36
 311 static float PIINV =  0.318309886183790671538;
 312
 313 /* Logarithm of gamma function */
 314
 315
 316 #ifdef ANSIC
 317 float lgamf( float xx )
 318 #else
 319 float lgamf(xx)
 320 double xx;
 321 #endif
 322 {
 323 float p, q, w, z, x;
 324 float nx, tx;
 325 int i, direction;
 326
 327 sgngamf = 1;
 328
 329 x = xx;
 330 if( x < 0.0 )
 331         {
 332         q = -x;
 333         w = lgamf(q); /* note this modifies sgngam! */
 334         p = floorf(q);
 335         if( p == q )
 336                 goto loverf;
 337         i = p;
 338         if( (i & 1) == 0 )
 339                 sgngamf = -1;
 340         else
 341                 sgngamf = 1;
 342         z = q - p;
 343         if( z > 0.5 )
 344                 {
 345                 p += 1.0;
 346                 z = p - q;
 347                 }
 348         z = q * sinf( PIF * z );
 349         if( z == 0.0 )
 350                 goto loverf;
 351         z = -logf( PIINV*z ) - w;
 352         return( z );
 353         }
 354
 355 if( x < 6.5 )
 356         {
 357         direction = 0;
 358         z = 1.0;
 359         tx = x;
 360         nx = 0.0;
 361         if( x >= 1.5 )
 362                 {
 363                 while( tx > 2.5 )
 364                         {
 365                         nx -= 1.0;
 366                         tx = x + nx;
 367                         z *=tx;
 368                         }
 369                 x += nx - 2.0;
 370 iv1r5:
 371                 p = x * polevlf( x, B, 7 );
 372                 goto cont;
 373                 }
 374         if( x >= 1.25 )
 375                 {
 376                 z *= x;
 377                 x -= 1.0; /* x + 1 - 2 */
 378                 direction = 1;
 379                 goto iv1r5;
 380                 }
 381         if( x >= 0.75 )
 382                 {
 383                 x -= 1.0;
 384                 p = x * polevlf( x, C, 7 );
 385                 q = 0.0;
 386                 goto contz;
 387                 }
 388         while( tx < 1.5 )
 389                 {
 390                 if( tx == 0.0 )
 391                         goto loverf;
 392                 z *=tx;
 393                 nx += 1.0;
 394                 tx = x + nx;
 395                 }
 396         direction = 1;
 397         x += nx - 2.0;
 398         p = x * polevlf( x, B, 7 );
 399
 400 cont:
 401         if( z < 0.0 )
 402                 {
 403                 sgngamf = -1;
 404                 z = -z;
 405                 }
 406         else
 407                 {
 408                 sgngamf = 1;
 409                 }
 410         q = logf(z);
 411         if( direction )
 412                 q = -q;
 413 contz:
 414         return( p + q );
 415         }
 416
 417 if( x > MAXLGM )
 418         {
 419 loverf:
 420         mtherr( "lgamf", OVERFLOW );
 421         return( sgngamf * MAXNUMF );
 422         }
 423
 424 /* Note, though an asymptotic formula could be used for x >= 3,
 425  * there is cancellation error in the following if x < 6.5.  */
 426 q = LS2PI - x;
 427 q += ( x - 0.5 ) * logf(x);
 428
 429 if( x <= 1.0e4 )
 430         {
 431         z = 1.0/x;
 432         p = z * z;
 433         q += ((    6.789774945028216E-004 * p
 434                  - 2.769887652139868E-003 ) * p
 435                 +  8.333316229807355E-002 ) * z;
 436         }
 437 return( q );
 438 }