[Haskell-cafe] Debugging Newton's method for square roots

[ajb at spamcop.net]

G'day all.

Quoting Tamas K Papp <tpapp at Princeton.EDU>:

> 2. Newton's method is not guaranteed to converge.

For computing square roots, it is.  The square root function is
exceedingly well-behaved.

But you can make things better by choosing a smart initial estimate.
The initial estimate in this version is so good that you only need a
fixed number of iterations, so there are no loops here.

mysqrt :: Double -> Double
mysqrt x
    | x <= 0    = if x == 0 then 0 else error "mysqrt domain error"
    | r == 1    = sqrt2 * sqrtx
    | otherwise = sqrtx
    where
        (q,r) = exponent x `divMod` 2

        sqrtx = scaleFloat q (sqrtSignificand (significand x))

        -- Feel free to add more digits if this is insufficient
        sqrt2 = 1.41421356237309504880168872420969807856967


-- Compute the square root of a number in the range [0.5,1)
sqrtSignificand :: Double -> Double
sqrtSignificand x
    = iter . iter . iter $ d0
    where
        -- d0 is a third-order Chebyshev approximation to sqrt x
        x' = x * 4.0 - 3.0
        d2 = -6.177921038278929e-3
        d1 = 2 * x' * d2 + 0.14589875298243973
        d0 = x' * d1 - d2 + 0.5 * 1.7196949654923188

        iter sx = 0.5 * (sx + x / sx)

Cheers,
Andrew Bromage