Apéry's constant

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

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

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)

ALGOL 68

Works with: ALGOL 68G version 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.

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

Output:

1.2020564036593442854830714115115999903483212709031775135036540966118572571921400836130084123260473111
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581

F#

// Apéry's constant. Nigel Galloway: March 3rd., 2023
open MathNet.Numerics
let fact=let g=Seq.unfold(fun(n,g)->Some(n,(n*g,g+1N)))(1N,2N)|>Seq.cache in (fun n->Seq.item (n-1) g)
let fN g=let g=BigRational.FromInt g in 126392N*g**5+412708N*g**4+531578N*g**3+336367N*g**2+104000N*g+12463N
let fG n g l=let i=n/g in (int i,Seq.unfold(fun(n,i)->if i=0 then None else let l=n/g in Some(int l,(10I*(n-l*g),i-1)))(10I*(n-i*g),l))
let r3=Seq.initInfinite(fun g->BigRational.PowN(((+)1>>BigRational.FromInt>>BigRational.Reciprocal)g,3))|>Seq.take 1000|>Seq.sum
let ma=(5N/2N)*(Seq.unfold(fun(n,g,l)->Some(n*g,(-n,(fact l)*(fact l)/(fact(2*l)*BigRational.FromInt(pown l 3)),l+1)))(1N,1N/2N,2)|>Seq.take 158|>Seq.sum)
let sw=(1N/24N)*(Seq.unfold(fun(n,g,l)->Some(n*g,(-n,(fact(2*l+1)**3*fact(2*l)**3*(fact l)**3*(fN l))/(fact(3*l+2)*(fact(4*l+3)**3)),l+1)))(1N,12463N/432N,1)|>Seq.take 20|>Seq.sum)
[r3;ma;sw]|>List.iter(fun n->let n,g=fG (n.Numerator) (n.Denominator) 100 in printf $"%d{n}."; g|>Seq.iter(printf "%d"); printfn "")

Output:

1.2020564036593442854830714115115999903483212709031775135036540966118572571921400836130084123260473111
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581

Bonus 4th. way as a continued fraction

I extend Continued fraction#Apéry's_constant to provide the 101 digits required here. This method requires 31 iterations to provide the 101 digits, not quite as good as Wedeniwski but much better than the other 2. This method has historic interest as Apéry used it to prove that this number is irrational.

let cf2br α β=let n0,g1,n1,g2=β(),α(),β(),β()
              seq{let (Π:BigRational)=g1/n1 in yield n0+Π; yield! Seq.unfold(fun(n,g,Π)->let a,b=α(),β() in let Π=Π*g/n in Some(n0+Π,(b+a/n,b+a/g,Π)))(g2+α()/n1,g2,Π)}
let aπ()=let mutable n=0N in (fun ()->n<-n+1N; -(BigRational.Pow(n,6)))
let bπ()=let mutable n=0N in (fun ()->n<-n+1N; (2N*n-1N)*(17N*n*n-17N*n+5N))
cf2br (aπ()) (bπ())|>Seq.skip 31|>Seq.take 1|>Seq.iter(fun n->let n=6N/n in let n,g=fG (n.Numerator) (n.Denominator) 100 in printf $"%d{n}."; g|>Seq.iter(printf "%d"); printfn "")

Output:

1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581

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

Java

import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;
import java.math.RoundingMode;
import java.util.HashMap;
import java.util.Map;

public final class AperysConstant {

	public static void main(String[] args) {
		System.out.println("Apéry's constant = ζ(3), truncated to 100 digits");
		System.out.println();
		apéry(1_000);
		markov(158);
		wedeniwski(20);
		System.out.println("True value: 1.20205690315959428539973816151144999" +
			"07649862923404988817922715553418382057863130901864558736093352581");		
	}
	
	private static void apéry(int termCount) {
		BigRational sum = BigRational.ZERO;
		for ( int k = 1; k <= termCount; k++ ) {
			sum = sum.add( new BigRational(BigInteger.ONE, BigInteger.valueOf(k * k * k)) );
		}
		System.out.println("Sum cubes : " + sum.truncate(100));
	}
	
