Pi Continued Fraction

Author

Eric W. Weisstein
January 4, 2019

This notebook downloaded from http://mathworld.wolfram.com/notebooks/Constants/PiContinuedFraction.nb.

For more information, see Eric's MathWorld entry http://mathworld.wolfram.com/PiContinuedFraction.html.

©2019 Wolfram Research, Inc. except for portions noted otherwise

Plots

cf=ArrayPlot[Transpose@PadLeft@IntegerDigits[ContinuedFraction[Pi,256],2],PixelConstrained->3,Frame->False]

SetDirectory["/Users/eww/Sites/mathworld.wolfram.com/htdocs/images/gifs"]

Export["PiContinuedFractionBinaryPlot.gif",cf]

Continued Fraction

p[-1]:=1
q[-1]:=0
p[0]:=b[0]
q[0]:=1
p[j_]:=p[j]=b[j] p[j-1]+a[j]p[j-2]
q[j_]:=q[j]=b[j]q[j-1]+a[j]q[j-2];

Simple

cf=ContinuedFraction[Pi,20]

N[Pi,20]

conv=FromContinuedFraction/@Reverse[NestList[Most,cf,Length[cf]-1]]

N[conv,20]

CorrectDigits[{Pi,#}]&/@conv

%46+1

Transpose[Through[{Numerator,Denominator}[#]]&/@conv]

Compute

Weisstein: 2010–2013

Update terms computed in A001203

(10^5)

Length[cf=ContinuedFraction[N[Pi,10^5]]]//Timing

(10^6)

Length[cf=ContinuedFraction[N[Pi,10^6]]]//Timing

Max[cf]

2 (10^7)

Length[cf=ContinuedFraction[N[Pi,20 10^6]]]//Timing

Max[cf]

Put[cf,"PiContinuedFraction20M.txt"]//Timing

5 (10^7)

Length[cf=ContinuedFraction[N[Pi,50 10^6]]]//Timing

Max[cf]

Put[cf,"PiContinuedFraction50M.txt"]//Timing

(10^8)

Length[cf=ContinuedFraction[N[Pi,10^8]]]//Timing

5 (10^8)

Length[cf=ContinuedFraction[N[Pi,5*^8]]]//Timing

Max[cf]

(10^9)

Length[cf=ContinuedFraction[N[Pi,10^9]]]//Timing

Max[cf]

2×(10^9)

{Length[cf=ContinuedFraction[N[Pi,2 10^9]]]//Timing,MaxMemoryUsed[]/10.^9GB}

Max[cf]

3×(10^9)

$Version

{Length[cf=ContinuedFraction[N[Pi,3 10^9]]]//Timing,MaxMemoryUsed[]/10.^9GB}

Date[]

Max[cf]

ByteCount[cf]/10^9.GB

HighWaterMarks[cf]

MaxMemoryUsed[]/10.^9GB

(pos=Flatten[Position[cf,#,1,1]&/@Range[10]])//Timing

...can't do the following; efficient but uses >140 GB of RAM...

p=SplitBy[SortBy[MapIndexed[{#1,#2[[1]]}&,Take[cf,100]],First],First][[All,1]]

p=SplitBy[SortBy[MapIndexed[{#1,#2[[1]]}&,cf],First],First][[All,1]];

5×(10^9)

$Version

{Date[],Length[cf=ContinuedFraction[N[Pi,5 10^9]]]//Timing,MaxMemoryUsed[]/10.^9GB,ByteCount[cf]/10^9.GB}

Max[cf]

Developer`PackedArrayQ[cf]

6×(10^9)

$Version

{Date[],Length[cf=ContinuedFraction[N[Pi,6 10^9]]]//Timing,MaxMemoryUsed[]/10.^9GB,ByteCount[cf]/10^9.GB}

Max[cf]//Timing

1.1×(10^10)

{$MachineName,$TopDirectory,$Version,DateString[]}

Precision[pi]

Length[cf=ContinuedFraction[pi]]//Timing

{ByteCount[cf],MaxMemoryUsed[]}/10^9. GB

DumpSave["/Developer/eww/ContinuedFractions/Pi_10672905501.mx",cf];

Precomputed

SetDirectory["/Users/eww/Data/Math"]

Length[cf=Get["PiContinuedFraction20M.txt"]]//Timing

Length[cf=Get["PiContinuedFraction50M.txt"]]//Timing

Length[cf=Get["PiContinuedFraction100M.txt"]]//Timing

Length[cf=ReadList["PiContinuedFraction500M.txt"]]//Timing

Length[cf=ReadList["PiContinuedFraction500M.txt"]]//Timing

Get["/Developer/eww/ContinuedFractions/Pi_10672905501.mx"];//Timing//AbsoluteTiming

{Length[cf],ByteCount[cf]/10^9. GB}

Syed Fahad: Apr 25, 2021

On 4/26/21 12:40 AM, Bill Gosper wrote:
> ---------- Forwarded message ---------
> From: Fahad Syed <sydfhd@gmail.com>
> Date: Sun, Apr 25, 2021 at 9:32 PM
> Subject: 30B terms of regular continued fraction of Pi
> To: <billgosper@gmail.com>
> 
> 
> Hello,
> 
> I am Syed Fahad, an EE undergrad at National Institute of Technology,
> Srinagar. Last year in March, I contacted you regarding my interest in
> computing and studying the regular continued fraction of mathematical
> constants. Well, it's been a long year since, but I've written this program
> called "Zuuv" (GitHub <https://github.com/sinandredemption/zuuv>) which
> computes the terms of continued fraction of a number from it's hex/dec
> floating-point representation. I wanted to inform you about the progress
> I've made since, and some issues and future plans.
> 
> I computed 30,113,021,586 terms of the simple continued fraction of Pi
> (from 30B hex digits). I have uploaded all the terms here: link to Google
> Drive folder
> <https://drive.google.com/drive/folders/1--Qh9Xxq1i6oeHnTXzKrQ9FoguHreBKy>
> (each
> zipped file contains 200M-400M terms separated by a newline). This
> computation adds a *new term* to A033089 <https://oeis.org/A033089> and
> A033090 <https://oeis.org/A033090>: we have a new term *52662113289* at
> position *15621034283*.
> 
> Algorithms
> The algorithm I used to calculate regular continued fraction terms from a
> num-den pair was basically this:
> 1. Recursively compute regular continued fraction of the upper half of the
> fraction.
> 2. Calculate correction term based on calculation done in step 1.
> 3. Compute regular continued fraction of the correction term.
> 4. Merge the two computed lists and return.
> The algorithm achieves log-linear runtime. I can discuss the algorithm more
> rigorously with the mathematics behind it if that interests you.
> 
> Issues
> On my poor laptop, I could do 30B terms with 8GB of RAM in less than a
> week. However, Zuuv barely uses multithreading at all, because of the
> sequential nature of the algorithm. So the program would really benefit
> from a parallel version of this algorithm, however I wasn't able to derive
> something as such.
> 
> Future plans and verifying the computation
> I am in talks with people with capable hardware and willingness to dedicate
> it for months. However, the issue with large computations would be that
> verification of the computation would be necessary. The idea I currently
> have is to use the theorem that if | r / s - a | < 1 / (2 * s^2) then r, s
> are convergents to a. However, I'm open to more efficient ideas.
> 
> I am currently in the process of calculating some 15B terms of regular
> continued fraction of Euler-Mascheroni constant. If there is some kind of
> analysis I could perform using the results of my computations (besides the
> standard of checking compliance with Gauss-Kuzmin distribution), then
> please do provide me some ideas/pointers.
> 
> I have also requested an account on OEIS.
> 
> Best,
> Syed Fahad

Length[cf=ReadList["/Users/eww/Documents/Math/ContinuedFractions/Pi/Fahad/iteration1.txt"]]//Timing//AbsoluteTiming

Length[cf=ReadList["/Users/eww/Documents/Math/ContinuedFractions/Pi/Fahad/iteration1.txt",Number]]//Timing//AbsoluteTiming

Length[cf=ReadList["/Users/eww/Documents/Math/ContinuedFractions/Pi/Fahad/iteration1.txt",String,NullRecords->True]]//Timing//AbsoluteTiming

>13 GB

Length[cf=Import["/Users/eww/Documents/Math/ContinuedFractions/Pi/Fahad/iteration1.txt","List"]]//Timing//AbsoluteTiming

Length[cf=Import["/Users/eww/Documents/Math/ContinuedFractions/Pi/Fahad/iteration1.txt","String"]]//Timing//AbsoluteTiming

StringTake[cf,4]//FullForm

Length[cfsplit=StringSplit[cf,"\r\n"]]//Timing//AbsoluteTiming

>11 GB

(cnt=StringCount[cf,"\r\n"])//Timing//AbsoluteTiming

ByteCount[cf=ConstantArray[0,cnt]]//Timing//AbsoluteTiming

Do[

?ReadList

High water marks

terms: A033089
indices: A033090

5×(10^9)

{HighWaterMarks[cf,Method->"LowMemory"],MaxMemoryUsed[]/10.^9GB}//Timing

6×(10^9)

{HighWaterMarks[cf,Method->"LowMemory"],MaxMemoryUsed[]/10.^9GB}//Timing

1.1×(10^10)

Length[cf]

(hwm=HighWaterMarks[cf])//Timing

(hwm=HighWaterMarks[cf,Method->"LowMemory"])//Timing

First Occurrences

A032523

5×(10^9)

(l=FirstOccurrences[Take[cf,10^6,Method->"LowMemory"],10^3]);//Timing

pos=Position[l,0,1,1]

(l=FirstOccurrences[cf,10^6,Method->"LowMemory"]);//Timing

ByteCount[l]/10.^9GB

pos=Position[l,0,1,1]

6×(10^9)

(l=FirstOccurrences[cf,10^6,Method->"LowMemory"]);//Timing

ByteCount[l]/10.^9GB

pos=Position[l,0,1,1]

1.1×(10^10)

Read

Get["/Developer/eww/ContinuedFractions/Pi_10672905501.mx"];//Timing//AbsoluteTiming

Length[cf]

MaxMemoryUsed[]/10^9. GB

Loop

fo=FirstOccurrences[cf];//Timing

p=Take[Position[fo,0]//Flatten,10]

Split

fo=FirstOccurrences[cf,47086,Method->"Split"];//Timing

Out of memory

Position through first missing

fo=FirstOccurrences[cf,47086];//Timing

Length[fo]

(p=Position[cf,47087,1,1])//Timing

BFile["~/b032523.txt",fo]

ListPlot[Most[fo],AxesLabel->{Subscript[a,n],n},PlotRange->All]

Clear[fo]

n-Digit Terms

1.1×(10^11)

Get["/Developer/eww/ContinuedFractions/Pi_10672905501.mx"];//Timing

Length[cf]

ByteCount[cf]

digs=IntegerLength[cf];//Timing

Developer`PackedArrayQ[digs]

Max[digs]//Timing

positions=Table[Flatten[Position[digs,n,1,20]],{n,10}];//Timing

data=MapIndexed[{#2[[1]],{#,cf[[#]]}}&,positions]

Grid[{#[[1]],Column@#[[2]]}&/@data,Dividers->All,Alignment->Left]

Khinchin and Lévy Constant

KhinchinList[x_,n_:200]:=MapIndexed[#1^(1/#2[[1]])&,Rest@FoldList[Times,1,Rest[ContinuedFraction[Mod[x,1],n+1]]]]
LevyList[x_,n_:200]:=MapIndexed[#1^(1/#2[[1]])&,Denominator/@FromContinuedFraction/@Reverse[NestList[Drop[#,-1]&,Rest[ContinuedFraction[Mod[x,1],n+1]],n]]]

<<Utilities`Typesetting`

With[{k=5000,c=Pi},
g1=ListPlot[KhinchinList[c,k],Joined->True,PlotStyle->Red,AxesLabel->{n,((Subscript[a, 1])(Subscript[⋯a, n]))^InlineFraction[1,n]},
PlotLabel->"Khinchin's constant",
Epilog->{Blue,Dashing[{.02}],Line[{{0,#},{k,#}}]&[2.685452001]}];
g2=ListPlot[LevyList[c,k],Joined->True,PlotStyle->Red,AxesLabel->{n,(Subscript[q, n])^InlineFraction[1,n]},
PlotLabel->"Lévy constant",
Epilog->{Blue,Dashing[{.02}],Line[{{0,#},{k,#}}]&[Exp[Pi^2/12/Log[2]]]}
];
]//Timing

GraphicsRow[{g1,g2},Spacings->-50,ImageSize->400]