ガンマ関数の近似式

初等関数の近似
を参考にしましたが、近似式でうまく出力できなかったので『C言語による最新アルゴリズム事典』を参考にしてガンマ関数について書きました。以下にコードとその近似式について書きます。

ガンマ関数
\Gamma (x) = \int_{0}^{\infty}e^{-t}t^{x-1} dt

漸化式
\Gamma (x+1)= x\Gamma(x)

対数log\Gamma(x)は漸近展開すると
log \Gamma(x) = (x - \frac{1}{2}) - x + \frac{1}{2}log 2 \pi +\sum \frac{B_{2n}}{2n(2n-1)x^{2n-1}}

また、x<0のときは
\Gamma(x) \Gamma(1-x) = \frac{\pi}{sin \pi x}
とした。

と参考文献には書かれているんだけどコード的にところどころ違う。

#-*- coding:utf-8 -*-
from numpy import *
PI=3.14159265358979324
LOG_2PI=1.837877066400934548
N=8
B0=1
B1=(-1.0/2.0)
B2=(1.0/6.0)
B4=(1.0/30.0)
B6=(1.0/42.0)
B8=(-1.0/30.0)
B10=(5.0/66.0)
B12=(-691.0/2730.0)
B14=(7.0/6.0)
B16=(-3617.0/510.0)
B18=(43867.0/798.0)
B20=(-174611.0/330.0)

def log_gamma(x):
	v=1
	while x < N:
		v*=x
		x+=1
	w=1/(x*x)
	return 0.5*LOG_2PI - log(v) -x + (x-0.5) *log(x) + (((((((( (B16/(16*15))*w + (B14/(14*13))) *w + (B12/(12*11))) *w + (B10*9)))*w +(B8/(8*7)))*w + (B6/(6*5)))*w +(B4/(4*3))) *w + (B2/(2*1)))/x

def gamma(x):
	if x < 0:
		return PI/(sin(PI*x)*exp(log_gamma(1-x)))
	return exp(log_gamma(x))