Haskell - Recursion Stack Overflow

4
votes

I am trying to sum all n from 1 to a very large number (10**9 for now) but it gives stack overflow. Also I don't think putting a stop at 1 and doing the sum n in different rows is the most efficient way but the code below is all my knowledge for Haskell. I really don't know much about functional programming and I would like as much explanation as possible. (Also I have tried putting $! strict in the last line which was told in other places but it changed nothing. I would be glad if you explained the most efficient way I could do this recursive function.)

main :: IO()

summ 1 = 1
summ n = 1/(n**2) + summ (n-1)

expected_value = pi*pi/6
errorPercent n = n / expected_value * 100

main = do
    putStr "%"
    print (errorPercent (summ $! (10**9)))
2
haskellrecursionstack-overflow
Oh by the way you can write summ as summ n = sum [ 1 / i^2 | i <- [1..n] ]. Nice, right? (Make sure to compile with optimizations for this one) – luqui
@luqui - I think you'd want sum [ 1 / i^2 | i <- reverse [1..n] ] to avoid loss of precision. – Omnifarious
@Omnifarious interesting! I had never thought about issues like this. I can't quite see why it happens though... – luqui
@luqui - Happens due to rounding. Floating point numbers work best when the numbers being added together are similar in magnitude. Otherwise, the bits being added from the small number sort of rotate off the end of the large number and don't even change it. So, even if you add several hundred thousand tiny numbers like that they don't change the original number at all. – Omnifarious
@luqui - As a rule of thumb, yes. I can also imagine that in some specific cases different orderings may be better. There are other rules for other operations as well. I don't know what they all are, but I do know that making the most of floating point precision is something to be thought through carefully. – Omnifarious

2 Answers

10
votes

chi has answered one bit of the question, which I think is the main problem, but there is something else that is bugging me. When you say 10**9, you get a floating point number (because ** is "fractional" exponentiation). And then you are using floating point equality to check for the base case of your recursion. 

summ 1 = ...

The problem with this is that it is possible, and as the argument gets larger, likely, that because of numerical error you will just barely miss the base case and descend into negative values forever.

summ 4 =        ... summ 3
summ 3 =        ... summ 2.000001
summ 2.000001 = ... summ 1.000001 
summ 1.000001 = ... summ 0.000001  -- BASE CASE MISSED!
summ 0.000001 = ... summ (-1.000001)
summ (-1.000001) = ... summ (-2.000001)

and so on. If you didn't get a stack overflow from 109 calls, you surely will with infinitely many.

You should either define your function on integers so there is no rounding error

summ :: Int -> Double
summ 1 = 1
summ n = 1 / (fromIntegral n ** 2) + summ (n - 1)
--            ^^^^^^^^^^^^
-- conversion necessary to go from Int to Double

main = ... print (summ (10 ^ 9))
--                      ^^^^^^
--      use integral exponentiation (^) instead of (**)

or use a more forgiving base case

summ :: Double -> Double
summ n | n <= 1 = 1
summ n = 1 / (n ** 2) + summ (n - 1)

In either case, you should definitely take chi's suggestion to do this in accumulator style, and you should also definitely put a type signature.

Here's more on how you get stack overflows in Haskell if you are curious.

7
votes

The problem here is that sums can not start being computed until the whole 10^9 recursion calls are over. Essentially, you are computing

1/(n**2) + ( 1/((n-1)**2) + ( 1/((n-2)**2) + ....

and the parentheses prevent to start summing. Instead, we would like to have

(( 1/(n**2) + 1/((n-1)**2) ) + 1/((n-2)**2) ) + ....

The easiest way is to use an "accumulator" additional argument:

summ 1 acc = 1 + acc
summ n acc = summ (n-1) $! acc + 1/(n**2)

main = do
    putStr "%"
    print (errorPercent (summ (10^9) 0))  -- set acc to 0 at the beginning

For improving the performance, I'd recommend to add a type signature to summ e.g. summ :: Int -> Double -> Double.

Full program below. This runs in 12s here (ghc -O2).

summ :: Int -> Double -> Double
summ 1 acc = 1 + acc
summ n acc = summ (n-1) $! acc + 1 / (fromIntegral n**2)

main :: IO ()
main = do
    putStr "%"
    print (summ (10^9) 0)  -- set acc to 0 at the beginning