home *** CD-ROM | disk | FTP | other *** search
/ Total Destruction / Total_Destruction.iso / addons / Lccwin32.exe / Lccwin32 / lccpub / lib / src / gamma.c < prev    next >
Encoding:
C/C++ Source or Header  |  1997-08-14  |  2.7 KB  |  109 lines
  1. #include <math.h>
  2. #include <stdio.h>
  3. #include <errno.h>
  4. /*@ Ensemble de fonctions pour calculer la fonction gamma
  5. Calcul du logarithme du valeur absolue de la fonction gamma pour \'eviter
  6. les d\'ebordements
  7. Les coefficients pour l'expansion au tour de z\'ero sont ceux de Hart and
  8. Cheney (\#5243) et pour l'expansion au tour de $\inf$ (\#5404)
  9.  
  10. */
  11. #define M 6
  12. #define N 8
  13. static double p1[] = {
  14.         0.83333333333333101837e-1, /* 1 / 12   */
  15.         -.277777777735865004e-2,   /* 1 / 360  */
  16.         0.793650576493454e-3,      /* 1 / 1260 */
  17.         -.5951896861197e-3,     /* 1 / 1680 */
  18.         0.83645878922e-3,
  19.         -.1633436431e-2
  20. };
  21. static double p2[] = {
  22.         -.42353689509744089647e5,
  23.         -.20886861789269887364e5,
  24.         -.87627102978521489560e4,
  25.         -.20085274013072791214e4,
  26.         -.43933044406002567613e3,
  27.         -.50108693752970953015e2,
  28.         -.67449507245925289918e1,
  29.         0.0
  30. };
  31. static double q2[] = {
  32.         -.42353689509744090010e5,
  33.         -.29803853309256649932e4,
  34.         0.99403074150827709015e4,
  35.         -.15286072737795220248e4,
  36.         -.49902852662143904834e3,
  37.         0.18949823415702801641e3,
  38.         -.23081551524580124562e2,
  39.         0.10000000000000000000e1
  40. };
  41. static double factor  = 0.9189385332046727417803297;
  42. int signgam = 0;
  43.  
  44. static double asym(double arg);
  45. static double neg(double arg);
  46. static double pos(double arg);
  47. double lgamma(double arg)
  48. {
  49.         signgam = 1;
  50.         if(arg <= 0.0) return(neg(arg));
  51.         if(arg > 8.0) return(asym(arg));
  52.         return((double)signgam * log(pos(arg)));
  53. }
  54.  
  55. static double asym(double arg)
  56. {
  57.         double n, argsq;
  58.         register int i;
  59.  
  60.         argsq = 1.0/(arg*arg);
  61.         for(n = 0.0, i = M-1; i>=0; i--){
  62.                 n = n*argsq + p1[i];
  63.         }
  64.         return((arg - 0.5)*log(arg) - arg + factor + n/arg);
  65. }
  66.  
  67. static double neg(double arg)
  68. {
  69.         double temp;
  70.  
  71.         arg = -arg;
  72.         temp = sin(M_PI*arg);
  73.         if(temp == 0.0) {
  74.                 return(HUGE_VAL);
  75.         }
  76.         if(temp < 0.0) temp = -temp;
  77.         else signgam = -1;
  78.         return(-log(arg*pos(arg)*temp/M_PI));
  79. }
  80.  
  81. static double pos(double arg)
  82. {
  83.         double n, d, s;
  84.         register int i;
  85.  
  86.         if(arg < 2.0) return(pos(arg+1.0)/arg);
  87.         if(arg > 3.0) return((arg-1.0)*pos(arg-1.0));
  88.  
  89.         s = arg - 2.;
  90.         for(n=0.0,d=0.0,i = N-1; i>=0; i--){
  91.                 n = n*s + p2[i];
  92.                 d = d*s + q2[i];
  93.         }
  94.         return(n/d);
  95. }
  96. #ifdef STANDALONE
  97. int main()
  98. {
  99.         int i=1;
  100.         double r;
  101.  
  102.         while (i < 10) {
  103.                 r= gamma((double)i);
  104.                 printf("[%3d] %10.15f\n",i,r);
  105.                 i++;
  106.         }
  107. } 
  108. #endif
  109.