I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

# Euler's constant 0.5772...

Euler's constant 0.5772...
You are encouraged to solve this task according to the task description, using any language you may know.

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)

## Burlesque

Subject to double rounding error, evaluates 33 partial terms in a series

`33ro{JroJJJL[\/JL[\/nr\/{-1}\/.*FLJL[ro?^?*\/JL[{2}\/.*FL\/?^?*\/{1.+}m[J{lg}m[\/?/?*++\/1 3.0./\/?^.*}m[++-2 3 2lg.*./.*2lg2./.+`
Output:
`0.5772156649015326`

## C

### Single precision

`/*********************************************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, 1997a = 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 indicesdouble 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 numberh = 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");}`
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

`/**************************************************Subject: Computation of Euler's constant 0.5772...         with the Brent-McMillan algorithm B1,         Math. Comp. 34 (1980), 305-312tested : tcc-0.9.27 with gmp 6.2.0-------------------------------------------------*/#include <gmp.h>#include <stdio.h>#include <stdlib.h>#include <time.h> //multi-precision float pointersmpf_ptr u, v, k2; //precision parametersunsigned 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 integersunsigned 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 //examplest = 1;switch (t) {case 1 :   n = 60;   r = 41;   s = 30;   t = 18;//100 digitsbreak;case 2 :   n = 4800;   r = 85;   s = 62;   t = 37;//8000 digits, 0.6 sbreak;case 3 :   n = 9375;   r = 91;   s = 68;   t = 40;//15625 digits, 2.5 sbreak;default :   n = 18750;   r = 98;   s = 73;   t = 43;//31250 digits, 12 s. @2.00GHz} //decimal precisione10 = n / .6;//binary precisione2 = (1 + e10) / .30103; //initialize mpf'smpf_set_default_prec (e2);mpf_inits (a, b, u, v, k2, (mpf_ptr)0); //Compute log terms ln (b, 16, 15); //a = r * bmpf_mul_ui (a, b, r); ln (b, 25, 24); //a += s * bmpf_mul_ui (u, b, s);mpf_add (a, a, u); ln (b, 81, 80); //a += t * bmpf_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 = 1mpf_neg (a, a);mpf_set_ui (b, 1);mpf_set (u, a);mpf_set (v, b); k = 0;n2 = n * n;//k2 = k * kmpf_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);}`
output:
```gamma 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495 (maxerr. 1e-100)
k = 255
```

### The easy way

`/*******************************************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 precisione10 = 100; //binary precisione2 = (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);}`

## FreeBASIC

### Single precision

`'**********************************************'Subject: Comparing five methods for'         computing Euler's constant 0.5772...'tested : FreeBasic 1.08.1'----------------------------------------------const eps = 1e-6dim as double a, b, h, n2, r, u, vdim as integer k, k2, m, n ? "From the definition, err. 3e-10" n = 400 h = 1for k = 2 to n   h += 1 / knext k'faster convergence: Negoi, 1997a = 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 = nk = 1do   k += 1   r *= n / k   s(k and 1) += r / kloop until r < eps ? "gamma"; s(1) - s(0) - log(n); !"\nk ="; k?  ? "Bailey, 1988" n = 5 a = 1h = 1n2 = 2^nr = 1k = 1do   k += 1   r *= n2 / k   h += 1 / k   b = a: a += r * hloop until abs(b - a) < epsa *= n2 / exp(n2) ? "gamma"; a - n * log(2); !"\nk ="; k?  ? "Brent-McMillan, 1980" n = 13 a = -log(n)b = 1u = av = bn2 = n * nk2 = 0k = 0do   k2 += 2*k + 1   k += 1   a *= n2 / k   b *= n2 / k2   a = (a + b) / k   u += a   v += bloop until abs(a) < eps ? "gamma"; u / v; !"\nk ="; k?  ? "How Euler did it in 1735"'Bernoulli numbers with even indicesdim as double B2(9) = {1,1/6,-1/30,1/42,_ -1/30,5/66,-691/2730,7/6,-3617/510,43867/798}m = 7if m > 9 then end n = 10 'n-th harmonic numberh = 1for k = 2 to n   h += 1 / knext k? "Hn   "; h h -= log(n)? "  -ln"; h 'expansion C = -digamma(1)a = -1 / (2*n)n2 = n * nr = 1for 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`
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

