view src/samer/functions/Gamma.java @ 5:b67a33c44de7

Remove some crap, etc

// 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  );
    }
    
}