projects / unix-history / blame
| Commit | Line | Data |
|---|---|---|
| cf89f98e SL |
1 | /* @(#)gamma.c 4.1 %G% */ |
| 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 | } |