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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
; 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!
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX