/*
 *  Draco Econometrics
 *  Copyright (C) 2006-2007 Approximatrix, LLC
 *  http://www.approximatrix.com
 *  
 *  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 3 of the License, or
 *  (at your option) 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, see <http://www.gnu.org/licenses/>.
 */

package com.approximatrix.draco.math.functions;

/** Approximations of Gamma
 * 
 * @author Jeffrey Armstrong
 * @see http://www.luschny.de/math/factorial/approx/SimpleCases.html
 */
public class Gamma {
	
	public static final int RAMANUJAN5 = 0;
	public static final int WEHMEIER4 = 1;
	public static final int GOSPER3 = 2;
	public static final int JTEM_LIB = 3;
	
	private static int gamma_technique = JTEM_LIB;
	private static int incomplete_gamma_technique = JTEM_LIB;
	
	public static void SetGammaTechnique(int i) {
		gamma_technique = i;
	}
	
	public static void SetIncompleteGammaTechnique(int i) {
		incomplete_gamma_technique = i;
	}
	
	public static double Gamma(double x) {
		switch(gamma_technique) {
			case RAMANUJAN5:
				return ramanujan5(x);
			case WEHMEIER4:
				return wehmeier4(x);
			case GOSPER3:
				return gosper3(x);
			case JTEM_LIB:
				return JTEMGamma.gamma(x);
		}
		return 0.0;
	}
	
	public static double IncompleteGamma(double a, double x) {
		switch(incomplete_gamma_technique) {
			case JTEM_LIB:
				return JTEMGamma.gamma(a,x,Gamma(a));
		}
		return 0.0;
	}
	
	private static double ramanujan5(double x) {
		double N = 2.0*x;
		return kern2(x)*Math.pow((1.0+1.0/N*(1.0+1.0/(2.0*N)*(1.0+1.0/(15.0*N)*(1.0-1.0/(4.0*N)*(11.0-79.0/(7.0*N)))))),(1.0/6.0));
	}
	
	private static double wehmeier4(double x) {
		return kern1(x)*Math.sqrt(x+1.0/6.0+(1.0/72.0)/x-(31.0/6480.0)/(x*x)-(139.0/155520.0)/(x*x*x));
	}
	
	private static double gosper3(double x) {
		return kern1(x)*Math.sqrt(x+1.0/6.0)*(1.0+(1.0/144.0)/(x*x)-(23.0/6480.0)/(x*x*x));
	}
	
	private static double kern1(double x) {
		return Math.sqrt(2.0*Math.PI)*Math.exp(x*(Math.log(x)-1.0));
	}
	
	private static double kern2(double x) {
		return kern1(x)*Math.sqrt(x);
	}
}
