Apéry's constant: Difference between revisions
m
→{{header|Phix}}: use pygments
(Created Nim solution.) |
m (→{{header|Phix}}: use pygments) |
||
| (9 intermediate revisions by 5 users not shown) | |||
Line 41:
;* [[oeis:A002117|OEIS:A002117 - Decimal expansion of zeta(3)]]
=={{header|ALGOL 68}}==
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}}
Uses Algol 68G's LONG LONG INT and LONG LONG REAL which have programmer specified precision. The factorials can get quuite large so more than 101 digits are needed.
<syntaxhighlight lang="algol68">
BEGIN # find Apéry's constant: the sum of the positive cubes' reciprocals #
# (this is the value of the Riemann zeta function applied to 3) #
PR precision 1000 PR # set precision of LONG LONG REAL #
# returns a string representation of z, truncated to 100 decimals #
PROC truncate100 = ( LONG LONG REAL z )STRING:
BEGIN
STRING result = fixed( z, 0, 101 ); # format with 101 decimals #
result[ : UPB result - 1 ] # remove the final digit #
END # truncate100 # ;
# methind 1 - sum the reciprocols of the cubes - 1000 terms from 1 #
BEGIN
LONG LONG REAL zeta3 := 0;
FOR k TO 1000 DO
LONG LONG INT llk = LENG LENG k;
zeta3 +:= 1 / ( llk * llk * llk )
OD;
print( ( truncate100( zeta3 ), newline ) )
END;
# method 2 - Markov's alternative representaton, 158 terms from 1 #
# 5/2 * sum [ (-1)^(k-1)( k!^2 / (2k!)k^2 ) ] from 1 #
BEGIN
LONG LONG INT fk := 1, f2k := 1;
LONG LONG REAL zeta3 := 0;
FOR k TO 158 DO
LONG LONG INT llk = k;
LONG LONG INT ll2k = llk * 2;
fk *:= llk;
f2k *:= ( ll2k - 1 ) *:= ll2k;
LONG LONG REAL term = ( fk * fk ) / ( f2k * llk * llk * llk );
IF ODD k THEN
zeta3 +:= term
ELSE
zeta3 -:= term
FI
OD;
zeta3 *:= 5 / 2;
print( ( truncate100( zeta3 ), newline ) )
END;
# method 3 - Wedeniwski representation - 20 terms from 0 #
# 1/24 * sum [ (-1)^k (2k+1)!^3(2k)!^3(k!)^3 #
# * ( 126392k^5 + 412708k^4 + 531578k^3 #
# + 336367k^2 + 104000k + 12463 #
# ) / (3k+2)! (4k + 3)!^3 #
# ] from 0 #
BEGIN
[]INT w coefficients = ( 126392, 412708, 531578, 336367, 104000, 12463 );
LONG LONG INT fk := 1, f2k := 1, f3k := 1, f4k := 1;
# ensure the divisor is a LONG LONG INT so the LHS is calculated as a #
# LONG LONG REAL value and not a REAL which is then widened #
LONG LONG REAL zeta3 := w coefficients[ UPB w coefficients ]
/ LENG LENG ( 2 * 6 * 6 * 6 );
FOR k TO 19 DO
LONG LONG INT llk = k;
LONG LONG INT ll2k = llk + llk;
LONG LONG INT ll3k = ll2k + llk;
LONG LONG INT ll4k = ll3k + llk;
fk *:= llk;
f2k *:= ( ll2k - 1 ) * ll2k;
f3k *:= ( ll3k - 2 ) * ( ll3k - 1 ) * ll3k;
f4k *:= ( ll4k - 3 ) * ( ll4k - 2 ) * ( ll4k - 1 ) * ll4k;
LONG LONG INT f2k1 = f2k * ( ll2k + 1 );
LONG LONG INT f3k2 = f3k * ( ll3k + 1 ) * ( ll3k + 2 );
LONG LONG INT f4k3 = f4k * ( ll4k + 1 ) * ( ll4k + 2 ) * ( ll4k + 3 );
LONG LONG REAL term := 0;
FOR c TO UPB w coefficients DO
term *:= llk +:= w coefficients[ c ]
OD;
LONG LONG INT fp = f2k1 * f2k * fk;
term *:= fp * fp * fp /:= f3k2 * f4k3 * f4k3 * f4k3;
IF ODD k THEN
zeta3 -:= term
ELSE
zeta3 +:= term
FI
OD;
zeta3 /:= 24;
print( ( truncate100( zeta3 ), newline ) )
END
END
</syntaxhighlight>
{{out}}
<pre>
1.2020564036593442854830714115115999903483212709031775135036540966118572571921400836130084123260473111
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
</pre>
=={{header|F_Sharp|F#}}==
Line 135 ⟶ 232:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
</pre>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight lang="Mathematica">
ClearAll["Global`*"];
TruancateTo100DecimalDigits = N[#, 100 + 1] &;
MyShowApéryConstant[expr_, caption_String] :=
Print[caption <>
ToString@Activate@TruancateTo100DecimalDigits[expr]];
MyShowApéryConstant[
Zeta[3], "Apéry's constant via Mathematica's Zeta:\n"]
MyShowApéryConstant[
Sum[1/(k^3), {k, 1,
1000}], "Apéry's constant via reciprocal cubes:\n"]
MyShowApéryConstant[(5/2*
Sum[(-1)^(k - 1)*(k!)^2/((2 k)!*k^3), {k, 1,
158}]), "Apéry's constant via Markov's summation:\n"]
MyShowApéryConstant[
1/24*Sum[(-1)^
k*((2 k + 1)!)^3*((2 k)!)^3*(k!)^3*(126392 k^5 + 412708 k^4 +
531578 k^3 + 336367 k^2 + 104000 k +
12463)/(((3 k + 2)!)*((4 k + 3)!)^3), {k, 0,
19}], "Apéry's constant via Wedeniwski's summation:\n"]
</syntaxhighlight>
{{out}}
<div style=background-color:#f8f9fa;padding:1em;white-space:pre-wrap;font-family:monospace;line-height:1.2em;>Apéry's constant via Mathematica's Zeta:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
Apéry's constant via reciprocal cubes (accurate to 6 decimal places):
1.202056<span style=color:red;>4036593442854830714115115999903483212709031775135036540966118572571921400836130084123260473112</span>
Apéry's constant via Markov's summation:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
Apéry's constant via Wedeniwski's summation:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
</div>
=={{header|Nim}}==
Line 196 ⟶ 330:
echo "Actual value to 100 decimal places:"
echo
echo()
apery(1000)
Line 215 ⟶ 349:
First 20 terms of Wedeniwski representation truncated to 100 decimal places:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
</pre>
=={{header|PARI/GP}}==
<syntaxhighlight lang="PARI/GP">
\\ Set the precision to, say, 110 digits
default(realprecision, 110);
\\ Function to display Apéry's constant
my_show_apery_constant(expr, caption) = printf("%s %.101g\n\n", caption, expr);
\\ Apéry's constant via PARI/GP's zeta function
my_show_apery_constant(zeta(3), "Apéry's constant via PARI/GP's Zeta:\n");
\\ Apéry's constant via reciprocal cubes
my_show_apery_constant(sum(k = 1, 1000, 1/k^3), "Apéry's constant via reciprocal cubes:\n");
\\ Apéry's constant via Markov's summation
my_show_apery_constant((5/2) * sum(k = 1, 158, (-1)^(k - 1) * (k!)^2 / ((2*k)! * k^3)), "Apéry's constant via Markov's summation:\n");
\\ Apéry's constant via Wedeniwski's summation
my_show_apery_constant(1/24 * sum(k = 0, 19, (-1)^k * ((2*k + 1)!)^3 * ((2*k)!)^3 * (k!)^3 * (126392*k^5 + 412708*k^4 + 531578*k^3 + 336367*k^2 + 104000*k + 12463) / (((3*k + 2)!) * ((4*k + 3)!)^3)), "Apéry's constant via Wedeniwski's summation:\n");
</syntaxhighlight>
{{out}}
<pre>
Apéry's constant via PARI/GP's Zeta:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
Apéry's constant via reciprocal cubes:
1.2020564036593442854830714115115999903483212709031775135036540966118572571921400836130084123260473112
Apéry's constant via Markov's summation:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
Apéry's constant via Wedeniwski's summation:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
</pre>
Line 259 ⟶ 430:
Ugh. If you ran this on the James Webb, you might just be able to pick out a faint small print outline of the word "elegant".<br>
Still, at least it is not like you do this sort of thing every day... and I got to fix a couple of bugs in my mpfr.js code.
<!--
<syntaxhighlight lang="phix">
with javascript_semantics
requires("1.0.2") -- (missing mpfr_ui_pow_ui() and bug in mpfr_mul_d(), both in mpfr.js)
include mpfr.e
mpfr_set_default_precision(-100)
mpfr {d,a,w,t} = mpfr_inits(4)
mpz {z,pk} = mpz_inits(2)
for k=1 to 1000 do
mpfr_ui_pow_ui(t,k,3)
mpfr_si_div(t,1,t)
mpfr_add(d,d,t)
end for
for k=1 to 158 do
mpz_fac_ui(z,k)
mpz_mul(z,z,z)
mpfr_set_z(t,z)
mpz_fac_ui(z,2*k)
mpfr_div_z(t,t,z)
mpz_ui_pow_ui(z,k,3)
mpfr_div_z(t,t,z)
if even(k) then
mpfr_sub(a,a,t)
else
mpfr_add(a,a,t)
end if
end for
mpfr_mul_d(a,a,5/2)
for k=0 to 19 do
mpz_ui_pow_ui(z,k,5)
mpz_mul_si(z,z,126392)
mpz_ui_pow_ui(pk,k,4)
mpz_mul_si(pk,pk,412708)
mpz_add(z,z,pk)
mpz_ui_pow_ui(pk,k,3)
mpz_mul_si(pk,pk,531578)
mpz_add(z,z,pk)
mpz_add_si(z,z,k*k*336367)
mpz_add_si(z,z,k*104000)
mpz_add_si(z,z,12463)
mpfr_set_z(t,z)
mpz_fac_ui(z,2*k+1)
mpz_pow_ui(z,z,3)
mpfr_mul_z(t,t,z)
mpz_fac_ui(z,2*k)
mpz_pow_ui(z,z,3)
mpfr_mul_z(t,t,z)
mpz_fac_ui(z,k)
mpz_pow_ui(z,z,3)
mpfr_mul_z(t,t,z)
mpz_fac_ui(z,3*k+2)
mpfr_div_z(t,t,z)
mpz_fac_ui(z,4*k+3)
mpz_pow_ui(z,z,3)
mpfr_div_z(t,t,z)
if odd(k) then
mpfr_sub(w,w,t)
else
mpfr_add(w,w,t)
end if
end for
mpfr_div_si(w,w,24)
constant fmt = """
Actual value to 100 decimal places:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
First 1000 terms of zeta(3) truncated to 100 decimal places. (accurate to 6 decimal places):
%s
First 158 terms of Markov / Apery representation truncated to 100 decimal places:
%s
First 20 terms of Wedeniwski representation truncated to 100 decimal places:
%s
"""
string direct = mpfr_get_fixed(d,100),
mapery = mpfr_get_fixed(a,100),
wdnski = mpfr_get_fixed(w,100)
printf(1,fmt,{direct,mapery,wdnski})
</syntaxhighlight>
{{out}}
<pre>
Line 358 ⟶ 530:
<small>Last digit of the 1000 terms line is 2 under pwa/p2js...</small><br>
As per Wren, you can verify or completely replace all this with mpfr_zeta_ui(w,3) <small>[on desktop/Phix only, not supported under pwa/p2js]</small>
=={{header|Python}}==
<syntaxhighlight lang="Python">
from sympy import zeta, factorial
from decimal import Decimal, getcontext
# Set the desired precision
getcontext().prec = 120
def my_sympy_format_to_decimal(sympy_result):
return Decimal(str(sympy_result.evalf(getcontext().prec)))
def print_apery_constant(description, value):
print(f"{description}:\n{str(value)[:102]}")
# Apéry's constant via SymPy's zeta function
zeta_3_str = str(zeta(3).evalf(getcontext().prec))
zeta_3_decimal = Decimal(zeta_3_str)
print_apery_constant("Apéry's constant via SymPy's zeta", zeta_3_decimal)
# Apéry's constant via Riemann summation of 1/(k cubed)
def apery_r(nterms=1_000):
total = sum(Decimal('1') / Decimal(k) ** 3 for k in range(1, nterms + 1))
return total
print_apery_constant("Apéry's constant via reciprocal cubes", apery_r())
# Apéry's constant via Markov's summation
def apery_m(nterms=158):
total = Decimal(2.5) * sum(
(Decimal(1) if k % 2 != 0 else Decimal(-1)) *
my_sympy_format_to_decimal(factorial(k) ** 2) /
my_sympy_format_to_decimal(factorial(2*k) * (k ** 3) )
for k in range(1, nterms + 1)
)
return total
print_apery_constant("Apéry's constant via Markov's summation", apery_m())
# Apéry's constant via Wedeniwski's summation
def apery_w(nterms=20):
total = Decimal('1') / Decimal('24') * sum(
(Decimal('1') if k % 2 == 0 else Decimal('-1')) *
my_sympy_format_to_decimal(factorial(2 * k + 1)) ** 3 *
my_sympy_format_to_decimal(factorial(2 * k)) ** 3 *
my_sympy_format_to_decimal(factorial(k)) ** 3 *
(Decimal('126392') * Decimal(k) ** 5 +
Decimal('412708') * Decimal(k) ** 4 +
Decimal('531578') * Decimal(k) ** 3 +
Decimal('336367') * Decimal(k) ** 2 +
Decimal('104000') * Decimal(k) +
Decimal('12463')) /
(my_sympy_format_to_decimal(factorial(3 * k + 2)) * my_sympy_format_to_decimal(factorial(4 * k + 3)) ** 3)
for k in range(0, nterms + 1)
)
return total
print_apery_constant("Apéry's constant via Wedeniwski's summation", apery_w())
</syntaxhighlight>
{{out}}
<pre>
Apéry's constant via SymPy's zeta:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
Apéry's constant via reciprocal cubes:
1.2020564036593442854830714115115999903483212709031775135036540966118572571921400836130084123260473111
Apéry's constant via Markov's summation:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
Apéry's constant via Wedeniwski's summation:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
</pre>
=={{header|Raku}}==
Line 389 ⟶ 628:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
</div>
=={{header|Sidef}}==
<syntaxhighlight lang="ruby">local Num!PREC = 4*101
say "Actual value to 100 decimal places:\n#{zeta(3)}"
say "\nFirst 1000 terms of ζ(3) truncated to 100 decimal places. (accurate to 6 decimal places):";
say sum(1..1000, {|k| 1/k**3 }).as_float(101)
say "\nFirst 158 terms of Markov / Apéry representation truncated to 100 decimal places:";
say ((5/2)*sum(1..158, {|k|
(-1)**(k-1) * (k!**2 / ((2*k)! * k**3))
}) -> as_float(101))
say "\nFirst 20 terms of Wedeniwski representation truncated to 100 decimal places:";
say ((1/24)*sum(^20, {|k|
(-1)**k * (2*k + 1)!**3 * (2*k)!**3 * k!**3 * (
126392*k**5 + 412708*k**4 + 531578*k**3 + 336367*k**2 + 104000*k + 12463
) / ((3*k + 2)! * (4*k + 3)!**3)
}) -> as_float(101))</syntaxhighlight>
{{out}}
<pre>
Actual value to 100 decimal places:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
First 1000 terms of ζ(3) truncated to 100 decimal places. (accurate to 6 decimal places):
1.2020564036593442854830714115115999903483212709031775135036540966118572571921400836130084123260473112
First 158 terms of Markov / Apéry representation truncated to 100 decimal places:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
First 20 terms of Wedeniwski representation truncated to 100 decimal places:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
</pre>
=={{header|Wren}}==
{{libheader|Wren-big}}
<syntaxhighlight lang="
var apery = Fn.new { |n|
Line 476 ⟶ 748:
We can also verify the actual value of Apéry's constant to 100 decimal places using MPFR which has a zeta function built in. A precision of 324 bits is needed.
{{libheader|Wren-gmp}}
<syntaxhighlight lang="
var x = Mpf.new(324)
| |||