Index: src/java/org/apache/commons/math/special/Gamma.java
===================================================================
--- src/java/org/apache/commons/math/special/Gamma.java	(revision 762121)
+++ src/java/org/apache/commons/math/special/Gamma.java	Sat May 23 21:35:09 PDT 2009
@@ -33,6 +33,8 @@
     /** Serializable version identifier */
     private static final long serialVersionUID = -6587513359895466954L;
 
+    public static final double GAMMA = 0.577215664901532860606512090082;
+
     /** Maximum allowed numerical error. */
     private static final double DEFAULT_EPSILON = 10e-15;
 
@@ -59,7 +61,7 @@
     /** Avoid repeated computation of log of 2 PI in logGamma */
     private static final double HALF_LOG_2_PI = 0.5 * Math.log(2.0 * Math.PI);
 
-    
+
     /**
      * Default constructor.  Prohibit instantiation.
      */
@@ -261,4 +263,28 @@
 
         return ret;
     }
+
+
+    // limits for switching algorithm in digamma
+    private static final double C_LIMIT = 49;
+    private static final double S_LIMIT = 1e-5;
+
+    public static double digamma(double x) {
+        if (x > 0 && x <= S_LIMIT) {
+            // use method 5 from Bernardo AS103
+            // accurate to O(x)
+            return -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 Math.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252));
+        }
+
+        return digamma(x + 1) - 1 / x;
+    }
+}
Index: src/test/org/apache/commons/math/special/GammaTest.java
===================================================================
--- src/test/org/apache/commons/math/special/GammaTest.java	(revision 670469)
+++ src/test/org/apache/commons/math/special/GammaTest.java	Sat May 23 21:36:32 PDT 2009
@@ -20,6 +20,7 @@
 import org.apache.commons.math.TestUtils;
 
 import junit.framework.TestCase;
+import junit.framework.Assert;
 
 /**
  * @version $Revision: 670469 $ $Date: 2008-06-23 01:01:38 -0700 (Mon, 23 Jun 2008) $
@@ -27,7 +28,7 @@
 public class GammaTest extends TestCase {
     /**
      * Constructor for BetaTest.
-     * @param name
+     * @param name  Name of the test.
      */
     public GammaTest(String name) {
         super(name);
@@ -92,4 +93,39 @@
     public void testLogGammaPositive() {
         testLogGamma(0.6931471805599457, 3.0);
     }
+
+    public void testDigammaLargeArgs() {
+        double eps = 1e-8;
+        Assert.assertEquals(4.6001618527380874002, Gamma.digamma(100), eps);
+        Assert.assertEquals(3.9019896734278921970, Gamma.digamma(50), eps);
+        Assert.assertEquals(2.9705239922421490509, Gamma.digamma(20), eps);
+        Assert.assertEquals(2.9958363947076465821, Gamma.digamma(20.5), eps);
+        Assert.assertEquals(2.2622143570941481605, Gamma.digamma(10.1), eps);
+        Assert.assertEquals(2.1168588189004379233, Gamma.digamma(8.8), eps);
+        Assert.assertEquals(1.8727843350984671394, Gamma.digamma(7), eps);
+        Assert.assertEquals(0.42278433509846713939, Gamma.digamma(2), eps);
+        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);
-}
+    }
+
+    public void testDigammaSmallArgs() {
+        // values for negative powers of 10 from 1 to 30 as computed by webMathematica with 20 digits
+        // see functions.wolfram.com
+        double[] expected = {-10.423754940411076795, -100.56088545786867450, -1000.5755719318103005,
+                -10000.577051183514335, -100000.57719921568107, -1.0000005772140199687e6, -1.0000000577215500408e7,
+                -1.0000000057721564845e8, -1.0000000005772156633e9, -1.0000000000577215665e10, -1.0000000000057721566e11,
+                -1.0000000000005772157e12, -1.0000000000000577216e13, -1.0000000000000057722e14, -1.0000000000000005772e15, -1e+16,
+                -1e+17, -1e+18, -1e+19, -1e+20, -1e+21, -1e+22, -1e+23, -1e+24, -1e+25, -1e+26,
+                -1e+27, -1e+28, -1e+29, -1e+30};
+        for (double n = 1; n < 30; n++) {
+            System.out.printf("10^-n = %.10f\n", Math.pow(10, -n));
+            checkRelativeError(String.format("Test %.0f: ", n), expected[(int) (n - 1)], Gamma.digamma(Math.pow(10.0, -n)), 1e-8);
+        }
+    }
+
+    private void checkRelativeError(String msg, double expected, double actual, double tolerance) {
+        System.out.printf("%s = %.2f (%s)\n", msg, Math.abs(expected - actual) / actual, Math.abs(expected - actual) > Math.abs(tolerance * actual));
+        Assert.assertEquals(msg, expected, actual, Math.abs(tolerance * actual));
+    }
+}