Max Alekseyev on Thu, 03 Aug 2017 19:35:28 +0200 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: Continued fractions for square roots |
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 2003cf(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 formedfor(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 da=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