HSPで実数の階乗を求める(スターリングの近似とオイラーの反射式) | Codetter（こーどったー）

; gamma.hsp
; bmpとしてデータ列を開く
; 2007.12.20, 石橋祥太

#module "Gamma"
;kzfunc(p1)
;p1=0～(0)
;p1の階乗を整数で返します。
;gfunc(p1)
;p1のΓ(ガンマ)関数を実数で返します。
;*スターリングの近似とオイラーの反射式を使用して
;いますが、100を超えたあたりから精度が著しく低下
;します。
#define M_E	2.7182818284590452354
#define PI	3.14159265358979323846
#defcfunc kzfunc int u1;階乗
	an=1.0
	if u1<0 : return -logf(0);0以下
	repeat u1,1
	an*=cnt
	loop
	return an
#defcfunc gfuncx double x1
	x2=x1-1
	if x2<0|x2>Pi/2 : return -logf(0)
	return sqrt(Pi*2*x2)*expf(logf(x2/M_E)*(x2/M_E))
#defcfunc gfunc double u1;ガンマ関数
	if u1<1 : u2=u1 : else : u2=-u1+1
	repeat
		if u2>Pi/2 : bnt=cnt : break
		u2+
	loop
	an1=gfuncx(u2)
	repeat bnt
		u2-
		if u2=0 : break
		an1/=u2
	loop
	if u2=0 : return 1.0*kzfunc(-1+u1)
	if u1>=1 : if an1*sin(Pi*u1)=0 : return logf(0) : else : an1=Pi/(an1*sin(Pi*u1))
	if u1<1 : an2=an1/M_E*Pi/2 : else : an2=an1*M_E/Pi*2
	return an2
#undef gfuncx
#global

; TEST
mes gfunc(0+1)	; 0!
mes gfunc(1.5)	; 0.5!
mes gfunc(1+1)	; 1!
mes gfunc(2.5)	; 1.5!
mes gfunc(2+1)	; 2!
mes gfunc(3.5)	; 2.5!
mes gfunc(3+1)	; 3!
mes gfunc(4.5)	; 3.5!
mes gfunc(4+1)	; 4!
mes gfunc(5.5)	; 4.5!
mes gfunc(5+1)	; 5!
mes gfunc(6.5)	; 5.5!