Apéry's constant

Apéry's constant is the sum of the reciprocals of the positive cubes.

Apéry's constant is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

That is, it is defined as the number $\zeta (3)=\sum _{k=1}^{\infty }{\frac {1}{k^{3}}},$ where ζ is the Riemann zeta function.

Approximately equal to:

1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581

This constant was known and studied by many early mathematicians, but was not named until 1978, after Roger Apéry, who was first to prove it was an irrational number.

$\zeta (3)$ by summing reciprocal cubes is easy to calculate, but converges very slowly. The first 1000 terms are only accurate to 6 decimal places.

There have been many fast convergence representations developed / discovered that generate correct digits much more quickly.

One of the earliest, discovered in the late 1800s by A. Markov and later widely published by Apéry is:

$\zeta (3)={\frac {5}{2}}\sum _{k=1}^{\infty }(-1)^{k-1}{\frac {k!^{2}}{(2k)!k^{3}}}.$

Much better than direct calculation of $\zeta (3)$, but still only yielding about .63 correct digits per iteration.

Several even faster converging representions are available. The fastest known to date, yielding about 5.04 correct digits per term, is by Sebastian Wedeniwski.

$\zeta (3)={\frac {1}{24}}\sum _{k=0}^{\infty }(-1)^{k}{\frac {(2k+1)!^{3}(2k)!^{3}k!^{3}(126392k^{5}+412708k^{4}+531578k^{3}+336367k^{2}+104000k+12463)}{(3k+2)!(4k+3)!^{3}}}.$

Task

• Show the value of Apéry's constant calculated at least three different ways.

1. Show the value of at least the first 1000 terms of $\zeta (3)$ by direct summing of reciprocal cubes, truncated to 100 decimal digits.
2. Show the value of the first 158 terms of Markov / Apéry representation truncated to 100 decimal digits.
3. Show the value of the first 20 terms of Wedeniwski representation truncated to 100 decimal digits.

See also

• Wikipedia: Apéry's constant
• OEIS:A002117 - Decimal expansion of zeta(3)

J

app3=: {{ +/_3x^~1+i.y }}
app3m=: {{ _5r2 * +/ {{ (_1^y)*(2^~!y)%(!2*y)*y^3}} 1x+i.y }}
app3sm=: {{ 1r24* +/ {{
  (_1^y)*(3^~*/!y,0 1+/2*y)*(12463 104000 336367 531578 412708 126392 p. y)%(!2 3 p.y)*(!3 4 p.y)^3
 }} x:i.y
}}

   0j100 ": app3 1000
1.2020564036593442854830714115115999903483212709031775135036540966118572571921400836130084123260473112
   0j100 ": app3m 158
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
   0j100 ": app3sm 20
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581

Julia

using SpecialFunctions

setprecision(120, base=10)

println("Apéry's constant via Julia's zeta:\n$(string(zeta(big"3"))[1:102])")

""" zeta(3) via Riemann summation of 1/(k cubed) """
Apéry_r(nterms = 1_000_000) = sum(big"1" / k^big"3" for k in 1:nterms)

println("\nApéry's constant via reciprocal cubes:\n$(string(Apéry_r())[1:102])")

""" zeta(3) via Markov's summation """
function Apéry_m(nterms = 158)
   return big"2.5" * sum((isodd(k) ? 1 : -1) * factorial(big(k))^2 /
   (factorial(big"2" * k) * k^big"3") for k in 1:nterms)
end

println("\nApéry's constant via Markov's summation:\n$(string(Apéry_m())[1:102])")

""" zeta(3) via Wedeniwski's summation """
function Apéry_w(nterms = 20)
   return big"1"/24 * sum((iseven(k) ? 1 : -1) * factorial(big"2" * k + 1)^3 *
      factorial(big"2" * k)^3 * factorial(big(k))^3 * 
      (126392 * k^big"5" + 412708 * k^big"4" + 531578 * k^big"3" + 336367 * k^big"2"
         + big"104000" * k + 12463) / (factorial(big"3" * k + 2) * factorial(big"4" * k+3)^3)
         for k in 0:nterms)
end

println("\nApéry's constant via Wedeniwski's summation:\n$(string(Apéry_w())[1:102])")

Output:

Apéry's constant via Julia's zeta:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581

Apéry's constant via reciprocal cubes:
1.2020569031590942858997379115114499908483196256737488817922717053418382053696464235214344450378979367

Apéry's constant via Markov's summation:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581

Apéry's constant via Wedeniwski's summation:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581

Perl

use v5.36;
use bigrat try => 'GMP';

sub f { my $r = 1; $r *= $_ for 1..shift; $r }

say 'Actual value to 100 decimal places:';
say '1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581';

say "\nFirst 1000 terms of ζ(3) truncated to 100 decimal places. (accurate to 6 decimal places):";
my $z3;
$z3 += 1/$_**3 for 1..1000;
say $z3->as_float(101);

say "\nFirst 158 terms of Markov / Apéry representation truncated to 100 decimal places:";
$z3 = 0;
$z3 += (-1)**($_-1) * (f($_)**2 / (f(2*$_) * $_**3)) for 1..158;
$z3 *= 5/2;
say $z3->as_float(101);

