I'm working on a project where I compute the periodic part of the continued fraction of the square root of a (nonsquare) positive integer.
The obvious approach is to take the numerical square root, compute the continued fraction (dropping the last, say, two terms), and see if it ends in a periodic part with enough members to feel confident in the pattern, maybe three. If not, increase precision and try again.
But I was concerned about the possibility of numerical error, so I rolled my own implementation which avoids t_REALs. I think the algorithm is classic; I modified it slightly from Beceanu 2003
if(D<0 || issquare(D), error("D must be a positive rational nonsquare"));
\\ Set a[n], b[n+1], and c[n+1]
\\ Check if a loop is formed
if(a[i]==a[n] && b[i+1]==b[n+1] && c[i+1]==c[n+1],
\\ Success! See if the loop can be split into a smaller loop -- this
\\ should only be possible if the parameters b and c have a cycle
\\ length which is a nontrivial multuple of the cycle length of a,
\\ which I'm not sure is possible.
for(j=1,k-d, if(loop[j]!=loop[j+d], next(2)));
\\ Least period is of length d
return([,[a[1..d]]]) \\ Unlikely to happen, but possible I suppose
addhelp(cf, "cf(D): Returns [u, v], where v is the periodic part of the continued fraction for sqrt(D) and u is the preperiodic part, where D is a positive rational nonsquare.");
This works, but it's slow. Is there a better way to do this? Are there tight enough bounds so I could prove that the approximation is close enough that the apparent period is in fact the period? Any other comments on my idea or code?