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})