	private static void markov(int termCount) {
		final BigInteger terms = BigInteger.valueOf(termCount);
	    BigInteger sign = BigInteger.ONE.negate();
	    BigRational sum = BigRational.ZERO;
	    for ( BigInteger k = BigInteger.ONE; k.compareTo(terms) <= 0; k = k.add(BigInteger.ONE) ) {
	        sign = sign.negate();
	        BigInteger numerator = sign.multiply(factorial(k).pow(2));
	        BigInteger denominator = factorial(BigInteger.TWO.multiply(k)).multiply(k.pow(3));
	        sum = sum.add( new BigRational(numerator, denominator) );	        
	    }
	    sum = sum.multiply( new BigRational(BigInteger.valueOf(5), BigInteger.TWO) );
	    System.out.println("Markov    : " + sum.truncate(100));
	}
	
	private static void wedeniwski(int termCount) {
		final BigInteger terms = BigInteger.valueOf(termCount);
	    BigInteger sign = BigInteger.ONE.negate();
	    BigRational sum = BigRational.ZERO;
	    for ( BigInteger k = BigInteger.ZERO; k.compareTo(terms) < 0; k = k.add(BigInteger.ONE) ) {
	        sign = sign.negate();	
	        BigInteger numerator = sign.multiply(factorial(BigInteger.TWO.multiply(k).add(BigInteger.ONE))
	        	.multiply(factorial(BigInteger.TWO.multiply(k))).multiply(factorial(k)).pow(3));
	        BigInteger term = ((((BigInteger.valueOf(126392).multiply(k)
	        	.add(BigInteger.valueOf(412708))).multiply(k)
	        	.add(BigInteger.valueOf(531578))).multiply(k)
	        	.add(BigInteger.valueOf(336367))).multiply(k)
	        	.add(BigInteger.valueOf(104000))).multiply(k)
	        	.add(BigInteger.valueOf(12463));	        		
	        numerator = numerator.multiply(term);
	        BigInteger denominator = factorial(BigInteger.valueOf(3).multiply(k).add(BigInteger.TWO))
	        	.multiply(factorial(BigInteger.valueOf(4).multiply(k).add(BigInteger.valueOf(3))).pow(3));
	        sum = sum.add( new BigRational(numerator, denominator) );
	    }
	    sum = sum.multiply( new BigRational(BigInteger.ONE, BigInteger.valueOf(24)) );
	    System.out.println("Wedeniwski: " + sum.truncate(100));
	}
	
	private static BigInteger factorial(BigInteger n) {
		if ( ! cache.containsKey(n) ) {
			return n.multiply(factorial(n.subtract(BigInteger.ONE)));
		}
		return cache.get(n);
	}
	
	private static Map<BigInteger, BigInteger> cache =
		new HashMap<BigInteger, BigInteger>(Map.of( BigInteger.ZERO, BigInteger.ONE ));

}

final class BigRational {
	
	public BigRational(BigInteger aNumerator, BigInteger aDenominator) {
		numerator = aNumerator;
		denominator = aDenominator;
		
    	BigInteger gcd = numerator.gcd(denominator);	    	
    	numerator = numerator.divide(gcd);
    	denominator = denominator.divide(gcd);
    }
	
	public BigRational add(BigRational other) {
		BigInteger numer = numerator.multiply(other.denominator).add(denominator.multiply(other.numerator));
		BigInteger denom = denominator.multiply(other.denominator);		
		return new BigRational(numer, denom);
	}	

	public BigRational multiply(BigRational other) {
		return new BigRational(numerator.multiply(other.numerator), denominator.multiply(other.denominator));
	}
	
	public BigDecimal truncate(int decimalPlaces) {
		BigDecimal numer = new BigDecimal(numerator);
		BigDecimal denom = new BigDecimal(denominator);
		return numer.divide(denom, new MathContext(decimalPlaces + 1, RoundingMode.FLOOR));
	}	
	
	public static final BigRational ZERO = new BigRational(BigInteger.ZERO, BigInteger.ONE);
	
	private BigInteger numerator;
	private BigInteger denominator;
	
}

Output:

Apéry's constant = ζ(3), truncated to 100 digits

Sum cubes : 1.2020564036593442854830714115115999903483212709031775135036540966118572571921400836130084123260473111
Markov    : 1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
Wedeniwski: 1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581
True value: 1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581

jq

Adapted from Wren

Works with gojq, the Go implementation of jq

In this entry, gojq is used as it supports unbounded-precision integer arithmetic.

A copy of the "rational" module may be found at Category:jq/rational.jq.

include "rational" {search: "."};  # see comment above

def factorial:
  reduce range(2; .+1) as $i (1; . * $i);

def apery($n):
  reduce range (1; $n+1) as $k ( r(0;1);
      radd(.; r(1; $k | . * . * .)))
  | "First \($n) terms of ζ(3) truncated to 100 decimal places (accurate to 6 decimal places):",
     r_to_decimal(100);

