ガンマ関数

<?php

/****************************/
/* ガンマ関数の計算         */
/*      coded by Y.Suganuma */
/****************************/

	printf("data? ");
	fscanf(STDIN, "%lf", $x);

	$y = gamma($x, $ind);

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

/****************************************/
/* Γ（ｘ）の計算（ガンマ関数，近似式） */
/*      ier : =0 : normal               */
/*            =-1 : x=-n (n=0,1,2,･･･)  */
/*      return : 結果                   */
/****************************************/
function gamma($x, &$ier)
{
	$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 = intval($x);
				$y = floatval($k) - $x;
				if (abs($y) < $err || abs(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;
}

?>