146 lines
2.5 KiB
Plaintext
146 lines
2.5 KiB
Plaintext
/*********************************************
|
|
Subject: Comparing five methods for
|
|
computing Euler's constant 0.5772...
|
|
// https://rosettacode.org/wiki/Euler%27s_constant_0.5772...
|
|
--------------------------------------------*/
|
|
double a, b, h, n2, r, u, v;
|
|
float floatA, floatB, floatN2;
|
|
int k, k2, m, n;
|
|
double eps = 1e-6;
|
|
|
|
void setup() {
|
|
size(100, 100);
|
|
noLoop();
|
|
}
|
|
|
|
void draw() {
|
|
println("From the definition, err. 3e-10\n");
|
|
|
|
n = 400;
|
|
|
|
h = 1;
|
|
|
|
for (int k = 2; k <= n; k++) {
|
|
h += 1.0 / k;
|
|
}
|
|
//faster convergence: Negoi, 1997
|
|
a = log(n +.5 + 1.0 / (24*n));
|
|
|
|
println("Hn ", h);
|
|
println("gamma ", h - a);
|
|
println("k = ", n);
|
|
println("");
|
|
|
|
|
|
println("Sweeney, 1963, err. idem");
|
|
n = 21;
|
|
|
|
double s[] = {0, n};
|
|
r = n;
|
|
k = 1;
|
|
while (r > eps) {
|
|
k ++;
|
|
r *= (double) n / k;
|
|
s[k & 1] = s[k & 1] + r / k;
|
|
}
|
|
|
|
// println("gamma %.16f\nk = %d\n\n", s[1] - s[0] - log(n), k);
|
|
|
|
println("Hn ", h);
|
|
println("gamma ", s[1] - s[0] - log(n));
|
|
println("k = ", k);
|
|
println("");
|
|
|
|
println("Bailey, 1988");
|
|
n = 5;
|
|
floatA = 1;
|
|
h = 1;
|
|
floatN2 = pow(2, n);
|
|
r = 1;
|
|
k = 1;
|
|
while (abs(floatB - floatA) > eps) {
|
|
k += 1;
|
|
r *= floatN2 / k;
|
|
h += 1.0 / k;
|
|
floatB = floatA;
|
|
floatA += r * h;
|
|
}
|
|
floatA *= floatN2 / exp(floatN2);
|
|
|
|
println("gamma ", floatA - n * log(2));
|
|
println("k = ", k);
|
|
println("");
|
|
|
|
println("Brent-McMillan, 1980");
|
|
|
|
n = 13;
|
|
|
|
floatA = -log(n);
|
|
floatB = 1;
|
|
u = a;
|
|
v = b;
|
|
n2 = n * n;
|
|
k2 = 0;
|
|
k = 0;
|
|
|
|
|
|
while (abs(floatA) > eps) {
|
|
k2 += 2*k + 1;
|
|
k += 1;
|
|
floatA *= n2 / k;
|
|
floatB *= n2 / k2;
|
|
floatA = (floatA + floatB) / k;
|
|
u += floatA;
|
|
v += floatB;
|
|
}
|
|
println("gamma ", u / v);
|
|
println("k ", k);
|
|
|
|
|
|
|
|
println("How Euler did it in 1735\n");
|
|
//Bernoulli numbers with even indices
|
|
|
|
double[] B2 = new double[11];
|
|
|
|
B2[1] = 1.0;
|
|
B2[2] = 1.0/6;
|
|
B2[3] = -1.0/30;
|
|
B2[4] = 1.0/42;
|
|
B2[5] = -1.0/30;
|
|
B2[6] = 5.0/66;
|
|
B2[7] = -691.0/2730;
|
|
B2[8] = 7.0/6;
|
|
B2[9] = -3617.0/510;
|
|
B2[10]= 43867.0/798;
|
|
|
|
m = 7;
|
|
n = 10;
|
|
|
|
//n-th harmonic number
|
|
h = 1;
|
|
for (k = 2; k <= n; k++) {
|
|
h += 1.0 / k;
|
|
}
|
|
println("Hn ", h);
|
|
|
|
h -= log(n);
|
|
println(" -ln ", h);
|
|
|
|
//expansion C = -digamma(1)
|
|
a = -1.0 / (2*n);
|
|
n2 = n * n;
|
|
r = 1;
|
|
for (k = 1; k <= m; k++) {
|
|
r *= n2;
|
|
a += B2[k] / (2*k * r);
|
|
}
|
|
|
|
println("");
|
|
println("err ", a);
|
|
println("gamma ", h + a );
|
|
println("k = ", n + m);
|
|
println("");
|
|
println("C = 0.57721566490153286...\n");
|
|
}
|