@Test public void testLogGamma() { System.out.println("lgamma"); assertTrue(Double.isInfinite(Gamma.lgamma(0))); assertEquals(0.0, Gamma.lgamma(1), 1E-7); assertEquals(0, Gamma.lgamma(2), 1E-7); assertEquals(Math.log(2.0), Gamma.lgamma(3), 1E-7); assertEquals(Math.log(6.0), Gamma.lgamma(4), 1E-7); assertEquals(-0.1207822, Gamma.lgamma(1.5), 1E-7); assertEquals(0.2846829, Gamma.lgamma(2.5), 1E-7); assertEquals(1.200974, Gamma.lgamma(3.5), 1E-6); assertEquals(2.453737, Gamma.lgamma(4.5), 1E-6); }
public static double lgamma(double x) { double xcopy = x; double fg = 0.0; double first = x + LANCZOS_SMALL_GAMMA + 0.5; double second = LANCZOS_COEFF[0]; if (x >= 0.0) { if (x >= 1.0 && x - (int) x == 0.0) { fg = lfactorial((int) x - 1); } else { first -= (x + 0.5) * log(first); for (int i = 1; i <= LANCZOS_N; i++) { second += LANCZOS_COEFF[i] / ++xcopy; } fg = log(sqrt(2.0 * PI) * second / x) - first; } } else { fg = PI / (gamma(1.0 - x) * sin(PI * x)); if (Double.isFinite(fg)) { if (fg < 0) { throw new IllegalArgumentException("The gamma function is negative: " + fg); } else { fg = log(fg); } } } return fg; }
Gamma { public static double lgamma(double x) { double xcopy = x; double fg = 0.0; double first = x + LANCZOS_SMALL_GAMMA + 0.5; double second = LANCZOS_COEFF[0]; if (x >= 0.0) { if (x >= 1.0 && x - (int) x == 0.0) { fg = lfactorial((int) x - 1); } else { first -= (x + 0.5) * log(first); for (int i = 1; i <= LANCZOS_N; i++) { second += LANCZOS_COEFF[i] / ++xcopy; } fg = log(sqrt(2.0 * PI) * second / x) - first; } } else { fg = PI / (gamma(1.0 - x) * sin(PI * x)); if (Double.isFinite(fg)) { if (fg < 0) { throw new IllegalArgumentException("The gamma function is negative: " + fg); } else { fg = log(fg); } } } return fg; } }
Gamma { public static double lgamma(double x) { double xcopy = x; double fg = 0.0; double first = x + LANCZOS_SMALL_GAMMA + 0.5; double second = LANCZOS_COEFF[0]; if (x >= 0.0) { if (x >= 1.0 && x - (int) x == 0.0) { fg = lfactorial((int) x - 1); } else { first -= (x + 0.5) * log(first); for (int i = 1; i <= LANCZOS_N; i++) { second += LANCZOS_COEFF[i] / ++xcopy; } fg = log(sqrt(2.0 * PI) * second / x) - first; } } else { fg = PI / (gamma(1.0 - x) * sin(PI * x)); if (Double.isFinite(fg)) { if (fg < 0) { throw new IllegalArgumentException("The gamma function is negative: " + fg); } else { fg = log(fg); } } } return fg; } private  Gamma(); }
Gamma { public static double lgamma(double x) { double xcopy = x; double fg = 0.0; double first = x + LANCZOS_SMALL_GAMMA + 0.5; double second = LANCZOS_COEFF[0]; if (x >= 0.0) { if (x >= 1.0 && x - (int) x == 0.0) { fg = lfactorial((int) x - 1); } else { first -= (x + 0.5) * log(first); for (int i = 1; i <= LANCZOS_N; i++) { second += LANCZOS_COEFF[i] / ++xcopy; } fg = log(sqrt(2.0 * PI) * second / x) - first; } } else { fg = PI / (gamma(1.0 - x) * sin(PI * x)); if (Double.isFinite(fg)) { if (fg < 0) { throw new IllegalArgumentException("The gamma function is negative: " + fg); } else { fg = log(fg); } } } return fg; } private  Gamma(); static double gamma(double x); static double lgamma(double x); static double regularizedIncompleteGamma(double s, double x); static double regularizedUpperIncompleteGamma(double s, double x); static double digamma(double x); static double inverseRegularizedIncompleteGamma(double a, double p); }
Gamma { public static double lgamma(double x) { double xcopy = x; double fg = 0.0; double first = x + LANCZOS_SMALL_GAMMA + 0.5; double second = LANCZOS_COEFF[0]; if (x >= 0.0) { if (x >= 1.0 && x - (int) x == 0.0) { fg = lfactorial((int) x - 1); } else { first -= (x + 0.5) * log(first); for (int i = 1; i <= LANCZOS_N; i++) { second += LANCZOS_COEFF[i] / ++xcopy; } fg = log(sqrt(2.0 * PI) * second / x) - first; } } else { fg = PI / (gamma(1.0 - x) * sin(PI * x)); if (Double.isFinite(fg)) { if (fg < 0) { throw new IllegalArgumentException("The gamma function is negative: " + fg); } else { fg = log(fg); } } } return fg; } private  Gamma(); static double gamma(double x); static double lgamma(double x); static double regularizedIncompleteGamma(double s, double x); static double regularizedUpperIncompleteGamma(double s, double x); static double digamma(double x); static double inverseRegularizedIncompleteGamma(double a, double p);  }