Index: src/java/org/apache/commons/math/special/Gamma.java
===================================================================
--- src/java/org/apache/commons/math/special/Gamma.java	(revision 778092)
+++ src/java/org/apache/commons/math/special/Gamma.java	Sat May 23 23:34:34 PDT 2009
@@ -271,26 +271,27 @@
 
     // limits for switching algorithm in digamma
     /** C limit */
-    private static final double C_LIMIT = 49;
-    /** S limit */
-    private static final double S_LIMIT = 1e-5;
+     private static final double C_LIMIT = 49;
+     /** S limit */
+     private static final double S_LIMIT = 1e-5;
 
     /**
-     * <p>Computes the <a href="http://en.wikipedia.org/wiki/Digamma_function">digamma function</a>
-     * using the algorithm defined in <br/>
-     * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.</p>
-     * 
-     * <p>Some of the constants have been changed to increase accuracy at the moderate expense
-     * of run-time performance.  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
+     * Computes the digamma function of x.
+     * This is an independently written implementation of the algorithm described in
+     * <p/>
+     * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976
+     * <p/>
+     * I have changed some of the constants 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/>
+     * 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.
-     * 
-     * @param x argument
-     * @return value of the digamma function
-     * @since 2.0
+     * @param x  The argument.
+     * @return   digamma(x) to within 10-8 relative or absolute error whichever is smaller.
+     * @see <a href="http://en.wikipedia.org/wiki/Digamma_function"> Digamma at wikipedia </a>
+     * @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf"> Bernardo's original article </a>
      */
     public static double digamma(double x) {
         if (x > 0 && x <= S_LIMIT) {
@@ -303,11 +304,39 @@
             // use method 4 (accurate to O(1/x^8)
             double inv = 1 / (x * x);
             //            1       1        1         1
-            // log(x) -  --- - ------ - ------- - -------
+            // 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;
     }
+
+    /**
+     * Computes the trigamma function of x.  This function is derived by taking the derivative of
+     * the implementation of digamma.
+     * <p/>
+     * 
+     * <p/>
+     * @param x  The argument.
+     * @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 at wikipedia </a>
+     * @see Gamma#digamma(double)
+     */
+    public static double trigamma(double x) {
+        if (x > 0 && x <= S_LIMIT) {
+            return 1 / (x * x);
-}
+        }
+
+        if (x >= C_LIMIT) {
+            double inv = 1 / (x * x);
+            //  1    1      1       1       1
+            //  - + ---- + ---- - ----- + -----
+            //  x      2      3       5       7
+            //      2 x    6 x    30 x    42 x
+            return 1 / x + inv / 2 + inv / x * (1.0 / 6 - inv * (1.0 / 30 + inv / 42));
+        }
+
+        return trigamma(x + 1) + 1 / (x * x);
+    }
+}
Index: src/test/org/apache/commons/math/special/GammaTest.java
===================================================================
--- src/test/org/apache/commons/math/special/GammaTest.java	(revision 778091)
+++ src/test/org/apache/commons/math/special/GammaTest.java	Sat May 23 23:36:24 PDT 2009
@@ -119,6 +119,31 @@
         }
     }
 
+    public void testTrigamma() {
+        double eps = 1e-8;
+        // computed using webMathematica.  For example, to compute trigamma($i) = Polygamma(1, $i), use
+        //
+        // http://functions.wolfram.com/webMathematica/Evaluated.jsp?name=PolyGamma2&plottype=0&vars={%221%22,%22$i%22}&digits=20
+        double[] data = {
+                1e-4, 1.0000000164469368793e8,
+                1e-3, 1.0000016425331958690e6,
+                1e-2, 10001.621213528313220,
+                1e-1, 101.43329915079275882,
+                1, 1.6449340668482264365,
+                2, 0.64493406684822643647,
+                3, 0.39493406684822643647,
+                4, 0.28382295573711532536,
+                5, 0.22132295573711532536,
+                10, 0.10516633568168574612,
+                20, 0.051270822935203119832,
+                50, 0.020201333226697125806,
+                100, 0.010050166663333571395
+        };
+        for (int i = data.length - 2; i >= 0; i -= 2) {
+            assertEquals(String.format("trigamma %.0f", data[i]), data[i + 1], Gamma.trigamma(data[i]), eps);
+        }
+    }
+
     private void checkRelativeError(String msg, double expected, double actual, double tolerance) {
         assertEquals(msg, expected, actual, Math.abs(tolerance * actual));
     }