def markov($n):
  {fact: 1,
   fact2: 1,
   sign: -1,
   sum: r(0;1)
  }
  | reduce range (1; $n + 1) as $k (.;
      .sign *= -1
      | .fact *= $k
      | (.fact * .fact * .sign) as $num
      | .mult = 2 * $k * (2*$k - 1)
      | .fact2 *= .mult
      | ($k | .*.*.) as $cube
      | (.fact2 * $cube) as $den
      | .sum |= radd(. ; r($num; $den) ) )
  | .sum |= rmult(.; r(5; 2))
  | "First \($n) terms of Markov / Apéry representation truncated to 100 decimal places:",
    (.sum | r_to_decimal(100));

def wedeniwski($n):
  def cube: .*.*.;
  {fact: 1,
   fact2: 1,
   sign: -1,
   sum: r(0;1),
   mult: 1
  }
  | reduce range (0; $n) as $k (.;
      .sign *= -1
      | if $k > 0
        then .fact *= $k
        | .mult = 2 * $k * (2*$k - 1)
        | .fact2 *= .mult
        end
      | (.fact2 * (2*$k + 1)) as $fact3
      | (.sign * ($fact3|cube) * (.fact2|cube) * (.fact|cube)) as $num
      | ($k * $k)  as $k2
      | ($k2 * $k) as $k3
      | ($k3 * $k) as $k4
      | ($k4 * $k) as $k5
      | (126392*$k5 + 412708*$k4 + 531578*$k3 + 336367*$k2 + 104000*$k + 12463) as $tmp
      | ($num * $tmp) as $num
      | ( (3*$k + 2 | factorial) * (4*$k + 3 | factorial | cube ) ) as $d
      | .sum |= radd(.; r($num; $d)) )
  | .sum |= r( .n; 24 * .d )
  | "First \($n) terms of Wedeniwski representation truncated to 100 decimal places:",
    (.sum | r_to_decimal(100));

