ref: 1b44caa0c40d7b722a7d9edd0a6eaefe1616c4b0
dir: /sys/src/ape/lib/ap/math/gamma.c/
#include <math.h> #include <errno.h> /* C program for floating point log gamma function gamma(x) computes the log of the absolute value of the gamma function. The sign of the gamma function is returned in the external quantity signgam. The coefficients for expansion around zero are #5243 from Hart & Cheney; for expansion around infinity they are #5404. Calls log and sin. */ int signgam; static double goobie = 0.9189385332046727417803297; static double pi = 3.1415926535897932384626434; #define M 6 #define N 8 static double p1[] = { 0.83333333333333101837e-1, -.277777777735865004e-2, 0.793650576493454e-3, -.5951896861197e-3, 0.83645878922e-3, -.1633436431e-2, }; static double p2[] = { -.42353689509744089647e5, -.20886861789269887364e5, -.87627102978521489560e4, -.20085274013072791214e4, -.43933044406002567613e3, -.50108693752970953015e2, -.67449507245925289918e1, 0.0, }; static double q2[] = { -.42353689509744090010e5, -.29803853309256649932e4, 0.99403074150827709015e4, -.15286072737795220248e4, -.49902852662143904834e3, 0.18949823415702801641e3, -.23081551524580124562e2, 0.10000000000000000000e1, }; static double asym(double arg) { double n, argsq; int i; argsq = 1 / (arg*arg); n = 0; for(i=M-1; i>=0; i--) n = n*argsq + p1[i]; return (arg-.5)*log(arg) - arg + goobie + n/arg; } static double pos(double arg) { double n, d, s; int i; if(arg < 2) return pos(arg+1)/arg; if(arg > 3) return (arg-1)*pos(arg-1); s = arg - 2; n = 0; d = 0; for(i=N-1; i>=0; i--){ n = n*s + p2[i]; d = d*s + q2[i]; } return n/d; } static double neg(double arg) { double temp; arg = -arg; temp = sin(pi*arg); if(temp == 0) { errno = EDOM; return HUGE_VAL; } if(temp < 0) temp = -temp; else signgam = -1; return -log(arg*pos(arg)*temp/pi); } double gamma(double arg) { signgam = 1; if(arg <= 0) return neg(arg); if(arg > 8) return asym(arg); return log(pos(arg)); }