4
\$\begingroup\$

I'm new to Haskell, and am trying to learn the basics with a simple function that will calculate the list of primes up to a given value (using the Sieve of Eratosthenes). Here is my code so far, which has the problem of becoming very slow for large inputs nmax:

primes :: Integer -> [Integer]
primes nmax
  | nmax <= 2 = [2]
  | otherwise = sieve (2:[3,5..nmax])
  where
    sieve :: [Integer] -> [Integer]
    sieve [] = []
    sieve [a] = [a]
    sieve (p:xs) = p : sieve [multiple*p+offset | multiple <- [1..(last xs) `quot` p], offset <-[1..(p-1)], multiple*p+offset `elem` xs]

I am aware of this resource which gives implementations of Eratosthenes' Sieve, however, their methods are beyond me as I'm only just beginning. I'm also aware of this this question, but they have used entirely different methods again. I would like to know what it is about my code that is so inefficient.

Example timings:

  • primes 100 takes 0.01s and uses 0 bytes of RAM.
  • primes 1000 takes 0.41s and uses ~38.5MB of RAM.
  • primes 2000 takes 3.15s and uses ~153MB of RAM.

I suspect the problem lies in checking if multiple*p+offset is an element of the list xs, but I'm not sure how to circumvent this check.

Improve this question
\$\endgroup\$
5
  • \$\begingroup\$ You could avoid that check by adding a function that intersects the sorted lists [multiple*p+offset | multiple <- [1..quot (last xs) p], offset <- [1..(p-1)] and xs in linear time of their length. For a constant factor, perhaps use nmax instead of last xs? \$\endgroup\$
    – Gurkenglas
    Sep 10, 2016 at 15:42
  • \$\begingroup\$ @Gurkenglas - Thanks for your suggestions. I will implement the latter. As for the former, isn't checking if multiple*p+offset is in xs already in linear time? \$\endgroup\$
    – Myridium
    Sep 10, 2016 at 15:49
  • 1
    \$\begingroup\$ @Gurkenglas why not write that as an answer? Seems good to me \$\endgroup\$
    – Vogel612
    Sep 10, 2016 at 15:52
  • \$\begingroup\$ You do that check linearly often, intersecting is linear only once. @Vogel612, my comment seems too insignificant to me to raise to answer status, you can do it :P \$\endgroup\$
    – Gurkenglas
    Sep 10, 2016 at 17:15
  • \$\begingroup\$ @Gurkenglas - Well I implemented a sorted list intersection function, but it is still dreadfully slow. I must be doing something else horribly wrong, because a Java program I made can generate primes up to 1 million in milliseconds. \$\endgroup\$
    – Myridium
    Sep 10, 2016 at 17:22

1 Answer 1

Reset to default
7
\$\begingroup\$

The underlying complexity

I suspect the problem lies …

Ah, but we don't like to suspect. Just as our compiler tells us that we're missing a semicolon, our profiler will tell us where we lose time. Let's add a main first:

main :: IO ()
main = readLn >>= print . length . primes

Now, it's easier to measure:

$ ghc -prof -auto-all -O2 Primes.hs
$ echo 1000 | ./Primes +RTS -s -p
168
      23,717,912 bytes allocated in the heap
       1,235,856 bytes copied during GC
         111,776 bytes maximum residency (2 sample(s))
          26,880 bytes maximum slop
               1 MB total memory in use (0 MB lost due to fragmentation)

                                     Tot time (elapsed)  Avg pause  Max pause
  Gen  0        44 colls,     0 par    0.000s   0.002s     0.0000s    0.0002s
  Gen  1         2 colls,     0 par    0.000s   0.001s     0.0005s    0.0008s

  INIT    time    0.000s  (  0.002s elapsed)
  MUT     time    0.375s  (  0.372s elapsed)
  GC      time    0.000s  (  0.003s elapsed)
  RP      time    0.000s  (  0.000s elapsed)
  PROF    time    0.000s  (  0.000s elapsed)
  EXIT    time    0.000s  (  0.000s elapsed)
  Total   time    0.375s  (  0.377s elapsed)

  %GC     time       0.0%  (0.8% elapsed)

  Alloc rate    63,247,765 bytes per MUT second

  Productivity 100.0% of total user, 99.5% of total elapsed

$ echo 2000 | ./Primes +RTS -s -p
303
      84,381,928 bytes allocated in the heap
       9,612,184 bytes copied during GC
         137,048 bytes maximum residency (2 sample(s))
          28,592 bytes maximum slop
               1 MB total memory in use (0 MB lost due to fragmentation)

                                     Tot time (elapsed)  Avg pause  Max pause
  Gen  0       161 colls,     0 par    0.062s   0.015s     0.0001s    0.0005s
  Gen  1         2 colls,     0 par    0.000s   0.001s     0.0003s    0.0004s

  INIT    time    0.000s  (  0.001s elapsed)
  MUT     time    2.297s  (  2.399s elapsed)
  GC      time    0.062s  (  0.015s elapsed)
  RP      time    0.000s  (  0.000s elapsed)
  PROF    time    0.000s  (  0.000s elapsed)
  EXIT    time    0.000s  (  0.000s elapsed)
  Total   time    2.359s  (  2.416s elapsed)

  %GC     time       2.6%  (0.6% elapsed)

  Alloc rate    36,737,710 bytes per MUT second

  Productivity  97.4% of total user, 95.1% of total elapsed

Before we have a look at the profile, we already see that something is amiss. For twice the maximum, your code takes almost 7 times the time: 0.375s for numbers up to 1000, and 2.297 seconds for 2000. That indicates that your algorithm doesn't have the right complexity to begin with. Let's check 4000:

$ echo 4000 | ./Primes +RTS -s -p
550
     302,844,912 bytes allocated in the heap
      70,455,392 bytes copied during GC
         431,144 bytes maximum residency (37 sample(s))
          36,976 bytes maximum slop
               2 MB total memory in use (0 MB lost due to fragmentation)

                                     Tot time (elapsed)  Avg pause  Max pause
  Gen  0       546 colls,     0 par    0.016s   0.093s     0.0002s    0.0006s
  Gen  1        37 colls,     0 par    0.031s   0.012s     0.0003s    0.0013s

  INIT    time    0.000s  (  0.001s elapsed)
  MUT     time   16.766s  ( 17.134s elapsed)
  GC      time    0.047s  (  0.105s elapsed)
  RP      time    0.000s  (  0.000s elapsed)
  PROF    time    0.000s  (  0.000s elapsed)
  EXIT    time    0.000s  (  0.000s elapsed)
  Total   time   16.812s  ( 17.240s elapsed)

  %GC     time       0.3%  (0.6% elapsed)

  Alloc rate    18,063,443 bytes per MUT second

  Productivity  99.7% of total user, 97.2% of total elapsed

For twice the elements, we have 8 times the time. This indicates that you have \$\mathcal O(n^3)\$ complexity.

An unhelpful profile

Let's have a look at the generated profile. In case you don't know GHC's profiling feature yet: -prof enables profiling during compilation, but you still need to tell GHC what to profile (in your code). -auto-all enables profiling of your functions without changing your code. And +RTS -p enables profiling at runtime.

This will create a .prof file. However, it's not very helpful, since most of your logic is in a list comprehension, and that's not profiled:

Sun Sep 11 08:24 2016 Time and Allocation Profiling Report  (Final)

       Primes +RTS -s -p -RTS

    total time  =       15.19 secs   (15193 ticks @ 1000 us, 1 processor)
    total alloc = 185,477,264 bytes  (excludes profiling overheads)

COST CENTRE  MODULE  %time %alloc

primes.sieve Main    100.0   99.9


                                                             individual     inherited
COST CENTRE     MODULE                     no.     entries  %time %alloc   %time %alloc

MAIN            MAIN                        44           0    0.0    0.0   100.0  100.0
 main           Main                        89           0    0.0    0.0   100.0  100.0
  primes        Main                        91           1    0.0    0.1   100.0  100.0
   primes.sieve Main                        92         550  100.0   99.9   100.0   99.9
 CAF            GHC.Read                    76           0    0.0    0.0     0.0    0.0
 CAF            GHC.IO.Encoding.CodePage    73           0    0.0    0.0     0.0    0.0
 CAF            GHC.IO.Encoding             68           0    0.0    0.0     0.0    0.0
 CAF            Text.Read.Lex               61           0    0.0    0.0     0.0    0.0
 CAF            GHC.IO.Handle.Text          57           0    0.0    0.0     0.0    0.0
 CAF            GHC.IO.Handle.FD            56           0    0.0    0.0     0.0    0.0
 CAF:main1      Main                        52           0    0.0    0.0     0.0    0.0
  main          Main                        88           1    0.0    0.0     0.0    0.0
 CAF:lvl5_r3TO  Main                        51           0    0.0    0.0     0.0    0.0
  main          Main                        90           0    0.0    0.0     0.0    0.0

We see that most of the work is done in primes.sieve. That's nothing new, though. We would like to see where the additional work happens. Remember that I said that we usually have to add profiling annotations (also called cost centers)? Now it's time to add them:

    sieve (p:xs) = p : sieve [multiple*p+offset 
                             | multiple <- {-# SCC multiple #-} [1..(last xs) `quot` p]
                             , offset <- {-# SCC offset #-} [1..(p-1)]
                             , {-# SCC elem #-} multiple*p+offset `elem` xs]

Now let's see what takes all that time (and space):

    Sun Sep 11 08:35 2016 Time and Allocation Profiling Report  (Final)

       Primes +RTS -s -p -RTS

    total time  =       14.27 secs   (14270 ticks @ 1000 us, 1 processor)
    total alloc = 186,588,744 bytes  (excludes profiling overheads)

COST CENTRE  MODULE  %time %alloc

elem         Main     98.9   46.7
primes.sieve Main      0.7    9.4
offset       Main      0.4   43.3


                                                             individual     inherited
COST CENTRE     MODULE                     no.     entries  %time %alloc   %time %alloc

MAIN            MAIN                        44           0    0.0    0.0   100.0  100.0
 main           Main                        89           0    0.0    0.0   100.0  100.0
  primes        Main                        91           1    0.0    0.1   100.0  100.0
   primes.sieve Main                        92         550    0.7    9.4   100.0   99.9
    elem        Main                        95     1813869   98.9   46.7    98.9   46.7
    offset      Main                        94        9263    0.4   43.3     0.4   43.3
    multiple    Main                        93         549    0.0    0.4     0.0    0.4
 CAF            GHC.Read                    76           0    0.0    0.0     0.0    0.0
 CAF            GHC.IO.Encoding.CodePage    73           0    0.0    0.0     0.0    0.0
 CAF            GHC.IO.Encoding             68           0    0.0    0.0     0.0    0.0
 CAF            Text.Read.Lex               61           0    0.0    0.0     0.0    0.0
 CAF            GHC.IO.Handle.Text          57           0    0.0    0.0     0.0    0.0
 CAF            GHC.IO.Handle.FD            56           0    0.0    0.0     0.0    0.0
 CAF:main1      Main                        52           0    0.0    0.0     0.0    0.0
  main          Main                        88           1    0.0    0.0     0.0    0.0
 CAF:lvl5_r3UH  Main                        51           0    0.0    0.0     0.0    0.0
  main          Main                        90           0    0.0    0.0     0.0    0.0

Let's focus on the new centers:

    elem        Main                        95     1813869   98.9   46.7    98.9   46.7
    offset      Main                        94        9263    0.4   43.3     0.4   43.3
    multiple    Main                        93         549    0.0

Wow. For 4000 elements, we have 1813869 elem checks. That's almost \$4000^{1.8}\$. Each of those elem checks is linear. While the size of the list decreases, at worst case it has size \$\mathcal O(n) \$. We end up with \$\mathcal O(n^{2.8})\$ complexity from our empirical point of view, which is really close to our real behaviour.

Prime time for efficiency

Now that we now that the elem checks destroy our runtime, and the creation of offset eats the other half of your RAM, let's get back to the drawing board and remember what Eratosthenes' Sieve is supposed to do:

  1. You write down the numbers from 2 to \$n\$.
  2. You start with the next number that you haven't looked at yet
  3. If it's crossed out, take the next one
  4. If it's not crossed out, cross all it's multiple out
  5. Go to step 2, until you've run out of numbers

However, here's what you do:

  1. You write down the numbers from 2 to \$n\$ (only the odds).
  2. You start with the next number that you haven't looked at and that isn't crossed out yet
  3. Put "maybe prime?" markers on all numbers that aren't multiples of it
  4. Cross out all multiples of your number (all that don't have a marker; that step is implicit)
  5. Go to step 2, until you've run out of numbers

The difference in step 3 is huge. For any number \$p\$, there are at most \$\lfloor \frac{n}{p} \rfloor\$ multiples in the range \$[1,n]\$. How many aren't? Well, \$n - \lfloor \frac{n}{p} \rfloor\$. And since \$p > 2\$, we will always check too many elements.

Now, let's change your sieve so that it actually implements the first algorithm:

    sieve :: [Integer] -> [Integer]
    sieve [] = []
    sieve [a] = [a]
    sieve (p:xs) = p : sieve [prime | prime <- xs
                                    , prime `rem` p /= 0]

Before I post the result, let us verify that this really does the same as yours: for a single \$p\$, you created a list of all non-divisible numbers that were also in xs. Therefore, you ended up with:

$$ \{x \in xs \} \cap \{n \in \mathbb N : n \neq 0 \mod p \}$$

Which is almost literally in the code above. Now for the result:

$ echo 4000 | ./Primes +RTS -s -p
2000
     339,984,336 bytes allocated in the heap
      24,834,416 bytes copied during GC
          89,728 bytes maximum residency (2 sample(s))
          27,736 bytes maximum slop
               2 MB total memory in use (0 MB lost due to fragmentation)

                                     Tot time (elapsed)  Avg pause  Max pause
  Gen  0       459 colls,     0 par    0.047s   0.040s     0.0001s    0.0003s
  Gen  1         2 colls,     0 par    0.000s   0.001s     0.0004s    0.0007s

  INIT    time    0.000s  (  0.001s elapsed)
  MUT     time    0.250s  (  0.290s elapsed)
  GC      time    0.047s  (  0.041s elapsed)
  RP      time    0.000s  (  0.000s elapsed)
  PROF    time    0.000s  (  0.000s elapsed)
  EXIT    time    0.000s  (  0.000s elapsed)
  Total   time    0.297s  (  0.332s elapsed)

  %GC     time      15.8%  (12.2% elapsed)

  Alloc rate    1,359,937,344 bytes per MUT second

  Productivity  84.2% of total user, 75.2% of total elapsed

Still too much garbage collection, but that's fine for a start.

Additional discussion

Try to avoid last. Instead, you could have used 

multiples <- [i * p | i <- [1..nmax`quot`p]]

or similar. last has to traverse all the list.

Next, it's great that you've added type signatures, but make sure that you can actually need Integer. On your platform, Int will probably go up to \$2^{63}-1\$, which is more than enough. It will also make your program a lot faster, since Integer has an overhead.

And last but not least, make the code slightly more readable by adding whitespace. If you get a vertical scrollbar on CR, it's usually a sign that your lines are too long.

Oh, and by the way: this variant will still be too slow compared to Java, but you're probably using another data structure with random access abilities in Java, whereas we're creating a "new" list in Haskell all the time.

TL;DR

Verify that you actually implemented the algorithm that you want to implement.

Improve this answer
\$\endgroup\$
4
  • \$\begingroup\$ Thank you; this is tremendously helpful. You know what they say about teaching a man to fish... I'm concerned by your code prime 'quot' p /= 0 however. This resource stresses that this is in fact not the Sieve of Eratosthenes and is needlessly inefficient, so that is why I avoided taking this approach. As for the Java: I am comparing with the timing of a Java program that creates an List of primes from scratch. I don't understand why Haskell cannot achieve the same performance. \$\endgroup\$
    – Myridium
    Sep 11, 2016 at 8:38
  • \$\begingroup\$ I figured by quot you meant mod. \$\endgroup\$
    – Myridium
    Sep 11, 2016 at 8:43
  • \$\begingroup\$ Just thought I would add to this answer that it's a good idea to enable profiling by default when installing packages via cabal. You can find out here. \$\endgroup\$
    – Myridium
    Sep 15, 2016 at 15:30
  • \$\begingroup\$ @Myridium: I recommend you to use stack instead of cabal for your projects (unless you want to build a library). That way you can use stack build --profile (see also) \$\endgroup\$
    – Zeta
    Sep 15, 2016 at 16:15

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.