121 lines
1.6 KiB
Plaintext
121 lines
1.6 KiB
Plaintext
'**********************************************
|
|
'Subject: Comparing five methods for
|
|
' computing Euler's constant 0.5772...
|
|
'tested : FreeBasic 1.08.1
|
|
'----------------------------------------------
|
|
const eps = 1e-6
|
|
dim as double a, b, h, n2, r, u, v
|
|
dim as integer k, k2, m, n
|
|
|
|
? "From the definition, err. 3e-10"
|
|
|
|
n = 400
|
|
|
|
h = 1
|
|
for k = 2 to n
|
|
h += 1 / k
|
|
next k
|
|
'faster convergence: Negoi, 1997
|
|
a = log(n +.5 + 1 / (24*n))
|
|
|
|
? "Hn "; h
|
|
? "gamma"; h - a; !"\nk ="; n
|
|
?
|
|
|
|
|
|
? "Sweeney, 1963, err. idem"
|
|
|
|
n = 21
|
|
|
|
dim as double s(1) = {0, n}
|
|
r = n
|
|
k = 1
|
|
do
|
|
k += 1
|
|
r *= n / k
|
|
s(k and 1) += r / k
|
|
loop until r < eps
|
|
|
|
? "gamma"; s(1) - s(0) - log(n); !"\nk ="; k
|
|
?
|
|
|
|
|
|
? "Bailey, 1988"
|
|
|
|
n = 5
|
|
|
|
a = 1
|
|
h = 1
|
|
n2 = 2^n
|
|
r = 1
|
|
k = 1
|
|
do
|
|
k += 1
|
|
r *= n2 / k
|
|
h += 1 / k
|
|
b = a: a += r * h
|
|
loop until abs(b - a) < eps
|
|
a *= n2 / exp(n2)
|
|
|
|
? "gamma"; a - n * log(2); !"\nk ="; k
|
|
?
|
|
|
|
|
|
? "Brent-McMillan, 1980"
|
|
|
|
n = 13
|
|
|
|
a = -log(n)
|
|
b = 1
|
|
u = a
|
|
v = b
|
|
n2 = n * n
|
|
k2 = 0
|
|
k = 0
|
|
do
|
|
k2 += 2*k + 1
|
|
k += 1
|
|
a *= n2 / k
|
|
b *= n2 / k2
|
|
a = (a + b) / k
|
|
u += a
|
|
v += b
|
|
loop until abs(a) < eps
|
|
|
|
? "gamma"; u / v; !"\nk ="; k
|
|
?
|
|
|
|
|
|
? "How Euler did it in 1735"
|
|
'Bernoulli numbers with even indices
|
|
dim as double B2(9) = {1,1/6,-1/30,1/42,_
|
|
-1/30,5/66,-691/2730,7/6,-3617/510,43867/798}
|
|
m = 7
|
|
if m > 9 then end
|
|
|
|
n = 10
|
|
|
|
'n-th harmonic number
|
|
h = 1
|
|
for k = 2 to n
|
|
h += 1 / k
|
|
next k
|
|
? "Hn "; h
|
|
|
|
h -= log(n)
|
|
? " -ln"; h
|
|
|
|
'expansion C = -digamma(1)
|
|
a = -1 / (2*n)
|
|
n2 = n * n
|
|
r = 1
|
|
for k = 1 to m
|
|
r *= n2
|
|
a += B2(k) / (2*k * r)
|
|
next k
|
|
|
|
? "err "; a; !"\ngamma"; h + a; !"\nk ="; n + m
|
|
?
|
|
? "C = 0.57721566490153286..."
|
|
end
|