home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / cprob / gamma.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  12.4 KB  |  608 lines
  1. /*                            gamma.c
  2.  *
  3.  *    Gamma function
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, gamma();
  10.  * extern int sgngam;
  11.  *
  12.  * y = gamma( 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 sgngam.
  21.  * This variable is also filled in by the logarithmic gamma
  22.  * function lgam().
  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. /*                            lgam()
  45.  *
  46.  *    Natural logarithm of gamma function
  47.  *
  48.  *
  49.  *
  50.  * SYNOPSIS:
  51.  *
  52.  * double x, y, lgam();
  53.  * extern int sgngam;
  54.  *
  55.  * y = lgam( x );
  56.  *
  57.  *
  58.  *
  59.  * DESCRIPTION:
  60.  *
  61.  * Returns the base e (2.718...) logarithm of the absolute
  62.  * value of the gamma function of the argument.
  63.  * The sign (+1 or -1) of the gamma function is returned in a
  64.  * global (extern) variable named sgngam.
  65.  *
  66.  * For arguments greater than 13, the logarithm of the gamma
  67.  * function is approximated by the logarithmic version of
  68.  * Stirling's formula using a polynomial approximation of
  69.  * degree 4. Arguments between -33 and +33 are reduced by
  70.  * recurrence to the interval [2,3] of a rational approximation.
  71.  * The cosecant reflection formula is employed for arguments
  72.  * less than -33.
  73.  *
  74.  * Arguments greater than MAXLGM return MAXNUM and an error
  75.  * message.  MAXLGM = 2.035093e36 for DEC
  76.  * arithmetic or 2.556348e305 for IEEE arithmetic.
  77.  *
  78.  *
  79.  *
  80.  * ACCURACY:
  81.  *
  82.  *
  83.  * arithmetic      domain        # trials     peak         rms
  84.  *    DEC     0, 3                  7000     5.2e-17     1.3e-17
  85.  *    DEC     2.718, 2.035e36       5000     3.9e-17     9.9e-18
  86.  *    IEEE    0, 3                 28000     5.4e-16     1.1e-16
  87.  *    IEEE    2.718, 2.556e305     40000     3.5e-16     8.3e-17
  88.  * The error criterion was relative when the function magnitude
  89.  * was greater than one but absolute when it was less than one.
  90.  *
  91.  * The following test used the relative error criterion, though
  92.  * at certain points the relative error could be much higher than
  93.  * indicated.
  94.  *    IEEE    -200, -4             10000     4.8e-16     1.3e-16
  95.  *
  96.  */
  97.  
  98. /*                            gamma.c    */
  99. /*    gamma function    */
  100.  
  101. /*
  102. Cephes Math Library Release 2.2:  July, 1992
  103. Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
  104. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  105. */
  106.  
  107.  
  108. #include "mconf.h"
  109.  
  110. #ifdef UNK
  111. static double P[] = {
  112.   1.60119522476751861407E-4,
  113.   1.19135147006586384913E-3,
  114.   1.04213797561761569935E-2,
  115.   4.76367800457137231464E-2,
  116.   2.07448227648435975150E-1,
  117.   4.94214826801497100753E-1,
  118.   9.99999999999999996796E-1
  119. };
  120. static double Q[] = {
  121. -2.31581873324120129819E-5,
  122.  5.39605580493303397842E-4,
  123. -4.45641913851797240494E-3,
  124.  1.18139785222060435552E-2,
  125.  3.58236398605498653373E-2,
  126. -2.34591795718243348568E-1,
  127.  7.14304917030273074085E-2,
  128.  1.00000000000000000320E0
  129. };
  130. #define MAXGAM 171.624376956302725
  131. static double LOGPI = 1.14472988584940017414;
  132. #endif
  133.  
  134. #ifdef DEC
  135. static short P[] = {
  136. 0035047,0162701,0146301,0005234,
  137. 0035634,0023437,0032065,0176530,
  138. 0036452,0137157,0047330,0122574,
  139. 0037103,0017310,0143041,0017232,
  140. 0037524,0066516,0162563,0164605,
  141. 0037775,0004671,0146237,0014222,
  142. 0040200,0000000,0000000,0000000
  143. };
  144. static short Q[] = {
  145. 0134302,0041724,0020006,0116565,
  146. 0035415,0072121,0044251,0025634,
  147. 0136222,0003447,0035205,0121114,
  148. 0036501,0107552,0154335,0104271,
  149. 0037022,0135717,0014776,0171471,
  150. 0137560,0034324,0165024,0037021,
  151. 0037222,0045046,0047151,0161213,
  152. 0040200,0000000,0000000,0000000
  153. };
  154. #define MAXGAM 34.84425627277176174
  155. static short LPI[4] = {
  156. 0040222,0103202,0043475,0006750,
  157. };
  158. #define LOGPI *(double *)LPI
  159. #endif
  160.  
  161. #ifdef IBMPC
  162. static short P[] = {
  163. 0x2153,0x3998,0xfcb8,0x3f24,
  164. 0xbfab,0xe686,0x84e3,0x3f53,
  165. 0x14b0,0xe9db,0x57cd,0x3f85,
  166. 0x23d3,0x18c4,0x63d9,0x3fa8,
  167. 0x7d31,0xdcae,0x8da9,0x3fca,
  168. 0xe312,0x3993,0xa137,0x3fdf,
  169. 0x0000,0x0000,0x0000,0x3ff0
  170. };
  171. static short Q[] = {
  172. 0xd3af,0x8400,0x487a,0xbef8,
  173. 0x2573,0x2915,0xae8a,0x3f41,
  174. 0xb44a,0xe750,0x40e4,0xbf72,
  175. 0xb117,0x5b1b,0x31ed,0x3f88,
  176. 0xde67,0xe33f,0x5779,0x3fa2,
  177. 0x87c2,0x9d42,0x071a,0xbfce,
  178. 0x3c51,0xc9cd,0x4944,0x3fb2,
  179. 0x0000,0x0000,0x0000,0x3ff0
  180. };
  181. #define MAXGAM 171.624376956302725
  182. static short LPI[4] = {
  183. 0xa1bd,0x48e7,0x50d0,0x3ff2,
  184. };
  185. #define LOGPI *(double *)LPI
  186. #endif 
  187.  
  188. #ifdef MIEEE
  189. static short P[] = {
  190. 0x3f24,0xfcb8,0x3998,0x2153,
  191. 0x3f53,0x84e3,0xe686,0xbfab,
  192. 0x3f85,0x57cd,0xe9db,0x14b0,
  193. 0x3fa8,0x63d9,0x18c4,0x23d3,
  194. 0x3fca,0x8da9,0xdcae,0x7d31,
  195. 0x3fdf,0xa137,0x3993,0xe312,
  196. 0x3ff0,0x0000,0x0000,0x0000
  197. };
  198. static short Q[] = {
  199. 0xbef8,0x487a,0x8400,0xd3af,
  200. 0x3f41,0xae8a,0x2915,0x2573,
  201. 0xbf72,0x40e4,0xe750,0xb44a,
  202. 0x3f88,0x31ed,0x5b1b,0xb117,
  203. 0x3fa2,0x5779,0xe33f,0xde67,
  204. 0xbfce,0x071a,0x9d42,0x87c2,
  205. 0x3fb2,0x4944,0xc9cd,0x3c51,
  206. 0x3ff0,0x0000,0x0000,0x0000
  207. };
  208. #define MAXGAM 171.624376956302725
  209. static short LPI[4] = {
  210. 0x3ff2,0x50d0,0x48e7,0xa1bd,
  211. };
  212. #define LOGPI *(double *)LPI
  213. #endif 
  214.  
  215. /* Stirling's formula for the gamma function */
  216. #if UNK
  217. static double STIR[5] = {
  218.  7.87311395793093628397E-4,
  219. -2.29549961613378126380E-4,
  220. -2.68132617805781232825E-3,
  221.  3.47222221605458667310E-3,
  222.  8.33333333333482257126E-2,
  223. };
  224. #define MAXSTIR 143.01608
  225. static double SQTPI = 2.50662827463100050242E0;
  226. #endif
  227. #if DEC
  228. static short STIR[20] = {
  229. 0035516,0061622,0144553,0112224,
  230. 0135160,0131531,0037460,0165740,
  231. 0136057,0134460,0037242,0077270,
  232. 0036143,0107070,0156306,0027751,
  233. 0037252,0125252,0125252,0146064,
  234. };
  235. #define MAXSTIR 26.77
  236. static short SQT[4] = {
  237. 0040440,0066230,0177661,0034055,
  238. };
  239. #define SQTPI *(double *)SQT
  240. #endif
  241. #if IBMPC
  242. static short STIR[20] = {
  243. 0x7293,0x592d,0xcc72,0x3f49,
  244. 0x1d7c,0x27e6,0x166b,0xbf2e,
  245. 0x4fd7,0x07d4,0xf726,0xbf65,
  246. 0xc5fd,0x1b98,0x71c7,0x3f6c,
  247. 0x5986,0x5555,0x5555,0x3fb5,
  248. };
  249. #define MAXSTIR 143.01608
  250. static short SQT[4] = {
  251. 0x2706,0x1ff6,0x0d93,0x4004,
  252. };
  253. #define SQTPI *(double *)SQT
  254. #endif
  255. #if MIEEE
  256. static short STIR[20] = {
  257. 0x3f49,0xcc72,0x592d,0x7293,
  258. 0xbf2e,0x166b,0x27e6,0x1d7c,
  259. 0xbf65,0xf726,0x07d4,0x4fd7,
  260. 0x3f6c,0x71c7,0x1b98,0xc5fd,
  261. 0x3fb5,0x5555,0x5555,0x5986,
  262. };
  263. #define MAXSTIR 143.01608
  264. static short SQT[4] = {
  265. 0x4004,0x0d93,0x1ff6,0x2706,
  266. };
  267. #define SQTPI *(double *)SQT
  268. #endif
  269.  
  270. int sgngam = 0;
  271. extern int sgngam;
  272. extern double MAXLOG, MAXNUM, PI;
  273.  
  274.  
  275. /* Gamma function computed by Stirling's formula.
  276.  * The polynomial STIR is valid for 33 <= x <= 172.
  277.  */
  278. static double stirf(x)
  279. double x;
  280. {
  281. double y, w, v;
  282. double pow(), exp(), polevl();
  283. double log();
  284.  
  285. w = 1.0/x;
  286. w = 1.0 + w * polevl( w, STIR, 4 );
  287. y = exp(x);
  288. if( x > MAXSTIR )
  289.     { /* Avoid overflow in pow() */
  290.     v = pow( x, 0.5 * x - 0.25 );
  291.     y = v * (v / y);
  292.     }
  293. else
  294.     {
  295.     y = pow( x, x - 0.5 ) / y;
  296.     }
  297. y = SQTPI * y * w;
  298. return( y );
  299. }
  300.  
  301.  
  302.  
  303. double gamma(x)
  304. double x;
  305. {
  306. double p, q, z;
  307. int i;
  308. double fabs(), lgam(), log(), exp(), gamma(), sin();
  309. double floor(), polevl(), p1evl(), stirf();
  310.  
  311. sgngam = 1;
  312. q = fabs(x);
  313.  
  314. if( q > 33.0 )
  315.     {
  316.     if( x < 0.0 )
  317.         {
  318.         p = floor(q);
  319.         if( p == q )
  320.             goto goverf;
  321.         i = p;
  322.         if( (i & 1) == 0 )
  323.             sgngam = -1;
  324.         z = q - p;
  325.         if( z > 0.5 )
  326.             {
  327.             p += 1.0;
  328.             z = q - p;
  329.             }
  330.         z = q * sin( PI * z );
  331.         if( z == 0.0 )
  332.             {
  333. goverf:
  334.             mtherr( "gamma", OVERFLOW );
  335.             return( sgngam * MAXNUM);
  336.             }
  337.         z = fabs(z);
  338.         z = PI/(z * stirf(q) );
  339.         }
  340.     else
  341.         {
  342.         z = stirf(x);
  343.         }
  344.     return( sgngam * z );
  345.     }
  346.  
  347. z = 1.0;
  348. while( x >= 3.0 )
  349.     {
  350.     x -= 1.0;
  351.     z *= x;
  352.     }
  353.  
  354. while( x < 0.0 )
  355.     {
  356.     if( x > -1.E-9 )
  357.         goto small;
  358.     z /=x;
  359.     x += 1.0;
  360.     }
  361.  
  362. while( x < 2.0 )
  363.     {
  364.     if( x < 1.e-9 )
  365.         goto small;
  366.     z /=x;
  367.     x += 1.0;
  368.     }
  369.  
  370. if( (x == 2.0) || (x == 3.0) )
  371.     return(z);
  372.  
  373. x -= 2.0;
  374. p = polevl( x, P, 6 );
  375. q = polevl( x, Q, 7 );
  376. return( z * p / q );
  377.  
  378. small:
  379. if( x == 0.0 )
  380.     {
  381.     mtherr( "gamma", SING );
  382.     return( MAXNUM );
  383.     }
  384. else
  385.     return( z/((1.0 + 0.5772156649015329 * x) * x) );
  386. }
  387.  
  388.  
  389.  
  390. /* A[]: Stirling's formula expansion of log gamma
  391.  * B[], C[]: log gamma function between 2 and 3
  392.  */
  393. #ifdef UNK
  394. static double A[] = {
  395.  8.11614167470508450300E-4,
  396. -5.95061904284301438324E-4,
  397.  7.93650340457716943945E-4,
  398. -2.77777777730099687205E-3,
  399.  8.33333333333331927722E-2
  400. };
  401. static double B[] = {
  402. -1.37825152569120859100E3,
  403. -3.88016315134637840924E4,
  404. -3.31612992738871184744E5,
  405. -1.16237097492762307383E6,
  406. -1.72173700820839662146E6,
  407. -8.53555664245765465627E5
  408. };
  409. static double C[] = {
  410.  1.00000000000000000000E0,
  411. -3.51815701436523470549E2,
  412. -1.70642106651881159223E4,
  413. -2.20528590553854454839E5,
  414. -1.13933444367982507207E6,
  415. -2.53252307177582951285E6,
  416. -2.01889141433532773231E6
  417. };
  418. /* log( sqrt( 2*pi ) ) */
  419. static double LS2PI  =  0.91893853320467274178;
  420. #define MAXLGM 2.556348e305
  421. #endif
  422.  
  423. #ifdef DEC
  424. static short A[] = {
  425. 0035524,0141201,0034633,0031405,
  426. 0135433,0176755,0126007,0045030,
  427. 0035520,0006371,0003342,0172730,
  428. 0136066,0005540,0132605,0026407,
  429. 0037252,0125252,0125252,0125132
  430. };
  431. static short B[] = {
  432. 0142654,0044014,0077633,0035410,
  433. 0144027,0110641,0125335,0144760,
  434. 0144641,0165637,0142204,0047447,
  435. 0145215,0162027,0146246,0155211,
  436. 0145322,0026110,0010317,0110130,
  437. 0145120,0061472,0120300,0025363
  438. };
  439. static short C[] = {
  440. /*0040200,0000000,0000000,0000000*/
  441. 0142257,0164150,0163630,0112622,
  442. 0143605,0050153,0156116,0135272,
  443. 0144527,0056045,0145642,0062332,
  444. 0145213,0012063,0106250,0001025,
  445. 0145432,0111254,0044577,0115142,
  446. 0145366,0071133,0050217,0005122
  447. };
  448. /* log( sqrt( 2*pi ) ) */
  449. static short LS2P[] = {040153,037616,041445,0172645,};
  450. #define LS2PI *(double *)LS2P
  451. #define MAXLGM 2.035093e36
  452. #endif
  453.  
  454. #ifdef IBMPC
  455. static short A[] = {
  456. 0x6661,0x2733,0x9850,0x3f4a,
  457. 0xe943,0xb580,0x7fbd,0xbf43,
  458. 0x5ebb,0x20dc,0x019f,0x3f4a,
  459. 0xa5a1,0x16b0,0xc16c,0xbf66,
  460. 0x554b,0x5555,0x5555,0x3fb5
  461. };
  462. static short B[] = {
  463. 0x6761,0x8ff3,0x8901,0xc095,
  464. 0xb93e,0x355b,0xf234,0xc0e2,
  465. 0x89e5,0xf890,0x3d73,0xc114,
  466. 0xdb51,0xf994,0xbc82,0xc131,
  467. 0xf20b,0x0219,0x4589,0xc13a,
  468. 0x055e,0x5418,0x0c67,0xc12a
  469. };
  470. static short C[] = {
  471. /*0x0000,0x0000,0x0000,0x3ff0,*/
  472. 0x12b2,0x1cf3,0xfd0d,0xc075,
  473. 0xd757,0x7b89,0xaa0d,0xc0d0,
  474. 0x4c9b,0xb974,0xeb84,0xc10a,
  475. 0x0043,0x7195,0x6286,0xc131,
  476. 0xf34c,0x892f,0x5255,0xc143,
  477. 0xe14a,0x6a11,0xce4b,0xc13e
  478. };
  479. /* log( sqrt( 2*pi ) ) */
  480. static short LS2P[] = {
  481. 0xbeb5,0xc864,0x67f1,0x3fed
  482. };
  483. #define LS2PI *(double *)LS2P
  484. #define MAXLGM 2.556348e305
  485. #endif
  486.  
  487. #ifdef MIEEE
  488. static short A[] = {
  489. 0x3f4a,0x9850,0x2733,0x6661,
  490. 0xbf43,0x7fbd,0xb580,0xe943,
  491. 0x3f4a,0x019f,0x20dc,0x5ebb,
  492. 0xbf66,0xc16c,0x16b0,0xa5a1,
  493. 0x3fb5,0x5555,0x5555,0x554b
  494. };
  495. static short B[] = {
  496. 0xc095,0x8901,0x8ff3,0x6761,
  497. 0xc0e2,0xf234,0x355b,0xb93e,
  498. 0xc114,0x3d73,0xf890,0x89e5,
  499. 0xc131,0xbc82,0xf994,0xdb51,
  500. 0xc13a,0x4589,0x0219,0xf20b,
  501. 0xc12a,0x0c67,0x5418,0x055e
  502. };
  503. static short C[] = {
  504. 0xc075,0xfd0d,0x1cf3,0x12b2,
  505. 0xc0d0,0xaa0d,0x7b89,0xd757,
  506. 0xc10a,0xeb84,0xb974,0x4c9b,
  507. 0xc131,0x6286,0x7195,0x0043,
  508. 0xc143,0x5255,0x892f,0xf34c,
  509. 0xc13e,0xce4b,0x6a11,0xe14a
  510. };
  511. /* log( sqrt( 2*pi ) ) */
  512. static short LS2P[] = {
  513. 0x3fed,0x67f1,0xc864,0xbeb5
  514. };
  515. #define LS2PI *(double *)LS2P
  516. #define MAXLGM 2.556348e305
  517. #endif
  518.  
  519.  
  520. /* Logarithm of gamma function */
  521.  
  522.  
  523. double lgam(x)
  524. double x;
  525. {
  526. double p, q, w, z;
  527. int i;
  528. double fabs(), sin(), log(), gamma(), polevl();
  529. double p1evl(), floor();
  530.  
  531. sgngam = 1;
  532.  
  533. if( x < -34.0 )
  534.     {
  535.     q = -x;
  536.     w = lgam(q); /* note this modifies sgngam! */
  537.     p = floor(q);
  538.     if( p == q )
  539.         goto loverf;
  540.     i = p;
  541.     if( (i & 1) == 0 )
  542.         sgngam = -1;
  543.     else
  544.         sgngam = 1;
  545.     z = q - p;
  546.     if( z > 0.5 )
  547.         {
  548.         p += 1.0;
  549.         z = p - q;
  550.         }
  551.     z = q * sin( PI * z );
  552.     if( z == 0.0 )
  553.         goto loverf;
  554. /*    z = log(PI) - log( z ) - w;*/
  555.     z = LOGPI - log( z ) - w;
  556.     return( z );
  557.     }
  558.  
  559. if( x < 13.0 )
  560.     {
  561.     z = 1.0;
  562.     while( x >= 3.0 )
  563.         {
  564.         x -= 1.0;
  565.         z *= x;
  566.         }
  567.     while( x < 2.0 )
  568.         {
  569.         if( x == 0.0 )
  570.             goto loverf;
  571.         z /=x;
  572.         x += 1.0;
  573.         }
  574.     if( z < 0.0 )
  575.         {
  576.         sgngam = -1;
  577.         z = -z;
  578.         }
  579.     else
  580.         sgngam = 1;
  581.     if( x == 2.0 )
  582.         return( log(z) );
  583.     x -= 2.0;
  584.     p = x * polevl( x, B, 5 ) / p1evl( x, C, 6);
  585.     return( log(z) + p );
  586.     }
  587.  
  588. if( x > MAXLGM )
  589.     {
  590. loverf:
  591.     mtherr( "lgam", OVERFLOW );
  592.     return( sgngam * MAXNUM );
  593.     }
  594.  
  595. q = ( x - 0.5 ) * log(x) - x + LS2PI;
  596. if( x > 1.0e8 )
  597.     return( q );
  598.  
  599. p = 1.0/(x*x);
  600. if( x >= 1000.0 )
  601.     q += ((   7.9365079365079365079365e-4 * p
  602.         - 2.7777777777777777777778e-3) *p
  603.         + 0.0833333333333333333333) / x;
  604. else
  605.     q += polevl( p, A, 4 ) / x;
  606. return( q );
  607. }
  608.