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.