Euler's constant 0.5772...: Difference between revisions
Added Easylang
TechAscent (talk | contribs) No edit summary |
(Added Easylang) |
||
| (24 intermediate revisions by 12 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}}==
Line 501 ⟶ 860:
</pre>
=={{header|
{{works with|Delphi|6.0}}
{{libheader|SysUtils,StdCtrls}}
This programs demonstrates the basic method of calculating the Euler number. It illustrates how the number of iterations increases the accuracy.
<syntaxhighlight lang="Delphi">
function ComputeEuler(N: int64): double;
{Compute Eurler number with N-number of iterations}
var I: integer;
var A: double;
begin
Result:=0;
for I:=1 to N-1 do
Result:=Result + 1 / I;
A:=Ln(N + 0.5 + 1/(24.0*N));
Result:=Result-A;
end;
procedure ShowEulersNumber(Memo: TMemo);
{Show Euler numbers at various levels of precision}
var Euler,G,A,Error: double;
var N: integer;
const Correct =0.57721566490153286060651209008240243104215933593992;
procedure ShowEulerError(N: int64);
{Show Euler number and Error}
begin
Euler:=ComputeEuler(N);
Error:=Correct-Euler;
Memo.Lines.Add('N = '+FloatToStrF(N,ffNumber,18,0));
Memo.Lines.Add('Euler='+FloatToStrF(Euler,ffFixed,18,18));
Memo.Lines.Add('Error='+FloatToStrF(Error,ffFixed,18,18));
Memo.Lines.Add('');
end;
begin
{Compute Euler number with iterations ranging 10 to 10^9}
for N:=1 to 9 do
begin
ShowEulerError(Trunc(Power(10,N)));
end;
end;
</syntaxhighlight>
{{out}}
<pre>
N = 10
Euler=0.477196250122325250
Error=0.100019414779207616
N = 100
Euler=0.567215644212103243
Error=0.010000020689429618
Euler=0.576215664880711742
Error=0.001000000020821119
Euler=0.577115664901475256
Error=0.000100000000057605
N = 100,000
Euler=0.577205664901386584
Error=0.000010000000146277
N = 1,000,000
Euler=0.577214664900591146
Error=0.000001000000941715
N = 10,000,000
Euler=0.577215564898391875
Error=0.000000100003140985
N = 100,000,000
Euler=0.577215654895336883
Error=0.000000010006195978
N = 1,000,000,000
Euler=0.577215664060567235
Error=0.000000000840965626
Elapsed Time: 4.716 Sec.
</pre>
=={{header|EasyLang}}==
<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}}==
<syntaxhighlight lang="futurebasic">
void local fn Euler
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
t = fn CACurrentMediaTime
fn Euler
printf @"\vCompute time: %.3f ms",(fn CACurrentMediaTime-t)*1000
HandleEvents
</syntaxhighlight>
{{output}}
<pre>
From the definition, err. 3e-10
Hn
gamma 0.
k =
C 0.577215664901533
Compute time: 67.300 ms
</pre>
=={{header|J}}==
We can approximate euler's constant using the difference between the reciprocal and the gamma function of a small number:
<syntaxhighlight lang=
0.577216</syntaxhighlight>
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]:
<syntaxhighlight lang=J>Euler=: {{
A=.B=. ^.1r13 1x1
r=. j=. 0
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:
<syntaxhighlight lang=J> Euler
0.5772156649015329</syntaxhighlight>
=={{header|Java}}==
<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) {
System.out.println(gamma(1_000_000));
}
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 931 ⟶ 1,126:
</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}}==
Line 945 ⟶ 1,154:
8328690001558193988042707411542227819716523011073565833967348717650491941812300040654693142999297779
569303100503086303418569803231083691640025892970890985486825777364288253954925873629596133298574739302</pre>
=={{header|Maxima}}==
<syntaxhighlight lang="maxima">
%gamma,numer;
</syntaxhighlight>
{{out}}
<pre>
0.5772156649015329
</pre>
=={{header|Nim}}==
Line 1,366 ⟶ 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}}==
Line 1,412 ⟶ 1,653:
<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>
Line 1,448 ⟶ 1,718:
{{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}}
Line 1,472 ⟶ 1,767:
Euler's gamma constant: 0.5772156649015328
</pre>
=={{header|Sidef}}==
Line 1,656 ⟶ 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.
<syntaxhighlight lang="
var eps = 1e-6
Line 1,776 ⟶ 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.
<syntaxhighlight lang="
var euler = Fn.new { |n, r, s, t|
Line 1,835 ⟶ 2,132:
</pre>
===The easy way (Embedded)===
<syntaxhighlight lang="
var prec = (101/0.30103).round
Line 1,844 ⟶ 2,141:
<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>
| |||