4

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)))
Improve this question
7
  • 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) Commented Feb 14, 2019 at 21:01
  • @luqui - I think you'd want sum [ 1 / i^2 | i <- reverse [1..n] ] to avoid loss of precision. Commented Apr 30, 2019 at 1:13
  • @Omnifarious interesting! I had never thought about issues like this. I can't quite see why it happens though... Commented Apr 30, 2019 at 21:47
  • @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. Commented Apr 30, 2019 at 22:37
  • 1
    @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. Commented Apr 30, 2019 at 23:20

2 Answers 2

Reset to default
10

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.

Improve this answer
Sign up to request clarification or add additional context in comments.

Comments

7

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
Improve this answer

8 Comments

It has worked and I could get output from 10^9 but it is really slow. I am currently waiting for 10^10 and 10^9 took like 20 minutes. Is this speed normal or is there ways to make it faster?
@Terobero Are you running the code after compilation with optimization e.g. -O2? If you are running it using GHCi instead, that is usually quite slow.
@Terobero I have no idea about VSCode. I guess it does compile code with GHC, since it runs main instead of giving you a REPL prompt (right?). There should be some place in its options where you can add -O2 if that's not already on. From the command line instead I would use stack ghc -- File.hs -O2
@Terobero See the last edit. Without errorPercent, and with -O2 on, the program with 10^9 here runs in 12s. Your 20 minutes is way too much.
@Terobero For integer variables, it's faster and safer to use Int or some other integral type. Using Double can cause subtle issues since e.g. x == x+1 when x is a large enough Double, and rounding errors can easily break the program. I'd tend to avoid floating point except when an inexact approximate value is good enough (e.g. not when we need precise comparisons such as if x == 1 then ...). luqui in his answer shows why using Double does not work here.
|

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.