`'***************************************************'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 pointersDim as mpf_ptr a, bDim shared as mpf_ptr k2, u, v'unsigned long integersDim as ulong k, n, n2, r, s, t'precision parametersDim shared as ulong e10, e2Dim shared e as clongDim shared f as doubleDim as double tim = TIMERCLS 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 = vDim 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 'examplest = 1select case tcase 1   n = 60   r = 41   s = 30   t = 18'100 digitscase 2   n = 4800   r = 85   s = 62   t = 37'8000 digits, 0.6 scase 3   n = 9375   r = 91   s = 68   t = 40'15625 digits, 2.5 scase else   n = 18750   r = 98   s = 73   t = 43'31250 digits, 12 s. @2.00GHzend select 'decimal precisione10 = n / .6'binary precisione2 = (1 + e10) / .30103 'initialize mpf'smpf_set_default_prec (e2)mpf_inits (a, b, u, v, k2, Cptr(mpf_ptr, 0)) 'Compute log terms ln b, 16, 15 'a = r * bmpf_mul_ui (a, b, r) ln b, 25, 24 'a += s * bmpf_mul_ui (u, b, s)mpf_add (a, a, u) ln b, 81, 80 'a += t * bmpf_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 = 1mpf_neg (a, a)mpf_set_ui (b, 1)mpf_set (u, a)mpf_set (v, b) k = 0n2 = n * n'k2 = k * kmpf_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`
output:
```gamma 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495 (maxerr. 1e-100)
k = 255
```

### The easy way

`' ******************************************'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, e10dim as double tim = TIMER 'decimal precisione10 = 100 'binary precisione2 = (1 + e10) / .30103mpfr_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`

## jq

Translated from C (Bailey, 1988)

Works with: jq

Works with gojq, the Go implementation of jq

`# Bailey, 1988def 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)) ;`
Output:
```bailey(5; 1E-6)
0.5772156649015341
```

## Julia

Translation of: PARI/GP
`display(MathConstants.γ)  # γ = 0.5772156649015... `

## Nim

`import std/math const n = 1e6var result = 1.0 for i in 2..int(n):  result += 1/i echo result - ln(n) `
Output:
```0.5772161649007153
```

## PARI/GP

built-in:

`\l "euler_const.log"\p 100print("gamma ", Euler);\q`

## Perl

`#!/usr/bin/perl use strict; # https://en.wikipedia.org/wiki/Euler%27s_constantuse warnings;use List::Util qw( sum ); print sum( map 1 / \$_, 1 .. 1e6) - log 1e6, "\n";`
Output:
```0.577216164900715
```

## Phix

### part 1 (max 12dp)

Translation of Perl, with the same accuracy limitation

```-- demo\rosetta\Eulers_constant.exw
with javascript_semantics
constant C = sum(sq_div(1,tagset(1e6)))-log(1e6)
printf(1,"gamma %.12f  (max 12d.p. of accuracy)\n",C)
```
Output:
```gamma 0.577216164901  (max 12d.p. of accuracy)
```

### part 2 (first principles)

Translation of C, from first principles.

```without js  -- no mpfr_get_d_2exp() in mpfr.js as yet
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)
mpfr_set_si(s, x+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_addmul_si(a, b, s)  // a += s * b
ln(b, 81, 80)  mpfr_addmul_si(a, b, t)  // 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;
mpfr_set_si(k2, 0)      // k2 = k * k (as below)
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_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})
```
Output:
```gamma 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467494 (maxerr. 1e-100)
```

### part 3 (the easy way)

```without js  -- no mpfr_const_euler() in mpfr.js as yet
requires("1.0.1") -- mpfr_const_euler()
include mpfr.e
mpfr gamma = mpfr_init(0,-100)
mpfr_const_euler(gamma)
printf(1,"gamma %s (mpfr_const_euler)\n",{mpfr_get_fixed(gamma,100)})
```

(matches part 2 except for the 100th decimal place)

Output:
```gamma 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495 (mpfr_const_euler)
```

## Processing

