Rosetta Code

Euler's constant 0.5772...: Difference between revisions

Page Discussion

Euler's constant 0.5772... (view source)

Revision as of 18:38, 20 January 2024

23,085 bytes added ,  3 months ago
Added Easylang
(Added Easylang)
 
(38 intermediate revisions by 24 users not shown)
Line 31:
__TOC__
 
=={{header|BASIC}}==
 
==={{header|FreeBASIC}}===
====Single precision====
<syntaxhighlight 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</syntaxhighlight>
 
{{out|output}}
<pre>
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...
</pre>
 
====Multi precision====
From first principles
<syntaxhighlight 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</syntaxhighlight>
 
{{out|output}}
<pre>
gamma 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495 (maxerr. 1e-100)
k = 255
</pre>
 
====The easy way====
<syntaxhighlight 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</syntaxhighlight>
 
=={{header|Burlesque}}==
Subject to double rounding error, evaluates 33 partial terms in a series
<langsyntaxhighlight lang="burlesque">33ro{JroJJJL[\/JL[\/nr\/{-1}\/.*FLJL[ro?^?*\/JL[{2}\/.*FL\/?^?*\/{1.+}m[J{lg}m[\/?/?*++\/1 3.0./\/?^.*}m[++-2 3 2lg.*./.*2lg2./.+</langsyntaxhighlight>
 
{{out}}
Line 42 ⟶ 401:
=={{header|C}}==
===Single precision===
<langsyntaxhighlight lang="c">/*********************************************
Subject: Comparing five methods for
computing Euler's constant 0.5772...
Line 162 ⟶ 521:
 
printf("\n\nC = 0.57721566490153286...\n");
}</langsyntaxhighlight>
 
{{out|output}}
Line 196 ⟶ 555:
From first principles
 
<langsyntaxhighlight lang="c">/**************************************************
Subject: Computation of Euler's constant 0.5772...
with the Brent-McMillan algorithm B1,
Line 373 ⟶ 732:
tim = clock() - tim;
printf("time: %.7f s\n",((double)tim)/CLOCKS_PER_SEC);
}</langsyntaxhighlight>
 
{{out|output}}
Line 382 ⟶ 741:
 
===The easy way===
<langsyntaxhighlight lang="c">/*******************************************
Subject: Euler's constant 0.5772...
tested : tcc-0.9.27 with mpfr 4.1.0
Line 409 ⟶ 768:
tim = clock() - tim;
gmp_printf ("time: %.7f s\n",((double)tim)/CLOCKS_PER_SEC);
}</langsyntaxhighlight>
 
=={{header|C++}}==
<syntaxhighlight lang="cpp">#include <array>
#include <cmath>
#include <iomanip>
#include <iostream>
 
double ByVaccaSeries(int numTerms)
=={{header|FreeBASIC}}==
{
===Single precision===
// this method is simple but converges slowly
<lang freebasic>'**********************************************
// calculate gamma by:
'Subject: Comparing five methods for
// 1 * (1/2 - 1/3) +
' computing Euler's constant 0.5772...
// 2 * (1/4 - 1/5 + 1/6 - 1/7) +
'tested : FreeBasic 1.08.1
// 3 * (1/8 - 1/9 + 1/10 - 1/11 + 1/12 - 1/13 + 1/14 - 1/15) +
'----------------------------------------------
// 4 * ( . . . ) +
const eps = 1e-6
dim as double a, b,// h,. n2,. r, u, v.
double gamma = 0;
dim as integer k, k2, m, n
size_t next = 4;
for(double numerator = 1; numerator < numTerms; ++numerator)
{
double delta = 0;
for(size_t denominator = next/2; denominator < next; denominator+=2)
{
// calculate terms two at a time
delta += 1.0/denominator - 1.0/(denominator + 1);
}
 
gamma += numerator * delta;
? "From the definition, err. 3e-10"
next *= 2;
}
return gamma;
}
 
// based on the C entry
n = 400
double ByEulersMethod()
{
//Bernoulli numbers with even indices
const std::array<double, 8> B2 {1.0, 1.0/6, -1.0/30, 1.0/42, -1.0/30,
5.0/66, -691.0/2730, 7.0/6};
const int n = 10;
 
//n-th harmonic number
h = 1
const double h = [] // immediately invoked lambda
for k = 2 to n
h += 1 / k{
double sum = 1;
next k
for (int k = 2; k <= n; k++) { sum += 1.0 / k; }
'faster convergence: Negoi, 1997
a = log(n +.5 + 1 / return sum - log(24*n));
}();
 
//expansion C = -digamma(1)
? "Hn "; h
double a = -1.0 / (2*n);
? "gamma"; h - a; !"\nk ="; n
double r = 1;
?
for (int k = 1; k < ssize(B2); k++)
{
r *= n * n;
a += B2[k] / (2*k * r);
}
 
return h + a;
}
 
? "Sweeney, 1963, err. idem"
 
int main()
n = 21
{
 
std::cout << std::setprecision(16) << "Vacca series: " << ByVaccaSeries(32);
dim as double s(1) = {0, n}
std::cout << std::setprecision(16) << "\nEulers method: " << ByEulersMethod();
r = n
k} = 1
</syntaxhighlight>
do
{{out}}
k += 1
<pre>
r *= n / k
Vacca series: 0.5772156610598274
s(k and 1) += r / k
Eulers method: 0.5772156649015325
loop until r < eps
</pre>
 
? "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
 
=={{header|Clojure}}==
? "err "; a; !"\ngamma"; h + a; !"\nk ="; n + m
<syntaxhighlight lang="clojure">
?
(let [n 1e10]
? "C = 0.57721566490153286..."
(loop [i 1
end</lang>
out (- (Math/log n))]
(if (<= i n)
(recur (inc i) (+ out (/ 1.0 i)))
out)))</syntaxhighlight>
 
{{out|output}}
<pre>
Gives: 0.5772156649484047
From the definition, err. 3e-10
Compared to the true value: 0.5772156649015328...
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...
</pre>
 
=={{header|Delphi}}==
===Multi precision===
{{works with|Delphi|6.0}}
From first principles
{{libheader|SysUtils,StdCtrls}}
<lang freebasic>'***************************************************
This programs demonstrates the basic method of calculating the Euler number. It illustrates how the number of iterations increases the accuracy.
'Subject: Computation of Euler's constant 0.5772...
<syntaxhighlight lang="Delphi">
' 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"
 
function ComputeEuler(N: int64): double;
'multi-precision float pointers
{Compute Eurler number with N-number of iterations}
Dim as mpf_ptr a, b
var I: integer;
Dim shared as mpf_ptr k2, u, v
var A: double;
'unsigned long integers
begin
Dim as ulong k, n, n2, r, s, t
Result:=0;
'precision parameters
for I:=1 to N-1 do
Dim shared as ulong e10, e2
Result:=Result + 1 / I;
Dim shared e as clong
A:=Ln(N + 0.5 + 1/(24.0*N));
Dim shared f as double
Result:=Result-A;
Dim as double tim = TIMER
end;
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
 
procedure ShowEulersNumber(Memo: TMemo);
If x <> 1 Then
{Show Euler numbers at various levels of precision}
Print "ln: illegal argument x - y <> 1"
var Euler,G,A,Error: double;
End
var N: integer;
End If
const Correct =0.57721566490153286060651209008240243104215933593992;
 
procedure ShowEulerError(N: int64);
's = 1 / (x + y)
{Show Euler number and Error}
mpf_set_ui (s, y)
begin
mpf_ui_div (s, 1, s)
Euler:=ComputeEuler(N);
'k2 = s * s
Error:=Correct-Euler;
mpf_mul (k2, s, s)
Memo.Lines.Add('N = '+FloatToStrF(N,ffNumber,18,0));
mpf_set (d, s)
Memo.Lines.Add('Euler='+FloatToStrF(Euler,ffFixed,18,18));
Memo.Lines.Add('Error='+FloatToStrF(Error,ffFixed,18,18));
Memo.Lines.Add('');
end;
 
begin
k = 1
{Compute Euler number with iterations ranging 10 to 10^9}
Do
for N:=1 to 9 do
k += 2
begin
'd *= k2
ShowEulerError(Trunc(Power(10,N)));
mpf_mul (d, d, k2)
end;
'q = d / k
mpf_div_ui (q, d, k)
's += q
mpf_add (s, s, q)
 
end;
f = mpf_get_d_2exp (@e, q)
Loop until abs(e) > e2
 
's *= 2
mpf_mul_2exp (s, s, 1)
End Sub
 
'Main
 
</syntaxhighlight>
'n = 2^i * 3^j * 5^k
{{out}}
<pre>
N = 10
Euler=0.477196250122325250
Error=0.100019414779207616
 
N = 100
'log(n) = r * log(16/15) + s * log(25/24) + t * log(81/80)
Euler=0.567215644212103243
Error=0.010000020689429618
 
N = 1,000
'solve linear system for r, s, t
Euler=0.576215664880711742
' 4 -3 -4| i
Error=0.001000000020821119
'-1 -1 4| j
'-1 2 -1| k
 
N = 10,000
'examples
Euler=0.577115664901475256
t = 1
Error=0.000100000000057605
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
 
N = 100,000
'decimal precision
Euler=0.577205664901386584
e10 = n / .6
Error=0.000010000000146277
'binary precision
e2 = (1 + e10) / .30103
 
N = 1,000,000
'initialize mpf's
Euler=0.577214664900591146
mpf_set_default_prec (e2)
Error=0.000001000000941715
mpf_inits (a, b, u, v, k2, Cptr(mpf_ptr, 0))
 
N = 10,000,000
'Compute log terms
Euler=0.577215564898391875
Error=0.000000100003140985
 
N = 100,000,000
ln b, 16, 15
Euler=0.577215654895336883
Error=0.000000010006195978
 
'aN = r * b1,000,000,000
Euler=0.577215664060567235
mpf_mul_ui (a, b, r)
Error=0.000000000840965626
 
Elapsed Time: 4.716 Sec.
ln b, 25, 24
</pre>
 
'a += s * b
mpf_mul_ui (u, b, s)
mpf_add (a, a, u)
 
=={{header|EasyLang}}==
ln b, 81, 80
<syntaxhighlight>
fastfunc gethn n .
i = 1
while i <= n
hn += 1 / i
i += 1
.
return hn
.
e = 2.718281828459045235
n = 10e8
numfmt 9 0
print gethn n - log10 n / log10 e
</syntaxhighlight>
 
=={{header|FutureBasic}}==
'a += t * b
<syntaxhighlight lang="futurebasic">
mpf_mul_ui (u, b, t)
void local fn Euler
mpf_add (a, a, u)
long n = 10000000, k
double a, h = 1.0
for k = 2 to n
h += 1 / k
next
a = log( n + 0.5 + 1 / ( 24 * n ) )
printf @"From the definition, err. 3e-10"
printf @"Hn %.15f", h
printf @"gamma %.15f", h - a
printf @"k = %ld\n", n
printf @"C %.15f", 0.5772156649015328
end fn
 
CFTimeInterval t
''gmp_printf (!"log(%lu) %.*Ff\n", n, e10, a)
t = fn CACurrentMediaTime
fn Euler
printf @"\vCompute time: %.3f ms",(fn CACurrentMediaTime-t)*1000
 
HandleEvents
'B&M, algorithm B1
</syntaxhighlight>
{{output}}
<pre>
From the definition, err. 3e-10
Hn 16.695311365857272
gamma 0.577215664898954
k = 10000000
 
C 0.577215664901533
'a = -a, b = 1
Compute time: 67.300 ms
mpf_neg (a, a)
</pre>
mpf_set_ui (b, 1)
mpf_set (u, a)
mpf_set (v, b)
 
=={{header|J}}==
k = 0
We can approximate euler's constant using the difference between the reciprocal and the gamma function of a small number:
n2 = n * n
<syntaxhighlight lang=J> (% - !@<:) 2^_27
'k2 = k * k
0.577216</syntaxhighlight>
mpf_set_ui (k2, 0)
do
'k2 += 2k + 1
mpf_add_ui (k2, k2, (k shl 1) + 1)
k += 1
 
For higher accuracy we can use a variation on the [[#C|C]] implementation which was attributed to [https://www.ams.org/journals/mcom/1980-34-149/S0025-5718-1980-0551307-4/S0025-5718-1980-0551307-4.pdf Richard P. Brent and Edwin M. McMillan]:
'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)
 
<syntaxhighlight lang=J>Euler=: {{
'u += a, v += b
A=.B=. ^.1r13 1x1
mpf_add (u, u, a)
r=. j=. 0
mpf_add (v, v, b)
whilst. (r=.%/B)~:!.0(r) do.
B=. B+A=. (j,1)%~+/\.A*169%(1,j)*(j=.j+1)
end.
r
}}0
</syntaxhighlight>
 
With J configured for 16 digits of display precision (<code>9!:11(16)</code>) we can see that this gave us:
f = mpf_get_d_2exp (@e, a)
Loop until abs(e) > e2
 
<syntaxhighlight lang=J> Euler
mpf_div (u, u, v)
0.5772156649015329</syntaxhighlight>
gmp_printf (!"gamma %.*Ff (maxerr. 1e-%lu)\n", e10, u, e10)
 
=={{header|Java}}==
gmp_printf (!"k = %lu\n\n", k)
<syntaxhighlight lang="java">
/**
* Using a simple formula derived from Hurwitz zeta function,
* as described on https://en.wikipedia.org/wiki/Euler%27s_constant,
* gives a result accurate to 12 decimal places.
*/
public class EulerConstant {
 
public static void main(String[] args) {
gmp_printf (!"time: %.7f s\n", TIMER - tim)
System.out.println(gamma(1_000_000));
end</lang>
}
 
{{out|output}}
<pre>
gamma 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495 (maxerr. 1e-100)
k = 255
</pre>
 
===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>
 
private static double gamma(int N) {
double gamma = 0.0;
for ( int n = 1; n <= N; n++ ) {
gamma += 1.0 / n;
}
gamma -= Math.log(N) + 1.0 / ( 2 * N );
return gamma;
}
}
</syntaxhighlight>
{{ out }}
<pre>0.5772156649007147</pre>
 
=={{header|jq}}==
Line 776 ⟶ 1,060:
'''Works with gojq, the Go implementation of jq'''
 
<langsyntaxhighlight lang="jq"># Bailey, 1988
def bailey($n; $eps):
pow(2; $n) as $n2
Line 786 ⟶ 1,070:
| .b = .a
| .a += (.r * .h) )
| (.a * $n2 / ($n2|exp) ) - ($n * (2|log)) ;</langsyntaxhighlight>
{{out}}
<pre>
Line 792 ⟶ 1,076:
0.5772156649015341
</pre>
 
 
=={{header|Julia}}==
{{trans|PARI/GP}}
<langsyntaxhighlight lang="julia">display(MathConstants.γ) # γ = 0.5772156649015...
</syntaxhighlight>
</lang>
 
=={{header|Lambdatalk}}==
Following the definition with an improvment from Negoi.
<syntaxhighlight lang="scheme">
{def negoi
{lambda {:n}
{let { {:n :n}
{:h {+ {S.map {lambda {:k} {/ 1 :k}} {S.serie 1 :n}}} }
{:a {log {+ :n 0.5 {/ 1 {* 24 :n}}}}} // Negoi, 1997
} {div}-> Hn :h
{div}gamma {- :h :a}
{div}k :n
}}}
-> negoi
 
{negoi 400}
-> Hn 6.5699296911765055
gamma 0.5772156645765731 with k = 400
(0.57721566457657 target)
</syntaxhighlight>
Following Sweeney
<syntaxhighlight lang="scheme">
{def sweeney
{def sweeney.set!
{lambda {:s :r :k :i}
{A.set! :i {+ {A.get :i :s} {/ :r :k}} :s}
}}
{def sweeney.loop
{lambda {:n :s :r :k}
{if {<= :r 1.e-10}
then gamma = {- {A.get 1 :s} {A.get 0 :s} {log :n}} with k=:k
else {sweeney.loop :n
{sweeney.set! :s {* :r {/ :n :k}} :k {% :k 2}}
{* :r {/ :n :k}}
{+ :k 1} }
}}}
{lambda {:n}
{sweeney.loop :n {A.new 0 :n} :n 2} }}
-> sweeney
 
{sweeney 21}
-> gamma = 0.577215664563631 with k=76
(0.57721566456363 target)
</syntaxhighlight>
 
=={{header|Lua}}==
A value using 100,000,000 iterations of the harmonic series computes in less than a second and is correct to eight decimal places.
<syntaxhighlight lang="lua">function computeGamma (iterations, decimalPlaces)
local Hn = 1
for i = 2, iterations do
Hn = Hn + (1/i)
end
local gamma = tostring(Hn - math.log(iterations))
return tonumber(gamma:sub(1, decimalPlaces + 2))
end
 
print(computeGamma(10^8, 8))</syntaxhighlight>
{{out}}
<pre>0.57721566</pre>
 
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">N[EulerGamma, 1000]</syntaxhighlight>
{{out}}
<pre>0.57721566490153286060651209008240243104215933593992359880576723488486772677766467093694706329174674
9514631447249807082480960504014486542836224173997644923536253500333742937337737673942792595258247094
9160087352039481656708532331517766115286211995015079847937450857057400299213547861466940296043254215
1905877553526733139925401296742051375413954911168510280798423487758720503843109399736137255306088933
1267600172479537836759271351577226102734929139407984301034177717780881549570661075010161916633401522
7893586796549725203621287922655595366962817638879272680132431010476505963703947394957638906572967929
6010090151251959509222435014093498712282479497471956469763185066761290638110518241974448678363808617
4945516989279230187739107294578155431600500218284409605377243420328547836701517739439870030237033951
8328690001558193988042707411542227819716523011073565833967348717650491941812300040654693142999297779
569303100503086303418569803231083691640025892970890985486825777364288253954925873629596133298574739302</pre>
 
=={{header|Maxima}}==
<syntaxhighlight lang="maxima">
%gamma,numer;
</syntaxhighlight>
{{out}}
<pre>
0.5772156649015329
</pre>
 
=={{header|Nim}}==
<syntaxhighlight lang="nim">import std/math
 
const n = 1e6
var result = 1.0
 
for i in 2..int(n):
result += 1/i
 
echo result - ln(n)
</syntaxhighlight>
{{out}}
<pre>
0.5772161649007153
</pre>
 
=={{header|PARI/GP}}==
built-in:
<langsyntaxhighlight lang="parigp">\l "euler_const.log"
\p 100
print("gamma ", Euler);
\q</langsyntaxhighlight>
 
=={{header|Perl}}==
<langsyntaxhighlight lang="perl">#!/usr/bin/perl
 
use strict; # https://en.wikipedia.org/wiki/Euler%27s_constant
Line 813 ⟶ 1,194:
use List::Util qw( sum );
 
print sum( map 1 / $_, 1 .. 1e6) - log 1e6, "\n";</langsyntaxhighlight>
{{out}}
<pre>
Line 822 ⟶ 1,203:
===part 1 (max 12dp)===
Translation of Perl, with the same accuracy limitation
<!--<langsyntaxhighlight Phixlang="phix">(phixonline)-->
<span style="color: #000080;font-style:italic;">-- demo\rosetta\Eulers_constant.exw</span>
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">C</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">tagset</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1e6</span><span style="color: #0000FF;">)))-</span><span style="color: #7060A8;">log</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1e6</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"gamma %.12f (max 12d.p. of accuracy)\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">C</span><span style="color: #0000FF;">)</span>
<!--</langsyntaxhighlight>-->
{{out}}
<pre>
Line 834 ⟶ 1,215:
===part 2 (first principles)===
Translation of C, from first principles.
<!--<langsyntaxhighlight Phixlang="phix">-->
<span style="color: #008080;">without</span> <span style="color: #008080;">js</span> <span style="color: #000080;font-style:italic;">-- no mpfr_get_d_2exp() in mpfr.js as yet</span>
<span style="color: #7060A8;">requires</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"1.0.2"</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- mpfr_get_d_2exp(), mpfr_addmul_si()</span>
Line 922 ⟶ 1,303:
<span style="color: #004080;">string</span> <span style="color: #000000;">su</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_get_fixed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">u</span><span style="color: #0000FF;">,</span><span style="color: #000000;">e10</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"gamma %s (maxerr. 1e-%d)\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">su</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">e10</span><span style="color: #0000FF;">})</span>
<!--</langsyntaxhighlight>-->
{{out}}
<pre>
Line 928 ⟶ 1,309:
</pre>
===part 3 (the easy way)===
<!--<langsyntaxhighlight Phixlang="phix">-->
<span style="color: #008080;">without</span> <span style="color: #008080;">js</span> <span style="color: #000080;font-style:italic;">-- no mpfr_const_euler() in mpfr.js as yet</span>
<span style="color: #7060A8;">requires</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"1.0.1"</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- mpfr_const_euler()</span>
Line 935 ⟶ 1,316:
<span style="color: #000000;">mpfr_const_euler</span><span style="color: #0000FF;">(</span><span style="color: #000000;">gamma</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"gamma %s (mpfr_const_euler)\n"</span><span style="color: #0000FF;">,{</span><span style="color: #7060A8;">mpfr_get_fixed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">gamma</span><span style="color: #0000FF;">,</span><span style="color: #000000;">100</span><span style="color: #0000FF;">)})</span>
<!--</langsyntaxhighlight>-->
<small>(matches part 2 except for the 100th decimal place)</small>
{{out}}
Line 941 ⟶ 1,322:
gamma 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495 (mpfr_const_euler)
</pre>
 
=={{header|Picat}}==
===List comprehension===
<syntaxhighlight lang="picat">main =>
Gamma = 0.57721566490153286060651209008240,
println(Gamma),
foreach(N in 1..8)
G = e(10**N),
println([n=N,g=G,diff=G-Gamma])
end.
 
e(N) = [1.0/I : I in 1..N].sum-log(N).</syntaxhighlight>
 
{{out}}
<pre>0.577215664901533
[n = 1,g = 0.626383160974208,diff = 0.049167496072675]
[n = 2,g = 0.582207331651529,diff = 0.004991666749996]
[n = 3,g = 0.577715581568206,diff = 0.000499916666674]
[n = 4,g = 0.577265664068165,diff = 0.000049999166632]
[n = 5,g = 0.577220664893106,diff = 0.000004999991574]
[n = 6,g = 0.577216164900715,diff = 0.000000499999182]
[n = 7,g = 0.577215714898951,diff = 0.000000049997419]
[n = 8,g = 0.577215669900188,diff = 0.000000004998655]</pre>
 
===Loop===
<syntaxhighlight lang="picat">e2(N) = E-log(N) =>
E = 1,
foreach(I in 2..N)
E := E + 1/I
end.</syntaxhighlight>
 
 
{{trans|Rust}}
<syntaxhighlight lang="picat">main =>
Gamma = 0.577215664901532860606512090082402,
println(gamma=Gamma),
member(N, 1..23),
G = gamma(N),
println([n=N,g=G,diff=G-Gamma]),
fail,
nl.
 
gamma(N) = Gamma =>
Gamma = 1/2 - 1/3,
foreach(I in 2..N)
Power = 2**I,
Sign = -1,
Term = 0,
foreach(Denominator in Power..(2*Power-1))
Sign := Sign * -1,
Term := Term + Sign / Denominator
end,
Gamma := Gamma + I*Term
end.</syntaxhighlight>
 
{{out}}
<pre>[n = 1,g = 0.166666666666667,diff = -0.410548998234866]
[n = 2,g = 0.314285714285714,diff = -0.262929950615819]
[n = 3,g = 0.416741591741592,diff = -0.160474073159941]
[n = 4,g = 0.482164184398886,diff = -0.095051480502647]
[n = 5,g = 0.522141654090275,diff = -0.055074010811257]
[n = 6,g = 0.545853770405349,diff = -0.031361894496184]
[n = 7,g = 0.559605750992416,diff = -0.017609913909116]
[n = 8,g = 0.567441138957738,diff = -0.009774525943794]
[n = 9,g = 0.571842107494027,diff = -0.005373557407506]
[n = 10,g = 0.574285301882304,diff = -0.002930363019229]
[n = 11,g = 0.57562856705805,diff = -0.001587097843483]
[n = 12,g = 0.576361123043496,diff = -0.000854541858037]
[n = 13,g = 0.576757887880701,diff = -0.000457777020832]
[n = 14,g = 0.576971520706463,diff = -0.00024414419507]
[n = 15,g = 0.577085964243776,diff = -0.000129700657756]
[n = 16,g = 0.577147000098518,diff = -0.000068664803014]
[n = 17,g = 0.577179425210813,diff = -0.00003623969072]
[n = 18,g = 0.577196591397621,diff = -0.000019073503912]
[n = 19,g = 0.577205651316587,diff = -0.000010013584946]
[n = 20,g = 0.57721041969158,diff = -0.000005245209953]
[n = 21,g = 0.577212923087556,diff = -0.000002741813977]
[n = 22,g = 0.577214234389975,diff = -0.000001430511558]
[n = 23,g = 0.577214919843452,diff = -0.000000745058081]</pre>
 
 
=={{header|Processing}}==
{{trans|C}}
<syntaxhighlight lang="processing">
<lang Processing>
/*********************************************
Subject: Comparing five methods for
Line 1,090 ⟶ 1,551:
println("C = 0.57721566490153286...\n");
}
</syntaxhighlight>
</lang>
{{out}}
<pre>
Line 1,123 ⟶ 1,584:
C = 0.57721566490153286...
</pre>
=={{header|Python}}==
<syntaxhighlight lang="python">
# /**************************************************
# Subject: Computation of Euler's constant 0.5772...
# with Euler's Zeta Series.
# tested : Python 3.11
# -------------------------------------------------*/
 
from scipy import special as s
 
def eulers_constant(n):
k = 2
euler = 0
while k <= n:
euler += (s.zeta(k) - 1)/k
k += 1
return 1 - euler
 
print(eulers_constant(47))
</syntaxhighlight>
{{out}}
<pre>
0.577215664901533
</pre>
=={{header|Racket}}==
 
<langsyntaxhighlight lang="racket">#lang racket/base
 
(require math/number-theory
Line 1,144 ⟶ 1,628:
(/ (bernoulli-number 2k) (* (expt n 2k) 2k)))))
 
(g)</langsyntaxhighlight>
 
{{out}}
Line 1,151 ⟶ 1,635:
 
=={{header|Raku}}==
<syntaxhighlight lang="raku" perl6line># 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;
# 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 ;</langsyntaxhighlight>
{{out}}
<pre>
0.5772149198434515
</pre>
 
=={{header|RPL}}==
{{trans|FreeBASIC}}
'''From the definition, with Negoi's convergence improvement'''
≪ → n
≪ 1
2 n '''FOR''' k k INV + '''NEXT'''
n .5 + 24 n * INV + LN -
≫ ≫ ‘<span style="color:blue">GAMMA1</span>’ STO
'''Sweeney method'''
≪ 0 OVER R→C 1E-6 → n s eps
≪ n 1
'''DO'''
1 + SWAP
OVER / n * SWAP
DUP2 /
IF OVER 2 MOD THEN (0,1) * END 's' STO+
'''UNTIL''' OVER eps < '''END'''
DROP2
s C→R SWAP - n LN -
≫ ≫ ‘<span style="color:blue">GAMMA2</span>’ STO
 
400 <span style="color:blue">GAMMA1</span>
21 <span style="color:blue">GAMMA2</span>
{{out}}
<pre>
2: .57721566456
1: .57717756228
</pre>
 
=={{header|Ruby}}==
<langsyntaxhighlight lang="ruby">n = 1e6
p (1..n).sum{ 1.0/_1 } - Math.log(n)</langsyntaxhighlight>
{{out}}
<pre>0.5772161649014507 #12 digits accurate
</pre>
 
=={{header|Rust}}==
{{trans|Raku}}
<syntaxhighlight lang="rust">// 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));
}</syntaxhighlight>
{{out}}
<pre>0.5772149198434514</pre>
 
=={{header|Scala}}==
<syntaxhighlight lang="Scala">
/**
* Using a simple formula derived from Hurwitz zeta function,
* as described on https://en.wikipedia.org/wiki/Euler%27s_constant,
* gives a result accurate to 11 decimal places: 0.57721566490...
*/
 
object EulerConstant extends App {
 
println(gamma(1_000_000))
 
private def gamma(N: Int): Double = {
val sumOverN = (1 to N).map(1.0 / _).sum
sumOverN - Math.log(N) - 1.0 / (2 * N)
}
 
}
</syntaxhighlight>
{{out}}
<pre>
0.5772156649007153
</pre>
 
=={{header|Scheme}}==
{{works with|Chez Scheme}}
<syntaxhighlight lang="scheme">; Procedure to compute factorial.
(define fact
(lambda (n)
(if (<= n 0)
1
(* n (fact (1- n))))))
 
; Compute Euler's gamma constant as the difference of log(n) from a sum.
; See section 2.3 of <http://numbers.computation.free.fr/Constants/Gamma/gamma.html>.
(define gamma
(lambda (n)
(let ((sum 0))
(do ((k 1 (1+ k)))
((> k (* 3.5911 n)) (- sum (log n)))
(set! sum (+ sum (/ (* (expt -1 (1- k)) (expt n k)) (* k (fact k)))))))))
 
; Show Euler's gamma constant computed at log(100).
(printf "Euler's gamma constant: ~a~%" (gamma 100))</syntaxhighlight>
{{out}}
<pre>
Euler's gamma constant: 0.5772156649015328
</pre>
 
 
 
=={{header|Sidef}}==
Built-in:
<syntaxhighlight lang="ruby"># 100 decimals of precision
local Num!PREC = 4*100
say Num.EulerGamma</syntaxhighlight>
{{out}}
<pre>
0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495
</pre>
 
Several formulas:
<syntaxhighlight lang="ruby">const n = (ARGV ? Num(ARGV[0]) : 50) # number of iterations
 
define ℯ = Num.e
define π = Num.pi
define γ = 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} = gamma
display(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 values
display(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
}), γ)</syntaxhighlight>
{{out}}
<pre>
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
</pre>
 
