Subversion Repositories Kolibri OS

Rev

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

  1. /*                                                      gammal.c
  2.  *
  3.  *      Gamma function
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * long double x, y, __tgammal_r();
  10.  * int* sgngaml;
  11.  * y = __tgammal_r( x, sgngaml );
  12.  *
  13.  * long double x, y, tgammal();
  14.  * y = tgammal( 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 |x| <= 13 are reduced by recurrence and the function
  24.  * approximated by a rational function of degree 7/8 in the
  25.  * interval (2,3).  Large arguments are handled by Stirling's
  26.  * formula. Large 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     -40,+40      10000       3.6e-19     7.9e-20
  35.  *    IEEE    -1755,+1755   10000       4.8e-18     6.5e-19
  36.  *
  37.  * Accuracy for large arguments is dominated by error in powl().
  38.  *
  39.  */
  40.  
  41. /*
  42. Copyright 1994 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. /*
  59. gamma(x+2)  = gamma(x+2) P(x)/Q(x)
  60. 0 <= x <= 1
  61. Relative error
  62. n=7, d=8
  63. Peak error =  1.83e-20
  64. Relative error spread =  8.4e-23
  65. */
  66.  
  67. #if UNK
  68. static const long double P[8] = {
  69.  4.212760487471622013093E-5L,
  70.  4.542931960608009155600E-4L,
  71.  4.092666828394035500949E-3L,
  72.  2.385363243461108252554E-2L,
  73.  1.113062816019361559013E-1L,
  74.  3.629515436640239168939E-1L,
  75.  8.378004301573126728826E-1L,
  76.  1.000000000000000000009E0L,
  77. };
  78. static const long double Q[9] = {
  79. -1.397148517476170440917E-5L,
  80.  2.346584059160635244282E-4L,
  81. -1.237799246653152231188E-3L,
  82. -7.955933682494738320586E-4L,
  83.  2.773706565840072979165E-2L,
  84. -4.633887671244534213831E-2L,
  85. -2.243510905670329164562E-1L,
  86.  4.150160950588455434583E-1L,
  87.  9.999999999999999999908E-1L,
  88. };
  89. #endif
  90. #if IBMPC
  91. static const unsigned short P[] = {
  92. 0x434a,0x3f22,0x2bda,0xb0b2,0x3ff0, XPD
  93. 0xf5aa,0xe82f,0x335b,0xee2e,0x3ff3, XPD
  94. 0xbe6c,0x3757,0xc717,0x861b,0x3ff7, XPD
  95. 0x7f43,0x5196,0xb166,0xc368,0x3ff9, XPD
  96. 0x9549,0x8eb5,0x8c3a,0xe3f4,0x3ffb, XPD
  97. 0x8d75,0x23af,0xc8e4,0xb9d4,0x3ffd, XPD
  98. 0x29cf,0x19b3,0x16c8,0xd67a,0x3ffe, XPD
  99. 0x0000,0x0000,0x0000,0x8000,0x3fff, XPD
  100. };
  101. static const unsigned short Q[] = {
  102. 0x5473,0x2de8,0x1268,0xea67,0xbfee, XPD
  103. 0x334b,0xc2f0,0xa2dd,0xf60e,0x3ff2, XPD
  104. 0xbeed,0x1853,0xa691,0xa23d,0xbff5, XPD
  105. 0x296e,0x7cb1,0x5dfd,0xd08f,0xbff4, XPD
  106. 0x0417,0x7989,0xd7bc,0xe338,0x3ff9, XPD
  107. 0x3295,0x3698,0xd580,0xbdcd,0xbffa, XPD
  108. 0x75ef,0x3ab7,0x4ad3,0xe5bc,0xbffc, XPD
  109. 0xe458,0x2ec7,0xfd57,0xd47c,0x3ffd, XPD
  110. 0x0000,0x0000,0x0000,0x8000,0x3fff, XPD
  111. };
  112. #endif
  113. #if MIEEE
  114. static const long P[24] = {
  115. 0x3ff00000,0xb0b22bda,0x3f22434a,
  116. 0x3ff30000,0xee2e335b,0xe82ff5aa,
  117. 0x3ff70000,0x861bc717,0x3757be6c,
  118. 0x3ff90000,0xc368b166,0x51967f43,
  119. 0x3ffb0000,0xe3f48c3a,0x8eb59549,
  120. 0x3ffd0000,0xb9d4c8e4,0x23af8d75,
  121. 0x3ffe0000,0xd67a16c8,0x19b329cf,
  122. 0x3fff0000,0x80000000,0x00000000,
  123. };
  124. static const long Q[27] = {
  125. 0xbfee0000,0xea671268,0x2de85473,
  126. 0x3ff20000,0xf60ea2dd,0xc2f0334b,
  127. 0xbff50000,0xa23da691,0x1853beed,
  128. 0xbff40000,0xd08f5dfd,0x7cb1296e,
  129. 0x3ff90000,0xe338d7bc,0x79890417,
  130. 0xbffa0000,0xbdcdd580,0x36983295,
  131. 0xbffc0000,0xe5bc4ad3,0x3ab775ef,
  132. 0x3ffd0000,0xd47cfd57,0x2ec7e458,
  133. 0x3fff0000,0x80000000,0x00000000,
  134. };
  135. #endif
  136. /*
  137. static const long double P[] = {
  138. -3.01525602666895735709e0L,
  139. -3.25157411956062339893e1L,
  140. -2.92929976820724030353e2L,
  141. -1.70730828800510297666e3L,
  142. -7.96667499622741999770e3L,
  143. -2.59780216007146401957e4L,
  144. -5.99650230220855581642e4L,
  145. -7.15743521530849602425e4L
  146. };
  147. static const long double Q[] = {
  148.  1.00000000000000000000e0L,
  149. -1.67955233807178858919e1L,
  150.  8.85946791747759881659e1L,
  151.  5.69440799097468430177e1L,
  152. -1.98526250512761318471e3L,
  153.  3.31667508019495079814e3L,
  154.  1.60577839621734713377e4L,
  155. -2.97045081369399940529e4L,
  156. -7.15743521530849602412e4L
  157. };
  158. */
  159. #define MAXGAML 1755.455L
  160. /*static const long double LOGPI = 1.14472988584940017414L;*/
  161.  
  162. /* Stirling's formula for the gamma function
  163. gamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x))
  164. z(x) = x
  165. 13 <= x <= 1024
  166. Relative error
  167. n=8, d=0
  168. Peak error =  9.44e-21
  169. Relative error spread =  8.8e-4
  170. */
  171. #if UNK
  172. static const long double STIR[9] = {
  173.  7.147391378143610789273E-4L,
  174. -2.363848809501759061727E-5L,
  175. -5.950237554056330156018E-4L,
  176.  6.989332260623193171870E-5L,
  177.  7.840334842744753003862E-4L,
  178. -2.294719747873185405699E-4L,
  179. -2.681327161876304418288E-3L,
  180.  3.472222222230075327854E-3L,
  181.  8.333333333333331800504E-2L,
  182. };
  183. #endif
  184. #if IBMPC
  185. static const unsigned short STIR[] = {
  186. 0x6ede,0x69f7,0x54e3,0xbb5d,0x3ff4, XPD
  187. 0xc395,0x0295,0x4443,0xc64b,0xbfef, XPD
  188. 0xba6f,0x7c59,0x5e47,0x9bfb,0xbff4, XPD
  189. 0x5704,0x1a39,0xb11d,0x9293,0x3ff1, XPD
  190. 0x30b7,0x1a21,0x98b2,0xcd87,0x3ff4, XPD
  191. 0xbef3,0x7023,0x6a08,0xf09e,0xbff2, XPD
  192. 0x3a1c,0x5ac8,0x3478,0xafb9,0xbff6, XPD
  193. 0xc3c9,0x906e,0x38e3,0xe38e,0x3ff6, XPD
  194. 0xa1d5,0xaaaa,0xaaaa,0xaaaa,0x3ffb, XPD
  195. };
  196. #endif
  197. #if MIEEE
  198. static const long STIR[27] = {
  199. 0x3ff40000,0xbb5d54e3,0x69f76ede,
  200. 0xbfef0000,0xc64b4443,0x0295c395,
  201. 0xbff40000,0x9bfb5e47,0x7c59ba6f,
  202. 0x3ff10000,0x9293b11d,0x1a395704,
  203. 0x3ff40000,0xcd8798b2,0x1a2130b7,
  204. 0xbff20000,0xf09e6a08,0x7023bef3,
  205. 0xbff60000,0xafb93478,0x5ac83a1c,
  206. 0x3ff60000,0xe38e38e3,0x906ec3c9,
  207. 0x3ffb0000,0xaaaaaaaa,0xaaaaa1d5,
  208. };
  209. #endif
  210. #define MAXSTIR 1024.0L
  211. static const long double SQTPI = 2.50662827463100050242E0L;
  212.  
  213. /* 1/gamma(x) = z P(z)
  214.  * z(x) = 1/x
  215.  * 0 < x < 0.03125
  216.  * Peak relative error 4.2e-23
  217.  */
  218. #if UNK
  219. static const long double S[9] = {
  220. -1.193945051381510095614E-3L,
  221.  7.220599478036909672331E-3L,
  222. -9.622023360406271645744E-3L,
  223. -4.219773360705915470089E-2L,
  224.  1.665386113720805206758E-1L,
  225. -4.200263503403344054473E-2L,
  226. -6.558780715202540684668E-1L,
  227.  5.772156649015328608253E-1L,
  228.  1.000000000000000000000E0L,
  229. };
  230. #endif
  231. #if IBMPC
  232. static const unsigned short S[] = {
  233. 0xbaeb,0xd6d3,0x25e5,0x9c7e,0xbff5, XPD
  234. 0xfe9a,0xceb4,0xc74e,0xec9a,0x3ff7, XPD
  235. 0x9225,0xdfef,0xb0e9,0x9da5,0xbff8, XPD
  236. 0x10b0,0xec17,0x87dc,0xacd7,0xbffa, XPD
  237. 0x6b8d,0x7515,0x1905,0xaa89,0x3ffc, XPD
  238. 0xf183,0x126b,0xf47d,0xac0a,0xbffa, XPD
  239. 0x7bf6,0x57d1,0xa013,0xa7e7,0xbffe, XPD
  240. 0xc7a9,0x7db0,0x67e3,0x93c4,0x3ffe, XPD
  241. 0x0000,0x0000,0x0000,0x8000,0x3fff, XPD
  242. };
  243. #endif
  244. #if MIEEE
  245. static const long S[27] = {
  246. 0xbff50000,0x9c7e25e5,0xd6d3baeb,
  247. 0x3ff70000,0xec9ac74e,0xceb4fe9a,
  248. 0xbff80000,0x9da5b0e9,0xdfef9225,
  249. 0xbffa0000,0xacd787dc,0xec1710b0,
  250. 0x3ffc0000,0xaa891905,0x75156b8d,
  251. 0xbffa0000,0xac0af47d,0x126bf183,
  252. 0xbffe0000,0xa7e7a013,0x57d17bf6,
  253. 0x3ffe0000,0x93c467e3,0x7db0c7a9,
  254. 0x3fff0000,0x80000000,0x00000000,
  255. };
  256. #endif
  257. /* 1/gamma(-x) = z P(z)
  258.  * z(x) = 1/x
  259.  * 0 < x < 0.03125
  260.  * Peak relative error 5.16e-23
  261.  * Relative error spread =  2.5e-24
  262.  */
  263. #if UNK
  264. static const long double SN[9] = {
  265.  1.133374167243894382010E-3L,
  266.  7.220837261893170325704E-3L,
  267.  9.621911155035976733706E-3L,
  268. -4.219773343731191721664E-2L,
  269. -1.665386113944413519335E-1L,
  270. -4.200263503402112910504E-2L,
  271.  6.558780715202536547116E-1L,
  272.  5.772156649015328608727E-1L,
  273. -1.000000000000000000000E0L,
  274. };
  275. #endif
  276. #if IBMPC
  277. static const unsigned short SN[] = {
  278. 0x5dd1,0x02de,0xb9f7,0x948d,0x3ff5, XPD
  279. 0x989b,0xdd68,0xc5f1,0xec9c,0x3ff7, XPD
  280. 0x2ca1,0x18f0,0x386f,0x9da5,0x3ff8, XPD
  281. 0x783f,0x41dd,0x87d1,0xacd7,0xbffa, XPD
  282. 0x7a5b,0xd76d,0x1905,0xaa89,0xbffc, XPD
  283. 0x7f64,0x1234,0xf47d,0xac0a,0xbffa, XPD
  284. 0x5e26,0x57d1,0xa013,0xa7e7,0x3ffe, XPD
  285. 0xc7aa,0x7db0,0x67e3,0x93c4,0x3ffe, XPD
  286. 0x0000,0x0000,0x0000,0x8000,0xbfff, XPD
  287. };
  288. #endif
  289. #if MIEEE
  290. static const long SN[27] = {
  291. 0x3ff50000,0x948db9f7,0x02de5dd1,
  292. 0x3ff70000,0xec9cc5f1,0xdd68989b,
  293. 0x3ff80000,0x9da5386f,0x18f02ca1,
  294. 0xbffa0000,0xacd787d1,0x41dd783f,
  295. 0xbffc0000,0xaa891905,0xd76d7a5b,
  296. 0xbffa0000,0xac0af47d,0x12347f64,
  297. 0x3ffe0000,0xa7e7a013,0x57d15e26,
  298. 0x3ffe0000,0x93c467e3,0x7db0c7aa,
  299. 0xbfff0000,0x80000000,0x00000000,
  300. };
  301. #endif
  302.  
  303. #ifndef __MINGW32__
  304. extern long double MAXLOGL, MAXNUML, PIL;
  305. /* #define PIL 3.14159265358979323846L */
  306. /* #define MAXNUML 1.189731495357231765021263853E4932L */
  307.  
  308. #ifdef ANSIPROT
  309. extern long double fabsl ( long double );
  310. extern long double lgaml ( long double );
  311. extern long double logl ( long double );
  312. extern long double expl ( long double );
  313. extern long double gammal ( long double );
  314. extern long double sinl ( long double );
  315. extern long double floorl ( long double );
  316. extern long double powl ( long double, long double );
  317. extern long double polevll ( long double, void *, int );
  318. extern long double p1evll ( long double, void *, int );
  319. extern int isnanl ( long double );
  320. extern int isfinitel ( long double );
  321. static long double stirf ( long double );
  322. #else
  323. long double fabsl(), lgaml(), logl(), expl(), gammal(), sinl();
  324. long double floorl(), powl(), polevll(), p1evll(), isnanl(), isfinitel();
  325. static long double stirf();
  326. #endif
  327. #ifdef INFINITIES
  328. extern long double INFINITYL;
  329. #endif
  330. #ifdef NANS
  331. extern long double NANL;
  332. #endif
  333.  
  334. #else /* __MINGW32__ */
  335. static long double stirf ( long double );
  336. #endif
  337.  
  338.  
  339. /* Gamma function computed by Stirling's formula.  */
  340.  
  341. static long double stirf(x)
  342. long double x;
  343. {
  344. long double y, w, v;
  345.  
  346. w = 1.0L/x;
  347. /* For large x, use rational coefficients from the analytical expansion.  */
  348. if( x > 1024.0L )
  349.         w = (((((6.97281375836585777429E-5L * w
  350.                 + 7.84039221720066627474E-4L) * w
  351.                 - 2.29472093621399176955E-4L) * w
  352.                 - 2.68132716049382716049E-3L) * w
  353.                 + 3.47222222222222222222E-3L) * w
  354.                 + 8.33333333333333333333E-2L) * w
  355.                 + 1.0L;
  356. else
  357.         w = 1.0L + w * polevll( w, STIR, 8 );
  358. y = expl(x);
  359. if( x > MAXSTIR )
  360.         { /* Avoid overflow in pow() */
  361.         v = powl( x, 0.5L * x - 0.25L );
  362.         y = v * (v / y);
  363.         }
  364. else
  365.         {
  366.         y = powl( x, x - 0.5L ) / y;
  367.         }
  368. y = SQTPI * y * w;
  369. return( y );
  370. }
  371.  
  372.  
  373. long double __tgammal_r(long double x, int* sgngaml)
  374. {
  375. long double p, q, z;
  376. int i;
  377.  
  378. *sgngaml = 1;
  379. #ifdef NANS
  380. if( isnanl(x) )
  381.         return(NANL);
  382. #endif
  383. #ifdef INFINITIES
  384. #ifdef NANS
  385. if( x == INFINITYL )
  386.         return(x);
  387. if( x == -INFINITYL )
  388.         return(NANL);
  389. #else
  390. if( !isfinite(x) )
  391.         return(x);
  392. #endif
  393. #endif
  394. q = fabsl(x);
  395.  
  396. if( q > 13.0L )
  397.         {
  398.         if( q > MAXGAML )
  399.                 goto goverf;
  400.         if( x < 0.0L )
  401.                 {
  402.                 p = floorl(q);
  403.                 if( p == q )
  404.                         {
  405. gsing:
  406.                         _SET_ERRNO(EDOM);
  407.                         mtherr( "tgammal", SING );
  408. #ifdef INFINITIES
  409.                         return (INFINITYL);
  410. #else
  411.                         return( *sgngaml * MAXNUML);
  412. #endif
  413.                         }
  414.                 i = p;
  415.                 if( (i & 1) == 0 )
  416.                         *sgngaml = -1;
  417.                 z = q - p;
  418.                 if( z > 0.5L )
  419.                         {
  420.                         p += 1.0L;
  421.                         z = q - p;
  422.                         }
  423.                 z = q * sinl( PIL * z );
  424.                 z = fabsl(z) * stirf(q);
  425.                 if( z <= PIL/MAXNUML )
  426.                         {
  427. goverf:
  428.                         _SET_ERRNO(ERANGE);
  429.                         mtherr( "tgammal", OVERFLOW );
  430. #ifdef INFINITIES
  431.                         return( *sgngaml * INFINITYL);
  432. #else
  433.                         return( *sgngaml * MAXNUML);
  434. #endif
  435.                         }
  436.                 z = PIL/z;
  437.                 }
  438.         else
  439.                 {
  440.                 z = stirf(x);
  441.                 }
  442.         return( *sgngaml * z );
  443.         }
  444.  
  445. z = 1.0L;
  446. while( x >= 3.0L )
  447.         {
  448.         x -= 1.0L;
  449.         z *= x;
  450.         }
  451.  
  452. while( x < -0.03125L )
  453.         {
  454.         z /= x;
  455.         x += 1.0L;
  456.         }
  457.  
  458. if( x <= 0.03125L )
  459.         goto Small;
  460.  
  461. while( x < 2.0L )
  462.         {
  463.         z /= x;
  464.         x += 1.0L;
  465.         }
  466.  
  467. if( x == 2.0L )
  468.         return(z);
  469.  
  470. x -= 2.0L;
  471. p = polevll( x, P, 7 );
  472. q = polevll( x, Q, 8 );
  473. return( z * p / q );
  474.  
  475. Small:
  476. if( x == 0.0L )
  477.         {
  478.         goto gsing;
  479.         }
  480. else
  481.         {
  482.         if( x < 0.0L )
  483.                 {
  484.                 x = -x;
  485.                 q = z / (x * polevll( x, SN, 8 ));
  486.                 }
  487.         else
  488.                 q = z / (x * polevll( x, S, 8 ));
  489.         }
  490. return q;
  491. }
  492.  
  493.  
  494. /* This is the C99 version. */
  495.  
  496. long double tgammal(long double x)
  497. {
  498.   int local_sgngaml=0;
  499.   return (__tgammal_r(x, &local_sgngaml));
  500. }
  501.  
  502.