Subversion Repositories Kolibri OS

Rev

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

  1. /*                                                      lgam()
  2.  *
  3.  *      Natural logarithm of gamma function
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, __lgamma_r();
  10.  * int* sgngam;
  11.  * y = __lgamma_r( x, sgngam );
  12.  *
  13.  * double x, y, lgamma();
  14.  * y = lgamma( x);
  15.  *
  16.  *
  17.  *
  18.  * DESCRIPTION:
  19.  *
  20.  * Returns the base e (2.718...) logarithm of the absolute
  21.  * value of the gamma function of the argument. In the reentrant
  22.  * version, the sign (+1 or -1) of the gamma function is returned
  23.  * in the variable referenced by sgngam.
  24.  *
  25.  * For arguments greater than 13, the logarithm of the gamma
  26.  * function is approximated by the logarithmic version of
  27.  * Stirling's formula using a polynomial approximation of
  28.  * degree 4. Arguments between -33 and +33 are reduced by
  29.  * recurrence to the interval [2,3] of a rational approximation.
  30.  * The cosecant reflection formula is employed for arguments
  31.  * less than -33.
  32.  *
  33.  * Arguments greater than MAXLGM return MAXNUM and an error
  34.  * message.  MAXLGM = 2.035093e36 for DEC
  35.  * arithmetic or 2.556348e305 for IEEE arithmetic.
  36.  *
  37.  *
  38.  *
  39.  * ACCURACY:
  40.  *
  41.  *
  42.  * arithmetic      domain        # trials     peak         rms
  43.  *    DEC     0, 3                  7000     5.2e-17     1.3e-17
  44.  *    DEC     2.718, 2.035e36       5000     3.9e-17     9.9e-18
  45.  *    IEEE    0, 3                 28000     5.4e-16     1.1e-16
  46.  *    IEEE    2.718, 2.556e305     40000     3.5e-16     8.3e-17
  47.  * The error criterion was relative when the function magnitude
  48.  * was greater than one but absolute when it was less than one.
  49.  *
  50.  * The following test used the relative error criterion, though
  51.  * at certain points the relative error could be much higher than
  52.  * indicated.
  53.  *    IEEE    -200, -4             10000     4.8e-16     1.3e-16
  54.  *
  55.  */
  56.  
  57. /*
  58.  * Cephes Math Library Release 2.8:  June, 2000
  59.  * Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
  60.  */
  61.  
  62. /*
  63.  * 26-11-2002 Modified for mingw.
  64.  * Danny Smith <dannysmith@users.sourceforge.net>
  65.  */
  66.  
  67.  
  68. #ifndef __MINGW32__
  69. #include "mconf.h"
  70. #ifdef ANSIPROT
  71. extern double pow ( double, double );
  72. extern double log ( double );
  73. extern double exp ( double );
  74. extern double sin ( double );
  75. extern double polevl ( double, void *, int );
  76. extern double p1evl ( double, void *, int );
  77. extern double floor ( double );
  78. extern double fabs ( double );
  79. extern int isnan ( double );
  80. extern int isfinite ( double );
  81. #else
  82. double pow(), log(), exp(), sin(), polevl(), p1evl(), floor(), fabs();
  83. int isnan(), isfinite();
  84. #endif
  85. #ifdef INFINITIES
  86. extern double INFINITY;
  87. #endif
  88. #ifdef NANS
  89. extern double NAN;
  90. #endif
  91. #else /* __MINGW32__ */
  92. #include "cephes_mconf.h"
  93. #endif /* __MINGW32__ */
  94.  
  95.  
  96. /* A[]: Stirling's formula expansion of log gamma
  97.  * B[], C[]: log gamma function between 2 and 3
  98.  */
  99. #ifdef UNK
  100. static double A[] = {
  101.  8.11614167470508450300E-4,
  102. -5.95061904284301438324E-4,
  103.  7.93650340457716943945E-4,
  104. -2.77777777730099687205E-3,
  105.  8.33333333333331927722E-2
  106. };
  107. static double B[] = {
  108. -1.37825152569120859100E3,
  109. -3.88016315134637840924E4,
  110. -3.31612992738871184744E5,
  111. -1.16237097492762307383E6,
  112. -1.72173700820839662146E6,
  113. -8.53555664245765465627E5
  114. };
  115. static double C[] = {
  116. /* 1.00000000000000000000E0, */
  117. -3.51815701436523470549E2,
  118. -1.70642106651881159223E4,
  119. -2.20528590553854454839E5,
  120. -1.13933444367982507207E6,
  121. -2.53252307177582951285E6,
  122. -2.01889141433532773231E6
  123. };
  124. /* log( sqrt( 2*pi ) ) */
  125. static double LS2PI  =  0.91893853320467274178;
  126. #define MAXLGM 2.556348e305
  127. static double LOGPI = 1.14472988584940017414;
  128. #endif
  129.  
  130. #ifdef DEC
  131. static const unsigned short A[] = {
  132. 0035524,0141201,0034633,0031405,
  133. 0135433,0176755,0126007,0045030,
  134. 0035520,0006371,0003342,0172730,
  135. 0136066,0005540,0132605,0026407,
  136. 0037252,0125252,0125252,0125132
  137. };
  138. static const unsigned short B[] = {
  139. 0142654,0044014,0077633,0035410,
  140. 0144027,0110641,0125335,0144760,
  141. 0144641,0165637,0142204,0047447,
  142. 0145215,0162027,0146246,0155211,
  143. 0145322,0026110,0010317,0110130,
  144. 0145120,0061472,0120300,0025363
  145. };
  146. static const unsigned short C[] = {
  147. /*0040200,0000000,0000000,0000000*/
  148. 0142257,0164150,0163630,0112622,
  149. 0143605,0050153,0156116,0135272,
  150. 0144527,0056045,0145642,0062332,
  151. 0145213,0012063,0106250,0001025,
  152. 0145432,0111254,0044577,0115142,
  153. 0145366,0071133,0050217,0005122
  154. };
  155. /* log( sqrt( 2*pi ) ) */
  156. static const unsigned short LS2P[] = {040153,037616,041445,0172645,};
  157. #define LS2PI *(double *)LS2P
  158. #define MAXLGM 2.035093e36
  159. static const unsigned short LPI[4] = {
  160. 0040222,0103202,0043475,0006750,
  161. };
  162. #define LOGPI *(double *)LPI
  163.  
  164. #endif
  165.  
  166. #ifdef IBMPC
  167. static const unsigned short A[] = {
  168. 0x6661,0x2733,0x9850,0x3f4a,
  169. 0xe943,0xb580,0x7fbd,0xbf43,
  170. 0x5ebb,0x20dc,0x019f,0x3f4a,
  171. 0xa5a1,0x16b0,0xc16c,0xbf66,
  172. 0x554b,0x5555,0x5555,0x3fb5
  173. };
  174. static const unsigned short B[] = {
  175. 0x6761,0x8ff3,0x8901,0xc095,
  176. 0xb93e,0x355b,0xf234,0xc0e2,
  177. 0x89e5,0xf890,0x3d73,0xc114,
  178. 0xdb51,0xf994,0xbc82,0xc131,
  179. 0xf20b,0x0219,0x4589,0xc13a,
  180. 0x055e,0x5418,0x0c67,0xc12a
  181. };
  182. static const unsigned short C[] = {
  183. /*0x0000,0x0000,0x0000,0x3ff0,*/
  184. 0x12b2,0x1cf3,0xfd0d,0xc075,
  185. 0xd757,0x7b89,0xaa0d,0xc0d0,
  186. 0x4c9b,0xb974,0xeb84,0xc10a,
  187. 0x0043,0x7195,0x6286,0xc131,
  188. 0xf34c,0x892f,0x5255,0xc143,
  189. 0xe14a,0x6a11,0xce4b,0xc13e
  190. };
  191. /* log( sqrt( 2*pi ) ) */
  192. static const union
  193. {
  194.   unsigned short  s[4];
  195.   double d;
  196. } ls2p  =  {{0xbeb5,0xc864,0x67f1,0x3fed}};
  197. #define LS2PI   (ls2p.d)
  198. #define MAXLGM 2.556348e305
  199. /* log (pi) */
  200. static const union
  201. {
  202.   unsigned short s[4];
  203.   double d;
  204. } lpi  =  {{0xa1bd,0x48e7,0x50d0,0x3ff2}};
  205. #define LOGPI (lpi.d)
  206. #endif
  207.  
  208. #ifdef MIEEE
  209. static const unsigned short A[] = {
  210. 0x3f4a,0x9850,0x2733,0x6661,
  211. 0xbf43,0x7fbd,0xb580,0xe943,
  212. 0x3f4a,0x019f,0x20dc,0x5ebb,
  213. 0xbf66,0xc16c,0x16b0,0xa5a1,
  214. 0x3fb5,0x5555,0x5555,0x554b
  215. };
  216. static const unsigned short B[] = {
  217. 0xc095,0x8901,0x8ff3,0x6761,
  218. 0xc0e2,0xf234,0x355b,0xb93e,
  219. 0xc114,0x3d73,0xf890,0x89e5,
  220. 0xc131,0xbc82,0xf994,0xdb51,
  221. 0xc13a,0x4589,0x0219,0xf20b,
  222. 0xc12a,0x0c67,0x5418,0x055e
  223. };
  224. static const unsigned short C[] = {
  225. 0xc075,0xfd0d,0x1cf3,0x12b2,
  226. 0xc0d0,0xaa0d,0x7b89,0xd757,
  227. 0xc10a,0xeb84,0xb974,0x4c9b,
  228. 0xc131,0x6286,0x7195,0x0043,
  229. 0xc143,0x5255,0x892f,0xf34c,
  230. 0xc13e,0xce4b,0x6a11,0xe14a
  231. };
  232. /* log( sqrt( 2*pi ) ) */
  233. static const union
  234. {
  235.   unsigned short  s[4];
  236.   double d;
  237. } ls2p  =  {{0x3fed,0x67f1,0xc864,0xbeb5}};
  238. #define LS2PI  ls2p.d
  239. #define MAXLGM 2.556348e305
  240. /* log (pi) */
  241. static const union
  242. {
  243.   unsigned short s[4];
  244.   double d;
  245. } lpi  =  {{0x3ff2, 0x50d0, 0x48e7, 0xa1bd}};
  246. #define LOGPI (lpi.d)
  247. #endif
  248.  
  249.  
  250. /* Logarithm of gamma function */
  251. /* Reentrant version */
  252.  
  253. double __lgamma_r(double x, int* sgngam)
  254. {
  255. double p, q, u, w, z;
  256. int i;
  257.  
  258. *sgngam = 1;
  259. #ifdef NANS
  260. if( isnan(x) )
  261.         return(x);
  262. #endif
  263.  
  264. #ifdef INFINITIES
  265. if( !isfinite(x) )
  266.         return(INFINITY);
  267. #endif
  268.  
  269. if( x < -34.0 )
  270.         {
  271.         q = -x;
  272.         w = __lgamma_r(q, sgngam); /* note this modifies sgngam! */
  273.         p = floor(q);
  274.         if( p == q )
  275.                 {
  276. lgsing:
  277.                 _SET_ERRNO(EDOM);
  278.                 mtherr( "lgam", SING );
  279. #ifdef INFINITIES
  280.                 return (INFINITY);
  281. #else
  282.                 return (MAXNUM);
  283. #endif
  284.                 }
  285.         i = p;
  286.         if( (i & 1) == 0 )
  287.                 *sgngam = -1;
  288.         else
  289.                 *sgngam = 1;
  290.         z = q - p;
  291.         if( z > 0.5 )
  292.                 {
  293.                 p += 1.0;
  294.                 z = p - q;
  295.                 }
  296.         z = q * sin( PI * z );
  297.         if( z == 0.0 )
  298.                 goto lgsing;
  299. /*      z = log(PI) - log( z ) - w;*/
  300.         z = LOGPI - log( z ) - w;
  301.         return( z );
  302.         }
  303.  
  304. if( x < 13.0 )
  305.         {
  306.         z = 1.0;
  307.         p = 0.0;
  308.         u = x;
  309.         while( u >= 3.0 )
  310.                 {
  311.                 p -= 1.0;
  312.                 u = x + p;
  313.                 z *= u;
  314.                 }
  315.         while( u < 2.0 )
  316.                 {
  317.                 if( u == 0.0 )
  318.                         goto lgsing;
  319.                 z /= u;
  320.                 p += 1.0;
  321.                 u = x + p;
  322.                 }
  323.         if( z < 0.0 )
  324.                 {
  325.                 *sgngam = -1;
  326.                 z = -z;
  327.                 }
  328.         else
  329.                 *sgngam = 1;
  330.         if( u == 2.0 )
  331.                 return( log(z) );
  332.         p -= 2.0;
  333.         x = x + p;
  334.         p = x * polevl( x, B, 5 ) / p1evl( x, C, 6);
  335.         return( log(z) + p );
  336.         }
  337.  
  338. if( x > MAXLGM )
  339.         {
  340.         _SET_ERRNO(ERANGE);
  341.         mtherr( "lgamma", OVERFLOW );
  342. #ifdef INFINITIES
  343.         return( *sgngam * INFINITY );
  344. #else
  345.         return( *sgngam * MAXNUM );
  346. #endif
  347.         }
  348.  
  349. q = ( x - 0.5 ) * log(x) - x + LS2PI;
  350. if( x > 1.0e8 )
  351.         return( q );
  352.  
  353. p = 1.0/(x*x);
  354. if( x >= 1000.0 )
  355.         q += ((   7.9365079365079365079365e-4 * p
  356.                 - 2.7777777777777777777778e-3) *p
  357.                 + 0.0833333333333333333333) / x;
  358. else
  359.         q += polevl( p, A, 4 ) / x;
  360. return( q );
  361. }
  362.  
  363. /* This is the C99 version */
  364.  
  365. double lgamma(double x)
  366. {
  367.   int local_sgngam=0;
  368.   return (__lgamma_r(x, &local_sgngam));
  369. }
  370.