Subversion Repositories Kolibri OS

Rev

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

  1. /*                                                      gamma.c
  2.  *
  3.  *      Gamma function
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, __tgamma_r();
  10.  * int* sgngam;
  11.  * y = __tgamma_r( x, sgngam );
  12.  *
  13.  * double x, y, tgamma();
  14.  * y = tgamma( x)
  15.  *
  16.  *
  17.  *
  18.  * DESCRIPTION:
  19.  *
  20.  * Returns gamma function of the argument.  The result is
  21.  * correctly signed. In the reentrant version the sign (+1 or -1)
  22.  * is returned in the variable referenced by sgngam.
  23.  *
  24.  * Arguments |x| <= 34 are reduced by recurrence and the function
  25.  * approximated by a rational function of degree 6/7 in the
  26.  * interval (2,3).  Large arguments are handled by Stirling's
  27.  * formula. Large 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.  *    DEC      -34, 34      10000       1.3e-16     2.5e-17
  36.  *    IEEE    -170,-33      20000       2.3e-15     3.3e-16
  37.  *    IEEE     -33,  33     20000       9.4e-16     2.2e-16
  38.  *    IEEE      33, 171.6   20000       2.3e-15     3.2e-16
  39.  *
  40.  * Error for arguments outside the test range will be larger
  41.  * owing to error amplification by the exponential function.
  42.  *
  43.  */
  44.  
  45. /*
  46. Cephes Math Library Release 2.8:  June, 2000
  47. Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
  48. */
  49.  
  50.  
  51. /*
  52.  * 26-11-2002 Modified for mingw.
  53.  * Danny Smith <dannysmith@users.sourceforge.net>
  54.  */
  55.  
  56.  
  57. #ifndef __MINGW32__
  58. #include "mconf.h"
  59. #else
  60. #include "cephes_mconf.h"
  61. #endif
  62.  
  63. #ifdef UNK
  64. static const double P[] = {
  65.   1.60119522476751861407E-4,
  66.   1.19135147006586384913E-3,
  67.   1.04213797561761569935E-2,
  68.   4.76367800457137231464E-2,
  69.   2.07448227648435975150E-1,
  70.   4.94214826801497100753E-1,
  71.   9.99999999999999996796E-1
  72. };
  73. static const double Q[] = {
  74. -2.31581873324120129819E-5,
  75.  5.39605580493303397842E-4,
  76. -4.45641913851797240494E-3,
  77.  1.18139785222060435552E-2,
  78.  3.58236398605498653373E-2,
  79. -2.34591795718243348568E-1,
  80.  7.14304917030273074085E-2,
  81.  1.00000000000000000320E0
  82. };
  83. #define MAXGAM 171.624376956302725
  84. static const double LOGPI = 1.14472988584940017414;
  85. #endif
  86.  
  87. #ifdef DEC
  88. static const unsigned short P[] = {
  89. 0035047,0162701,0146301,0005234,
  90. 0035634,0023437,0032065,0176530,
  91. 0036452,0137157,0047330,0122574,
  92. 0037103,0017310,0143041,0017232,
  93. 0037524,0066516,0162563,0164605,
  94. 0037775,0004671,0146237,0014222,
  95. 0040200,0000000,0000000,0000000
  96. };
  97. static const unsigned short Q[] = {
  98. 0134302,0041724,0020006,0116565,
  99. 0035415,0072121,0044251,0025634,
  100. 0136222,0003447,0035205,0121114,
  101. 0036501,0107552,0154335,0104271,
  102. 0037022,0135717,0014776,0171471,
  103. 0137560,0034324,0165024,0037021,
  104. 0037222,0045046,0047151,0161213,
  105. 0040200,0000000,0000000,0000000
  106. };
  107. #define MAXGAM 34.84425627277176174
  108. #endif
  109.  
  110. #ifdef IBMPC
  111. static const unsigned short P[] = {
  112. 0x2153,0x3998,0xfcb8,0x3f24,
  113. 0xbfab,0xe686,0x84e3,0x3f53,
  114. 0x14b0,0xe9db,0x57cd,0x3f85,
  115. 0x23d3,0x18c4,0x63d9,0x3fa8,
  116. 0x7d31,0xdcae,0x8da9,0x3fca,
  117. 0xe312,0x3993,0xa137,0x3fdf,
  118. 0x0000,0x0000,0x0000,0x3ff0
  119. };
  120. static const unsigned short Q[] = {
  121. 0xd3af,0x8400,0x487a,0xbef8,
  122. 0x2573,0x2915,0xae8a,0x3f41,
  123. 0xb44a,0xe750,0x40e4,0xbf72,
  124. 0xb117,0x5b1b,0x31ed,0x3f88,
  125. 0xde67,0xe33f,0x5779,0x3fa2,
  126. 0x87c2,0x9d42,0x071a,0xbfce,
  127. 0x3c51,0xc9cd,0x4944,0x3fb2,
  128. 0x0000,0x0000,0x0000,0x3ff0
  129. };
  130. #define MAXGAM 171.624376956302725
  131. #endif
  132.  
  133. #ifdef MIEEE
  134. static const unsigned short P[] = {
  135. 0x3f24,0xfcb8,0x3998,0x2153,
  136. 0x3f53,0x84e3,0xe686,0xbfab,
  137. 0x3f85,0x57cd,0xe9db,0x14b0,
  138. 0x3fa8,0x63d9,0x18c4,0x23d3,
  139. 0x3fca,0x8da9,0xdcae,0x7d31,
  140. 0x3fdf,0xa137,0x3993,0xe312,
  141. 0x3ff0,0x0000,0x0000,0x0000
  142. };
  143. static const unsigned short Q[] = {
  144. 0xbef8,0x487a,0x8400,0xd3af,
  145. 0x3f41,0xae8a,0x2915,0x2573,
  146. 0xbf72,0x40e4,0xe750,0xb44a,
  147. 0x3f88,0x31ed,0x5b1b,0xb117,
  148. 0x3fa2,0x5779,0xe33f,0xde67,
  149. 0xbfce,0x071a,0x9d42,0x87c2,
  150. 0x3fb2,0x4944,0xc9cd,0x3c51,
  151. 0x3ff0,0x0000,0x0000,0x0000
  152. };
  153. #define MAXGAM 171.624376956302725
  154. #endif
  155.  
  156. /* Stirling's formula for the gamma function */
  157. #if UNK
  158. static const double STIR[5] = {
  159.  7.87311395793093628397E-4,
  160. -2.29549961613378126380E-4,
  161. -2.68132617805781232825E-3,
  162.  3.47222221605458667310E-3,
  163.  8.33333333333482257126E-2,
  164. };
  165. #define MAXSTIR 143.01608
  166. static const double SQTPI = 2.50662827463100050242E0;
  167. #endif
  168. #if DEC
  169. static const unsigned short STIR[20] = {
  170. 0035516,0061622,0144553,0112224,
  171. 0135160,0131531,0037460,0165740,
  172. 0136057,0134460,0037242,0077270,
  173. 0036143,0107070,0156306,0027751,
  174. 0037252,0125252,0125252,0146064,
  175. };
  176. #define MAXSTIR 26.77
  177. static const unsigned short SQT[4] = {
  178. 0040440,0066230,0177661,0034055,
  179. };
  180. #define SQTPI *(double *)SQT
  181. #endif
  182. #if IBMPC
  183. static const unsigned short STIR[20] = {
  184. 0x7293,0x592d,0xcc72,0x3f49,
  185. 0x1d7c,0x27e6,0x166b,0xbf2e,
  186. 0x4fd7,0x07d4,0xf726,0xbf65,
  187. 0xc5fd,0x1b98,0x71c7,0x3f6c,
  188. 0x5986,0x5555,0x5555,0x3fb5,
  189. };
  190. #define MAXSTIR 143.01608
  191.  
  192. static const union
  193. {
  194.   unsigned short s[4];
  195.   double d;
  196. } sqt = {{0x2706,0x1ff6,0x0d93,0x4004}};
  197. #define SQTPI (sqt.d)
  198. #endif
  199. #if MIEEE
  200. static const unsigned short STIR[20] = {
  201. 0x3f49,0xcc72,0x592d,0x7293,
  202. 0xbf2e,0x166b,0x27e6,0x1d7c,
  203. 0xbf65,0xf726,0x07d4,0x4fd7,
  204. 0x3f6c,0x71c7,0x1b98,0xc5fd,
  205. 0x3fb5,0x5555,0x5555,0x5986,
  206. };
  207. #define MAXSTIR 143.01608
  208. static const unsigned short SQT[4] = {
  209. 0x4004,0x0d93,0x1ff6,0x2706,
  210. };
  211. #define SQTPI *(double *)SQT
  212. #endif
  213.  
  214. #ifndef __MINGW32__
  215. int sgngam = 0;
  216. extern int sgngam;
  217. extern double MAXLOG, MAXNUM, PI;
  218. #ifdef ANSIPROT
  219. extern double pow ( double, double );
  220. extern double log ( double );
  221. extern double exp ( double );
  222. extern double sin ( double );
  223. extern double polevl ( double, void *, int );
  224. extern double p1evl ( double, void *, int );
  225. extern double floor ( double );
  226. extern double fabs ( double );
  227. extern int isnan ( double );
  228. extern int isfinite ( double );
  229. static double stirf ( double );
  230. double lgam ( double );
  231. #else
  232. double pow(), log(), exp(), sin(), polevl(), p1evl(), floor(), fabs();
  233. int isnan(), isfinite();
  234. static double stirf();
  235. double lgam();
  236. #endif
  237. #ifdef INFINITIES
  238. extern double INFINITY;
  239. #endif
  240. #ifdef NANS
  241. extern double NAN;
  242. #endif
  243. #else /* __MINGW32__ */
  244. static double stirf ( double );
  245. #endif
  246.  
  247. /* Gamma function computed by Stirling's formula.
  248.  * The polynomial STIR is valid for 33 <= x <= 172.
  249.  */
  250. static double stirf(x)
  251. double x;
  252. {
  253. double y, w, v;
  254.  
  255. w = 1.0/x;
  256. w = 1.0 + w * polevl( w, STIR, 4 );
  257. y = exp(x);
  258. if( x > MAXSTIR )
  259.         { /* Avoid overflow in pow() */
  260.         v = pow( x, 0.5 * x - 0.25 );
  261.         y = v * (v / y);
  262.         }
  263. else
  264.         {
  265.         y = pow( x, x - 0.5 ) / y;
  266.         }
  267. y = SQTPI * y * w;
  268. return( y );
  269. }
  270.  
  271.  
  272.  
  273. double __tgamma_r(double x, int* sgngam)
  274. {
  275. double p, q, z;
  276. int i;
  277.  
  278. *sgngam = 1;
  279. #ifdef NANS
  280. if( isnan(x) )
  281.         return(x);
  282. #endif
  283. #ifdef INFINITIES
  284. #ifdef NANS
  285. if( x == INFINITY )
  286.         return(x);
  287. if( x == -INFINITY )
  288.         return(NAN);
  289. #else
  290. if( !isfinite(x) )
  291.         return(x);
  292. #endif
  293. #endif
  294. q = fabs(x);
  295.  
  296. if( q > 33.0 )
  297.         {
  298.         if( x < 0.0 )
  299.                 {
  300.                 p = floor(q);
  301.                 if( p == q )
  302.                         {
  303. gsing:
  304.                         _SET_ERRNO(EDOM);
  305.                         mtherr( "tgamma", SING );
  306. #ifdef INFINITIES
  307.                         return (INFINITY);
  308. #else
  309.                         return (MAXNUM);
  310. #endif
  311.                         }
  312.                 i = p;
  313.                 if( (i & 1) == 0 )
  314.                         *sgngam = -1;
  315.                 z = q - p;
  316.                 if( z > 0.5 )
  317.                         {
  318.                         p += 1.0;
  319.                         z = q - p;
  320.                         }
  321.                 z = q * sin( PI * z );
  322.                 if( z == 0.0 )
  323.                         {
  324.                         _SET_ERRNO(ERANGE);
  325.                         mtherr( "tgamma", OVERFLOW );
  326. #ifdef INFINITIES
  327.                         return( *sgngam * INFINITY);
  328. #else
  329.                         return( *sgngam * MAXNUM);
  330. #endif
  331.                         }
  332.                 z = fabs(z);
  333.                 z = PI/(z * stirf(q) );
  334.                 }
  335.         else
  336.                 {
  337.                 z = stirf(x);
  338.                 }
  339.         return( *sgngam * z );
  340.         }
  341.  
  342. z = 1.0;
  343. while( x >= 3.0 )
  344.         {
  345.         x -= 1.0;
  346.         z *= x;
  347.         }
  348.  
  349. while( x < 0.0 )
  350.         {
  351.         if( x > -1.E-9 )
  352.                 goto Small;
  353.         z /= x;
  354.         x += 1.0;
  355.         }
  356.  
  357. while( x < 2.0 )
  358.         {
  359.         if( x < 1.e-9 )
  360.                 goto Small;
  361.         z /= x;
  362.         x += 1.0;
  363.         }
  364.  
  365. if( x == 2.0 )
  366.         return(z);
  367.  
  368. x -= 2.0;
  369. p = polevl( x, P, 6 );
  370. q = polevl( x, Q, 7 );
  371. return( z * p / q );
  372.  
  373. Small:
  374. if( x == 0.0 )
  375.         {
  376.         goto gsing;
  377.         }
  378. else
  379.         return( z/((1.0 + 0.5772156649015329 * x) * x) );
  380. }
  381.  
  382. /* This is the C99 version */
  383.  
  384. double tgamma(double x)
  385. {
  386.   int local_sgngam=0;
  387.   return (__tgamma_r(x, &local_sgngam));
  388. }
  389.