double
digamma(double x) Stolen from Radford Neal's fbm package.
assert x > 0 : x;
double r, f, t;
r = 0;
while (x <= 5) {
r -= 1 / x;
x += 1;
f = 1 / (x * x);
...
double
digamma(double x) digamma
if (x > 6.0) {
double x2 = x * x;
double x4 = x2 * x2;
double x6 = x2 * x4;
double x8 = x4 * x4;
double x10 = x6 * x4;
double x12 = x6 * x6;
double x14 = x10 * x4;
...
double
digamma(double x) Compute the value of digamma function
double y, y2, sum = 0;
for (y = x; y < 10; y++) {
sum -= 1.0 / y;
y2 = 1.0 / (y * y);
sum += Math.log(y) + c1 / y + y2 * (c2 + y2 * (c4 + y2 * c6));
return sum;
double
digammaByDefinition(int d) Compute the digamma function from first principles
if (d < 1) {
return Double.NaN;
double result = 0;
for (int n = d; n > 1; n--) {
result += 1.0 / (double) (n - 1);
result += -EULER_MASCHERONI_CONSTANT;
...
double
digammaDiff(double x, int d) Compute digamma difference
double sum = 0;
int dcutoff = 16;
if (d > dcutoff) {
return (digamma(x + d) - digamma(x));
for (int i = 0; i < d; ++i) {
sum += 1 / (x + i);
return (sum);