home *** CD-ROM | disk | FTP | other *** search
/ Magazyn Amiga 13 / MA_Cover_13.bin / source / c / apl / gamma.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-11-23  |  2.1 KB  |  116 lines
  1. /*
  2.     C program for floating point log gamma function
  3.  
  4.     gamma(x) computes the log of the absolute
  5.     value of the gamma function.
  6.     The sign of the gamma function is returned in the
  7.     external quantity signgam.
  8.  
  9.     The coefficients for expansion around zero
  10.     are #5243 from Hart & Cheney; for expansion
  11.     around infinity they are #5404.
  12.  
  13.     Calls log and sin.
  14. */
  15.  
  16. #include <errno.h>
  17. #include <math.h>
  18.  
  19. int    errno;
  20. int    signgam = 0;
  21. static double goobie    = 0.9189385332046727417803297;
  22. static double pi    = 3.1415926535897932384626434;
  23. static double asym(), pos(), neg();
  24.  
  25. #define M 6
  26. #define N 8
  27. static double p1[] = {
  28.     0.83333333333333101837e-1,
  29.     -.277777777735865004e-2,
  30.     0.793650576493454e-3,
  31.     -.5951896861197e-3,
  32.     0.83645878922e-3,
  33.     -.1633436431e-2,
  34. };
  35. static double p2[] = {
  36.     -.42353689509744089647e5,
  37.     -.20886861789269887364e5,
  38.     -.87627102978521489560e4,
  39.     -.20085274013072791214e4,
  40.     -.43933044406002567613e3,
  41.     -.50108693752970953015e2,
  42.     -.67449507245925289918e1,
  43.     0.0,
  44. };
  45. static double q2[] = {
  46.     -.42353689509744090010e5,
  47.     -.29803853309256649932e4,
  48.     0.99403074150827709015e4,
  49.     -.15286072737795220248e4,
  50.     -.49902852662143904834e3,
  51.     0.18949823415702801641e3,
  52.     -.23081551524580124562e2,
  53.     0.10000000000000000000e1,
  54. };
  55.  
  56. double
  57. gamma(arg)
  58. double arg;
  59. {
  60.     double log(), pos(), neg(), asym();
  61.  
  62.     signgam = 1.;
  63.     if(arg <= 0.) return(neg(arg));
  64.     if(arg > 8.) return(asym(arg));
  65.     return(log(pos(arg)));
  66. }
  67.  
  68. static double
  69. asym(arg)
  70. double arg;
  71. {
  72.     double log();
  73.     double n, argsq;
  74.     int i;
  75.  
  76.     argsq = 1./(arg*arg);
  77.     for(n=0,i=M-1; i>=0; i--) n = n*argsq + p1[i];
  78.     return((arg-.5)*log(arg) - arg + goobie + n/arg);
  79. }
  80.  
  81. static double
  82. neg(arg)
  83. double arg;
  84. {
  85.     double temp;
  86.     double log(), sin(), pos();
  87.  
  88.     arg = -arg;
  89.     temp = sin(pi*arg);
  90.     if(temp == 0.) {
  91.         errno = EDOM;
  92.         return(HUGE_VAL);
  93.     }
  94.     if(temp < 0.) temp = -temp;
  95.     else signgam = -1;
  96.     return(-log(arg*pos(arg)*temp/pi));
  97. }
  98.  
  99. static double
  100. pos(arg)
  101. double arg;
  102. {
  103.     double n, d, s;
  104.     int i;
  105.  
  106.     if(arg < 2.) return(pos(arg+1.)/arg);
  107.     if(arg > 3.) return((arg-1.)*pos(arg-1.));
  108.  
  109.     s = arg - 2.;
  110.     for(n=0,d=0,i=N-1; i>=0; i--){
  111.         n = n*s + p2[i];
  112.         d = d*s + q2[i];
  113.     }
  114.     return(n/d);
  115. }
  116.