96 lines
2.3 KiB
Plaintext
96 lines
2.3 KiB
Plaintext
include mpfr.e
|
|
mpfr_set_default_prec(-87) -- 87 decimal places.
|
|
|
|
sequence c = mpfr_inits(40)
|
|
|
|
function gamma(atom z)
|
|
mpfr accm = c[1]
|
|
if mpfr_cmp_si(accm,0)=0 then
|
|
-- c[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(c) do
|
|
-- c[k] = exp(13-k)*power(13-k,k-1.5)/k1_factrl
|
|
mpfr_set_si(tmk,length(c)+1-k)
|
|
mpfr_exp(c[k],tmk)
|
|
mpfr_set_d(p,k-1.5)
|
|
mpfr_pow(p,tmk,p)
|
|
mpfr_div(p,p,k1_factrl)
|
|
mpfr_mul(c[k],c[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(c) do
|
|
-- accm += c[k]/(z+k-1)
|
|
mpfr ck = mpfr_init_set(c[k]),
|
|
zk = mpfr_init(z+k-1)
|
|
mpfr_div(ck,ck,zk)
|
|
mpfr_add(accm,accm,ck)
|
|
end for
|
|
atom zc = z+length(c)
|
|
-- accm *= exp(-zc)*power(zc,z+0.5) -- Gamma(z+1)
|
|
mpfr ez = mpfr_init(-zc),
|
|
p = mpfr_init(zc),
|
|
zh = mpfr_init(z+0.5)
|
|
mpfr_exp(ez,ez)
|
|
mpfr_pow(p,p,zh)
|
|
mpfr_mul(accm,accm,ez)
|
|
mpfr_mul(accm,accm,p)
|
|
-- return accm/z
|
|
mpfr_set_d(ez,z)
|
|
mpfr_div(accm,accm,ez)
|
|
return accm
|
|
end function
|
|
|
|
function gamma2(atom z)
|
|
mpfr r = mpfr_init(z)
|
|
mpfr_gamma(r,r)
|
|
return r
|
|
end function
|
|
|
|
constant FMT = "%43.40Rf"
|
|
|
|
procedure sq(mpfr x, integer n, d=1)
|
|
mpfr p = mpfr_init_set(x)
|
|
mpfr_mul_si(p,p,n)
|
|
mpfr_div_si(p,p,d)
|
|
mpfr_mul(p,p,p)
|
|
string xs = mpfr_sprintf(FMT,x),
|
|
ps = mpfr_sprintf(FMT,p)
|
|
printf(1,"%s,%s\n",{xs,ps})
|
|
end procedure
|
|
|
|
procedure si(mpfr x)
|
|
string xs = mpfr_sprintf(FMT,x)
|
|
printf(1,"%s\n",trim_tail(xs,".0"))
|
|
end procedure
|
|
|
|
sq(gamma(-3/2),3,4)
|
|
sq(gamma(-1/2),-1,2)
|
|
sq(gamma(1/2),1)
|
|
si(gamma(1))
|
|
sq(gamma(3/2),2)
|
|
si(gamma(2))
|
|
sq(gamma(5/2),4,3)
|
|
si(gamma(3))
|
|
sq(gamma(7/2),8,15)
|
|
si(gamma(4))
|
|
puts(1,"mpfr_gamma():\n")
|
|
sq(gamma2(-3/2),3,4)
|
|
sq(gamma2(-1/2),-1,2)
|
|
sq(gamma2(1/2),1)
|
|
si(gamma2(1))
|
|
sq(gamma2(3/2),2)
|
|
si(gamma2(2))
|
|
sq(gamma2(5/2),4,3)
|
|
si(gamma2(3))
|
|
sq(gamma2(7/2),8,15)
|
|
si(gamma2(4))
|