Re: Continued fractions for square roots

On Wed, Aug 2, 2017 at 4:28 PM, Charles Greathouse <charles.greathouse@case.edu> wrote:

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?

Charles Greathouse
Case Western Reserve University