ダイガンマ関数

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;
  }
}