Below is a description of an algorithm which, with
each iteration, will double the number of significant 
digits in the computation of a rational square root
approximation.

  I do not know if this algorithm is new, but I found
it interesting nonetheless.  Lisp code for implementing
this algorithm can be found here.



If by convention we say that:

  s = isqrt(C)

  d = C - s^2

  Then the square root of C can be expressed as as the infinite
continued fraction:

         d
  s + -----------------
              d
      2s + ------------
                   d
           2s + -------
                
                2s + .. 

We designate the tail of this continued fraction using

  e(n) = nth error term

and we say that the nth error term of the continued fraction
representation has the form:

              d
  e(n) =  ----------
          2s + e(n+1)

  Although we sometimes drop the subscript on the error term for
notational convenience.

  A pretty good first order approximation for a square root can
be computed as follows (even allowing e(1) to be zero):

                    d  
  sqrt(C) ~= s + --------
                 2s + e(1)

  Without proof, we claim that a generalized expression for a partial
evaluation of our continued fraction can be represented as:

                A + Be
  sqrt(C) = s + -------
                C + De

  It is easy to see that the first order approximation given above is
an instance of this expression when A = d, B = 0, C = 2s, D = 1.  The
generalized representation is useful for representing the result of
evaluating some number of successive terms in the continued fraction
representation.

  Assuming that the above representation is the result of evaluating n
terms of the continued fraction for sqrt(C), then the n+1 term would be
computed by substituting the next term into the error expression in
the representation.

      A + B(d/(2s + e))
  s + ---------------
      C + D(d/(2s + e))

      A(2s + e) + Bd
  s + --------------
      C(2s + e) + Dd
  
      (A2s + Bd) + Ae
  s + -----------------
      (C2s + Dd) + Ce

  Which, we observe, is once again in the general representational
form.  

If the evaluation of the first n terms produced

  A + Be
  ------
  C + De
  
and the evaluation of the next m terms were to produce

  W + Xe
  ------
  Y + Ze
  
Then the evaluation of the first n+m terms would be:

        W + Xe 
  A + B(------)
        Y + Ze      A(Y + Ze) + B(W + Xe)
  ------------  =   ---------------------
        W + Xe      C(Y + Ze) + D(W + Xe)
  C + D(------)
        Y + Ze
  
  (AY + BW) + (AZ + BX)e
  ----------------------
  (CY + DW) + (CZ + DX)e
  
  Because the continued fraction representation of the square root is
uniform, the evaluation of any n sucessive terms will always produce
the same result.

  We can take advantage of this fact to refine our first order
approximation by substituting A,B,C, and D in for W,X,Y and Z in the
above expression.  The will double the number of terms we have
evaluated.  Of course, this procedure can be repeated again and again,
with each iteration of the algorithm doubling the number of
significant digits in our representation.

 Here is an example run computing the sqrt of 1973 for 1 to 10
iterations of the algorithm.  Note that after iteration 6 we have
more than 100 significant digits.

 1 : 44.4184552114124148567022233646060917619843207813905667651972754144711476673949363834982650044981364863
 2 : 44.4184646288146721950566598272783899534327268025085249602382778806857857085449370627285274658896791048
 3 : 44.4184646290256187642752431232557204024098591484297090342230548605659661647016608600579504981095676558
 4 : 44.4184646290256187643810796574090605395683327294135496111246850857337544379108155113103260155775597834
 5 : 44.4184646290256187643810796574090605395949744270465990361024620576194006618026805283703947239542091268
 6 : 44.4184646290256187643810796574090605395949744270465990361024620576194006618043686917147360058911830087
 7 : 44.4184646290256187643810796574090605395949744270465990361024620576194006618043686917147360058911830087
 8 : 44.4184646290256187643810796574090605395949744270465990361024620576194006618043686917147360058911830087
 9 : 44.4184646290256187643810796574090605395949744270465990361024620576194006618043686917147360058911830087
10 : 44.4184646290256187643810796574090605395949744270465990361024620576194006618043686917147360058911830087

Lisp code for implementing this algorithm can be found here.