Translation of: C
` /********************************************* 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");} `
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...
```

## Racket

`#lang racket/base (require math/number-theory         math/base) gamma.0 ;; if you want to work it out the hard way...(define (H n)  (for/sum ((i n)) (/ (add1 i)))) (define (g #:n (n 10) #:k (k 7))  (+ (- (H n)        (log n)        (/ (* n 2)))     (for/sum ((2k (in-range 2 (* 2 (add1 k)) 2)))       (/ (bernoulli-number 2k) (* (expt n 2k) 2k))))) (g)`
Output:
```0.5772156649015329
0.5772156649015324```

## Raku

`# 20211124 Raku programming solution  sub gamma (\N where N > 1) { # Vacca series https://w.wiki/4ybp                             # convert terms to FatRat for arbitrary precision   return  (1/2 - 1/3) + [+] (2..N).race.map: -> \n {       my (\$power, \$sign, \$term) = 2**n, -1;       for (\$power..^2*\$power) { \$term += (\$sign = -\$sign) / \$_ }       n*\$term   }} say gamma 23 ;`
Output:
```0.5772149198434515
```

## Ruby

`n = 1e6p (1..n).sum{ 1.0/_1 } - Math.log(n)`
Output:
```0.5772161649014507 #12 digits accurate
```

## Rust

Translation of: Raku
`// 20220322 Rust programming solution fn gamma(N: u32) -> f64 { // Vacca series https://w.wiki/4ybp     return 1f64 / 2f64 - 1f64 / 3f64        + ((2..=N).map(|n| {            let power: u32 = 2u32.pow(n);            let mut sign: f64 = -1f64;            let mut term: f64 = 0f64;             for denominator in power..=(2 * power - 1) {                sign *= -1f64;                term += sign / f64::from(denominator);            }             return f64::from(n) * term;        }))        .sum::<f64>();} fn main() {    println!("{}", gamma(23));}`
Output:
`0.5772149198434514`

## Sidef

Built-in:

`# 100 decimals of precisionlocal Num!PREC = 4*100say Num.EulerGamma`
Output:
```0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495
```

Several formulas:

`const n = (ARGV ? Num(ARGV[0]) : 50)       # number of iterations define ℯ = Num.edefine π = Num.pidefine γ = Num.EulerGamma func display(r, t) {    say "#{r}\terror: #{ '%.0g' % abs(r - t) }"} # Original definition of the Euler-Mascheroni constant, due to Euler (1731)display(sum(1..n, {|n| 1/n }) - log(n), γ) # Formula due to Euler (best convergence)display(harmfrac(n) - log(n) - 1/(2*n) - sum(1..n, {|k|    -bernoulli(2*k) / (2*k) / n**(2*k)}), γ) # Formula derived from the above formula of Euler,# using approximations of Bernoulli numbers.display(harmfrac(n) - log(n) - 1/(2*n) - sum(1..n, {|k|    (-1)**k * 4 * sqrt(π*k) * (π * ℯ)**(-2*k) * k**(2*k) / (2*k) / n**(2*k)}), γ) # Euler-Mascheroni constant, involving zeta(n)display(1 - sum(2..(n+1), {|n|    (zeta(n) - 1) / n}), γ) # Limit_{n->Infinity} zeta((n+1)/n) - n} = gammadisplay(zeta((n+1)/n) - n, γ) # Series due to Euler (1731).display(sum(2..(n+1), {|n|    (-1)**n * zeta(n) / n}),  γ) # Formula due to Euler in terms of log(2) and the odd zeta valuesdisplay(3/4 - log(2)/2 + sum(1..n, {|n|    (1 - 1/(2*n + 1)) * (zeta(2*n + 1) - 1)}), γ) # Formula due to Euler in terms of log(2) and the odd zeta values (VII)display(log(2) - sum(1..n, {|n|    zeta(2*n + 1) / (2*n + 1) / 2**(2*n)}), γ) # Formula due to Vacca (1910)display(sum(1..n, {|n|    (-1)**n * floor(log2(n)) / n}), γ)`
Output:
```0.58718233290127899894172100505421724484389898107	error: 0.01
0.57721566490153286060651209008240243104215933594	error: 2e-58
0.577201775274185974649592917416402750312879661543	error: 1e-05
0.577215664901532868991217643213156967011395887777	error: 8e-18
0.57867004101560321761330253540741574969540128177	error: 0.001
0.567507841748305076772446986005418728718501189232	error: 0.01
0.577215664901532860606512090082272222693164663104	error: 1e-31
0.57721566490153286060651209008240496797924366087	error: 3e-33
0.611009392556929160631547597803236563977259982583	error: 0.03
```

