Euler's constant 0.5772...
You are encouraged to solve this task according to the task description, using any language you may know.
- Task.
Compute the Euler constant 0.5772...
Discovered by Leonhard Euler around 1730, it is the most ubiquitous mathematical constant after pi and e, but appears more arcane than these.
Denoted gamma (γ), it measures the amount by which the partial sums of the harmonic series (the simplest diverging series) differ from the logarithmic function (its approximating integral): lim n → ∞ (1 + 1/2 + 1/3 + … + 1/n − log(n)).
The definition of γ converges too slowly to be numerically useful, but in 1735 Euler himself applied his recently discovered summation formula to compute ‘the notable number’ accurate to 15 places. For a single-precision implementation this is still the most economic algorithm.
In 1961, the young Donald Knuth used Euler's method to evaluate γ to 1271 places. Knuth found that the computation of the Bernoulli numbers required in the Euler-Maclaurin formula was the most time-consuming part of the procedure.
The next year Dura Sweeney computed 3566 places, using a formula based on the expansion of the exponential integral which didn't need Bernoulli numbers. It's a bit-hungry method though: 2d digits of working precision obtain d correct places only.
This was remedied in 1988 by David Bailey; meanwhile Richard Brent and Ed McMillan had published an even more efficient algorithm based on Bessel function identities and found 30100 places in 20 hours time.
Nowadays the old records have far been exceeded: over 6·1011 decimal places are already known. These massive computations suggest that γ is neither rational nor algebraic, but this is yet to be proven.
- References.
[1] Gourdon and Sebah, The Euler constant γ. (for all formulas)
[2] Euler's original journal article translated from the latin (p. 9)
C
Single precision
<lang c>/********************************************* Subject: Comparing five methods for
computing Euler's constant 0.5772...
tested : tcc-0.9.27
*/
- include <math.h>
- include <stdio.h>
- define eps 1e-6
int main(void) { double a, b, h, n2, r, u, v; int k, k2, m, n;
printf("From the definition, err. 3e-10\n");
n = 400;
h = 1; for (k = 2; k <= n; k++) {
h += 1.0 / k;
} //faster convergence: Negoi, 1997 a = log(n +.5 + 1.0 / (24*n));
printf("Hn %.16f\n", h); printf("gamma %.16f\nk = %d\n\n", h - a, n);
printf("Sweeney, 1963, err. idem\n");
n = 21;
double s[] = {0, n}; r = n; k = 1; do {
k += 1; r *= (double) n / k; s[k & 1] += r / k;
} while (r > eps);
printf("gamma %.16f\nk = %d\n\n", s[1] - s[0] - log(n), k);
printf("Bailey, 1988\n");
n = 5;
a = 1; h = 1; n2 = pow(2,n); r = 1; k = 1; do {
k += 1; r *= n2 / k; h += 1.0 / k; b = a; a += r * h;
} while (fabs(b - a) > eps); a *= n2 / exp(n2);
printf("gamma %.16f\nk = %d\n\n", a - n * log(2), k);
printf("Brent-McMillan, 1980\n");
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;
} while (fabs(a) > eps);
printf("gamma %.16f\nk = %d\n\n", u / v, k);
printf("How Euler did it in 1735\n");
//Bernoulli numbers with even indices
double B2[] = {1.0,1.0/6,-1.0/30,1.0/42,-1.0/30,\
5.0/66,-691.0/2730,7.0/6,-3617.0/510,43867.0/798};
m = 7; if (m > 9) return(0);
n = 10;
//n-th harmonic number h = 1; for (k = 2; k <= n; k++) {
h += 1.0 / k;
} printf("Hn %.16f\n", h);
h -= log(n); printf(" -ln %.16f\n", 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);
}
printf("err %.16f\ngamma %.16f\nk = %d", a, h + a, n + m);
printf("\n\nC = 0.57721566490153286...\n"); }</lang>
- output:
From the definition, err. 3e-10 Hn 6.5699296911765055 gamma 0.5772156645765731 k = 400 Sweeney, 1963, err. idem gamma 0.5772156645636311 k = 68 Bailey, 1988 gamma 0.5772156649015341 k = 89 Brent-McMillan, 1980 gamma 0.5772156649015329 k = 40 How Euler did it in 1735 Hn 2.9289682539682538 -ln 0.6263831609742079 err -0.0491674960726754 gamma 0.5772156649015325 k = 17 C = 0.57721566490153286...
Multi precision
From first principles
<lang c>/************************************************** Subject: Computation of Euler's constant 0.5772...
with the Brent-McMillan algorithm B1,
Math. Comp. 34 (1980), 305-312
tested : tcc-0.9.27 with gmp 6.2.0
*/
- include <gmp.h>
- include <stdio.h>
- include <stdlib.h>
- include <time.h>
//multi-precision float pointers mpf_ptr u, v, k2;
//precision parameters unsigned long e10, e2; long e; double f;
//log(x/y) with the Taylor series for atanh(x-y/x+y) void ln (mpf_ptr s, unsigned long x, unsigned long y) { mpf_ptr d = u, q = v; unsigned long k;
//Möbius transformation k = x; x -= y; y += k;
if (x != 1) {
printf ("ln: illegal argument x - y != 1");
exit;
}
//s = 1 / (x + y) mpf_set_ui (s, y); mpf_ui_div (s, 1, s); //k2 = s * s mpf_mul (k2, s, s); mpf_set (d, s);
k = 1;
do {
k += 2;
//d *= k2
mpf_mul (d, d, k2);
//q = d / k
mpf_div_ui (q, d, k);
//s += q
mpf_add (s, s, q);
f = mpf_get_d_2exp (&e, q); } while (abs(e) < e2);
//s *= 2 mpf_mul_2exp (s, s, 1);
}
int main (void) { mpf_ptr a = malloc(sizeof(__mpf_struct)); mpf_ptr b = malloc(sizeof(__mpf_struct)); u = malloc(sizeof(__mpf_struct)); v = malloc(sizeof(__mpf_struct)); k2 = malloc(sizeof(__mpf_struct)); //unsigned long integers unsigned long k, n, n2, r, s, t;
clock_t tim = clock();
// n = 2^i * 3^j * 5^k
// log(n) = r * log(16/15) + s * log(25/24) + t * log(81/80)
// solve linear system for r, s, t // 4 -3 -4| i // -1 -1 4| j // -1 2 -1| k
//examples t = 1; switch (t) { case 1 :
n = 60; r = 41; s = 30; t = 18;
//100 digits break; case 2 :
n = 4800; r = 85; s = 62; t = 37;
//8000 digits, 0.6 s break; case 3 :
n = 9375; r = 91; s = 68; t = 40;
//15625 digits, 2.5 s break; default :
n = 18750; r = 98; s = 73; t = 43;
//31250 digits, 12 s. @2.00GHz }
//decimal precision e10 = n / .6; //binary precision e2 = (1 + e10) / .30103;
//initialize mpf's mpf_set_default_prec (e2); mpf_inits (a, b, u, v, k2, (mpf_ptr)0);
//Compute log terms
ln (b, 16, 15);
//a = r * b mpf_mul_ui (a, b, r);
ln (b, 25, 24);
//a += s * b mpf_mul_ui (u, b, s); mpf_add (a, a, u);
ln (b, 81, 80);
//a += t * b mpf_mul_ui (u, b, t); mpf_add (a, a, u);
//gmp_printf ("log(%lu) %.*Ff\n", n, e10, a);
//B&M, algorithm B1
//a = -a, b = 1 mpf_neg (a, a); mpf_set_ui (b, 1); mpf_set (u, a); mpf_set (v, b);
k = 0; n2 = n * n; //k2 = k * k mpf_set_ui (k2, 0); do {
//k2 += 2k + 1 mpf_add_ui (k2, k2, (k << 1) + 1); k += 1;
//b = b * n2 / k2 mpf_div (b, b, k2); mpf_mul_ui (b, b, n2); //a = (a * n2 / k + b) / k mpf_div_ui (a, a, k); mpf_mul_ui (a, a, n2); mpf_add (a, a, b); mpf_div_ui (a, a, k);
//u += a, v += b mpf_add (u, u, a); mpf_add (v, v, b);
f = mpf_get_d_2exp (&e, a);
} while (abs(e) < e2);
mpf_div (u, u, v); gmp_printf ("gamma %.*Ff (maxerr. 1e-%lu)\n", e10, u, e10);
gmp_printf ("k = %lu\n\n", k);
tim = clock() - tim; printf("time: %.7f s\n",((double)tim)/CLOCKS_PER_SEC); }</lang>
- output:
gamma 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495 (maxerr. 1e-100) k = 255
The easy way
<lang c>/******************************************* Subject: Euler's constant 0.5772... tested : tcc-0.9.27 with mpfr 4.1.0
*/
- include <gmp.h>
- include <mpfr.h>
- include <stdio.h>
- include <stdlib.h>
- include <time.h>
int main (void) { mpfr_ptr a = malloc(sizeof(__mpfr_struct)); unsigned long e2, e10; clock_t tim = clock();
//decimal precision e10 = 100;
//binary precision e2 = (1 + e10) / .30103; mpfr_init2 (a, e2);
mpfr_const_euler (a, MPFR_RNDN); mpfr_printf ("gamma %.*Rf\n\n", e10, a);
tim = clock() - tim; gmp_printf ("time: %.7f s\n",((double)tim)/CLOCKS_PER_SEC); }</lang>
FreeBASIC
Single precision
<lang freebasic>'********************************************** '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</lang>
- output:
From the definition, err. 3e-10 Hn 6.569929691176506 gamma 0.5772156645765731 k = 400 Sweeney, 1963, err. idem gamma 0.5772156645636311 k = 68 Bailey, 1988 gamma 0.5772156649015341 k = 89 Brent-McMillan, 1980 gamma 0.5772156649015329 k = 40 How Euler did it in 1735 Hn 2.928968253968254 -ln 0.6263831609742079 err -0.04916749607267539 gamma 0.5772156649015325 k = 17 C = 0.57721566490153286...
Multi precision
From first principles <lang freebasic>'*************************************************** 'Subject: Computation of Euler's constant 0.5772... ' with the Brent-McMillan algorithm B1, ' Math. Comp. 34 (1980), 305-312 'tested : FreeBasic 1.08.1 with gmp 6.2.0 '---------------------------------------------------
- include "gmp.bi"
'multi-precision float pointers Dim as mpf_ptr a, b Dim shared as mpf_ptr k2, u, v 'unsigned long integers Dim as ulong k, n, n2, r, s, t 'precision parameters Dim shared as ulong e10, e2 Dim shared e as clong Dim shared f as double Dim as double tim = TIMER CLS
a = allocate(len(__mpf_struct)) b = allocate(len(__mpf_struct)) u = allocate(len(__mpf_struct)) v = allocate(len(__mpf_struct)) k2 = allocate(len(__mpf_struct))
'log(x/y) with the Taylor series for atanh(x-y/x+y) Sub ln (byval s as mpf_ptr, byval x as ulong, byval y as ulong) Dim as mpf_ptr d = u, q = v Dim k as ulong
'Möbius transformation k = x: x -= y: y += k
If x <> 1 Then
Print "ln: illegal argument x - y <> 1"
End
End If
's = 1 / (x + y) mpf_set_ui (s, y) mpf_ui_div (s, 1, s) 'k2 = s * s mpf_mul (k2, s, s) mpf_set (d, s)
k = 1
Do
k += 2
'd *= k2
mpf_mul (d, d, k2)
'q = d / k
mpf_div_ui (q, d, k)
's += q
mpf_add (s, s, q)
f = mpf_get_d_2exp (@e, q) Loop until abs(e) > e2
's *= 2 mpf_mul_2exp (s, s, 1)
End Sub
'Main
'n = 2^i * 3^j * 5^k
'log(n) = r * log(16/15) + s * log(25/24) + t * log(81/80)
'solve linear system for r, s, t ' 4 -3 -4| i '-1 -1 4| j '-1 2 -1| k
'examples t = 1 select case t case 1
n = 60 r = 41 s = 30 t = 18
'100 digits case 2
n = 4800 r = 85 s = 62 t = 37
'8000 digits, 0.6 s case 3
n = 9375 r = 91 s = 68 t = 40
'15625 digits, 2.5 s case else
n = 18750 r = 98 s = 73 t = 43
'31250 digits, 12 s. @2.00GHz end select
'decimal precision e10 = n / .6 'binary precision e2 = (1 + e10) / .30103
'initialize mpf's mpf_set_default_prec (e2) mpf_inits (a, b, u, v, k2, Cptr(mpf_ptr, 0))
'Compute log terms
ln b, 16, 15
'a = r * b mpf_mul_ui (a, b, r)
ln b, 25, 24
'a += s * b mpf_mul_ui (u, b, s) mpf_add (a, a, u)
ln b, 81, 80
'a += t * b mpf_mul_ui (u, b, t) mpf_add (a, a, u)
gmp_printf (!"log(%lu) %.*Ff\n", n, e10, a)
'B&M, algorithm B1
'a = -a, b = 1 mpf_neg (a, a) mpf_set_ui (b, 1) mpf_set (u, a) mpf_set (v, b)
k = 0 n2 = n * n 'k2 = k * k mpf_set_ui (k2, 0) do
'k2 += 2k + 1 mpf_add_ui (k2, k2, (k shl 1) + 1) k += 1
'b = b * n2 / k2 mpf_div (b, b, k2) mpf_mul_ui (b, b, n2) 'a = (a * n2 / k + b) / k mpf_div_ui (a, a, k) mpf_mul_ui (a, a, n2) mpf_add (a, a, b) mpf_div_ui (a, a, k)
'u += a, v += b mpf_add (u, u, a) mpf_add (v, v, b)
f = mpf_get_d_2exp (@e, a)
Loop until abs(e) > e2
mpf_div (u, u, v) gmp_printf (!"gamma %.*Ff (maxerr. 1e-%lu)\n", e10, u, e10)
gmp_printf (!"k = %lu\n\n", k)
gmp_printf (!"time: %.7f s\n", TIMER - tim) end</lang>
- output:
gamma 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495 (maxerr. 1e-100) k = 255
The easy way
<lang freebasic>' ****************************************** 'Subject: Euler's constant 0.5772... 'tested : FreeBasic 1.08.1 with mpfr 4.1.0 '-------------------------------------------
- include "gmp.bi"
- include "mpfr.bi"
dim as mpfr_ptr a = allocate(len(__mpfr_struct)) dim as ulong e2, e10 dim as double tim = TIMER
'decimal precision e10 = 100
'binary precision e2 = (1 + e10) / .30103 mpfr_init2 (a, e2)
mpfr_const_euler (a, MPFR_RNDN) mpfr_printf (!"gamma %.*Rf\n\n", e10, a)
gmp_printf (!"time: %.7f s\n", TIMER - tim) end</lang>
jq
Translated from C (Bailey, 1988)
Works with gojq, the Go implementation of jq
<lang jq># Bailey, 1988 def bailey($n; $eps):
pow(2; $n) as $n2
| {a :1, b: 0, h: 1, r: 1, k: 1}
| until( (.b - .a)|fabs <= $eps;
.k += 1
| .r *= ($n2 / .k)
| .h += (1.0 / .k)
| .b = .a
| .a += (.r * .h) )
| (.a * $n2 / ($n2|exp) ) - ($n * (2|log)) ;</lang>
- Output:
bailey(5; 1E-6) 0.5772156649015341
Julia
<lang julia>display(MathConstants.γ) # γ = 0.5772156649015... </lang>
PARI/GP
built-in: <lang parigp>\l "euler_const.log" \p 100 print("gamma ", Euler); \q</lang>
Perl
<lang perl>#!/usr/bin/perl
use strict; # https://en.wikipedia.org/wiki/Euler%27s_constant use warnings; use List::Util qw( sum );
print sum( map 1 / $_, 1 .. 1e6) - log 1e6, "\n";</lang>
- Output:
0.577216164900715
Phix
-- demo\rosetta\Eulers_constant.exw without js -- no mpfr_get_d_2exp() or mpfr_const_euler() in mpfr.js as yet -- part 1, translation of Perl, with the same inaccuracy constant C = sum(sq_div(1,tagset(1e6)))-log(1e6) printf(1,"gamma %.12f (max 12d.p. of accuracy)\n",C) -- part 2, translation of C, from first principles. requires("1.0.1") -- mpfr_get_d_2exp(), mpfr_const_euler() include mpfr.e mpfr u, v, k2; integer e, e10, e2 atom f //log(x/y) with the Taylor series for atanh(x-y/x+y) procedure ln(mpfr s, integer x, y) mpfr d = u, q = v; assert((x-y)==1) y += x mpfr_set_si(s, y) mpfr_si_div(s, 1, s) // s = 1 / (x + y) mpfr_mul(k2, s, s) // k2 = s * s mpfr_set(d, s) integer k = 1 while true do k += 2; mpfr_mul(d, d, k2) // d *= k2 mpfr_div_si(q, d, k) // q = d / k mpfr_add(s, s, q) // s += q {f,e} = mpfr_get_d_2exp(q) if abs(e)>=e2 then exit end if end while mpfr_mul_si(s, s, 2) //s *= 2 end procedure mpfr a, b integer k, n = 60, -- (required precision in decimal dp *6/10) n2, r = 41, s = 30, t = 18; // n = 2^i * 3^j * 5^k // log(n) = r * log(16/15) + s * log(25/24) + t * log(81/80) // solve linear system for r, s, t // 4 -3 -4| i // -1 -1 4| j // -1 2 -1| k //decimal precision e10 = floor(n/0.6) //binary precision e2 = floor((1 + e10) / 0.30103) mpfr_set_default_precision(e2) {a, b, u, v, k2} = mpfr_inits(5) //Compute log terms ln(b, 16, 15) mpfr_mul_si(a, b, r) // a = r * b ln(b, 25, 24) mpfr_mul_si(u, b, s) mpfr_add (a, a, u) // a += s * b ln(b, 81, 80) mpfr_mul_si(u, b, t) mpfr_add (a, a, u) // a += t * b mpfr_neg(a, a) // a = -a mpfr_set_si(b, 1) // b = 1 mpfr_set (u, a) mpfr_set (v, b) k = 0; n2 = n * n; //k2 = k * k mpfr_set_si(k2, 0) while true do mpfr_add_si(k2, k2, k*2+1) // k2 += 2k + 1 k += 1; mpfr_div(b, b, k2) mpfr_mul_si(b, b, n2) // b = b * n2 / k2 mpfr_div_si(a, a, k) mpfr_mul_si(a, a, n2) mpfr_add (a, a, b) mpfr_div_si(a, a, k) // a = (a * n2 / k + b) / k mpfr_add(u, u, a) // u += a mpfr_add(v, v, b) // v += b {f,e} = mpfr_get_d_2exp (a) if abs(e)>=e2 then exit end if end while mpfr_div(u, u, v) string su = mpfr_get_fixed(u,e10) printf(1,"gamma %s (maxerr. 1e-%d)\n", {su, e10}) -- part 3, the easy way mpfr gamma = mpfr_init(0,-100) mpfr_const_euler(gamma) printf(1,"gamma %s (mpfr_const_euler)\n",{mpfr_get_fixed(gamma,100)})
- Output:
gamma 0.577216164901 (max 12d.p. of accuracy) gamma 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467494 (maxerr. 1e-100) gamma 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495 (mpfr_const_euler)
Processing
<lang Processing> /*********************************************
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");
} </lang>
- Output:
From the definition, err. 3e-10 Hn 6.569929751800373 gamma 0.577215823577717 k = 400 Sweeney, 1963, err. idem Hn 6.569929751800373 gamma 0.5772155784070492 k = 68 Bailey, 1988 gamma 0.5772176 k = 67 Brent-McMillan, 1980 gamma 0.5772157269353247 k 40 How Euler did it in 1735 Hn 2.9289682805538177 -ln 0.6263831555843353 err -0.044995839604388355 gamma 0.581387315979947 k = 17 C = 0.57721566490153286...
Raku
<lang perl6># 20211124 Raku programming solution
sub gamma (\N) {
# my $sum = (1/2 - 1/3).FatRat; # for arbitrary precision
my $sum = (1/2 - 1/3);
return $sum + [+] (2..N).race.map: -> \n {
my $power = 2**n;
my $sign = -1;
# my FatRat $term; # for arbitrary precision
# for ($power..^2*$power) { $term += (1 / (($sign *= -1)*$_)).FatRat }
my $term;
for ($power..^2*$power) { $term += ($sign = -$sign) / $_ }
n*$term }
}
say gamma 23 ;</lang>
- Output:
0.5772149198434515
Wren
Single precision (Cli)
Note that, whilst internally double arithmetic is carried out to the same precision as C (Wren is written in C), printing doubles is effectively limited to a maximum of 14 decimal places. <lang ecmascript>import "./fmt" for Fmt
var eps = 1e-6
System.print("From the definition, err. 3e-10") var n = 400 var h = 1 for (k in 2..n) h = h + 1/k //faster convergence: Negoi, 1997 var a = (n + 0.5 + 1/(24*n)).log
Fmt.print("Hn $0.14f", h) Fmt.print("gamma $0.14f\nk = $d\n", h - a, n)
System.print("Sweeney, 1963, err. idem") n = 21 var s = [0, n] var r = n var k = 1 while (true) {
k = k + 1 r = r * n / k s[k & 1] = s[k & 1] + r/k if (r <= eps) break
} Fmt.print("gamma $0.14f\nk = $d\n", s[1] - s[0] - n.log, k)
System.print("Bailey, 1988") n = 5 a = 1 h = 1 var n2 = 2.pow(n) r = 1 k = 1 while (true) {
k = k + 1 r = r * n2 / k h = h + 1/k var b = a a = a + r * h if ((b-a).abs <= eps) break
} a = a * n2 / n2.exp Fmt.print("gamma $0.14f\nk = $d\n", a - n * 2.log, k)
System.print("Brent-McMillan, 1980") n = 13 a = -n.log var b = 1 var u = a var v = b n2 = n * n var k2 = 0 k = 0 while (true) {
k2 = k2 + 2*k + 1 k = k + 1 a = a * n2 / k b = b * n2 / k2 a = (a + b)/k u = u + a v = v + b if (a.abs <= eps) break
} Fmt.print("gamma $0.14f\nk = $d\n", u / v, k)
System.print("How Euler did it in 1735") // Bernoulli numbers with even indices var b2 = [1, 1/6, -1/30, 1/42, -1/30, 5/66, -691/2730, 7/6, -3617/510, 43867/798] var m = 7 n = 10 // n'th harmonic number h = 1 for (k in 2..n) h = h + 1/k Fmt.print("Hn $0.14f", h) h = h - n.log Fmt.print(" -ln $0.14f", h) // expansion C = -digamma(1) a = -1 / (2*n) n2 = n * n r = 1 for (k in 1..m) {
r = r * n2 a = a + b2[k] / (2*k*r)
} Fmt.print("err $0.14f\ngamma $0.14f\nk = $d", a, h + a, n + m) System.print("\nC = 0.57721566490153286...")</lang>
- Output:
From the definition, err. 3e-10 Hn 6.56992969117651 gamma 0.57721566457657 k = 400 Sweeney, 1963, err. idem gamma 0.57721566456363 k = 68 Bailey, 1988 gamma 0.57721566490154 k = 89 Brent-McMillan, 1980 gamma 0.57721566490153 k = 40 How Euler did it in 1735 Hn 2.92896825396825 -ln 0.62638316097421 err -0.04916749607268 gamma 0.57721566490153 k = 17 C = 0.57721566490153286...
Multi precision (Embedded)
The display is limited to 100 digits for all four examples as I couldn't see much point in showing them all. <lang ecmascript>import "./gmp" for Mpf
var euler = Fn.new { |n, r, s, t|
// decimal precision var e10 = (n/0.6).floor
// binary precision var e2 = ((1 + n/0.6)/0.30103).round Mpf.defaultPrec = e2
var b = Mpf.new().log(Mpf.from(16).div(15))
var a = b.mul(r)
b = Mpf.new().log(Mpf.from(25).div(24))
a.add(b.mul(s))
b = Mpf.new().log(Mpf.from(81).div(80))
var u = b * t
a.add(u).neg
b.set(1)
u.set(a)
var v = Mpf.from(b)
var k = 0
var n2 = n * n
var k2 = Mpf.zero
while (true) {
k2.add((k << 1) + 1)
k = k + 1
b.mul(n2).div(k2)
a.mul(n2).div(k).add(b).div(k)
u.add(a)
v.add(b)
var e = Mpf.frexp(a)[1]
if (e.abs >= e2) break
}
u.div(v)
System.print("gamma %(u.toString(10, 100)) (maxerr. 1e-%(e10))")
System.print("k = %(k)")
}
var start = System.clock euler.call(60, 41, 30, 18) euler.call(4800, 85, 62, 37) euler.call(9375, 91, 68, 40) euler.call(18750, 98, 73, 43) System.print("\nTook %(System.clock - start) seconds.")</lang>
- Output:
gamma 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495 (maxerr. 1e-100) k = 255 gamma 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495 (maxerr. 1e-8000) k = 20462 gamma 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495 (maxerr. 1e-15625) k = 39967 gamma 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495 (maxerr. 1e-31250) k = 79936 Took 2.330538 seconds.
The easy way (Embedded)
<lang ecmascript>import "./gmp" for Mpf
var prec = (101/0.30103).round var gamma = Mpf.euler(prec) System.print(gamma.toString(10, 100))</lang>
- Output:
0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495