Return-Path: Delivered-To: apmail-commons-commits-archive@minotaur.apache.org Received: (qmail 83324 invoked from network); 25 May 2009 16:07:06 -0000 Received: from hermes.apache.org (HELO mail.apache.org) (140.211.11.3) by minotaur.apache.org with SMTP; 25 May 2009 16:07:06 -0000 Received: (qmail 34938 invoked by uid 500); 25 May 2009 13:20:39 -0000 Delivered-To: apmail-commons-commits-archive@commons.apache.org Received: (qmail 34853 invoked by uid 500); 25 May 2009 13:20:39 -0000 Mailing-List: contact commits-help@commons.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: dev@commons.apache.org Delivered-To: mailing list commits@commons.apache.org Received: (qmail 34837 invoked by uid 99); 25 May 2009 13:20:39 -0000 Received: from athena.apache.org (HELO athena.apache.org) (140.211.11.136) by apache.org (qpsmtpd/0.29) with ESMTP; Mon, 25 May 2009 13:20:39 +0000 X-ASF-Spam-Status: No, hits=-2000.0 required=10.0 tests=ALL_TRUSTED X-Spam-Check-By: apache.org Received: from [140.211.11.4] (HELO eris.apache.org) (140.211.11.4) by apache.org (qpsmtpd/0.29) with ESMTP; Mon, 25 May 2009 13:20:28 +0000 Received: by eris.apache.org (Postfix, from userid 65534) id 0FE3923888BD; Mon, 25 May 2009 13:20:08 +0000 (UTC) Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Subject: svn commit: r778416 - in /commons/proper/math/trunk/src: java/org/apache/commons/math/special/Gamma.java test/org/apache/commons/math/special/GammaTest.java Date: Mon, 25 May 2009 13:20:07 -0000 To: commits@commons.apache.org From: psteitz@apache.org X-Mailer: svnmailer-1.0.8 Message-Id: <20090525132008.0FE3923888BD@eris.apache.org> X-Virus-Checked: Checked by ClamAV on apache.org Author: psteitz Date: Mon May 25 13:20:07 2009 New Revision: 778416 URL: http://svn.apache.org/viewvc?rev=778416&view=rev Log: Added trigamma, javadoc fixes for digamma. JIRA: MATH-267. Patched by Ted Dunning. Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/special/Gamma.java commons/proper/math/trunk/src/test/org/apache/commons/math/special/GammaTest.java Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/special/Gamma.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/special/Gamma.java?rev=778416&r1=778415&r2=778416&view=diff ============================================================================== --- commons/proper/math/trunk/src/java/org/apache/commons/math/special/Gamma.java (original) +++ commons/proper/math/trunk/src/java/org/apache/commons/math/special/Gamma.java Mon May 25 13:20:07 2009 @@ -271,25 +271,28 @@ // 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; /** - *

Computes the digamma function - * using the algorithm defined in
+ *

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.

* *

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 + * 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 + *

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. + * less than 10^5 and 10^-8 relative for results larger than that.

* - * @param x argument - * @return value of the digamma function + * @param x the argument + * @return digamma(x) to within 10-8 relative or absolute error whichever is smaller + * @see Digamma at wikipedia + * @see Bernardo's original article * @since 2.0 */ public static double digamma(double x) { @@ -303,11 +306,38 @@ // 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.

+ * + * @param x the argument + * @return trigamma(x) to within 10-8 relative or absolute error whichever is smaller + * @see Trigamma at wikipedia + * @see Gamma#digamma(double) + * @since 2.0 + */ + 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); + } } Modified: commons/proper/math/trunk/src/test/org/apache/commons/math/special/GammaTest.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/special/GammaTest.java?rev=778416&r1=778415&r2=778416&view=diff ============================================================================== --- commons/proper/math/trunk/src/test/org/apache/commons/math/special/GammaTest.java (original) +++ commons/proper/math/trunk/src/test/org/apache/commons/math/special/GammaTest.java Mon May 25 13:20:07 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)); }