## Vlang

Translation of: C
`import mathconst eps = 1e-6fn main() {    //.5772    println("From the definition, err. 3e-10")    mut n := 400    mut h := 1.0    for k in 2..n+1 {        h += 1.0/f64(k)    }    //faster convergence: Negoi, 1997    mut a := math.log(f64(n) + 0.5 + 1.0/f64(24*n))     println("Hn    \${h:0.16f}")    println("gamma \${h-a:0.16f}\nk = \$n\n")     println("Sweeney, 1963, err. idem")    n = 21    mut s := [0.0, f64(n)]    mut r := f64(n)    mut k := 1    for {        k++        r *= f64(n) / f64(k)        s[k & 1] += r/f64(k)        if r <= eps {            break        }    }    println("gamma \${s[1] - s[0] - math.log(n):0.16f}\nk = \$k\n")     println("Bailey, 1988")    n = 5    a = 1.0    h = 1.0    mut n2 := math.pow(f64(2),f64(n))    r = 1    k = 1    for {        k++        r *= n2 / f64(k)        h += 1/f64(k)        b := a        a += r * h        if math.abs(b-a) <= eps {            break        }    }    a *= n2 / math.exp(n2)    println("gamma \${a - n * math.log(2):.16f}\nk = \$k\n")     println("Brent-McMillan, 1980")    n = 13    a = -math.log(n)    mut b := 1.0    mut u := a    mut v := b    n2 = n * n    mut k2 := 0    k = 0    for {        k2 += 2*k + 1        k++        a *= n2 / f64(k)        b *= n2 / f64(k2)        a = (a + b)/f64(k)        u += a        v += b        if math.abs(a) <= eps {            break        }    }    println("gamma \${u/v:0.16f}\nk = \$k\n")     println("How Euler did it in 1735")    // Bernoulli numbers with even indices    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    n = 10    // n'th harmonic number    h = 1.0    for kz in 2..n+1 {        h += 1.0/f64(kz)    }     println("Hn    \${h:0.16f}")    h -= math.log(n)    println("  -ln \${h:0.16f}")    // expansion C = -digamma(1)    a = -1.0 / (2.0*f64(n))    n2 = f64(n * n)    r = 1    for kq in 1..m+1 {        r *= n2        a += b2[kq] / (2.0*f64(kq)*r)    }    println("err  \${a:0.16f}\ngamma \${h+a:0.16f}\nk = \${n+m}")     println("\nC = 0.57721566490153286...")}`
Output:
`Same as C entry`

## Wren

Translation of: C

### Single precision (Cli)

Library: Wren-fmt

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.

`import "./fmt" for Fmt var eps = 1e-6 System.print("From the definition, err. 3e-10")var n = 400var h = 1for (k in 2..n) h = h + 1/k//faster convergence: Negoi, 1997var 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 = 21var s = [0, n]var r = nvar k = 1while (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 = 5a = 1h = 1var n2 = 2.pow(n)r = 1k = 1while (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.expFmt.print("gamma \$0.14f\nk = \$d\n", a - n * 2.log, k) System.print("Brent-McMillan, 1980")n = 13a = -n.logvar b = 1var u = avar v = bn2 = n * nvar k2 = 0k = 0while (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 indicesvar b2 = [1, 1/6, -1/30, 1/42, -1/30, 5/66, -691/2730, 7/6, -3617/510, 43867/798]var m = 7n = 10// n'th harmonic numberh = 1for (k in 2..n) h = h + 1/kFmt.print("Hn    \$0.14f", h)h = h - n.logFmt.print("  -ln \$0.14f", h)// expansion C = -digamma(1)a = -1 / (2*n)n2 = n * nr = 1for (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...")`
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)

Library: Wren-gmp

The display is limited to 100 digits for all four examples as I couldn't see much point in showing them all.

`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.clockeuler.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.")`
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)

`import "./gmp" for Mpf var prec = (101/0.30103).roundvar gamma = Mpf.euler(prec)System.print(gamma.toString(10, 100))`
Output:
```0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495
```