ガンマ関数

/****************************/
/* ガンマ関数の計算         */
/*      coded by Y.Suganuma */
/****************************/
import java.io.*;

public class Test {
	public static void main(String args[]) throws IOException
	{
		BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
		double x, y;
		int ind[] = new int [1];

		System.out.print("data? ");
		x = Double.parseDouble(in.readLine());

		y = gamma(x, ind);

		System.out.println("   x " + x + " gamma(x) " + y + " ind " + ind[0]);
	}

	/****************************************/
	/* Γ（ｘ）の計算（ガンマ関数，近似式） */
	/*      ier : =0 : normal               */
	/*            =-1 : x=-n (n=0,1,2,･･･)  */
	/*      return : 結果                   */
	/****************************************/
	static double gamma(double x, int ier[])
	{
		double err, g, s, t, v, w, y;
		int k;

		ier[0] = 0;

		if (x > 5.0) {
			v = 1.0 / x;
			s = ((((((-0.000592166437354 * v + 0.0000697281375837) * v +
                0.00078403922172) * v - 0.000229472093621) * v -
                0.00268132716049) * v + 0.00347222222222) * v +
                0.0833333333333) * v + 1.0;
			g = 2.506628274631001 * Math.exp(-x) * Math.pow(x,x-0.5) * s;
		}

		else {

			err = 1.0e-20;
			w   = x;
			t   = 1.0;

			if (x < 1.5) {

				if (x < err) {
					k = (int)x;
					y = (double)k - x;
					if (Math.abs(y) < err || Math.abs(1.0-y) < err)
						ier[0] = -1;
				}

				if (ier[0] == 0) {
					while (w < 1.5) {
						t /= w;
						w += 1.0;
					}
				}
			}

			else {
				if (w > 2.5) {
					while (w > 2.5) {
						w -= 1.0;
						t *= w;
					}
				}
			}

			w -= 2.0;
			g  = (((((((0.0021385778 * w - 0.0034961289) * w + 
                 0.0122995771) * w - 0.00012513767) * w + 0.0740648982) * w +
                 0.0815652323) * w + 0.411849671) * w + 0.422784604) * w +
                 0.999999926;
			g *= t;
		}

		return g;
	}
}