projects / unix-history / blame
| Commit | Line | Data |
|---|---|---|
| 9cc34517 BJ |
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 | ||
| 24 | #define M 6 | |
| 25 | #define N 8 | |
| 26 | static double p1[] = { | |
| 27 | 0.83333333333333101837e-1, | |
| 28 | -.277777777735865004e-2, | |
| 29 | 0.793650576493454e-3, | |
| 30 | -.5951896861197e-3, | |
| 31 | 0.83645878922e-3, | |
| 32 | -.1633436431e-2, | |
| 33 | }; | |
| 34 | static double p2[] = { | |
| 35 | -.42353689509744089647e5, | |
| 36 | -.20886861789269887364e5, | |
| 37 | -.87627102978521489560e4, | |
| 38 | -.20085274013072791214e4, | |
| 39 | -.43933044406002567613e3, | |
| 40 | -.50108693752970953015e2, | |
| 41 | -.67449507245925289918e1, | |
| 42 | 0.0, | |
| 43 | }; | |
| 44 | static double q2[] = { | |
| 45 | -.42353689509744090010e5, | |
| 46 | -.29803853309256649932e4, | |
| 47 | 0.99403074150827709015e4, | |
| 48 | -.15286072737795220248e4, | |
| 49 | -.49902852662143904834e3, | |
| 50 | 0.18949823415702801641e3, | |
| 51 | -.23081551524580124562e2, | |
| 52 | 0.10000000000000000000e1, | |
| 53 | }; | |
| 54 | ||
| 55 | double | |
| 56 | gamma(arg) | |
| 57 | double arg; | |
| 58 | { | |
| 59 | double log(), pos(), neg(), asym(); | |
| 60 | ||
| 61 | signgam = 1.; | |
| 62 | if(arg <= 0.) return(neg(arg)); | |
| 63 | if(arg > 8.) return(asym(arg)); | |
| 64 | return(log(pos(arg))); | |
| 65 | } | |
| 66 | ||
| 67 | static double | |
| 68 | asym(arg) | |
| 69 | double arg; | |
| 70 | { | |
| 71 | double log(); | |
| 72 | double n, argsq; | |
| 73 | int i; | |
| 74 | ||
| 75 | argsq = 1./(arg*arg); | |
| 76 | for(n=0,i=M-1; i>=0; i--){ | |
| 77 | n = n*argsq + p1[i]; | |
| 78 | } | |
| 79 | return((arg-.5)*log(arg) - arg + goobie + n/arg); | |
| 80 | } | |
| 81 | ||
| 82 | static double | |
| 83 | neg(arg) | |
| 84 | double arg; | |
| 85 | { | |
| 86 | double temp; | |
| 87 | double log(), sin(), pos(); | |
| 88 | ||
| 89 | arg = -arg; | |
| 90 | temp = sin(pi*arg); | |
| 91 | if(temp == 0.) { | |
| 92 | errno = EDOM; | |
| 93 | return(HUGE); | |
| 94 | } | |
| 95 | if(temp < 0.) temp = -temp; | |
| 96 | else signgam = -1; | |
| 97 | return(-log(arg*pos(arg)*temp/pi)); | |
| 98 | } | |
| 99 | ||
| 100 | static double | |
| 101 | pos(arg) | |
| 102 | double arg; | |
| 103 | { | |
| 104 | double n, d, s; | |
| 105 | register i; | |
| 106 | ||
| 107 | if(arg < 2.) return(pos(arg+1.)/arg); | |
| 108 | if(arg > 3.) return((arg-1.)*pos(arg-1.)); | |
| 109 | ||
| 110 | s = arg - 2.; | |
| 111 | for(n=0,d=0,i=N-1; i>=0; i--){ | |
| 112 | n = n*s + p2[i]; | |
| 113 | d = d*s + q2[i]; | |
| 114 | } | |
| 115 | return(n/d); | |
| 116 | } |