Re: Calculating SQRT in ADA

["Howard W. LUDWIG" <howard.w.ludwig@lmco.com>]

David C. Hoos, Sr. wrote:

> Mark A Biggar wrote in message <36FAA3DF.42C31CF@lmco.com>...
> >Hans Marqvardsen wrote:
> >>
> >> John Herro wrote:
> >
> >No,no,no,no!
> >
> >This algorithm has quadratic convergence.  This means that each time
> through the
> >loop doubles the number of correct bits in the answer.  So with IEEE double
> >precession float you need loop at most 6 times, and even with extended
> float
> >7 times, so don't test against a torerance at all, just loop 6 or 7 times
> >and be done with it.
> >
> No,no,no,no! (to coin a phrase)
>
> The number of correct bits doubles with every iteration, only if you are
> counting place holding bits as well as mantissa bits.  So,
> if you take the square root of 2.0e100, for example, 173 iterations are
> required if the argument is used as the initial estimate.
>
> similar number is required for 0.5e-100, but only six iterations are
> required for 2.0e00 or 0.5e00.

No, no, no, no! (to plagiarize a coined phrase)

It is a common oversimplified assumption that Newton[-Raphson] has quadratic
convergence.  There are many well-known cases where the method does not converge
at all.  Even when the method does converge, the convergence can be linear.
Quadratic convergence applies when the quadratic term is negligible compared to
the linear term in the Taylor series expansion of the function in question.
When the initial guess is far off, "negligible" does not describe the situation.

To find the actual  a = sqrt(c), where c is known, involves calculating
    x(n+1) = 0.5 * (x(n) + c/x(n)).
Let x(n) = a + d(n), where d(n) measures the deviation of the estimate x(n) from
the actual value a.  Then
    d(n+1) = x(n+1) - a
           = 0.5 * (a + d(n) + a**2/(a + d(n))) - a
           = 0.5 * d(n)**2/(a + d(n)).
When d(n) is small compared to a, the denominator behaves more and more like a,
so the deviation is indeed proportional to d(n)**2.  However, when d(n) is large
compared to a, the denominator for some iterations looks more like d(n),
canceling out one of the d(n) factors in the numerator, so the deviation is cut
each time by a factor of 2, thus linear convergence, until the iterations get to
the place of having the deviation be of the order of a, at which time the
convergence transitions to quadratic.

In the case of the example given [c = 2E100, a =1.4142...E50], starting with an
initial guess of c requires 166 iterations of halving the deviation each time to
get to the point of the deviation being smaller than the actual root a.  Then
quadratic convergence kicks in, apparently requiring only another 7 iterations
to achieve the ultimate desired result. If we simply counted place-holder bits,
roughly 167, plus the desired 53 bits of precision, then we have about 220 bits
to deal with.  _Assuming_ the first bit to be fine and doubling the number of
good bits each iteration, only 8 iterations would be needed (2^8 = 256), nowhere
close to the 173 iterations quoted above.  This shows that you cannot assume
that the first bit is good.  If you want to do a "back of the envelope"
calculation based on quadratic convergence, the 2E100 case using 2E100 as the
initial guess would say your first guess has less than 1E-50 bit good, in order
to double the number of good bits each iteration and to obtain 53 good bits
after 173 iterations.  A better way of looking at it is that the place-holder
bits are lopped off at the rate of one per iteration for the first 166
iterations.

HWL