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