comparison src/samer/functions/Gamma.java @ 0:bf79fb79ee13

Initial Mercurial check in.

+// JEM - Java tools for Experimental Mathematics
+// Copyright (C) 2001 JEM-Group
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+package samer.functions;
+
+/** 
+ * Implementation of the real (incomplete) Gamma function.
+ * <p>
+ * This class provides an implementaion of the real Gamma function
+ * <p align=center>
+ * <code>&Gamma;( z ) = &int;<sub>0</sub><sup>&infin;</sup> t<sup>z-1</sup>e<sup>-t</sup>dt</code>
+ * </p>
+ * and its sibling, the <i>incomplete</i> Gamma function:
+ * <p align=center>
+ * <code>&Gamma;( a, x ) = &int;<sub>0</sub><sup>x</sup>t<sup>a-1</sup>e<sup>-t</sup>dt</code>
+ * </p><p>
+ * The Implementation follows the ideas of C.Lanczos and their realizations
+ * presented in <em>Numerical Recipes in C</em>.
+ * </p><p> 
+ * Further reading:
+ * <table>
+ * <tr><td> - </td><td>C.Lanczos, 1964, SIAM Journal of Numerical Analysis, ser. B, vol. 1,
+ * pp. 86-96.</td></tr>
+ * <tr><td> - </td><td>W.H.Press, S.A.Teukolsky, W.T.Vetterling, B.P.Flannery.
+ * <em>Numerical Recipes in C</em>. Cambridge University Press, 1992. </td></tr>
+ * </table>
+ *</p>
+ */
+public class Gamma {
+
+    private static final double a = Math.sqrt( 2*Math.PI ) *    1.000000000190015;
+    private static final double b = Math.sqrt( 2*Math.PI ) *   76.18009172947146;
+    private static final double c = Math.sqrt( 2*Math.PI ) *  -86.50532032941677;
+    private static final double d = Math.sqrt( 2*Math.PI ) *   24.01409824083091;
+    private static final double e = Math.sqrt( 2*Math.PI ) *   -1.231739572450155;
+    private static final double f = Math.sqrt( 2*Math.PI ) *    0.1208650973866179e-2;
+    private static final double g = Math.sqrt( 2*Math.PI ) *   -0.5395239384953e-5;
+
+    private Gamma() {}
+
+    static double gamma_( double z ) {
+
+	if( z <= 0 )
+	    throw new IllegalArgumentException("argument must be positive");
+
+	double s = Math.exp( Math.log( z + 5.5 ) * ( z + 0.5 ) - ( z + 5.5 ) );
+	
+	return s * ( a + 
+		     b/(z+1) + c/(z+2) + d/(z+3) + 
+		     e/(z+4) + f/(z+5) + g/(z+6)  ) / z;
+    }
+
+    /**
+     * Retunrs the gamma function at <code>z</code>.
+     * @param z a real number
+     * @return <code> &Gamma;( z ) = &int;<sub>0</sub><sup>&infin;</sup> 
+     * t<sup>z-1</sup>e<sup>-t</sup>dt</code>
+     */
+    public static double gamma( double z ) {
+
+	return z <= 0.5 ? Math.PI / gamma_ ( 1 - z ) / Math.sin( Math.PI * z ) : gamma_ ( z );
+    }
+
+    static double logOfGamma_( double z ) {
+
+	if( z <= 0 )
+	    throw new IllegalArgumentException("argument must be positive");
+
+	double s = Math.log( z + 5.5 ) * ( z + 0.5 ) - ( z + 5.5 );
+	
+	return s + Math.log(  ( a + 
+				b/(z+1) + c/(z+2) + d/(z+3) + 
+				e/(z+4) + f/(z+5) + g/(z+6)  ) / z );
+    }
+
+
+    /**
+     * Retunrs the logarithm of the gamma function at <code>z</code>.
+     * @param z a real number
+     * @return <code>log( &Gamma;( z ) ) = log ( &int;<sub>0</sub><sup>&infin;</sup> 
+     * t<sup>z-1</sup>e<sup>-t</sup>dt )</code>
+     */
+    public static double logOfGamma( double z ) {
+
+	return z <= 0.5 
+	    ? Math.log( Math.PI / Math.sin( Math.PI * z ) ) - logOfGamma_ ( 1 - z )  
+	    : logOfGamma_ ( z );
+    
+    }
+
+
+    private static double gamma( double a, double x, double gammaOfA, boolean computeGammaOfA ) {
+
+	if( a<=0 )
+	    throw new IllegalArgumentException( "a is not positive" );
+
+	if( x<0 )
+	    throw new IllegalArgumentException( "x is negative" );
+
+	if( x < a + 1 ) {
+	    if( computeGammaOfA )
+		gammaOfA = gamma ( a );
+
+	    return computeViaSeries( a, x, gammaOfA );
+	
+	} else
+	    return computeViaContinuedFraction( a, x );
+    }
+
+    /**
+     * Returns the incomplete gamma function at <code>a, x</code>.
+     * The algorithms for the incomplete gamma function needs
+     * to compute <code>&Gamma;(a)</code>; providing this value
+     * will therefor lead to an optimization if you need
+     * this value also for a diffrent porpuse.
+     * @param a a positive real number
+     * @param x a real number greater than <code>a</code>
+     * @return <code>&Gamma;( a, x ) = &int;<sub>0</sub><sup>&infin;</sup> 
+     * t<sup>z-1</sup>e<sup>-t</sup>dt</code>
+     */
+    public static double gamma( double a, double x, double gammaOfA ) {
+	return gamma ( a, x, gammaOfA, false );
+    }
+
+
+    /**
+     * Returns the incomplete gamma function at <code>a, x</code>.
+     * @param a a positive real number
+     * @param x a real number greater than <code>a</code>
+     * @return <code>&Gamma;( a, x ) = &int;<sub>0</sub><sup>&infin;</sup> 
+     * t<sup>z-1</sup>e<sup>-t</sup>dt</code>
+     */
+    public static double gamma( double a, double x ) {
+	return gamma ( a, x, 123, true );
+    }
+    
+    static int MAX_NUMBER_OF_ITERATIONS = 200;
+
+    static double EPS = 1e-16;
+
+    static double computeViaSeries( double a, double x, double gammaOfA ) {
+
+	double factor = a;
+	double delta  = 1 / a;
+	double sum    = delta;
+	
+	int counter   = 0;
+
+	while( Math.abs(delta) > Math.abs( sum ) * EPS ) {
+	    
+	    delta *= x/++factor;
+
+	    sum   += delta;
+
+	    if( counter > MAX_NUMBER_OF_ITERATIONS )
+		throw new RuntimeException( "exeeded max number of iterations="+
+					    MAX_NUMBER_OF_ITERATIONS );
+	}
+
+	return gammaOfA  - sum * Math.exp( a*Math.log( x ) - x );
+    }
+
+    static double MIN_DP = 1e-30;
+
+    static double computeViaContinuedFraction( double a, double x ) {
+
+	    double b = x + 1 - a;
+	    double c = 1/MIN_DP;
+	    double d = 1/b;
+	    double h = d;
+
+	    for( int i=1; i<MAX_NUMBER_OF_ITERATIONS; i++ ) {
+
+		double an = -i * (i-a);
+
+		b += 2;
+
+		d  = an*d + b;
+		c  = an/c + b;
+
+		if( Math.abs( c ) < MIN_DP )
+		    c = MIN_DP;
+		if( Math.abs( d ) < MIN_DP )
+		    d = MIN_DP;
+
+		d = 1 / d;
+
+		double delta = c * d;
+
+		h *= delta;
+
+		if( Math.abs( delta - 1 ) < EPS )
+		    return h * Math.exp( a * Math.log( x ) - x );
+	    }
+
+	    throw new RuntimeException( "exeeded max number of iterations="+
+					MAX_NUMBER_OF_ITERATIONS  );
+    }
+    
+}
+