Search code examples
debugginghaskellnumericsicp

bug in nthroot using fixedpoint


Here is my attempted solution to exercise 1-45 from Structure and Interpretation of Computer Programs, implemented in Haskell:

fixedpoint :: (Double->Double) -> Double -> Double -> Double
fixedpoint f guess err
  | err > abs next-guess  = next
  | otherwise             = fixedpoint f next err
  where next = f guess

avg :: Double -> Double -> Double
avg a b = (a+b)/2

avgdamp :: (Double -> Double) -> (Double -> Double)
avgdamp f = \ x -> avg x (f x)

log2 :: Double -> Double
log2 = (/ log 2) . log

repeated :: (a->a) -> Int -> (a->a)
repeated f n = \ x -> (iterate f x) !! n

nthroot :: Int -> Double -> Double -> Double
nthroot n x err = fixedpoint (repeated avgdamp numdampsneeded $ eqn) 1.0 err where
  numdampsneeded = floor $ log2 $ fromIntegral n  -- experimentally determined
  eqn = \ y -> x/y^(n-1)

I tested the helper procedures one-by-one, they seemed fine. I looked at 3 blogs on which solutions are posted and mine seems to correspond to them with some minor exceptions like implementation language and that I use explicit error argument.

Nevertheless I get ridicolous results from nthroot, and I do not have a clue.

> nthroot 2 16 0.000000001
5.1911764705882355

> nthroot 3 8 0.000001
2.447530864197531

> nthroot 3 27 0.001
7.0688775510204085

> nthroot 5 32 0.001
6.563864764681383

or in this case I expect to get back (approximation of) [0..16] as in octave

 > [0:16].^(1/2).^2

ans =
    0.00000    1.00000    2.00000    3.00000    4.00000    5.00000    6.00000    7.00000    8.00000    9.00000   10.00000   11.00000   12.00000   13.00000   14.00000   15.00000   16.00000

but

> map (^2) $ map (\x -> nthroot 2 x 0.000001) [0..16]
[0.25,1.0,2.006944444444444,3.0625,4.2025,5.4444444444444455,6.79719387755102,8.265625,9.852623456790123,11.559999999999999,13.388946280991735,15.340277777777777,17.41457100591716,19.612244897959187,21.933611111111112,24.37890625,26.94831314878893]

Any guess of what could be behind these bad results is welcome.


Solution

  • Function application binds tighter than any infix operator. Therefore,

    abs next - guess  ==  (abs next) - guess
    

    which isn't really what you want. If you fix this issue, your algorithm works fine.