public static double gammaFunction(double x) {

		double xcopy = x;
		double first = x + lgfGamma + 0.5;
		double second = lgfCoeff[0];
		double fg = 0.0D;

		if (x >= 0.0) {
			first = Math.pow(first, x + 0.5) * Math.exp(-first);
			for (int i = 1; i <= lgfN; i++)
				second += lgfCoeff[i] / ++xcopy;
			fg = first * Math.sqrt(2.0 * Math.PI) * second / x;
		} else {
			fg = -Math.PI / (x * Stat.gamma(-x) * Math.sin(Math.PI * x));
		}
		return fg;
	}