http://mallet.cs.umass.edu/より
class Digamma {
public static final double EULER_MASCHERONI = -0.5772156649015328606065121;
public static final double DIGAMMA_COEF_1 = 1.0 / 12;
public static final double DIGAMMA_COEF_2 = 1.0 / 120;
public static final double DIGAMMA_COEF_3 = 1.0 / 252;
public static final double DIGAMMA_COEF_4 = 1.0 / 240;
public static final double DIGAMMA_COEF_5 = 1.0 / 132;
public static final double DIGAMMA_COEF_6 = 691.0 / 32760;
public static final double DIGAMMA_COEF_7 = 1.0 / 12;
public static final double DIGAMMA_LARGE = 9.5;
public static final double DIGAMMA_SMALL = .000001;
public static double digamma(double z) {
double psi = 0;
if (z < DIGAMMA_SMALL) {
psi = EULER_MASCHERONI - (1 / z);
return psi;
}
while (z < DIGAMMA_LARGE) {
psi -= 1 / z;
z++;
}
double invZ = 1 / z;
double invZSquared = invZ * invZ;
psi += Math.log(z)
- .5
* invZ
- invZSquared
* (DIGAMMA_COEF_1 - invZSquared
* (DIGAMMA_COEF_2 - invZSquared
* (DIGAMMA_COEF_3 - invZSquared
* (DIGAMMA_COEF_4 - invZSquared
* (DIGAMMA_COEF_5 - invZSquared
* (DIGAMMA_COEF_6 - invZSquared * DIGAMMA_COEF_7))))));
return psi;
}
}