Re: Continued fractions for square roots

[Karim Belabas <Karim.Belabas@math.u-bordeaux.fr>]

* Charles Greathouse [2017-08-02 22:29]:
> 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],

This last loop is going to kill you ( O(n^2), where n ~ D^(1/2) is expected ).
You should use a map (or a hashtable in C...)

Something like this [ completely untested :-) ]

cf(D)={
  my(a0=sqrtint(D\1), a=List(a0), b=List(a0), c=List(D-a0^2), n=1);
  if(D<0 || issquare(D), error("D must be a positive rational nonsquare"));
  my(M=Map());
  mapput(M,[a0,a0,c[1]], n);
  while(1,
    my (A,B,C); \\ Set a[n+1], b[n+1], and c[n+1]
    listput(a, A = (a0+b[n])\c[n]);
    listput(b, B = A*c[n]-b[n]);
    listput(c, C = (D-B^2)\c[n]);
    if (!mapisdefined(M, [A,B,C]), n++; mapput(M, [A,B,C], n); next);
    \\ a loop is formed
    i = mapget(M,[A,B,C]);
    my(k=n+1-i,loop=vector(k,j,a[n+1-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=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
    )
  );
}


Cheers,

    K.B.
--
Karim Belabas, IMB (UMR 5251)  Tel: (+33) (0)5 40 00 26 17
Universite de Bordeaux         Fax: (+33) (0)5 40 00 21 23
351, cours de la Liberation    http://www.math.u-bordeaux.fr/~kbelabas/
F-33405 Talence (France)       http://pari.math.u-bordeaux.fr/  [PARI/GP]
`