86 lines
1.6 KiB
Text
86 lines
1.6 KiB
Text
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...")
|