ガンマ関数

# -*- coding: UTF-8 -*-
from math import *
import numpy as np

############################################
# Γ（ｘ）の計算（ガンマ関数，近似式）
#      ier : =0 : normal
#            =-1 : x=-n (n=0,1,2,･･･)
#      return : 結果
#      coded by Y.Suganuma
############################################

def gamma(x, ier) :

	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 * 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 = int(x)
				y = float(k) - x
				if abs(y) < err or 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

----------------------------------

# -*- coding: UTF-8 -*-
import numpy as np
import sys
from math import *
from function import gamma

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

ind = np.empty(1, np.int)
s   = input("data? ")
x   = float(s)

y   = gamma(x, ind)

print("   x " + str(x) + " gamma(x) " + str(y) + " ind " + str(ind[0]))