local fmt = require "fmt"

local eps <const> = 1e-6

print("From the definition, err. 3e-10")
local n = 400
local h = 1
for k = 2, n do h += 1 / k end
-- Faster convergence: Negoi, 1997.
local a = math.log(n + 0.5 + 1 / (24 * n))

fmt.print("Hn    %0.14f", h)
fmt.print("gamma %0.14f\nk = %d\n", h - a, n)

print("Sweeney, 1963, err. idem")
n = 21
local s = {0, n}
local r = n
local k = 1
repeat
    k += 1
    r *= n / k
    s[(k & 1) + 1] += r / k
until r <= eps

fmt.print("gamma %0.14f\nk = %d\n", s[2] - s[1] - math.log(n), k)

print("Bailey, 1988")
n = 5
a = 1
h = 1
local n2 = 1 << n
r = 1
k = 1
repeat
    k += 1
    r *= n2 / k
    h += 1 / k
    local b = a
    a += r * h
until math.abs(b - a) <= eps
a *= n2 / math.exp(n2)
fmt.print("gamma %0.14f\nk = %d\n", a - n * math.log(2), k)

print("Brent-McMillan, 1980")
n = 13
a = -math.log(n)
local b = 1
local u = a
local v = b
n2 = n * n
local k2 = 0
k = 0
repeat
    k2 += 2 * k + 1
    k += 1
    a *= n2 / k
    b *= n2 / k2
    a = (a + b) / k
    u += a
    v += b
until math.abs(a) <= eps

fmt.print("gamma %0.14f\nk = %d\n", u / v, k)

print("How Euler did it in 1735")
-- Bernoulli numbers with even indices.
local b2 = {1, 1/6, -1/30, 1/42, -1/30, 5/66, -691/2730, 7/6, -3617/510, 43867/798}
local m = 7
n = 10
-- n'th harmonic number.
h = 1
for l = 2, n do h += 1 / l end
fmt.print("Hn    %0.14f", h)
h -= math.log(n)
fmt.print("  -ln %0.14f", h)
-- Expansion C = -digamma(1).
a = -1 / (2 * n)
n2 = n * n
r = 1
for l = 1, m do
    r *= n2
    a += b2[l + 1] / (2 * l * r)
end
fmt.print("err  %0.14f\ngamma %0.14f\nk = %d", a, h + a, n + m)
print("\nC = 0.57721566490153286...")