This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Numerics Volume, Chapter 2

Find integers (Subscript[a, k])∈{1,2,3,4}, such that the following identity holds within 100 digits of precision.

(Subscript[a, 0])+(1/(((Subscript[a, 1])/π)+(1/(((Subscript[a, 2])/π)+(1/(((Subscript[a, 3])/π)+(1/…)))))))=(Subscript[a, 0])+((Subscript[a, 1])/π)+((Subscript[a, 2])/(π^2))+((Subscript[a, 3])/(π^3))+…

nextBestApproximation[cA_, 𝒸_, n_, o_] :=
Module[{λ = Length[cA](*, minδ, prec, cASeriesPart, 𝒸NList, cfF, 
        newTerms, L, aux, bS*)},
 (* difference between continued fraction and power series *)
 minδ = Infinity; 
 (* precision needed for a reliable calculation *)
 prec = Ceiling[Abs[Log[10, N[𝒸^(λ + n)]]]] + 20; 𝒸N = N[𝒸, prec]; 
 (* existing power series part--avoid multiple evaluation *)
 cASeriesPart = cA.(𝒸N^Range[0, λ - 1]); 
 (* powers of the constant *)
 𝒸NList = 𝒸N^Range[Length[cA], λ + n - 1];
 (* factors of the continued fraction *)
 cfF = Table[If[j == 0, 1, 𝒸N], {j, 0, λ + n - 1}];
 (* new series terms *)
 newTerms = Table[t[k], {k, n}];
 (* find best possible of the o^n possible continuations *)
 Do[L = Join[cA, newTerms];
    (* difference between continued fraction and power series *)
    aux = Abs[FromContinuedFraction[cfF L] - cASeriesPart - newTerms.𝒸NList];
    If[aux < minδ, minδ = aux; bS = L],
    (* generate iterators *)
    Evaluate[Sequence @@ ({#, o}& /@ newTerms)]]; 
 (* return best sequence continuation and difference *)
   {bS, N @ minδ}]

If we do not take enough “look ahead” terms, the process of finding the next series terms gets stuck.

NestList[nextBestApproximation[Drop[#[[1]], 
             If[Length[#[[1]]] < 2, 0, -2]], 1/Pi, 3, 4]&, 
         {{}, Infinity}, 20]

Next, we use six “look ahead” terms and discard five of them after each call to nextBestApproximation. This time, we do not get stuck so early.

NestList[nextBestApproximation[Drop[#[[1]], 
             If[Length[#[[1]]] < 5, 0, -5]], 1/Pi, 6, 4]&, 
        {{}, Infinity}, 20]

Running now the following input for several hours, we get the sequence we are looking for.

NestList[((* see progress *) Print[#]; 
          nextBestApproximation[Drop[#[[1]], 
                  Max[0, Length[#[[1]] - 7]]], 1/Pi, 8, 4])&, 
         {{}, Infinity}, 12]

πCFPSSequence = {1, 2, 1, 1, 4, 1, 3, 1, 3, 3, 3, 4, 2, 3, 3, 2, 3, 4, 3, 2, 3, 
  4, 4, 3, 2, 3, 3, 1, 2, 3, 4, 2, 4, 2, 1, 4, 1, 2, 1, 2, 2, 2, 1, 4, 1, 2, 3, 
  3, 1, 3, 3, 3, 4, 4, 1, 2, 4, 3, 4, 3, 4, 4, 3, 4, 4, 4, 1, 4, 1, 3, 4, 4, 1, 
  4, 4, 4, 3, 4, 1, 2, 3, 1, 2, 1, 3, 3, 4, 1, 2, 3, 4, 3, 3, 3, 3, 2, 2, 4, 3, 
  4, 2, 4, 3, 4, 4, 1, 3, 4, 4, 2, 1, 4, 4, 3, 4, 1, 1, 2, 2, 3, 3, 1, 1, 2, 4, 
  1, 4, 1, 3, 1, 1, 3, 3, 4, 2, 3, 4, 3, 3, 2, 1, 1, 4, 1, 1, 4, 1, 3, 1, 2, 3, 
  1, 2, 1, 4, 4, 2, 4, 3, 2, 4, 1, 2, 3, 3, 2, 1, 3, 3, 2, 3, 1, 3, 2, 3, 4, 2, 
  3, 1, 3, 3, 2, 1, 1, 1, 4, 3, 3, 4, 4, 1, 4, 1, 3, 4, 4, 1, 4, 2, 1, 4, 3, 1, 
  2, 1, 2, 4, 2, 4, 4, 2, 1, 2, 2, 1, 2, 1, 3, 4, 1, 3, 2, 2, 3, 4, 1, 3, 3, 4, 
  2, 3, 1, 2, 2, 3, 3, 1, 4, 1, 4, 1, 1, 3, 1, 4, 2, 4, 3, 1, 1, 1, 4, 2, 4, 1, 
  4, 4, 3, 3, 2, 3, 3, 2, 3, 3, 1, 1, 1, 4, 1, 1, 2, 2, 2, 4, 2, 1, 1, 2, 2, 4,
  1, 3, 4, 4, 1, 2, 1, 1, 4, 1, 2, 2, 2, 3, 3, 3, 2, 4, 3, 1, 4, 3, 3, 2, 4, 1, 
  2, 3, 3, 4, 4, 3, 1, 1, 1, 3, 3, 2, 3, 1, 2, 4, 3, 4, 3, 2, 2, 1, 1, 2, 1, 3, 
  1, 3, 3, 3, 2, 4, 2, 2, 2, 3, 4, 3, 3, 2, 2};

The sequence has about 350 terms.

Length[πCFPSSequence]

The difference between the continued fraction and the power series is less than (10^(-100)).

(* value of the continued fraction *)
cfL = MapIndexed[If[#2[[1]] === 1, 1, 1/Pi]#1&, πCFPSSequence];

(* value of the power series *)
cfS = MapIndexed[#1 (1/Pi)^(#2[[1]] - 1)&, πCFPSSequence];

(* difference *)
(FromContinuedFraction[N[cfL, 300]] - Total[cfS]) // N

FromContinuedFraction[N[cfL, 300]] // N[#, 30]&