diff --git a/src/main/java/org/apache/commons/math4/special/Gamma.java b/src/main/java/org/apache/commons/math4/special/Gamma.java index 70aaec3cb..c59049cc9 100644 --- a/src/main/java/org/apache/commons/math4/special/Gamma.java +++ b/src/main/java/org/apache/commons/math4/special/Gamma.java @@ -427,18 +427,16 @@ public class Gamma { *
Computes the digamma function of x.
* *This is an independently written implementation of the algorithm described in - * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.
+ * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976. + * A reflection formula (https://en.wikipedia.org/wiki/Digamma_function#Reflection_formula) + * is incorporated to improve performance on negative values. * *Some of the constants have been changed to increase accuracy at the moderate expense - * of run-time. The result should be accurate to within 10^-8 absolute tolerance for - * x >= 10^-5 and within 10^-8 relative tolerance for x > 0.
- * - *Performance for large negative values of x will be quite expensive (proportional to - * |x|). Accuracy for negative values of x should be about 10^-8 absolute for results - * less than 10^5 and 10^-8 relative for results larger than that.
+ * of run-time. The result should be accurate to within 10^-8 relative tolerance for + * 0 < x < 10^-5 and within 10^-8 absolute tolerance otherwise. * * @param x Argument. - * @return digamma(x) to within 10-8 relative or absolute error whichever is smaller. + * @return digamma(x) to within 10-8 relative or absolute error whichever is larger. * @see Digamma * @see Bernardo's original article * @since 2.0 @@ -448,22 +446,32 @@ public class Gamma { return x; } + double digamma = 0.0; + if (x < 0) { + // use reflection formula to fall back into positive values + digamma -= FastMath.PI / FastMath.tan(FastMath.PI * x); + x = 1 - x; + } + if (x > 0 && x <= S_LIMIT) { // use method 5 from Bernardo AS103 // accurate to O(x) - return -GAMMA - 1 / x; + return digamma -GAMMA - 1 / x; } - if (x >= C_LIMIT) { - // use method 4 (accurate to O(1/x^8) - double inv = 1 / (x * x); - // 1 1 1 1 - // log(x) - --- - ------ + ------- - ------- - // 2 x 12 x^2 120 x^4 252 x^6 - return FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252)); + while (x < C_LIMIT) { + digamma -= 1.0 / x; + x += 1.0; } - return digamma(x + 1) - 1 / x; + // use method 4 (accurate to O(1/x^8) + double inv = 1 / (x * x); + // 1 1 1 1 + // log(x) - --- - ------ + ------- - ------- + // 2 x 12 x^2 120 x^4 252 x^6 + digamma += FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252)); + + return digamma; } /** diff --git a/src/test/java/org/apache/commons/math4/special/GammaTest.java b/src/test/java/org/apache/commons/math4/special/GammaTest.java index 51982a3a9..c3205267b 100644 --- a/src/test/java/org/apache/commons/math4/special/GammaTest.java +++ b/src/test/java/org/apache/commons/math4/special/GammaTest.java @@ -107,6 +107,7 @@ public class GammaTest { Assert.assertEquals(-100.56088545786867450, Gamma.digamma(0.01), eps); Assert.assertEquals(-4.0390398965921882955, Gamma.digamma(-0.8), eps); Assert.assertEquals(4.2003210041401844726, Gamma.digamma(-6.3), eps); + Assert.assertEquals(-3.110625123035E-5, Gamma.digamma(1.4616), eps); } @Test