[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
More information about the Haskell-Cafe
mailing list