rapidly converging rational sqrt

From: TheBean (TheBeaNerd_at_yahoo.com)
Date: 10/26/04 

Date: Tue, 26 Oct 2004 14:00:11 +0000 (UTC)

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 at:

  http://thegreves.com/david/sqrt/sqrt.html

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 sucessive 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 error 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

Dave

Relevant Pages

  • Re: No Set Contains Every Computable Natural (was Church-Turing compared to Zuse-Fredkin thesis)
    ... Russell is right about Turing saying this, ... You might not want to call this an algorithm which has ... decimal representation of x digit by digit. ...
    (comp.theory)
  • Re: printf "%.100f" the value 2 ^ (-127)
    ... > I was trying to imagine an algorithm that, for at least one finite FP ... > representation, would output an infinite series of digits. ... This range contains an infinite number of rational values. ... or as a periodic decimal fraction (e.g. ...
    (microsoft.public.vc.language)
  • Re: Modelling an OO data structure
    ... strike me as especially O-O representations. ... The choice of algorithm affects the representation. ... extending the programming language into the problem domain (a ...
    (comp.object)
  • Re: all paths between 2 nodes
    ... On 9 Apr 2004, Glenn C. Rhoads wrote: ... Under this representation, you ... one could easily implement this algorithm. ... The usual adjacency list can find a topsort in O ...
    (comp.theory)
  • Re: Cohens paper on byte order
    ... This was assumed in suggesting that the early evaluation of a bit stream ... that a fractional numeric representation did not have this limitation. ... >> string on a particular machine but this is irrelevant if the bit string ... semantics match those of the objects that they are manipulating. ...
    (sci.crypt)