Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

  1. /*                                                      gammaf.c
  2.  *
  3.  *      Gamma function
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * float x, y, __tgammaf_r();
  10.  * int* sgngamf;
  11.  * y = __tgammaf_r( x, sgngamf );
  12.  *
  13.  * float x, y, tgammaf();
  14.  * y = tgammaf( x);
  15.  *
  16.  *
  17.  * DESCRIPTION:
  18.  *
  19.  * Returns gamma function of the argument.  The result is
  20.  * correctly signed. In the reentrant version the sign (+1 or -1)
  21.  * is returned in the variable referenced by sgngamf.
  22.  *
  23.  * Arguments between 0 and 10 are reduced by recurrence and the
  24.  * function is approximated by a polynomial function covering
  25.  * the interval (2,3).  Large arguments are handled by Stirling's
  26.  * formula. Negative arguments are made positive using
  27.  * a reflection formula.  
  28.  *
  29.  *
  30.  * ACCURACY:
  31.  *
  32.  *                      Relative error:
  33.  * arithmetic   domain     # trials      peak         rms
  34.  *    IEEE       0,-33      100,000     5.7e-7      1.0e-7
  35.  *    IEEE       -33,0      100,000     6.1e-7      1.2e-7
  36.  *
  37.  *
  38.  */
  39.  
  40. /*
  41. Cephes Math Library Release 2.7:  July, 1998
  42. Copyright 1984, 1987, 1989, 1992, 1998 by Stephen L. Moshier
  43. */
  44.  
  45.  
  46. /*
  47.  * 26-11-2002 Modified for mingw.
  48.  * Danny Smith <dannysmith@users.sourceforge.net>
  49.  */
  50.  
  51.  
  52. #ifndef __MINGW32__
  53. #include "mconf.h"
  54. #else
  55. #include "cephes_mconf.h"
  56. #endif
  57.  
  58. /* define MAXGAM 34.84425627277176174 */
  59.  
  60. /* Stirling's formula for the gamma function
  61.  * gamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) ( 1 + 1/x P(1/x) )
  62.  * .028 < 1/x < .1
  63.  * relative error < 1.9e-11
  64.  */
  65. static const float STIR[] = {
  66. -2.705194986674176E-003,
  67.  3.473255786154910E-003,
  68.  8.333331788340907E-002,
  69. };
  70. static const float MAXSTIR = 26.77;
  71. static const float SQTPIF = 2.50662827463100050242; /* sqrt( 2 pi ) */
  72.  
  73. #ifndef __MINGW32__
  74.  
  75. extern float MAXLOGF, MAXNUMF, PIF;
  76.  
  77. #ifdef ANSIC
  78. float expf(float);
  79. float logf(float);
  80. float powf( float, float );
  81. float sinf(float);
  82. float gammaf(float);
  83. float floorf(float);
  84. static float stirf(float);
  85. float polevlf( float, float *, int );
  86. float p1evlf( float, float *, int );
  87. #else
  88. float expf(), logf(), powf(), sinf(), floorf();
  89. float polevlf(), p1evlf();
  90. static float stirf();
  91. #endif
  92.  
  93. #else /* __MINGW32__ */
  94. static float stirf(float);
  95. #endif
  96.  
  97. /* Gamma function computed by Stirling's formula,
  98.  * sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x))
  99.  * The polynomial STIR is valid for 33 <= x <= 172.
  100.  */
  101. static float stirf( float x )
  102. {
  103. float  y, w, v;
  104.  
  105. w = 1.0/x;
  106. w = 1.0 + w * polevlf( w, STIR, 2 );
  107. y = expf( -x );
  108. if( x > MAXSTIR )
  109.         { /* Avoid overflow in pow() */
  110.         v = powf( x, 0.5 * x - 0.25 );
  111.         y *= v;
  112.         y *= v;
  113.         }
  114. else
  115.         {
  116.         y = powf( x, x - 0.5 ) * y;
  117.         }
  118. y = SQTPIF * y * w;
  119. return( y );
  120. }
  121.  
  122.  
  123. /* gamma(x+2), 0 < x < 1 */
  124. static const float P[] = {
  125.  1.536830450601906E-003,
  126.  5.397581592950993E-003,
  127.  4.130370201859976E-003,
  128.  7.232307985516519E-002,
  129.  8.203960091619193E-002,
  130.  4.117857447645796E-001,
  131.  4.227867745131584E-001,
  132.  9.999999822945073E-001,
  133. };
  134.  
  135. float __tgammaf_r( float x, int* sgngamf)
  136. {
  137. float p, q, z, nz;
  138. int i, direction, negative;
  139.  
  140. #ifdef NANS
  141. if( isnan(x) )
  142.         return(x);
  143. #endif
  144. #ifdef INFINITIES
  145. #ifdef NANS
  146. if( x == INFINITYF )
  147.         return(x);
  148. if( x == -INFINITYF )
  149.         return(NANF);
  150. #else
  151. if( !isfinite(x) )
  152.         return(x);
  153. #endif
  154. #endif
  155.  
  156. *sgngamf = 1;
  157. negative = 0;
  158. nz = 0.0;
  159. if( x < 0.0 )
  160.         {
  161.         negative = 1;
  162.         q = -x;
  163.         p = floorf(q);
  164.         if( p == q )
  165.                 {
  166. gsing:
  167.                 _SET_ERRNO(EDOM);
  168.                 mtherr( "tgammaf", SING );
  169. #ifdef INFINITIES
  170.                 return (INFINITYF);
  171. #else
  172.                 return (MAXNUMF);
  173. #endif
  174.                         }
  175.         i = p;
  176.         if( (i & 1) == 0 )
  177.                 *sgngamf = -1;
  178.         nz = q - p;
  179.         if( nz > 0.5 )
  180.                 {
  181.                 p += 1.0;
  182.                 nz = q - p;
  183.                 }
  184.         nz = q * sinf( PIF * nz );
  185.         if( nz == 0.0 )
  186.                 {
  187.                 _SET_ERRNO(ERANGE);
  188.                 mtherr( "tgamma", OVERFLOW );
  189. #ifdef INFINITIES
  190.                 return( *sgngamf * INFINITYF);
  191. #else
  192.                 return( *sgngamf * MAXNUMF);
  193. #endif
  194.                 }
  195.         if( nz < 0 )
  196.                 nz = -nz;
  197.         x = q;
  198.         }
  199. if( x >= 10.0 )
  200.         {
  201.         z = stirf(x);
  202.         }
  203. if( x < 2.0 )
  204.         direction = 1;
  205. else
  206.         direction = 0;
  207. z = 1.0;
  208. while( x >= 3.0 )
  209.         {
  210.         x -= 1.0;
  211.         z *= x;
  212.         }
  213. /*
  214. while( x < 0.0 )
  215.         {
  216.         if( x > -1.E-4 )
  217.                 goto Small;
  218.         z *=x;
  219.         x += 1.0;
  220.         }
  221. */
  222. while( x < 2.0 )
  223.         {
  224.         if( x < 1.e-4 )
  225.                 goto Small;
  226.         z *=x;
  227.         x += 1.0;
  228.         }
  229.  
  230. if( direction )
  231.         z = 1.0/z;
  232.  
  233. if( x == 2.0 )
  234.         return(z);
  235.  
  236. x -= 2.0;
  237. p = z * polevlf( x, P, 7 );
  238.  
  239. gdone:
  240.  
  241. if( negative )
  242.         {
  243.         p = *sgngamf * PIF/(nz * p );
  244.         }
  245. return(p);
  246.  
  247. Small:
  248. if( x == 0.0 )
  249.         {
  250.         goto gsing;
  251.         }
  252. else
  253.         {
  254.         p = z / ((1.0 + 0.5772156649015329 * x) * x);
  255.         goto gdone;
  256.         }
  257. }
  258.  
  259. /* This is the C99 version */
  260.  
  261. float tgammaf(float x)
  262. {
  263.   int local_sgngamf=0;
  264.   return (__tgammaf_r(x, &local_sgngamf));
  265. }
  266.