/* digamma.c
*
* Mark Johnson, 2nd September 2007
*
* Computes the Ψ(x) or digamma function, i.e., the derivative of the
* log gamma function, using a series expansion.
*
* Warning: I'm not a numerical analyst, so I may have made errors here!
*
* The parameters of the series were computed using the Maple symbolic
* algebra program as follows:
*
* series(Psi(x+1/2), x=infinity, 21);
*
* which produces:
*
* ln(x)+1/(24*x^2)-7/960/x^4+31/8064/x^6-127/30720/x^8+511/67584/x^10-1414477/67092480/x^12+8191/98304/x^14-118518239/267386880/x^16+5749691557/1882718208/x^18-91546277357/3460300800/x^20+O(1/(x^21))
*
* It looks as if the terms in this expansion *diverge* as the powers
* get larger. However, for large x, the x^-n term will dominate.
*
* I used Maple to examine the difference between this series and
* Digamma(x+1/2) over the range 7 < x < 20, and discovered that the
* difference is less that 1e-8 if the terms up to x^-8 are included.
* This determined the power used in the code here. Of course,
* Maple uses some kind of approximation to calculate Digamma,
* so all I've really done here is find the smallest power that produces
* the Maple approximation; still, that should be good enough for our
* purposes.
*
* This expansion is accurate for x > 7; we use the recurrence
*
* digamma(x) = digamma(x+1) - 1/x
*
* to make x larger than 7.
*/
#include
#include
double digamma(double x) {
double result = 0, xx, xx2, xx4;
assert(x > 0);
for ( ; x < 7; ++x)
result -= 1/x;
x -= 1.0/2.0;
xx = 1.0/x;
xx2 = xx*xx;
xx4 = xx2*xx2;
result += log(x)+(1./24.)*xx2-(7.0/960.0)*xx4+(31.0/8064.0)*xx4*xx2-(127.0/30720.0)*xx4*xx4;
return result;
}
#ifdef TEST
#include
#include
int main(int argc, char** argv) {
double x;
for (x = 0.1; x < 10; x += 0.1)
printf("digamma(%g) = %g, exp(digamma(%g)) = %g\n", x, digamma(x), x, exp(digamma(x)));
exit(EXIT_SUCCESS);
}
#endif /* TEST */