"Actual value to 100 decimal places:",
"1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581",
apery(1000),
markov(158),
wedeniwski(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

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

Mathematica /Wolfram Language

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

Output:

Apéry's constant via Mathematica's Zeta:

1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581

Apéry's constant via reciprocal cubes (accurate to 6 decimal places):

1.2020564036593442854830714115115999903483212709031775135036540966118572571921400836130084123260473112

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

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

Nim

Translation of: Wren

Library: bignum

import std/strformat
import bignum

func toDecimal100(r: Rat): string =
  ## Return the representation of a rational up to 100 decimals.
  r *= newInt(10)^100
  result.setLen(102)
  result = ($r.toInt)[0..100]
  result.insert(".", 1)

proc apery(n: Positive) =
  var sum = newRat()
  for k in 1..n:
    sum += newRat(1, k^3)
  echo &"First {n} terms of ζ(3) truncated to 100 decimal places (accurate to 6 decimal places):"
  echo sum.toDecimal100
  echo()

proc markov(n: Positive) =
  var neg = true
  var fact1, fact2 = newInt(1)
  var sum = newRat()
  for k in 1..n:
    neg = not neg
    fact1 *= k
    var num = fact1 * fact1
    if neg: num = -num
    fact2 *= 2 * k * (2 * k - 1)
    let denom = fact2 * k^3
    sum += newRat(num, denom)
  sum *= newRat(5, 2)
  echo &"First {n} terms of Markov / Apéry representation truncated to 100 decimal places:"
  echo sum.toDecimal100
  echo()

proc wedeniwski(n: Positive) =
  var fact1, fact2 = newInt(1)
  var neg = true
  var sum = newRat()
  for k in 0..<n:
    neg = not neg
    if k > 0:
      fact1 *= k
      fact2 *= 2 * k * (2 * k - 1)
    let fact3 = fact2 * (2 * k + 1)
    var num = (fact1 * fact2 * fact3)^3
    num *= ((((126392 * k + 412708) * k + 531578) * k + 336367) * k + 104000) * k + 12463
    if neg: num = -num
    let denom = fac(4 * k + 3)^3 * fac(3 * k + 2)
    sum += newRat(num, denom)
  sum /= 24
  echo &"First {n} terms of Wedeniwski representation truncated to 100 decimal places:"
  echo sum.toDecimal100
  echo()


echo "Actual value to 100 decimal places:"
echo "1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581"
echo()
apery(1000)
markov(158)
wedeniwski(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

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");

Output:

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

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]

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

Output:

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

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

REXX

Libraries: How to use
Library: Numbers
Library: Functions

Version 1 Straightforward

Just apply the formulas.

arg n; if n = '' then n = 101; numeric digits n
parse version version; say version; glob. = ''
say 'Apery''s constant = zeta(3) to' n 'decimal places'
say 'Version 1 Three methods, using formulas'
say
call time('r'); a = Definition(); e = format(time('e'),,3)
say 'Definition' a '('e 'seconds)'
call time('r'); a = Markov(); e = format(time('e'),,3)
say 'Markov    ' a '('e 'seconds)'
call time('r'); a = Wedeniwski(); e = format(time('e'),,3)
say 'Wedeniwski' a '('e 'seconds)'
call time('r'); a = TrueValue(); e = format(time('e'),,3)
say 'True value' a '('e 'seconds)'
exit

Definition:
procedure
numeric digits digits()+2
y = 0
do k = 1 to 1000
   y = y + 1/(k**3)
end
numeric digits digits()-2
return y+0

Markov:
/* Markov Apery series */
procedure expose fact. glob.
fact. = 0
numeric digits digits()+2
y = 0
do k = 1 to 158
   y = y + (-1)**(k-1) * Fact(k)**2 / (Fact(2*k)*k**3)
end
y = y*2.5
numeric digits digits()-2
return y+0

Wedeniwski:
/* Wedeniwski series */
procedure expose fact. glob.
fact. = 0
numeric digits Digits()+2
y = 0
do k = 0 to 19
   y = y + ((-1)**k * Fact(2*k+1)**3 * Fact(2*k)**3 * Fact(k)**3,
         * (((((126392*k+412708)*k+531578)*k+336367)*k+104000)*k+12463)),
         / (Fact(3*k+2) * Fact(4*k+3)**3)
end
y = y/24
numeric digits Digits()-2
return y+0

TrueValue:
procedure
return 1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581+0

include Functions
include Numbers

Output:

You may specify the precision as parameter. Running with the default (101 digits = 100 decimal places) gives:

REXX-Regina_3.9.6(MT) 5.00 29 Apr 2024
Apery's constant = zeta(3) to 100 decimal places
Version 1 Formulas with factorials

Definition 1.2020564036593442854830714115115999903483212709031775135036540966118572571921400836130084123260473112 (0.020 seconds)
Markov     1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581 (0.071 seconds)
Wedeniwski 1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581 (0.007 seconds)
True value 1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581 (0.000 seconds)

Version 2 Wedeniwski optimized

All factorials replaced by recurring updates and a few other optimizations (compare ALGOL 68 and Nim). This version does not require additional routines. Much more verbose as Version 1, however slightly faster.

arg n; if n = '' then n = 101; numeric digits n
parse version version; say version; glob. = ''
say 'Apery''s constant = zeta(3) to' n 'decimal places'
say 'Version 2 Wedeniwski factorials replaced by recurring values'
say
call time('r'); a = Wedeniwski(); e = format(time('e'),,3)
say 'Wedeniwski' a '('e 'seconds)'
call time('r'); a = TrueValue(); e = format(time('e'),,3)
say 'True value' a '('e 'seconds)'
exit

Wedeniwski:
/* Wedeniwski */
procedure
numeric digits Digits()+2
/* Term for k = 0 */
f10 = 1; f20 = 1; f21 = 1; f32 = 2; f43 = 6; f00 = 12463; s = 1; y = 0
do k = 1 to 20
/* Add term to series */
   y = y + s * ((f10*f20*f21)**3*f00) / (f32*f43**3)
/* Prepare next term */
   k2 = k+k; k3 = k2+k; k4 = k3+k
/* Recurring factorials */
   f10 = f10 * k
   f20 = f20 * k2 * (k2-1)
   f21 = f21 * k2 * (k2+1)
   f32 = f32 * k3 * (k3+1) * (k3+2)
   f43 = f43 * k4 * (k4+1) * (k4+2) * (k4+3)
/* Polynomial */
   f00 = ((((126392*k+412708)*k+531578)*k+336367)*k+104000)*k+12463
/* Change sign */
   s = -s
end
y = y/24
numeric digits Digits()-2
return y+0

TrueValue:
procedure
return 1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581+0

Output:

Running with the default 101 digits, this program gives:

REXX-Regina_3.9.6(MT) 5.00 29 Apr 2024
Apery's constant = zeta(3) to 100 decimal places
Version 2 Factorials replaced by recurring values

Wedeniwski 1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581 (0.006 seconds)
True value 1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581 (0.000 seconds)

Sidef

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

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

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