-->
without javascript_semantics -- (no mpfr_exp(), mpfr_gamma() in pwa/p2js)
constant dp = 30
string fmt = "%5s: %33s, %33s, %32s\n"
requires("1.0.2") -- (mpfr_get_fixed(maxlen), mpfr_gamma)
include mpfr.e
mpfr_set_default_precision(-87) -- 87 decimal places.
sequence mc = mpfr_inits(40)
function mpfr_spouge_gamma(mpfr z)
mpfr accm = mc[1]
if mpfr_cmp_si(accm,0)=0 then
-- mc[1] := sqrt(2*PI)
mpfr_const_pi(accm)
mpfr_mul_si(accm,accm,2)
mpfr_sqrt(accm,accm)
-- k1_factrl = (k - 1)!*(-1)^k with 0!==1
mpfr k1_factrl = mpfr_init(1),
tmk = mpfr_init(),
p = mpfr_init()
for k=2 to length(mc) do
-- mc[k] = exp(13-k)*power(13-k,k-1.5)/k1_factrl
mpfr_set_si(tmk,length(mc)+1-k)
mpfr_exp(mc[k],tmk)
mpfr_set_d(p,k-1.5)
mpfr_pow(p,tmk,p)
mpfr_div(p,p,k1_factrl)
mpfr_mul(mc[k],mc[k],p)
-- k1_factrl *= -(k-1)
mpfr_mul_si(k1_factrl,k1_factrl,-(k-1))
end for
end if
accm = mpfr_init_set(accm)
for k=2 to length(mc) do
-- accm += mc[k]/(z+k-1)
mpfr ck = mpfr_init_set(mc[k]),
zk = mpfr_init_set(z)
mpfr_add_si(zk,zk,k-1)
mpfr_div(ck,ck,zk)
mpfr_add(accm,accm,ck)
end for
-- atom zc = z+length(mc)
-- accm *= exp(-zc)*power(zc,z+0.5) -- Gamma(z+1)
mpfr p = mpfr_init_set(z),
ez = mpfr_init(),
zh = mpfr_init(0.5)
mpfr_add_si(p,p,length(mc))
mpfr_neg(ez,p)
mpfr_exp(ez,ez)
mpfr_add(zh,zh,z)
mpfr_pow(p,p,zh)
mpfr_mul(accm,accm,ez)
mpfr_mul(accm,accm,p)
-- return accm/z
mpfr_div(accm,accm,z)
return accm
end function
constant mPI = mpfr_init(),
mqPI = mpfr_init()
mpfr_const_pi(mPI)
string pistr = mpfr_get_fixed(mPI,dp)
mpfr_sqrt(mqPI,mPI)
function makempfr(object x)
mpfr res
if string(x) then
x = split(x,'/')
res = mpfr_init(x[1])
mpfr_div_si(res,res,to_integer(x[2]))
elsif sequence(x) then
{mpfr x1, object x2} = x
res = mpfr_init_set(x1)
if string(x2) then
mpfr d = mpfr_init(x2)
mpfr_div(res,res,d)
else
mpfr_div_d(res,res,x2)
end if
else
res = mpfr_init(x)
end if
return res
end function
procedure mq(sequence zm, integer d=dp)
mpfr {z, mul} = apply(zm,makempfr)
mpfr s = mpfr_spouge_gamma(z)
string t = mpfr_get_fixed(s,d,10,maxlen:=dp+2)
mpfr e = mpfr_init()
mpfr_gamma(e,z)
mpfr p = mpfr_init_set(s)
mpfr_mul(p,p,mul)
mpfr_mul(p,p,p)
string zs = mpfr_get_fixed(z,3),
es = mpfr_get_fixed(e,d,10,maxlen:=dp+2),
ps = mpfr_get_fixed(p,dp,10)
printf(1,fmt,{zs,t,es,ps})
end procedure
printf(1," z %s %s %s\n",{pad(" spouge ",dp+2,"BOTH",'-'),pad(" expected ",dp+2,"BOTH",'-'),pistr})
papply({{-1.5,0.75},{-0.5,-0.5},{0.5,1},{1,{mqPI,1}},{1.5,2},{2,{mqPI,1}},{2.5,"4/3"},{3,{mqPI,2}},{3.5,"8/15"},{4,{mqPI,6}}},mq)
mq({"1/1000",{mqPI,"999.4237724845954661149822012996"}},28)
mq({"1/100",{mqPI,"99.43258511915060371353298887051"}},29)
mq({"1/10",{mqPI,"9.513507698668731836292487177265"}},30)
mq({10,{mqPI,362880}},25)
mq({100,{mqPI,"9.332621544394415268169923885627e155"}},0)
mq({150,{mqPI,"3.80892263763056972698595524350735e260"}},0)