Merge branch 'feature-MATH-1318' into develop

Completes issue MATH-1318 (see JIRA).

src/changes/changes.xml

@@ -57,6 +57,10 @@ If the output is not quite correct, check for invisible trailing spaces!
       <action dev="erans" type="add" issue="MATH-1015" due-to="Thomas Neidhart">
         Gauss-Laguerre quadrature.
       </action>
+      <action dev="erans" type="update" issue="MATH-1318" due-to="Eric Prescott-Gagnon">
+        "o.a.c.m.special.Gamma.digamma": Improved performance (through the use of
+        the reflection formula for negative arguments).
+      </action>
       <action dev="erans" type="add" issue="MATH-1350" due-to="Rob Tompkins">
         Improved code coverage (unit tests).
       </action>

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;
     }
 
     /**
@@ -472,7 +480,7 @@ public class Gamma {
      * of digamma.
      *
      * @param x Argument.
-     * @return trigamma(x) to within 10-8 relative or absolute error whichever is smaller
+     * @return trigamma(x) to within 10^-8 relative or absolute error whichever is smaller
      * @see <a href="http://en.wikipedia.org/wiki/Trigamma_function">Trigamma</a>
      * @see Gamma#digamma(double)
      * @since 2.0

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