87 lines
2.1 KiB
Text
87 lines
2.1 KiB
Text
without js -- no mpfr_get_d_2exp() in mpfr.js as yet
|
|
requires("1.0.2") -- mpfr_get_d_2exp(), mpfr_addmul_si()
|
|
include mpfr.e
|
|
mpfr u, v, k2;
|
|
integer e, e10, e2
|
|
atom f
|
|
|
|
//log(x/y) with the Taylor series for atanh(x-y/x+y)
|
|
procedure ln(mpfr s, integer x, y)
|
|
mpfr d = u, q = v;
|
|
assert((x-y)==1)
|
|
mpfr_set_si(s, x+y)
|
|
mpfr_si_div(s, 1, s) // s = 1 / (x + y)
|
|
mpfr_mul(k2, s, s) // k2 = s * s
|
|
mpfr_set(d, s)
|
|
integer k = 1
|
|
while true do
|
|
k += 2;
|
|
mpfr_mul(d, d, k2) // d *= k2
|
|
mpfr_div_si(q, d, k) // q = d / k
|
|
mpfr_add(s, s, q) // s += q
|
|
{f,e} = mpfr_get_d_2exp(q)
|
|
if abs(e)>=e2 then exit end if
|
|
end while
|
|
mpfr_mul_si(s, s, 2) //s *= 2
|
|
end procedure
|
|
|
|
mpfr a, b
|
|
integer k,
|
|
n = 60, -- (required precision in decimal dp *6/10)
|
|
n2,
|
|
r = 41,
|
|
s = 30,
|
|
t = 18;
|
|
|
|
// n = 2^i * 3^j * 5^k
|
|
|
|
// log(n) = r * log(16/15) + s * log(25/24) + t * log(81/80)
|
|
|
|
// solve linear system for r, s, t
|
|
// 4 -3 -4| i
|
|
// -1 -1 4| j
|
|
// -1 2 -1| k
|
|
|
|
//decimal precision
|
|
e10 = floor(n/0.6)
|
|
//binary precision
|
|
e2 = floor((1 + e10) / 0.30103)
|
|
|
|
mpfr_set_default_precision(e2)
|
|
{a, b, u, v, k2} = mpfr_inits(5)
|
|
|
|
//Compute log terms
|
|
ln(b, 16, 15) mpfr_mul_si(a, b, r) // a = r * b
|
|
ln(b, 25, 24) mpfr_addmul_si(a, b, s) // a += s * b
|
|
ln(b, 81, 80) mpfr_addmul_si(a, b, t) // a += t * b
|
|
|
|
mpfr_neg(a, a) // a = -a
|
|
mpfr_set_si(b, 1) // b = 1
|
|
mpfr_set (u, a)
|
|
mpfr_set (v, b)
|
|
|
|
k = 0;
|
|
n2 = n * n;
|
|
mpfr_set_si(k2, 0) // k2 = k * k (as below)
|
|
while true do
|
|
mpfr_add_si(k2, k2, k*2+1) // k2 += 2k + 1
|
|
k += 1;
|
|
|
|
mpfr_div(b, b, k2)
|
|
mpfr_mul_si(b, b, n2) // b = b * n2 / k2
|
|
|
|
mpfr_div_si(a, a, k)
|
|
mpfr_mul_si(a, a, n2)
|
|
mpfr_add (a, a, b)
|
|
mpfr_div_si(a, a, k) // a = (a * n2 / k + b) / k
|
|
|
|
mpfr_add(u, u, a) // u += a
|
|
mpfr_add(v, v, b) // v += b
|
|
|
|
{f,e} = mpfr_get_d_2exp (a)
|
|
if abs(e)>=e2 then exit end if
|
|
end while
|
|
|
|
mpfr_div(u, u, v)
|
|
string su = mpfr_get_fixed(u,e10)
|
|
printf(1,"gamma %s (maxerr. 1e-%d)\n", {su, e10})
|