From 5366158572bd88f18df7883ab823ea62b7701d77 Mon Sep 17 00:00:00 2001
From: Eric Prescott-Gagnon <eric.prescottgagnon@jda.com>
Date: Fri, 8 Jan 2016 14:18:30 -0500
Subject: [PATCH] Improved digamma function, especially for negative values.

---
 .../apache/commons/math4/special/Gamma.java   | 42 +++++++++++--------
 .../commons/math4/special/GammaTest.java      |  1 +
 2 files changed, 26 insertions(+), 17 deletions(-)

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 {
      * <p>Computes the digamma function of x.</p>
      *
      * <p>This is an independently written implementation of the algorithm described in
-     * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.</p>
+     * 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.</p>
      *
      * <p>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.</p>
-     *
-     * <p>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.</p>
+     * 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.</p>
      *
      * @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 <a href="http://en.wikipedia.org/wiki/Digamma_function">Digamma</a>
      * @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf">Bernardo&apos;s original article </a>
      * @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