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
cf(D)={
my(a0=sqrtint(D\1),a=List(),b=List([a0]),c=List([D-a0^2]),n=1);
if(D<0 || issquare(D), error("D must be a positive rational nonsquare"));
while(1,
\\ Set a[n], b[n+1], and c[n+1]
listput(a, (a0+b[n])\c[n]);
listput(b, a[n]*c[n]-b[n]);
listput(c, (D-b[n+1]^2)\c[n]);
\\ Check if a loop is formed
for(i=1,n-1,
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.
my(k=n-i,loop=vector(k,j,a[n-k+j]));
fordiv(k,d,
for(j=1,k-d, if(loop[j]!=loop[j+d], next(2)));
\\ Least period is of length d
a=concat(a0,Vec(a));
forstep(j=n-k,1,-1,
if(a[j]!=a[j+d],
return([a[1..j], a[j+1..j+d]])
)
);
return([[],[a[1..d]]]) \\ Unlikely to happen, but possible I suppose
)
)
);
n++;
);
}
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?