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

ajb at spamcop.net ajb at spamcop.net
Mon Oct 16 01:14:03 EDT 2006

```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
```