(phixonline)-->
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