say "\nFirst 20 terms of Wedeniwski representation truncated to 100 decimal places:";
$z3 = 0;
$z3 += (-1)**$_ *  f(2*$_+1)**3 * f(2*$_)**3 * f($_)**3 * (126392*$_**5 + 412708*$_**4 + 531578*$_**3 + 336367*$_**2 + 104000*$_ + 12463)
                 / ( f(3*$_+2) * f(4*$_+3)**3 )
       for 0..19;
$z3 *= 1/24;
say $z3->as_float(101);

Output:

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

Phix

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".
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.

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})

Output:

Actual value to 100 decimal places:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581

First 1000 terms of zeta(3) truncated to 100 decimal places. (accurate to 6 decimal places):
1.2020564036593442854830714115115999903483212709031775135036540966118572571921400836130084123260473111

First 158 terms of Markov / Apery representation truncated to 100 decimal places:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581

First 20 terms of Wedeniwski representation truncated to 100 decimal places:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581

Last digit of the 1000 terms line is 2 under pwa/p2js...
As per Wren, you can verify or completely replace all this with mpfr_zeta_ui(w,3) [on desktop/Phix only, not supported under pwa/p2js]

Raku

sub postfix:<!> (Int $n) { (constant f = 1, |[\×] 1..*)[$n] }

say 'Actual value to 100 decimal places:';
say '1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581';

say "\nFirst 1000 terms of ζ(3) truncated to 100 decimal places. (accurate to 6 decimal places):";
say (1..1000).map({FatRat.new: 1, .³}).sum.substr: 0, 102;

say "\nFirst 158 terms of Markov / Apéry representation truncated to 100 decimal places:";
say (5/2 × (1..158).map( -> \k { (-1)**(k-1) × FatRat.new: k!², ((2×k)! × k³) } ).sum).substr: 0, 102;

say "\nFirst 20 terms of Wedeniwski representation truncated to 100 decimal places:";
say (1/24 × ((^20).map: -> \k {
    (-1)**k × FatRat.new: (2×k+1)!³ × (2×k)!³ × k!³ × (126392×k⁵ + 412708×k⁴ + 531578×k³ + 336367×k² + 104000×k + 12463), (3×k+2)! × (4×k+3)!³
}).sum).substr: 0, 102;

Output:

Actual value to 100 decimal places:

1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581

First 1000 terms of ζ(3) truncated to 100 decimal places. (accurate to 6 decimal places): 1.2020564036593442854830714115115999903483212709031775135036540966118572571921400836130084123260473111

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

Wren

Library: Wren-big

import "./big" for BigInt, BigRat

var apery = Fn.new { |n|
    var sum = BigRat.zero
    for (k in 1..n) sum = sum + BigRat.new(1, k*k*k)
    System.print("First %(n) terms of ζ(3) truncated to 100 decimal places (accurate to 6 decimal places):")
    System.print(sum.toDecimal(100, false))
    System.print()
}

var markov = Fn.new { |n|
    var fact = BigInt.one
    var fact2 = BigInt.one
    var sign = BigInt.minusOne
    var sum = BigRat.zero
    for (k in 1..n) {
        sign = sign * BigInt.minusOne
        fact = fact * k
        var num = fact.square * sign
        var mult = 2 * k * (2*k - 1)
        fact2 = fact2 * mult
        var cube = k * k * k
        var den = fact2 * cube
        sum = sum + BigRat.new(num, den)
    }
    sum = sum * BigRat.new(5, 2)
    System.print("First %(n) terms of Markov / Apéry representation truncated to 100 decimal places:")
    System.print(sum.toDecimal(100, false))
    System.print()
}

var wedeniwski = Fn.new { |n|
    var fact = BigInt.one
    var fact2 = BigInt.one
    var sign = BigInt.minusOne
    var sum = BigRat.zero
    var mult = 1
    for (k in 0..n-1) {
        sign = sign * BigInt.minusOne
        if (k > 0) {
            fact = fact * k
            mult = 2 * k * (2*k - 1)
            fact2 = fact2 * mult
        }
        var fact3 = fact2 * (2*k + 1)
        var num = sign * fact3.cube * fact2.cube * fact.cube
        var cube = k * k * k
        var quad = cube * k
        var pent = quad * k
        var tmp =  126392*pent + 412708*quad + 531578*cube + 336367*k*k + 104000*k + 12463
        num = num * tmp
        var den = BigInt.factorial(3*k + 2) * BigInt.factorial(4*k + 3).cube
        sum = sum + BigRat.new(num, den)
    }
    sum = sum / 24
    System.print("First %(n) terms of Wedeniwski representation truncated to 100 decimal places:")
    System.print(sum.toDecimal(100, false))
    System.print()
}

System.print("Actual value to 100 decimal places:")
System.print("1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581")
System.print()
apery.call(1000)
markov.call(158)
wedeniwski.call(20)

Output:

Actual value to 100 decimal places:
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581

First 1000 terms of ζ(3) truncated to 100 decimal places (accurate to 6 decimal places):
1.2020564036593442854830714115115999903483212709031775135036540966118572571921400836130084123260473111

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

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.

Library: Wren-gmp

import "./gmp" for Mpf

var x = Mpf.new(324)
var zeta = x.zetaUi(3)
System.print(zeta.toString(101))

Output:

1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581