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