Vraj Mohan writes:

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.

Generally it's better to separate out the different parts of the algorithm. So sqrt_step x candidate = (candidate + x/candidate)/2 Now you can try take 20 $ iterate (sqrt_step 2) 1 and watch the convergence. When you're happy with that, you can use something like “head . dropWhile unconverged”.

(The equivalent code is well-behaved on MIT Scheme)

Is it? Is there equivalent code to “my_sqrt :: Float -> Float”? (that might be pertinent). -- Jón Fairbairn Jon.Fairbairn@cl.cam.ac.uk http://www.chaos.org.uk/~jf/Stuff-I-dont-want.html (updated 2006-09-13)

Vraj Mohan wrote:

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?

1) As Jon said, by seperating the iteration from the selection of the result. iterate, filter, head, etc. are your friends. Makes the code more readable and maintainable, too. 2) By learning more about floating point numbers. There is no Float y such that | y*y - 96 | < 0.00001. That's why such an iteration is better terminated not when the result is good enough, but when it stops getting better. It's also the reason why you code might work with -fexcess-precision or in an untyped language or with the next or previous compiler release or on rainy days or whatever else. Udo. -- Walk softly and carry a BFG-9000.

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: http://www.haskell.org/pipermail/haskell-cafe/2006-September/017858.html 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

Vraj Mohan wrote:

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)

There was some excellent advice in the other responses, but I thought it worth mentioning that your Haskell code converges if you step up from Float -> Float to Double -> Double.