1 /* @(#)lgamma.c 4.2 06/06/85 */ 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, drem and sin. 16 */ 17 18 #include <errno.h> 19 #include <math.h> 20 #ifdef VAX 21 static long NaN_[] = {0x8000, 0x0}; 22 #define NaN (*(double *) NaN_) 23 #endif 24 25 26 int errno; 27 int signgam = 0; 28 static double goobie = 0.9189385332046727417803297; /* log(2*pi)/2 */ 29 static double pi = 3.1415926535897932384626434; 30 31 #define M 6 32 #define N 8 33 static double p1[] = { 34 0.83333333333333101837e-1, 35 -.277777777735865004e-2, 36 0.793650576493454e-3, 37 -.5951896861197e-3, 38 0.83645878922e-3, 39 -.1633436431e-2, 40 }; 41 static double p2[] = { 42 -.42353689509744089647e5, 43 -.20886861789269887364e5, 44 -.87627102978521489560e4, 45 -.20085274013072791214e4, 46 -.43933044406002567613e3, 47 -.50108693752970953015e2, 48 -.67449507245925289918e1, 49 0.0, 50 }; 51 static double q2[] = { 52 -.42353689509744090010e5, 53 -.29803853309256649932e4, 54 0.99403074150827709015e4, 55 -.15286072737795220248e4, 56 -.49902852662143904834e3, 57 0.18949823415702801641e3, 58 -.23081551524580124562e2, 59 0.10000000000000000000e1, 60 }; 61 62 double 63 gamma(arg) 64 double arg; 65 { 66 double log(), pos(), neg(), asym(); 67 68 signgam = 1.; 69 if(arg <= 0.) return(neg(arg)); 70 if(arg > 8.) return(asym(arg)); 71 return(log(pos(arg))); 72 } 73 74 static double 75 asym(arg) 76 double arg; 77 { 78 double log(); 79 double n, argsq; 80 int i; 81 82 argsq = 1./(arg*arg); 83 for(n=0,i=M-1; i>=0; i--){ 84 n = n*argsq + p1[i]; 85 } 86 return((arg-.5)*log(arg) - arg + goobie + n/arg); 87 } 88 89 static double 90 neg(arg) 91 double arg; 92 { 93 double temp; 94 double log(), sin(), drem(), pos(); 95 96 arg = -arg; 97 /* 98 * to see if arg were a true integer, the old code used the 99 * mathematically correct observation: 100 * sin(n*pi) = 0 <=> n is an integer. 101 * but in finite precision arithmetic, sin(n*PI) will NEVER 102 * be zero simply because n*PI is a rational number. hence 103 * it failed to work with our newer, more accurate sin() 104 * which uses true pi to do the argument reduction... 105 * temp = sin(pi*arg); 106 */ 107 temp = drem(arg, 1.e0); 108 if(temp == 0.) { 109 errno = EDOM; 110 #ifdef VAX 111 return (NaN); 112 #else 113 return(HUGE); 114 #endif 115 } 116 117 temp = drem(arg, 2.e0); 118 if (temp < 0.) 119 temp = -temp; 120 else 121 signgam = -1; 122 return(-log(arg*pos(arg)*sin(pi*temp)/pi)); 123 } 124 125 static double 126 pos(arg) 127 double arg; 128 { 129 double n, d, s; 130 register i; 131 132 if(arg < 2.) return(pos(arg+1.)/arg); 133 if(arg > 3.) return((arg-1.)*pos(arg-1.)); 134 135 s = arg - 2.; 136 for(n=0,d=0,i=N-1; i>=0; i--){ 137 n = n*s + p2[i]; 138 d = d*s + q2[i]; 139 } 140 return(n/d); 141 } 142