In this notebook the efficient ways to compute Champernowne constant ζ and its continued fraction expansion are discussed.

ζ is defined as the decimal number having after the decimal point the digits of consecutive natural numbers beginning with 1: ζ = 0.12345678910111213...(ie we concatenate positive integers).

Firstly a fast way to compute ζ is needed. This can be done arithmetically: \(TraditionalForm
      \`ζ = 
        \(1\/10 + 2\/10\^2 + ...\) + 9\/10\^9 + 10\/10\^11 + 11\/10\^13 + ...
          \) or, after a transformation,  \(TraditionalForm
      \`ζ = 
        \(1\/10\) 
          \((1 + \(1\/10\) 
                \((\(\(2 + \)...\) \(1\/10\) 
                    \((9 + \(1\/100\) 
                          \((10 + \(1\/100\) \((\(11 + \)...)\))\))\))\))
            \)\).

Here both ways are implemented.

\(champ1[n_] := 
      Module[
		{i, k = 1, p = 1, res = 0}, 
		
        For[i = 1, i <= n, \(i++\), 
			
          If[Mod[i, 10^p] == 0, 
				\(p++\)]; 
			k *= 10^p; 

          			res += i/k]; 
		res]\)

\(champ2[n_] := 
      Module[
		{i, k, res = 0}, 
		k = 10^\((Floor[Log[10, n]] + 1)\); 
        
		For[i = n, i >= 1, \(i--\), 
			res = \((res + i)\)/k; 

          			If[k == 10  i, 
				k = i]]; 
		res]\)

champ1[13]

champ2[13] == %

Timing[champ1[500]; ]Timing[champ2[500]; ]

The second way is evidently much faster, mainly because on each step it deals with numbers of the same order as n, while champ1[n] uses numbers at least as large as (10^n). Another approach to computing ζ is to try to avoid arithmetic operations and build it from list of digits using function IntegerDigits (or using strings and ToExpression).

\(champ3[n_] := 
	
      \(ToExpression[#]/10^StringLength[#]&\)[
		
        StringJoin[ToString\ /@\ Range[n]\ ]\ ]\)

\(champ4[n_] := 
	
      FromDigits[
		{Flatten[IntegerDigits\ /@\ Range[n]], 0}]\)

champ2[13] == champ3[13] == champ4[13]

Timing[champ2[3000]; ]Timing[champ3[3000]; ]Timing[champ4[3000]; ]

It turns out that so far the arithmetical way is the fastest.

Now, let's observe some interesting properties of ζ. Note that champ[(10^n)] gives ζ with (Subscript[L, n]) = 9 + 90*2 + ... + 9*(10^(n - 1))*n correct digits after the decimal point (because there are 9 one-digit numbers, 90 two-digit numbers, and so on).

L /: (Subscript[L, n])_ = 9 (∑_(k = 1)^n) k (10^(k - 1))

ch = N[champ2[100], (Subscript[L, 2])]

9  ch

The first digits are 1.111..., so 9*ζ is approximately (10/9) (ζ denoting the exact value of Champernowne constant).

(10/9) - 9  ch

If we multiply the result again, by factor of 99 (or 11) this time, another pattern appears.

11 ((10/9) - 9  ch)

The result can be written as \(TraditionalForm
      \`\((10 + 0.121212...)\)*10\^\(-8\) ≈ 
        \((10 + 12\/99)\)*10\^\(-8\)\). To obtain the next digits more precision is required.

Timing[ch = N[champ2[(10^4)], (Subscript[L, 4])]; ]

Subtract and multiply by 111.

N[111 (11 ((10/9) - 9  ch) - (10^(-8)) (10 + (12/99))), 
      100]

Use Rationalize to find the approximation for the repeating part.

Rationalize[0.011102011102011102011102]

Repeat the same procedure again.

N[1111 
        (111 (11 ((10/9) - 9  ch) - (10^(-8)) (10 + (12/99)))
                 - (10^(-187)) (100 + (122/10989))), 100]

Rationalize[0.012111121022012111121022]

11111 
      (1111 (
              111 (11 ((10/9) - 9  ch) - 
                      (10^(-8)) (10 + (12/99))) - 
                (10^(-187)) (100 + (122/10989))) - 
          (10^(-2886)) (11000 + (1222/100899)))

After each iteration the number of digits that get cancelled out is close to (Subscript[L, n]), ie each successive approximation is correct to (( ≈ 10)^(-(Subscript[L, n]))). So in order to see the beginning of the next 'period' ζ should be calculated using (10^n) plus several more natural numbers.

The following iterations will take a considerable amount of time but you may skip the evaluation of those cells since all the output is given here.

For number of digits that large, it becomes inefficient to perform arithmetic operations, and now champ3 and champ4 work noticeably faster than champ2 (and champ3 slightly faster than champ4, probably because strings are more 'compact' than lists).

Timing[ch = champ3[(10^5) + 100]; ]

Convert the vulgar fraction ch to decimal.

$MaxPrecision = (10^6);

Timing[chf = N[ch, (Subscript[L, 5]) + 100*6]; ]

And find the next repeating parts.

Timing[
      N[11111 (
            1111 (111 
                      (11 ((10/9) - 9  chf) - 
                          (10^(-8)) (10 + (12/99))) - 
                    (10^(-187)) (100 + (122/10989))) - 
              (10^(-2886)) (11000 + (1222/100899))), 100]]

Rationalize[0.13432233332234312442]

Timing[
      111111 (
          11111 (
                1111 (
                      111 (11 ((10/9) - 9  chf) - 
                              (10^(-8)) (10 + (12/99))) - 
                        (10^(-187)) (100 + (122/10989))) - 
                  (10^(-2886)) (11000 + (1222/100899))) - 
            (10^(-38885)) (12210000 + (150738/1122211)))]

Rationalize[14.921978010988901097802186813076]

And so the last acquired term (Subscript[a, n]) can be written as \(TraditionalForm
      \`\(10\^\(-488884\)\) \((135653100000 + 135788642\/9099909)\)\). Now the general form of the sequence can be guessed. The negative exponents of 10:

{8, 187, 2886, 38885, 488884} == {(Subscript[L, 1]) - 1, (Subscript[L, 2]) - 2, (Subscript[L, 3]) - 3, 
        (Subscript[L, 4]) - 4, (Subscript[L, 5]) - 5}

The integers in parenthesis:

{100, 11000, 12210000, 135653100000} == {1*(10^2), 11*(10^3), 
        11*111*(10^4), 11*111*1111*(10^5)}

The fractions:

{(12/99), (122/10989), (1222/100899), (150738/1122211), (135788642/9099909)}
         == (1/9) {12  (1/(1*11)), 122  (1/(11*111)), 
          1222  (11/(111*1111)), 12222 (( 11*111)/(1111*11111)), 
          122222 (( 11*111*1111)/(11111*111111))}

Let \(TraditionalForm\`r\_k = \(\(111...\) 1\)\&_\), \(TraditionalForm\`s\_k = \(\(1222...\) 2\)\&_\), where overbar means the number denoted by such digits, and (Subscript[r, k]), (Subscript[s, k]) each has (k + 1) digits. 

Let also (Subscript[R, n]) = (Subscript[r, 1])*(Subscript[r, 2])*...*(Subscript[r, n]) (so, for instance, \(TraditionalForm
      \`122222*\(11*111*1111\)\/\(11111*111111\) = 
        \(s\_5\) R\_3\/\(\(r\_4\) r\_5\)\)). 

Then for each term (Subscript[a, n]) found above \(TraditionalForm
      \`a\_n = \(10\^\(\(-L\_n\) + n\)\) 
          \((\(10\^n\) R\_\(n - 2\) + 
              \(s\_n\) \(R\_\(n - 2\)\/\(\(r\_\(n - 1\)\) r\_n\)\) 1\/9)
            \)\).

\(r /: r\_n_ = ∑\+\(k = 0\)\%n 10\^k; 

    R /: R\_n_ = ∏\+\(k = 1\)\%n r\_k; 

    s /: s\_n_ = 10\^n + 2 \(∑\+\(k = 0\)\%\(n - 1\)10\^k\); 

    a /: a\_n_ = 
      \(10\^\(\(-L\_n\) + n\)\) 
        \((\(10\^n\) R\_\(n - 2\) + 
            \(s\_n\) \(R\_\(n - 2\)\/\(\(r\_\(n - 1\)\) r\_n\)\) 1\/9)\); \)

\({\(10\^\(-8\)\) \((10 + 12\/99)\), 
        \(10\^\(-187\)\) \((100 + 122\/10989)\), 
        \(10\^\(-2886\)\) \((11000 + 1222\/100899)\), 
        \(10\^\(-38885\)\) \((12210000 + 150738\/1122211)\)} == 
	
      Table[a\_n, {n, 4}]\)

Now we can find ζ explicitly. Let's solve for ζ' = (10/9) - 9  ζ.
We started with ζ' and on each step multiplied the current value by (Subscript[r, n]) (by 11, then by 111, etc) and subtracted (Subscript[a, n]).
Thus, on the first step
\(TraditionalForm
      \`11  ζ' - a\_1 ≈ 0 => 
        ζ' ≈ a\_1\/11\);
on the next step
\(TraditionalForm
      \`111 \((11  ζ' - a\_1)\) - a\_2 ≈ 0 => 
          ζ' ≈ \(a\_2/111 + a\_1\)\/11 = 
        a\_2\/\(11*111\) + a\_1\/11\);
then
\(TraditionalForm
      \`1111 \((111 \((11  ζ' - a\_1)\) - a\_2)\) - a\_3 ≈ 0
             => 
          ζ' ≈ \(\((a\_3/1111 + a\_2)\)/111 + a\_1\)\/11 = 
        a\_3\/\(11*111*1111\) + a\_2\/\(11*111\) + a\_1\/11\).

nth iteration gives \(TraditionalForm
      \`\(ζ' ≈ a\_1\/R\_1 + a\_2\/R\_2 + ...\) + 
          a\_n\/R\_n = ∑\+\(k = 1\)\%n μ\_k\) => \(TraditionalForm
      \`ζ = 10\/81 - \(1\/9\) \(∑\_\(k = 1\)\%n μ\_k\)\). 

This can be simplified because \(TraditionalForm
      \`R\_\(n - 2\)\/R\_n = 1\/\(\(r\_\(n - 1\)\) r\_n\)\), and  
      ((Subscript[a, n])/(Subscript[R, n])) = 
        (10^(-(Subscript[L, n]) + n)) 
          (((10^n)/((Subscript[r, n - 1]) (Subscript[r, n]))) + 
              ((Subscript[s, n])/(9 ((((Subscript[r, n - 1]) (Subscript[r, n])))^2)))).

\(μ /: μ\_n_ = 
      \(10\^\(\(-L\_n\) + n\)\) 
        \((10\^n\/\(\(r\_\(n - 1\)\) r\_n\) + 
            s\_n\/\(9 \((\(r\_\(n - 1\)\) r\_n)\)\^2\))\); 

    ζ /: ζ\_n_ := 
      10\/81 - \(1\/9\) \(∑\+\(k = 1\)\%n μ\_k\); \)

Of course, we didn't prove that this sum converges to ζ, but if our suppositions are correct then this is an explicit analytical expression for Champernowne constant:

Evaluate //@ (Subscript[ζ, ∞])

N[(Subscript[ζ, 1]), (Subscript[L, 2])]

As we saw, the error of nth iteration is approximately (10^(-(Subscript[L, n]))) (so (Subscript[ζ, 1]) gives almost (Subscript[L, 2]) correct digits of ζ, etc). Actually, the error appears in the next to last group of digits (99 instead of 98 above, 999 instead of 998 in (Subscript[ζ, 2]) and so on). So the absolute value of the remainder of the series \(TraditionalForm
      \`\(|\(∑\+\(k = n\)\%∞ μ\_k 
          | \)\) = 
        9 | 
          \(ζ - ζ\_\(n - 1\) | 
            \( < 9*10\^\(\(-L\_n\) + n + 1\)\)\)\). This allows us to compute ζ very efficiently.

(Subscript[ζ, 4]) gives the same (Subscript[L, 5]) = 488889 digits as computing champ3[100000] above but about 450 times faster.

Timing[(Subscript[ζ, 4]); ]

You need not evaluate the next cell, zeta5 is used here and below only to compare timings.

Timing[zeta5 = (Subscript[ζ, 5]); ]

So adding more terms still requires significant time. The main reason is that the fractions appearing in the computation are huge and most of the time is spent trying to reduce them. It can be avoided by storing numerator and denominator of each fraction separately and performing all operations on them without using Mathematica's Rational type.

Let's represent each such fraction (m/n) as Rat[m, n] and define the basic operations for them.

The arithmetic operations involving integers, rationals and Rat numbers:

Rat /: r_Rat + n_Integer = r + Rat[n, 1];
Rat /: r_Rat n_Integer = r Rat[n, 1];

Rat/:r_Rat + q_Rational:=r+Rat[Numerator[q],Denominator[q]];
Rat /: r_Rat q_Rational := r Rat[Numerator[q], Denominator[q]];

Rat /: Rat[a_, b_] + Rat[c_, d_] = Rat[a*d + b*c, b*d];
Rat /: Rat[a_, b_] Rat[c_, d_] = Rat[a*c, b*d];

Rat /: Rat[a_, b_]^n_Integer?NonNegative = Rat[a^n, b^n];
Rat /: Rat[a_, b_]^n_Integer?Negative = Rat[b^-n, a^-n];

Other functions that will be needed later:

\(Rat /: Floor[Rat[a_, b_]] = Quotient[a, b]; 

    Rat /: FractionalPart[Rat[a_, b_]] = Rat[Mod[a, b], b]; 


    Rat /: Equal[Rat[a_, b_], n_Integer] = a == b*n; 

    Rat /: Equal[n_Integer, Rat[a_, b_]] = a == b*n; 

    Rat /: IntegerQ[r_Rat] = r == Floor[r]; 


    Rat /: Normal[r_Rat] := Divide\ @@\ r; 

    Rat /: N[Rat[a_, b_], n_:  $MachinePrecision] := 
      N[a, n + 1]/N[b, n + 1]\)

In fact, the N function for Rat won't be used, since, as was already shown, it takes a long time to compute the required decimal digits, and moreover, using them to find the continued fractions would be much slower than using Rat numbers (the point of using which is that all operations are integer).

Here is the formula for (Subscript[ζ, n]) using Rat ((Subscript[μ, n]) is changed to (Subscript[η, n]) and (Subscript[ζ, n]) to (Subscript[ξ, n])).

\(η /: η\_n_ = 
      \(Rat[Numerator[#], Denominator[#]]&\)\ @
		
        Simplify[Together[μ\_n]]\)

ξ /: (Subscript[ξ, n])_ = 
      (10/81) - (1/9) (∑_(k = 1)^n) (Subscript[η, k]);

(Subscript[ξ, 1])

N[(Subscript[ξ, 1]), (Subscript[L, 2])]

Normal[(Subscript[ξ, 1])] == (Subscript[ζ, 1])

The computation of (Subscript[ξ, 5]) is about 45 times faster than (Subscript[ζ, 5]).

Timing[xi5 = (Subscript[ξ, 5]); ]

Each of the three functions below creates a list of the first n partial quotients of the number x (if the expansion doesn't terminate earlier). FixedPointList in cf is called with three positional parameters, the third one defining the maximal number of steps.

\(cf[x_, n_:  Infinity] := 
	
      Floor[
		
        FixedPointList[
			1/FractionalPart[#]&, 
			x, n - 1, 

          			SameTest -> \((IntegerQ[#2]&)\)]\ ]\)

\(\(cf2[x_, n_:  Infinity] := 
      Module[
		{i, a = x, b, res = {}}, 
		
        For[i = 1, i <= n, \(i++\), 
			res = {res, b = Floor[a]}; 

          			If[a == b, 
				Break[]\ ]; 
			a = 1/\((a - b)\)]; 
        
		Flatten[res]\ ]; \)\)

\(\(cf3[x_Rat, n_:  Infinity] := 
      Module[
		{i, a1, a2, b, res = {}}, 
		{a1, a2} = List\ @@\ x; 

        		For[i = 1, i <= n, \(i++\), 
			
          If[a2 == 0, 
				Break[]\ ]; 
			
          res = {res, b = Quotient[a1, a2]}; 

          			{a1, a2} = {a2, a1 - a2*b}\ ]; 
		Flatten[res]\ ]; \)\)

You need not evaluate the two following cells.

Timing[cf[zeta5]; ]Timing[cf2[zeta5]; ]

Timing[cf[xi5]; ]Timing[cf2[xi5]; ]

The timings above show that cf is the slowest version (it uses FractionalPart and IntegerQ on each step, and in the case of Rat argument it actually does double work). Also, cf stores Rat objects before taking Floor. Notice that owing to the fact that the necessary arithmetical operations for Rat are properly defined, cf and cf2 can work with Rat as well as with internal number types.

Timing[dg = cf3[xi5]; ]

cf3 (specifically written for Rat) turns out to be not much slower than cf2 with rational argument (even though working with Rat involves much larger integers). Thus the difference between computing the initial data (zeta5 versus xi5) is dominating.

This is the table of the large (>1000) partial quotients, showing their positions in the expansion and their decimal logarithms.

\(Transpose[

        	{Flatten[Position[dg, x_ /; Log[10., x] > 3]], 
		
          Cases[dg, x_ /; Log[10., x] > 3 -> Floor[Log[10., x]]\ ]}\ ] // 
	
      TableForm\)

This shows where expansions for (Subscript[ξ, n]) terminate (in fact, the expansion for ζ is infinite).

\(Append[Array[Length[cf3[ξ\_#]]&, 4], Length[dg]]\)

Interestingly, the expansions for (Subscript[ξ, n]) terminate just before next large partial quotient (positions 19, 41, 163, 527). In fact, the sequence of partial quotients of ζ doesn't have upper bound. We can make a guess that there will be a number having approximately 4888000 digits in the position 1711.

Another question is the minimal precision required to obtain the desired number of digits in partial quotients (ie the maximal error that can be made in the initial approximation to ζ). This can be estimated by changing the input slightly and comparing the results. It turns out that to get n total digits in the expansion we need to know approximately 2*n decimal digits (this rule doesn't work for (Subscript[ζ, n]) which is a 'special case').

For the first 160 partial quotients to be correct, 5787 decimal digits are required.

cf3[(Subscript[ξ, 3]) + (10^(-5787)), 160] == cf3[(Subscript[ξ, 3]), 160]cf3[(Subscript[ξ, 3]) + (10^(-5788)), 160] == cf3[(Subscript[ξ, 3]), 160]

(∑_k) lg((Subscript[a, k]))is close to the total number of digits in the expansion with partial quotients (Subscript[a, k]). Rest drops the leading zero of the expansion.

Plus @@ (Log[10., Rest[cf3[(Subscript[ξ, 3]), 160]]])5787/%

Copyright © 1999, Maxim Rytin