# [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,

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
where
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
myself.

Best,

Tamas
```