59 lines
1.4 KiB
Text
59 lines
1.4 KiB
Text
require "bignum"
|
|
|
|
local function euler(n, r, s, t)
|
|
-- Decimal precision.
|
|
local e10 = math.floor(n / 0.6)
|
|
|
|
-- Binary precision.
|
|
local e2 = math.round((1 + n / 0.6) / 0.30103)
|
|
|
|
-- Start with mpfr for the logs.
|
|
mpfr.init(e2)
|
|
local b = mpfr.log(new bigrat(16, 15, true, false), 2)
|
|
mpfr.mul(b, r)
|
|
local a = mpfr.setVar(1, b)
|
|
mpfr.log(new bigrat(25, 24, true, false), b)
|
|
mpfr.mul(b, s)
|
|
mpfr.addVar(a, b)
|
|
mpfr.log(new bigrat(81, 80, true, false), b)
|
|
mpfr.mul(b, t)
|
|
local u = mpfr.setVar(3, b)
|
|
mpfr.addVar(a, u)
|
|
mpfr.neg(a)
|
|
|
|
-- Switch to mpf for the basic arithmetic.
|
|
mpf.init(e2)
|
|
a = mpf.liftVar(1, a)
|
|
b = mpf.set(b, 1)
|
|
mpf.setVar(u, a)
|
|
local v = mpf.setVar(4, b)
|
|
local k = 0
|
|
local n2 = mpf.set(5, n * n)
|
|
local k2 = mpf.set(6, 0)
|
|
repeat
|
|
mpf.add(k2, (k << 1) + 1)
|
|
k += 1
|
|
mpf.mulVar(b, n2)
|
|
mpf.divVar(b, k2)
|
|
mpf.mulVar(a, n2)
|
|
mpf.div(a, k)
|
|
mpf.addVar(a, b)
|
|
mpf.div(a, k)
|
|
mpf.addVar(u, a)
|
|
mpf.addVar(v, b)
|
|
local e = mpf.frx(a, true)
|
|
until math.abs(e) >= e2
|
|
mpf.divVar(u, v)
|
|
local st = mpf.getStr(u, 10, 100)
|
|
print($"gamma {st} (maxerr. 1e-{e10})")
|
|
print($"k = {k}\n")
|
|
for i = 1, 6 do
|
|
mpfr.clear(i)
|
|
mpf.clear(i)
|
|
end
|
|
end
|
|
|
|
euler(60, 41, 30, 18)
|
|
euler(4800, 85, 62, 37)
|
|
euler(9375, 91, 68, 40)
|
|
euler(18750, 98, 73, 43)
|