ガンマ関数

/****************************/
/* ガンマ関数の計算         */
/*      coded by Y.Suganuma */
/****************************/
#include <stdio.h>
#include <math.h>

double gamma(double, int *);

int main()
{
	double x, y;
	int ind;

	printf("data? ");
	scanf("%lf", &x);

	y = gamma(x, &ind);

	printf("   x %f gamma(x) %f ind %d\n", x, y, ind);

	return 0;
}

/****************************************/
/* Γ（ｘ）の計算（ガンマ関数，近似式） */
/*      ier : =0 : normal               */
/*            =-1 : x=-n (n=0,1,2,･･･)  */
/*      return : 結果                   */
/****************************************/
#include <math.h>

double gamma(double x, int *ier)
{
	double err, g, s, t, v, w, y;
	long k;

	*ier = 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 * exp(-x) * pow(x,x-0.5) * s;
	}

	else {

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

		if (x < 1.5) {

			if (x < err) {
				k = (long)x;
				y = (double)k - x;
				if (fabs(y) < err || fabs(1.0-y) < err)
					*ier = -1;
			}

			if (*ier == 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;
}