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

Tamas K Papp tpapp at Princeton.EDU
Sun Oct 15 11:21:16 EDT 2006

On Sun, Oct 15, 2006 at 12:17:20PM +0000, Vraj Mohan wrote:
> I am new to Haskell and need help in debugging my code.
> I wrote the following function for calculating square roots using Newton's 
> method:
> my_sqrt :: Float -> Float
> my_sqrt x = improve 1 x
>          where improve y x = if abs (y * y - x) < epsilon 
>                                 then y 
>                                 else improve ((y + (x/y))/ 2) x
>                epsilon = 0.00001
> This works for several examples that I tried out but goes into an infinite loop
> for my_sqrt 96. How do I go about debugging this code in GHC or Hugs?
> (The equivalent code is well-behaved on MIT Scheme)

Hi Vraj,

1. First, generally about debugging: I asked this question a month
ago, on this mailing list.  See this thread:


2. Newton's method is not guaranteed to converge.  For examples, and a
nice introduction to uni- and multivariate rootfinding, see Numerical
Methods in Economics, Kenneth L Judd, Chapter 5.

I suggest that you use bisection, it is much more robust.  I have
written a bisection routine which I have been planning to post after a
bit of testing, maybe you will find it useful:

bisection tol reltol f a b | tol < 0 || reltol < 0 =
                               error "bisection needs positive tolerances"
                           | abs fa <= tol = Just a -- a is a root
                           | abs fb <= tol = Just b -- b is a root
                           | fa*fb > 0 = Nothing -- IVT not guaranteed
                           | a > b = bis b a fb fa -- exchange endpoints
                           | otherwise = bis a b fa fb -- standard case
      fa = f a                  -- function evaluated only once at each
      fb = f b                  -- endpoint, store them here
      bis a b fa fb = let m = (a+b)/2
                          fm = f m in
                      if (b-a) <= reltol*(1+abs(a)+abs(b)) || abs fm <= tol
                      then Just m
                      -- implicit assumption: fa * fb > 0
                      else if fm*fa < 0 -- new interval: [a,m]
                           then bis a m fa fm
                           else bis m b fm fb -- new interval: [m,b]

-- uniroot is a specialisation of bisection, with commonly used
-- tolerance values
uniroot = bisection 1e-20 1e-15

Untested code:

mysqrt x = uniroot f 0 x
       where f y = y**2-x

It will give you a Maybe a type value, ie Nothing for negative
numbers.  Comments about my code are of course welcome, I am a newbie



More information about the Haskell-Cafe mailing list