Subversion Repositories Kolibri OS

Rev

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

  1. /*                                                      lgaml()
  2.  *
  3.  *      Natural logarithm of gamma function
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * long double x, y, __lgammal_r();
  10.  * int* sgngaml;
  11.  * y = __lgammal_r( x, sgngaml );
  12.  *
  13.  * long double x, y, lgammal();
  14.  * y = lgammal( 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 sgngaml.
  24.  *
  25.  * For arguments greater than 33, 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 MAXLGML (10^4928) return MAXNUML.
  34.  *
  35.  *
  36.  *
  37.  * ACCURACY:
  38.  *
  39.  *
  40.  * arithmetic      domain        # trials     peak         rms
  41.  *    IEEE         -40, 40        100000     2.2e-19     4.6e-20
  42.  *    IEEE    10^-2000,10^+2000    20000     1.6e-19     3.3e-20
  43.  * The error criterion was relative when the function magnitude
  44.  * was greater than one but absolute when it was less than one.
  45.  *
  46.  */
  47.  
  48. /*
  49.  * Copyright 1994 by Stephen L. Moshier
  50.  */
  51.  
  52. /*
  53.  * 26-11-2002 Modified for mingw.
  54.  * Danny Smith <dannysmith@users.sourceforge.net>
  55.  */
  56.  
  57. #ifndef __MINGW32__
  58. #include "mconf.h"
  59. #ifdef ANSIPROT
  60. extern long double fabsl ( long double );
  61. extern long double lgaml ( long double );
  62. extern long double logl ( long double );
  63. extern long double expl ( long double );
  64. extern long double gammal ( long double );
  65. extern long double sinl ( long double );
  66. extern long double floorl ( long double );
  67. extern long double powl ( long double, long double );
  68. extern long double polevll ( long double, void *, int );
  69. extern long double p1evll ( long double, void *, int );
  70. extern int isnanl ( long double );
  71. extern int isfinitel ( long double );
  72. #else
  73. long double fabsl(), lgaml(), logl(), expl(), gammal(), sinl();
  74. long double floorl(), powl(), polevll(), p1evll(), isnanl(), isfinitel();
  75. #endif
  76. #ifdef INFINITIES
  77. extern long double INFINITYL;
  78. #endif
  79. #ifdef NANS
  80. extern long double NANL;
  81. #endif
  82. #else /* __MINGW32__ */
  83. #include "cephes_mconf.h"
  84. #endif /* __MINGW32__ */
  85.  
  86. #if UNK
  87. static long double S[9] = {
  88. -1.193945051381510095614E-3L,
  89.  7.220599478036909672331E-3L,
  90. -9.622023360406271645744E-3L,
  91. -4.219773360705915470089E-2L,
  92.  1.665386113720805206758E-1L,
  93. -4.200263503403344054473E-2L,
  94. -6.558780715202540684668E-1L,
  95.  5.772156649015328608253E-1L,
  96.  1.000000000000000000000E0L,
  97. };
  98. #endif
  99. #if IBMPC
  100. static const unsigned short S[] = {
  101. 0xbaeb,0xd6d3,0x25e5,0x9c7e,0xbff5, XPD
  102. 0xfe9a,0xceb4,0xc74e,0xec9a,0x3ff7, XPD
  103. 0x9225,0xdfef,0xb0e9,0x9da5,0xbff8, XPD
  104. 0x10b0,0xec17,0x87dc,0xacd7,0xbffa, XPD
  105. 0x6b8d,0x7515,0x1905,0xaa89,0x3ffc, XPD
  106. 0xf183,0x126b,0xf47d,0xac0a,0xbffa, XPD
  107. 0x7bf6,0x57d1,0xa013,0xa7e7,0xbffe, XPD
  108. 0xc7a9,0x7db0,0x67e3,0x93c4,0x3ffe, XPD
  109. 0x0000,0x0000,0x0000,0x8000,0x3fff, XPD
  110. };
  111. #endif
  112. #if MIEEE
  113. static long S[27] = {
  114. 0xbff50000,0x9c7e25e5,0xd6d3baeb,
  115. 0x3ff70000,0xec9ac74e,0xceb4fe9a,
  116. 0xbff80000,0x9da5b0e9,0xdfef9225,
  117. 0xbffa0000,0xacd787dc,0xec1710b0,
  118. 0x3ffc0000,0xaa891905,0x75156b8d,
  119. 0xbffa0000,0xac0af47d,0x126bf183,
  120. 0xbffe0000,0xa7e7a013,0x57d17bf6,
  121. 0x3ffe0000,0x93c467e3,0x7db0c7a9,
  122. 0x3fff0000,0x80000000,0x00000000,
  123. };
  124. #endif
  125.  
  126. #if UNK
  127. static long double SN[9] = {
  128.  1.133374167243894382010E-3L,
  129.  7.220837261893170325704E-3L,
  130.  9.621911155035976733706E-3L,
  131. -4.219773343731191721664E-2L,
  132. -1.665386113944413519335E-1L,
  133. -4.200263503402112910504E-2L,
  134.  6.558780715202536547116E-1L,
  135.  5.772156649015328608727E-1L,
  136. -1.000000000000000000000E0L,
  137. };
  138. #endif
  139. #if IBMPC
  140. static const unsigned SN[] = {
  141. 0x5dd1,0x02de,0xb9f7,0x948d,0x3ff5, XPD
  142. 0x989b,0xdd68,0xc5f1,0xec9c,0x3ff7, XPD
  143. 0x2ca1,0x18f0,0x386f,0x9da5,0x3ff8, XPD
  144. 0x783f,0x41dd,0x87d1,0xacd7,0xbffa, XPD
  145. 0x7a5b,0xd76d,0x1905,0xaa89,0xbffc, XPD
  146. 0x7f64,0x1234,0xf47d,0xac0a,0xbffa, XPD
  147. 0x5e26,0x57d1,0xa013,0xa7e7,0x3ffe, XPD
  148. 0xc7aa,0x7db0,0x67e3,0x93c4,0x3ffe, XPD
  149. 0x0000,0x0000,0x0000,0x8000,0xbfff, XPD
  150. };
  151. #endif
  152. #if MIEEE
  153. static long SN[27] = {
  154. 0x3ff50000,0x948db9f7,0x02de5dd1,
  155. 0x3ff70000,0xec9cc5f1,0xdd68989b,
  156. 0x3ff80000,0x9da5386f,0x18f02ca1,
  157. 0xbffa0000,0xacd787d1,0x41dd783f,
  158. 0xbffc0000,0xaa891905,0xd76d7a5b,
  159. 0xbffa0000,0xac0af47d,0x12347f64,
  160. 0x3ffe0000,0xa7e7a013,0x57d15e26,
  161. 0x3ffe0000,0x93c467e3,0x7db0c7aa,
  162. 0xbfff0000,0x80000000,0x00000000,
  163. };
  164. #endif
  165.  
  166.  
  167. /* A[]: Stirling's formula expansion of log gamma
  168.  * B[], C[]: log gamma function between 2 and 3
  169.  */
  170.  
  171.  
  172. /* log gamma(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x A(1/x^2)
  173.  * x >= 8
  174.  * Peak relative error 1.51e-21
  175.  * Relative spread of error peaks 5.67e-21
  176.  */
  177. #if UNK
  178. static long double A[7] = {
  179.  4.885026142432270781165E-3L,
  180. -1.880801938119376907179E-3L,
  181.  8.412723297322498080632E-4L,
  182. -5.952345851765688514613E-4L,
  183.  7.936507795855070755671E-4L,
  184. -2.777777777750349603440E-3L,
  185.  8.333333333333331447505E-2L,
  186. };
  187. #endif
  188. #if IBMPC
  189. static const unsigned short A[] = {
  190. 0xd984,0xcc08,0x91c2,0xa012,0x3ff7, XPD
  191. 0x3d91,0x0304,0x3da1,0xf685,0xbff5, XPD
  192. 0x3bdc,0xaad1,0xd492,0xdc88,0x3ff4, XPD
  193. 0x8b20,0x9fce,0x844e,0x9c09,0xbff4, XPD
  194. 0xf8f2,0x30e5,0x0092,0xd00d,0x3ff4, XPD
  195. 0x4d88,0x03a8,0x60b6,0xb60b,0xbff6, XPD
  196. 0x9fcc,0xaaaa,0xaaaa,0xaaaa,0x3ffb, XPD
  197. };
  198. #endif
  199. #if MIEEE
  200. static long A[21] = {
  201. 0x3ff70000,0xa01291c2,0xcc08d984,
  202. 0xbff50000,0xf6853da1,0x03043d91,
  203. 0x3ff40000,0xdc88d492,0xaad13bdc,
  204. 0xbff40000,0x9c09844e,0x9fce8b20,
  205. 0x3ff40000,0xd00d0092,0x30e5f8f2,
  206. 0xbff60000,0xb60b60b6,0x03a84d88,
  207. 0x3ffb0000,0xaaaaaaaa,0xaaaa9fcc,
  208. };
  209. #endif
  210.  
  211. /* log gamma(x+2) = x B(x)/C(x)
  212.  * 0 <= x <= 1
  213.  * Peak relative error 7.16e-22
  214.  * Relative spread of error peaks 4.78e-20
  215.  */
  216. #if UNK
  217. static long double B[7] = {
  218. -2.163690827643812857640E3L,
  219. -8.723871522843511459790E4L,
  220. -1.104326814691464261197E6L,
  221. -6.111225012005214299996E6L,
  222. -1.625568062543700591014E7L,
  223. -2.003937418103815175475E7L,
  224. -8.875666783650703802159E6L,
  225. };
  226. static long double C[7] = {
  227. /* 1.000000000000000000000E0L,*/
  228. -5.139481484435370143617E2L,
  229. -3.403570840534304670537E4L,
  230. -6.227441164066219501697E5L,
  231. -4.814940379411882186630E6L,
  232. -1.785433287045078156959E7L,
  233. -3.138646407656182662088E7L,
  234. -2.099336717757895876142E7L,
  235. };
  236. #endif
  237. #if IBMPC
  238. static const unsigned short B[] = {
  239. 0x9557,0x4995,0x0da1,0x873b,0xc00a, XPD
  240. 0xfe44,0x9af8,0x5b8c,0xaa63,0xc00f, XPD
  241. 0x5aa8,0x7cf5,0x3684,0x86ce,0xc013, XPD
  242. 0x259a,0x258c,0xf206,0xba7f,0xc015, XPD
  243. 0xbe18,0x1ca3,0xc0a0,0xf80a,0xc016, XPD
  244. 0x168f,0x2c42,0x6717,0x98e3,0xc017, XPD
  245. 0x2051,0x9d55,0x92c8,0x876e,0xc016, XPD
  246. };
  247. static const unsigned short C[] = {
  248. /*0x0000,0x0000,0x0000,0x8000,0x3fff, XPD*/
  249. 0xaa77,0xcf2f,0xae76,0x807c,0xc008, XPD
  250. 0xb280,0x0d74,0xb55a,0x84f3,0xc00e, XPD
  251. 0xa505,0xcd30,0x81dc,0x9809,0xc012, XPD
  252. 0x3369,0x4246,0xb8c2,0x92f0,0xc015, XPD
  253. 0x63cf,0x6aee,0xbe6f,0x8837,0xc017, XPD
  254. 0x26bb,0xccc7,0xb009,0xef75,0xc017, XPD
  255. 0x462b,0xbae8,0xab96,0xa02a,0xc017, XPD
  256. };
  257. #endif
  258. #if MIEEE
  259. static long B[21] = {
  260. 0xc00a0000,0x873b0da1,0x49959557,
  261. 0xc00f0000,0xaa635b8c,0x9af8fe44,
  262. 0xc0130000,0x86ce3684,0x7cf55aa8,
  263. 0xc0150000,0xba7ff206,0x258c259a,
  264. 0xc0160000,0xf80ac0a0,0x1ca3be18,
  265. 0xc0170000,0x98e36717,0x2c42168f,
  266. 0xc0160000,0x876e92c8,0x9d552051,
  267. };
  268. static long C[21] = {
  269. /*0x3fff0000,0x80000000,0x00000000,*/
  270. 0xc0080000,0x807cae76,0xcf2faa77,
  271. 0xc00e0000,0x84f3b55a,0x0d74b280,
  272. 0xc0120000,0x980981dc,0xcd30a505,
  273. 0xc0150000,0x92f0b8c2,0x42463369,
  274. 0xc0170000,0x8837be6f,0x6aee63cf,
  275. 0xc0170000,0xef75b009,0xccc726bb,
  276. 0xc0170000,0xa02aab96,0xbae8462b,
  277. };
  278. #endif
  279.  
  280. /* log( sqrt( 2*pi ) ) */
  281. static const long double LS2PI  =  0.91893853320467274178L;
  282. #define MAXLGM 1.04848146839019521116e+4928L
  283.  
  284.  
  285. /* Logarithm of gamma function */
  286. /* Reentrant version */
  287.  
  288. long double __lgammal_r(long double x, int* sgngaml)
  289. {
  290. long double p, q, w, z, f, nx;
  291. int i;
  292.  
  293. *sgngaml = 1;
  294. #ifdef NANS
  295. if( isnanl(x) )
  296.         return(NANL);
  297. #endif
  298. #ifdef INFINITIES
  299. if( !isfinitel(x) )
  300.         return(INFINITYL);
  301. #endif
  302. if( x < -34.0L )
  303.         {
  304.         q = -x;
  305.         w = __lgammal_r(q, sgngaml); /* note this modifies sgngam! */
  306.         p = floorl(q);
  307.         if( p == q )
  308.                 {
  309. lgsing:
  310.                 _SET_ERRNO(EDOM);
  311.                 mtherr( "lgammal", SING );
  312. #ifdef INFINITIES
  313.                 return (INFINITYL);
  314. #else
  315.                 return (MAXNUML);
  316. #endif
  317.                 }
  318.         i = p;
  319.         if( (i & 1) == 0 )
  320.                 *sgngaml = -1;
  321.         else
  322.                 *sgngaml = 1;
  323.         z = q - p;
  324.         if( z > 0.5L )
  325.                 {
  326.                 p += 1.0L;
  327.                 z = p - q;
  328.                 }
  329.         z = q * sinl( PIL * z );
  330.         if( z == 0.0L )
  331.                 goto lgsing;
  332. /*      z = LOGPI - logl( z ) - w; */
  333.         z = logl( PIL/z ) - w;
  334.         return( z );
  335.         }
  336.  
  337. if( x < 13.0L )
  338.         {
  339.         z = 1.0L;
  340.         nx = floorl( x +  0.5L );
  341.         f = x - nx;
  342.         while( x >= 3.0L )
  343.                 {
  344.                 nx -= 1.0L;
  345.                 x = nx + f;
  346.                 z *= x;
  347.                 }
  348.         while( x < 2.0L )
  349.                 {
  350.                 if( fabsl(x) <= 0.03125 )
  351.                         goto lsmall;
  352.                 z /= nx +  f;
  353.                 nx += 1.0L;
  354.                 x = nx + f;
  355.                 }
  356.         if( z < 0.0L )
  357.                 {
  358.                 *sgngaml = -1;
  359.                 z = -z;
  360.                 }
  361.         else
  362.                 *sgngaml = 1;
  363.         if( x == 2.0L )
  364.                 return( logl(z) );
  365.         x = (nx - 2.0L) + f;
  366.         p = x * polevll( x, B, 6 ) / p1evll( x, C, 7);
  367.         return( logl(z) + p );
  368.         }
  369.  
  370. if( x > MAXLGM )
  371.         {
  372.         _SET_ERRNO(ERANGE);
  373.         mtherr( "lgammal", OVERFLOW );
  374. #ifdef INFINITIES
  375.         return( *sgngaml * INFINITYL );
  376. #else
  377.         return( *sgngaml * MAXNUML );
  378. #endif
  379.         }
  380.  
  381. q = ( x - 0.5L ) * logl(x) - x + LS2PI;
  382. if( x > 1.0e10L )
  383.         return(q);
  384. p = 1.0L/(x*x);
  385. q += polevll( p, A, 6 ) / x;
  386. return( q );
  387.  
  388.  
  389. lsmall:
  390. if( x == 0.0L )
  391.         goto lgsing;
  392. if( x < 0.0L )
  393.         {
  394.         x = -x;
  395.         q = z / (x * polevll( x, SN, 8 ));
  396.         }
  397. else
  398.         q = z / (x * polevll( x, S, 8 ));
  399. if( q < 0.0L )
  400.         {
  401.         *sgngaml = -1;
  402.         q = -q;
  403.         }
  404. else
  405.         *sgngaml = 1;
  406. q = logl( q );
  407. return(q);
  408. }
  409.  
  410. /* This is the C99 version */
  411.  
  412. long double lgammal(long double x)
  413. {
  414.   int local_sgngaml=0;
  415.   return (__lgammal_r(x, &local_sgngaml));
  416. }
  417.