=={{header|V (Vlang)}}==
{{trans|C}}
<syntaxhighlight lang="v (vlang)">import math
const eps = 1e-6
fn 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...")
}</syntaxhighlight>
{{out}}
<pre>Same as C entry</pre>
 
=={{header|Wren}}==
{{trans|C}}
Line 1,189 ⟶ 1,953:
{{libheader|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.
<langsyntaxhighlight ecmascriptlang="wren">import "./fmt" for Fmt
 
var eps = 1e-6
Line 1,275 ⟶ 2,039:
}
Fmt.print("err $0.14f\ngamma $0.14f\nk = $d", a, h + a, n + m)
System.print("\nC = 0.57721566490153286...")</langsyntaxhighlight>
 
{{out}}
Line 1,309 ⟶ 2,073:
{{libheader|Wren-gmp}}
The display is limited to 100 digits for all four examples as I couldn't see much point in showing them all.
<langsyntaxhighlight ecmascriptlang="wren">import "./gmp" for Mpf
 
var euler = Fn.new { |n, r, s, t|
Line 1,352 ⟶ 2,116:
euler.call(9375, 91, 68, 40)
euler.call(18750, 98, 73, 43)
System.print("\nTook %(System.clock - start) seconds.")</langsyntaxhighlight>
 
{{out}}
Line 1,368 ⟶ 2,132:
</pre>
===The easy way (Embedded)===
<langsyntaxhighlight ecmascriptlang="wren">import "./gmp" for Mpf
 
var prec = (101/0.30103).round
var gamma = Mpf.euler(prec)
System.print(gamma.toString(10, 100))</langsyntaxhighlight>
 
{{out}}
<pre>
0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495
</pre>
 
=={{header|XPL0}}==
{{trans|C}}
<syntaxhighlight lang "XPL0">\*********************************************
\Subject: Comparing five methods for
\ computing Euler's constant 0.5772...
\---------------------------------------------
 
include xpllib; \for Print
define Epsilon = 1e-6;
 
real A, B, H, N2, R, U, V, S(2), B2;
int K, K2, M, N;
[Print("From the definition, error 3e-10\n");
N:= 400; H:= 1.;
for K:= 2 to N do
H:= H + 1.0/float(K);
\Faster convergence: Negoi, 1997
A:= Ln(float(N) + 0.5 + 1.0/(24.*float(N)));
Print("Hn %1.16f\n", H);
Print("gamma %1.16f\nK = %d\n\n", H-A, N);
 
Print("Sweeney, 1963, error 3e-10\n");
N:= 21; S(0):= 0.; S(1):= float(N);
R:= float(N); K:= 1;
repeat
K:= K+1;
R:= R * float(N) / float(K);
S(K&1):= S(K&1) + R/float(K);
until R <= Epsilon;
Print("gamma %1.16f\nK = %d\n\n", S(1)-S(0)-Ln(float(N)), K);
 
Print("Bailey, 1988\n");
N:= 5; A:= 1.; H:= 1.;
N2:= Pow(2., float(N));
R:= 1.; K:= 1;
repeat
K:= K+1;
R:= R * N2 / float(K);
H:= H + 1.0/float(K);
B:= A; A:= A + R*H;
until abs(B-A) <= Epsilon;
A:= A * N2 / Exp(N2);
Print("gamma %1.16f\nK = %d\n\n", A-float(N)*Ln(2.), K);
 
Print("Brent-McMillan, 1980\n");
N:= 13; A:= -Ln(float(N));
B:= 1.; U:= A; V:= B;
N2:= float(N*N); K2:= 0; K:= 0;
repeat
K2:= K2 + 2*K + 1;
K:= K+1;
A:= A * N2 / float(K);
B:= B * N2 / float(K2);
A:= (A + B) / float(K);
U:= U + A;
V:= V + B;
until abs(A) <= Epsilon;
Print("gamma %1.16f\nK = %d\n\n", U/V, K);
 
Print("How Euler did it in 1735\n");
\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;
\Nth harmonic number
H:= 1.;
for K:= 2 to N do
H:= H + 1.0/float(K);
Print("Hn %1.16f\n", H);
H:= H - Ln(float(N));
Print(" -ln %1.16f\n", H);
\Expansion C:= -digamma(1)
A:= -1.0 / (2.*float(N));
N2:= float(N*N);
R:= 1.;
for K:= 1 to M do [
R:= R * N2;
A:= A + B2(K)/(2.*float(K)*R);
];
Print("err %1.16f\ngamma %1.16f\nK = %d", A, H+A, N+M);
 
Print("\n\nC = 0.57721566490153286...\n");
]</syntaxhighlight>
{{out}}
<pre>
From the definition, error 3e-10
Hn 6.5699296911765100
gamma 0.5772156645765730
K = 400
 
Sweeney, 1963, error 3e-10
gamma 0.5772156645636310
K = 68
 
Bailey, 1988
gamma 0.5772156649015350
K = 89
 
Brent-McMillan, 1980
gamma 0.5772156649015330
K = 40
 
How Euler did it in 1735
Hn 2.9289682539682500
-ln 0.6263831609742080
err -0.0491674960726750
gamma 0.5772156649015330
K = 17
 
C = 0.57721566490153286...
</pre>
Chkas
1,983

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/317877...355238"