Continued fraction code

Routine

This routine behaves differently for a rational input, or a quadratic irrational input with an option set to seek the fundamental unit of the quadratic field, or a general irrational input.  Another option tells the routine that an input is irrational, for inputs such as Log[3/2].  Other options control which parts of the computation are shown: 
* Digits: how many digits
* NumSteps: how many Euclidean algorithm steps for a generic irrational input
* ShowAlpha: whether to show the sequence of α-values in the Euclidean algorithm
* ShowConvergents: whether to show the convergents
* SizeCheck: whether to show the check that the convergents are approaching the limit as 1/k^2 where k is the denominator of the convergent.

cf[α0_,opts___]:=Module[

{α={α0},q={Floor[α0]},
hk={{1,0},{Floor[α0],1}},
digits=Digits/.{opts}/.{Digits->10},
irrationalSwitch=IrrationalSwitch/.{opts}/.{IrrationalSwitch->False},
oFgeneratorSwitch=OFgeneratorSwitch/.{opts}/.{OFgeneratorSwitch->False},
oFn,
norm,
step=0,numSteps,
showAlpha=ShowAlpha/.{opts}/.{ShowAlpha->False},
showConvergents=ShowConvergents/.{opts}/.{ShowConvergents->False},
sizeCheck=SizeCheck/.{opts}/.{SizeCheck->False}},

If[Element[α0,Rationals],
Print["α0 is rational"],
irrationalSwitch=True
];
If[irrationalSwitch,
Print["α0 is irrational"];
numSteps=NumSteps/.{opts}/.{NumSteps->10}
];
If[oFgeneratorSwitch,
oFn=OFn/.{opts};
Print["Quadratic field is Q(sqrt ",oFn,")"];
norm=If[Mod[oFn,4]==1,
(hk[[2,1]]+hk[[2,2]]/2)^2-hk[[2,2]]^2/4oFn ,
hk[[2,1]]^2-hk[[2,2]]^2oFn 
];
While[Abs[norm]>1,AppendTo[α,1/(α[[-1]]-q[[-1]])];
AppendTo[q,Floor[α[[-1]]]];
AppendTo[hk,hk[[-1]]q[[-1]]+hk[[-2]]];
norm=If[Mod[oFn,4]==1,
(hk[[Length[hk],1]]+hk[[Length[hk],2]]/2)^2-hk[[Length[hk],2]]^2/4oFn ,
hk[[Length[hk],1]]^2-hk[[Length[hk],2]]^2oFn 
];
],
While[If[!irrationalSwitch,!q[[-1]]==α[[-1]],step<numSteps],
AppendTo[α,1/(α[[-1]]-q[[-1]])];
AppendTo[q,Floor[α[[-1]]]];
AppendTo[hk,hk[[-1]]q[[-1]]+hk[[-2]]];
step++
];
];

Print["q = ",q];
If[showAlpha&&!irrationalSwitch,Print["α = ",α]];
If[showAlpha,Print["N[α] = ",N[α]]];
Print["{h, k} = ",hk];
If[oFgeneratorSwitch,
If[Mod[oFn,4]==1,
Print["h^2 + hk - (",oFn,"-1)/4 k^2 = ",(#[[1]]^2+#[[1]]#[[2]]-(oFn-1)/4 #[[2]]^2)&/@hk],
Print["h^2 - ",oFn," k^2 = ",(#[[1]]^2-oFn #[[2]]^2 )&/@hk]
]
];
If[showConvergents,
Print["{α0, h/k} = ",{N[α0,digits],
TableForm[N[hk[[#,1]]/hk[[#,2]],digits]&/@Range[2,Length[hk]]]}]
];
If[sizeCheck,Print["{1/k, |α0-h/k|, 1/k^2, sizecheck} = ",TableForm[{N[1/hk[[#,2]],digits],
N[Abs[α0-hk[[#,1]]/hk[[#,2]]],digits],N[1/hk[[#,2]]^2,digits],
Abs[α0-hk[[#,1]]/hk[[#,2]]]<1/hk[[#,2]]^2}&/@Range[2,Length[hk]]]]
];
];

Calls

963/731

cf[963/731,ShowAlpha->True,ShowConvergents->True,SizeCheck->True]

(Sqrt[13])

cf[Sqrt[13]]

e

cf[E,NumSteps->20]

π

cf[Pi,NumSteps->5,ShowConvergents->True,SizeCheck->True]

(Subscript[Log, 2])(3/2)

cf[Log[2,3/2],IrrationalSwitch->True]

(Sqrt[n]) or (1+(Sqrt[n]))/2 for a range of squarefree n

{nmin,nmax}={2,21}
(If[SquareFreeQ[#],
r=If[Mod[#,4]==1,(-1+Sqrt[#])/2,Sqrt[#]];Print[Style[StringForm["Continued fraction of ``",r],Bold]];cf[r,OFgeneratorSwitch->True,OFn->#]])&/@Range[nmin,nmax]

(Sqrt[(4^n)-1]) for a range of  n

{nmin,nmax}={10,13}
(r=Sqrt[4^#-1];Print[Style[StringForm["Continued fraction of ``",r],Bold]];cf[r,OFgeneratorSwitch->True,OFn->4^#-1])&/@Range[nmin,nmax]

(Sqrt[(4^n)+1]) for a range of  n

{nmin,nmax}={10,13}
(r=(-1+Sqrt[4^#+1])/2;Print[Style[StringForm["Continued fraction of ``",r],Bold]];cf[r,OFgeneratorSwitch->True,OFn->4^#+1])&/@Range[nmin,nmax]