114 lines
3.9 KiB
Text
114 lines
3.9 KiB
Text
with javascript_semantics
|
|
sequence c = repeat(0,12)
|
|
function spouge_gamma(atom z)
|
|
atom accm = c[1]
|
|
if accm=0 then
|
|
accm = sqrt(2*PI)
|
|
c[1] = accm
|
|
atom k1_factrl = 1 -- (k - 1)!*(-1)^k with 0!==1
|
|
for k=2 to 12 do
|
|
c[k] = exp(13-k)*power(13-k,k-1.5)/k1_factrl
|
|
k1_factrl *= -(k-1)
|
|
end for
|
|
end if
|
|
for k=2 to 12 do
|
|
accm += c[k]/(z+k-1)
|
|
end for
|
|
accm *= exp(-(z+12))*power(z+12,z+0.5) -- Gamma(z+1)
|
|
return accm/z
|
|
end function
|
|
|
|
function taylor_gamma(atom x)
|
|
-- (good for values between 0 and 1, apparently)
|
|
constant t = { 1.00000_00000_00000_00000,
|
|
0.57721_56649_01532_86061,
|
|
-0.65587_80715_20253_88108,
|
|
-0.04200_26350_34095_23553,
|
|
0.16653_86113_82291_48950,
|
|
-0.04219_77345_55544_33675,
|
|
-0.00962_19715_27876_97356,
|
|
0.00721_89432_46663_09954,
|
|
-0.00116_51675_91859_06511,
|
|
-0.00021_52416_74114_95097,
|
|
0.00012_80502_82388_11619,
|
|
-0.00002_01348_54780_78824,
|
|
-0.00000_12504_93482_14267,
|
|
0.00000_11330_27231_98170,
|
|
-0.00000_02056_33841_69776,
|
|
0.00000_00061_16095_10448,
|
|
0.00000_00050_02007_64447,
|
|
-0.00000_00011_81274_57049,
|
|
0.00000_00001_04342_67117,
|
|
0.00000_00000_07782_26344,
|
|
-0.00000_00000_03696_80562,
|
|
0.00000_00000_00510_03703,
|
|
-0.00000_00000_00020_58326,
|
|
-0.00000_00000_00005_34812,
|
|
0.00000_00000_00001_22678,
|
|
-0.00000_00000_00000_11813,
|
|
0.00000_00000_00000_00119,
|
|
0.00000_00000_00000_00141,
|
|
-0.00000_00000_00000_00023,
|
|
0.00000_00000_00000_00002 }
|
|
atom y = x-1,
|
|
s = t[$]
|
|
for n=length(t)-1 to 1 by -1 do
|
|
s = s*y + t[n]
|
|
end for
|
|
return 1/s
|
|
end function
|
|
|
|
function lanczos_gamma(atom z)
|
|
if z<0.5 then
|
|
return PI / (sin(PI*z)*lanczos_gamma(1-z))
|
|
end if
|
|
-- use a lanczos approximation:
|
|
atom x = 0.99999999999980993,
|
|
t = z + 6.5;
|
|
sequence p = { 676.5203681218851,
|
|
-1259.1392167224028,
|
|
771.32342877765313,
|
|
-176.61502916214059,
|
|
12.507343278686905,
|
|
-0.13857109526572012,
|
|
9.9843695780195716e-6,
|
|
1.5056327351493116e-7 }
|
|
z -= 1
|
|
for i=1 to length(p) do
|
|
x += p[i] / (z + i)
|
|
end for
|
|
return sqrt(2*PI) * power(t,z+0.5) * exp(-t) * x
|
|
end function
|
|
|
|
constant sqPI = sqrt(PI)
|
|
|
|
procedure sq(sequence zm, string fmt="%19.16f")
|
|
atom {z, mul} = zm
|
|
atom e = sqPI/mul
|
|
sequence s = {spouge_gamma(z),
|
|
taylor_gamma(z),
|
|
lanczos_gamma(z)},
|
|
error = sq_abs(sq_sub(s,e))
|
|
string t = join(s,", ",fmt:=fmt)&", "
|
|
integer bdx = smallest(error,return_index:=true)
|
|
atom best = s[bdx],
|
|
p = s[bdx]*mul
|
|
for i=1 to length(s) do -- (potentially mark >1)
|
|
if s[i]=best then
|
|
t[i*22-2..i*22-1] = "*,"
|
|
end if
|
|
end for
|
|
string es = sprintf(fmt,e)
|
|
printf(1,"%5g: %s %s, %19.16f\n",{z,t,es,p*p})
|
|
end procedure
|
|
|
|
printf(1," z ------ spouge ----- ----- taylor ------ ----- lanczos ----- ---- expected ----- %19.16f\n",PI)
|
|
papply({{-3/2,3/4},{-1/2,-1/2},{1/2,1},{1,sqPI},{3/2,2},{2,sqPI},{5/2,4/3},{3,sqPI/2},{7/2,8/15},{4,sqPI/6}},sq)
|
|
sq({0.001,sqPI/999.4237725},"%19.15f")
|
|
sq({0.01,sqPI/99.43258512},"%19.16f")
|
|
sq({0.1,sqPI/9.513507699},"%19.16f")
|
|
sq({10,sqPI/362880},"%19.12f")
|
|
sq({100,sqPI/9.332621544e155},"%19.13g")
|
|
if machine_bits()=64 then
|
|
sq({150,sqPI/3.808922638e260},"%19.13g") -- (fatal power overflow error on 32 bits)
|
|
end if
|