home *** CD-ROM | disk | FTP | other *** search
/ Il CD di internet / CD.iso / SOURCE / D / LIBC / LIBC-4.6 / LIBC-4 / libc-linux / sysdeps / linux / i386 / math / lgammal.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-26  |  4.4 KB  |  192 lines

  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.  The Berkeley software License Agreement
  4.  * specifies the terms and conditions for redistribution.
  5.  */
  6.  
  7. /* Coefficients for long double precision by Stephen L. Moshier
  8.    (moshier@world.std.com).  This routine was spot chaecked at
  9.    10,000 random arguments between -40 and +40.  The largest error
  10.    observed was 2.2e-19.  Error criterion was absolute for function
  11.    value < 1, relative otherwise.  */
  12.  
  13. #if 0
  14. #ifndef lint
  15. static char sccsid[] = "@(#)lgamma.c    5.3 (Berkeley) 9/22/88";
  16. #endif /* not lint */
  17. #endif
  18.  
  19. /*
  20.     C program for floating point log Gamma function
  21.  
  22.     lgamma(x) computes the log of the absolute
  23.     value of the Gamma function.
  24.     The sign of the Gamma function is returned in the
  25.     external quantity signgaml.
  26.  
  27.     The coefficients for expansion around zero
  28.     are #5243 from Hart & Cheney; for expansion
  29.     around infinity they are #5404.
  30.  
  31.     Calls log, floor and sin.
  32. */
  33. #include "mathl.h"
  34. long double floorl(long double);
  35. long double logl(long double);
  36. long double sinl(long double);
  37.  
  38. #if 0
  39. #include "mathimpl.h"
  40. #else
  41. #include <math.h>
  42. #endif
  43. #if defined(vax)||defined(tahoe)
  44. #include <errno.h>
  45. #endif    /* defined(vax)||defined(tahoe) */
  46.  
  47. int signgaml;
  48.  
  49.   /* log(2*pi)/2 */
  50. static const long double goobie = 0.9189385332046727417803297L;
  51. static const long double pi   = 3.1415926535897932384626434L;
  52.  
  53. /*
  54. log gamma(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x P(1/x^2)
  55. Relative error
  56. n=6, d=0
  57. Peak error =  1.51e-21
  58. Relative error spread =  5.7e-21
  59. */
  60. #define M 7
  61. static const long double p1[] = {
  62.  8.33333333333333144750e-2L,
  63. -2.77777777775034960344e-3L,
  64.  7.93650779585507075567e-4L,
  65. -5.95234585176568851461e-4L,
  66.  8.41272329732249808063e-4L,
  67. -1.88080193811937690718e-3L,
  68.  4.88502614243227078116e-3L
  69. };
  70. /*
  71. gamma(x+2) = P(x)/Q(x)
  72. Relative error
  73. n=7, d=8
  74. Peak error =  1.84e-20
  75. Relative error spread =  8.4e-23
  76. */
  77. #define N 9
  78. static const long double p2[] = {
  79. -7.15743521530849602425e4L,
  80. -5.99650230220855581642e4L,
  81. -2.59780216007146401957e4L,
  82. -7.96667499622741999770e3L,
  83. -1.70730828800510297666e3L,
  84. -2.92929976820724030353e2L,
  85. -3.25157411956062339893e1L,
  86. -3.01525602666895735709e0L,
  87.  0.0L
  88. };
  89. static const long double q2[] = {
  90. -7.15743521530849602412e4L,
  91. -2.97045081369399940529e4L,
  92.  1.60577839621734713377e4L,
  93.  3.31667508019495079814e3L,
  94. -1.98526250512761318471e3L,
  95.  5.69440799097468430177e1L,
  96.  8.85946791747759881659e1L,
  97. -1.67955233807178858919e1L,
  98.  1.00000000000000000000e0L,
  99. };
  100.  
  101. static long double pos(long double);
  102. static long double neg(long double);
  103. static long double asym(long double);
  104.  
  105. long double
  106. lgammal(long double arg)
  107. {
  108.  
  109.     signgaml = 1;
  110.     if(arg <= 0.L) return(neg(arg));
  111.     if(arg > 8.L) return(asym(arg));
  112.     return(logl(pos(arg)));
  113. }
  114.  
  115. static long double
  116. asym(long double arg)
  117. {
  118.     long double n, argsq;
  119.     int i;
  120.  
  121.     n = 0.0L;
  122.  
  123.     /* Avoid overflow of arg*arg. */
  124.     if(arg > 1.0e10L)
  125.       goto noexpan;
  126.  
  127.     argsq = 1./(arg*arg);
  128.     for(i=M-1; i>=0; i--){
  129.         n = n*argsq + p1[i];
  130.     }
  131. noexpan:
  132.     return((arg-.5L)*logl(arg) - arg + goobie + n/arg);
  133. }
  134.  
  135. static long double
  136. neg(long double arg)
  137. {
  138.     long double t;
  139.  
  140.     arg = -arg;
  141.      /*
  142.       * to see if arg were a true integer, the old code used the
  143.       * mathematically correct observation:
  144.       * sin(n*pi) = 0 <=> n is an integer.
  145.       * but in finite precision arithmetic, sin(n*PI) will NEVER
  146.       * be zero simply because n*PI is a rational number.  hence
  147.       *    it failed to work with our newer, more accurate sin()
  148.       * which uses true pi to do the argument reduction...
  149.       *    temp = sin(pi*arg);
  150.       */
  151.     t = floorl(arg);
  152.     if (arg - t  > 0.5e0L)
  153.         t += 1.e0L;                /* t := integer nearest arg */
  154. #if defined(vax)||defined(tahoe)
  155.     if (arg == t) {
  156.         return(__infnan(ERANGE));        /* +INF */
  157.     }
  158. #endif    /* defined(vax)||defined(tahoe) */
  159.     signgaml = (int) (t - 2*floorl(t/2));    /* signgam =  1 if t was odd, */
  160.                         /*            0 if t was even */
  161.     signgaml = signgam - 1 + signgam;    /* signgam =  1 if t was odd, */
  162.                         /*           -1 if t was even */
  163.     t = arg - t;                /*  -0.5 <= t <= 0.5 */
  164.     if (t < 0.e0L) {
  165.         t = -t;
  166.         signgaml = -signgaml;
  167.     }
  168.     return(-logl(arg*pos(arg)*sinl(pi*t)/pi));
  169. }
  170.  
  171. static long double
  172. pos(long double arg)
  173. {
  174.     long double n, d, s;
  175.     register i;
  176.  
  177.     if(arg < 2.L) return(pos(arg+1.L)/arg);
  178.     if(arg > 3.L) return((arg-1.0L)*pos(arg-1.0L));
  179.  
  180.     s = arg - 2.0L;
  181. /*    for(n=0,d=0,i=N-1; i>=0; i--){*/
  182. /*
  183.     n = 0.0L;
  184.     d = 0.0L;
  185. */
  186.     for(n=0,d=0,i=N-1; i>=0; i--){
  187.         n = n*s + p2[i];
  188.         d = d*s + q2[i];
  189.     }
  190.     return(n/d);
  191. }
  192.