Real quadratic field fundamental unit

New version: case free
Compute the fundamental unit naively and by continued fractions

cfcasefree[OFn_,opts___]:=Module[

{disc,
a1,b1,offset,h,k,
α0,α,q0,q,hk,
norm,
cfOnly,
reportClassNumber
},

IsSquare[b_]:=Module[
{},
Return[b==1||Union[Mod[Flatten[FactorInteger[b][[All,2]]],2]]=={0} ]
];

NormFcn:=
(hk[[Length[hk],1]]
+disc hk[[Length[hk],2]]/2)^2
-disc hk[[Length[hk],2]]^2/4;

PrintResult[h_,k_]:=Module[
{},
If[Mod[disc,4]==1,
Print[StringForm[" `` + `` (``+Sqrt(``))/2",h,k,disc,disc]],
Print[StringForm[" `` + `` Sqrt(``)",h+disc/2 k,k,disc/4]]]
];

cfOnly=CfOnly/.{opts}/.{CfOnly->False};

reportClassNumber=ReportClassNumber/.{opts}/.{ReportClassNumber->True};

disc=If[Mod[OFn,4]==1,OFn,4OFn];
α0=(-disc+Sqrt[disc])/2(*+Floor[(disc+Sqrt[disc])/2]*);
α={α0};
q0=Floor[α0];
q={q0};
hk={{1,0},{q0,1}};

Print[StringForm["Q(sqrt ``)",OFn]];

If[!cfOnly,
Print["Fundamental unit by naive search"];
b1=1;
While[(!IsSquare[disc b1^2-4])&&(!IsSquare[disc b1^2+4]),
b1++];
Print[" Number of steps: ",b1];
If[IsSquare[disc b1^2-4],
a1=Sqrt[disc b1^2-4],
a1=Sqrt[disc b1^2+4]
];
{h,k}={(a1-disc b1)/2,b1};
PrintResult[h,k]
];

Print["Fundamental unit by continued fractions"];
norm=NormFcn;
While[(*!norm==1*)Abs[norm]>1,AppendTo[α,1/(α[[-1]]-q[[-1]])];
AppendTo[q,Floor[α[[-1]]]];
AppendTo[hk,hk[[-1]]*q[[-1]]+hk[[-2]]];
norm=NormFcn
];
Print[" Number of steps: ",Length[hk]-1];
PrintResult[hk[[-1,1]],hk[[-1,2]]];

If[!cfOnly,
If[{h,k}=={hk[[-1,1]],hk[[-1,2]]},
Print[Style["The two results agree",Darker[Green]]],
Print[Style["The two results disagree",Red]]
]
];

If[reportClassNumber,
Print["The class number formula gives h = ",
N[-Sum[JacobiSymbol[disc,r]Log[Sin[Pi r/disc]],{r,Floor[disc/2]}]/Log[hk[[-1,1]]+hk[[-1,2]](disc+Sqrt[disc])/2]]]
]
];

New version: case free
Calls

Q((Sqrt[5]))

cfcasefree[5]

Q((Sqrt[13]))

cfcasefree[13]

Q((Sqrt[31]))

cfcasefree[31]

Q((Sqrt[94]))

cfcasefree[94,CfOnly->True]

Q((Sqrt[n])) for a range of squarefree n

{nmin,nmax}={102,130};
For[n=nmin, n≤nmax,n++,
If[SquareFreeQ[n],
cfcasefree[n];
Print[""]
]
]

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

{nmin,nmax}={10,11};
For[n=nmin,n≤nmax,n++,
Print[StringForm["Q(sqrt(4^``-1))",n]];cfcasefree[4^n-1];
Print[""]
]

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

{nmin,nmax}={10,11};
For[n=nmin,n≤nmax,n++,
Print[StringForm["Q(sqrt(4^``+1))",n]];cfcasefree[4^n+1];
Print[""]
]

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

{nmin,nmax}={10,11};
For[n=nmin,n≤nmax,n++,
Print[StringForm["Q(sqrt(4^``+2))",n]];cfcasefree[4^n+2];
Print[""]
]

Q((Sqrt[1069]))

cfcasefree[1069,CfOnly->True]

Q((Sqrt[1141]))

cfcasefree[1141,CfOnly->True]

Old version: casewise
Compute the fundamental unit naively and by continued fractions

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

{a1,b1,offset,h,k,
α={α0},q={Floor[α0]},
hk={{1,0},{Floor[α0],1}},
norm
},

IsSquare[b_]:=Module[
{},
Return[b==1||Union[Mod[Flatten[FactorInteger[b][[All,2]]],2]]=={0} ]
];

NormFcn:=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 
];

PrintResult[h_,k_]:=Module[
{},
If[Mod[OFn,4]==1,
Print[StringForm[" `` + `` (1+Sqrt(``))/2",h,k,OFn]],
Print[StringForm[" `` + `` Sqrt(``)",h,k,OFn]]
];
];

Print[StringForm["Q(sqrt ``)",OFn]];

Print["Fundamental unit by naive search"];
b1=1;
offset=If[Mod[OFn,4]==1,4,1];
While[(!IsSquare[OFn b1^2-offset])&&(!IsSquare[OFn b1^2+offset]),b1++];
Print[" Number of steps: ",b1];
If[IsSquare[OFn b1^2-offset],
a1=Sqrt[OFn b1^2-offset],
a1=Sqrt[OFn b1^2+offset]
];
{h,k}=If[Mod[OFn,4]==1,{(a1-b1)/2,b1},{a1,b1}];
PrintResult[h,k];

Print["Fundamental unit by continued fractions"];
norm=NormFcn;
While[(*!norm==1*)Abs[norm]>1,AppendTo[α,1/(α[[-1]]-q[[-1]])];
AppendTo[q,Floor[α[[-1]]]];
AppendTo[hk,hk[[-1]]*q[[-1]]+hk[[-2]]];
norm=NormFcn
];
Print[" Number of steps: ",Length[hk]-1];
PrintResult[hk[[-1,1]],hk[[-1,2]]];

If[{h,k}=={hk[[-1,1]],hk[[-1,2]]},
Print[Style["The two results agree",Darker[Green]]],
Print[Style["The two results disagree",Red]]
]
];

Old version: casewise
Calls

(1+(Sqrt[5]))/2

cf[5,(-1+Sqrt[5])/2]

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

{nmin,nmax}={102,130};
For[n=nmin, n≤nmax,n++,
If[SquareFreeQ[n],
r=If[Mod[n,4]==1,(-1+Sqrt[n])/2,Sqrt[n]];
cf[n,r];
Print[""]
]
]

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

{nmin,nmax}={10,13};
For[n=nmin,n≤nmax,n++,
r=Sqrt[4^n-1];Print[StringForm["Q(sqrt(4^``-1))",n]];cf[4^n-1,r];
Print[""]
]

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

{nmin,nmax}={10,13};
For[n=nmin,n≤nmax,n++,
r=(-1+Sqrt[4^n+1])/2;Print[StringForm["Q(sqrt(4^``+1))",n]];cf[4^n+1,r];
Print[""]
]