4
\$\begingroup\$

I have this code, which is a pseudo-Sieve of Eratosthenes for generating primes:

module Prime where
  
import Data.List (unfoldr)
nth :: Int -> Maybe Integer
nth n
  | n < 1 = Nothing
  | otherwise = Just $ primes !! (n - 1)

primes :: [Integer]
primes = unfoldr (Just . f) [2..]
  where
    f :: [Integer] -> (Integer, [Integer])
    f cands = (next, filter (\n -> n `mod` next /= 0) (tail cands))
      where
        next = head cands

This code is correct, but it's quite slow - it takes about 6 seconds to generate the 10,001st prime. A simple algorithm, without fancy performance work, that I picked off the Haskell Wiki can find the same prime in a fraction of a second.

As per the limited profiling I've done, it doesn't seem there's any space leaks, but rather genuine algorithm slowness. What is the source of this slowness?

I'm looking for minimal changes that can be simply understood, to target the inefficiency and acheive a decent speedup.

EDIT

@arrowd suggested that the extra work involved in re-traversing the list using !! was wasted. But a version of the algorithm without this is even slower:

nth' :: Int -> Integer
nth' = loop [2..]
  where
    loop :: [Integer] -> Int -> Integer
    loop cands 1 = head cands
    loop cands n = loop (filter (\n' -> n' `mod` head cands /= 0) cands) (n - 1)
\$\endgroup\$
7
  • 4
    \$\begingroup\$ Data.List.(!!) is O(n). \$\endgroup\$
    – arrowd
    Aug 14, 2023 at 7:58
  • \$\begingroup\$ head is O(n) too \$\endgroup\$ Aug 14, 2023 at 11:21
  • \$\begingroup\$ Pardon? hackage.haskell.org/package/base-4.18.0.0/docs/… \$\endgroup\$ Aug 14, 2023 at 11:22
  • \$\begingroup\$ @AriFordsham my bad, was looking at DLists hackage.haskell.org/package/dlist-0.8.0.2/docs/Data-DList.html \$\endgroup\$ Aug 14, 2023 at 11:29
  • \$\begingroup\$ Intuitively, the first thing I’d look at for improving performance in real-world code would be replacing the list with a mutable Data.Vector. If you still want to do it with lists as a learning exercise, I would suggest traversing from back to front, appending each non-multiple from the current list to the head of the updated list with a tail-recursive-modulo-cons function. \$\endgroup\$
    – Davislor
    Aug 14, 2023 at 13:07

2 Answers 2

6
\$\begingroup\$

Looking at your code, one performance bottleneck is the repeated use of mod, which is the slowest integer operation on modern CPUs if it’s in the ISA at all. You also use unfoldr, which can be slower than scanl' or a tail-recursive-modulo-cons function that strictly evaluates its arguments (bang patterns or seq). The !! operator also works in linear rather than constant time, but since you only invoke it once (thanks, benrg), this would not be a large issue.

Intuitively, you’ll probably speed up this algorithm the most by using a mutable Data.Vector instead of a list. However, that’s not really idiomatic Haskell, so you might prefer lists for that reason.

To speed up this sieve using lists, here is one approach you might try: It’s more of a true sieve of Eratosthenes that finds all primes within a given range, using Int rather than Integer. It’s therefore less flexible than what you wrote.

  • Define a tail-recursive helper function sieve :: [Int] -> [Int].
  • Call sieve [2..n]. You could also optimize by calling sieve with only the odd numbers between 3 and n, and prepending 2 to the result.
  • If the argument to sieve is [], there are no more numbers in range to sieve. Return [].
  • Within each call to sieve n (p:potentials), n is the limit of the range, p is definitely prime and potentials is a list of potential primes in range. Splitting the second argument into a head and tail is a constant-time operation.
    • If p*p > n, you have already found all composites in range. All remaining numbers you were given are prime. Return (p:potentials).
    • Otherwise:
      • Generate a list of all multiples of p in range. Do this with a well-behaved list comprehension that uses only addition (much faster than mod) and avoids creating the list as a data structure in memory.
      • Take the set difference of potentials and the multiples, optimizing based on the fact that both are ordered lists with no duplicates. This function should be either tail-recursive-modulo-cons, or tail-recursive with an accumulator.
      • Return p:(sieve setDifference). This makes your function tail-recursive-modulo-cons.

You might convert this algorithm (and many other tail-recursive-modulo-cons algorithms) into a Data.List.scanl'. This will give you code equivalent to what an imperative language would compile a strictly-evaluated while loop that generates list from head to tail into. This is typically much faster than unfoldr.

\$\endgroup\$
1
  • 2
    \$\begingroup\$ As I understand it, there is no "repeated use of !!" in their benchmark. It's used only once to find primes !! 10000. \$\endgroup\$
    – benrg
    Aug 14, 2023 at 22:36
1
\$\begingroup\$

In your code, before a natural number n can appear in primes, it must pass a n `mod` p /= 0 filter-predicate test for each prime p less than it. This is similar to trial division, but slower for two reasons: first, trial division normally stops at √n, since a factor larger than that can't exist unless a smaller factor does, and second, since you are only testing against primes, you have to find all of the smaller primes first, by the same slow algorithm. As a result, where trial division needs O(√n) divisions to test n, you are doing at least O((n/log n)²) divisions to test only n, or O(n/log n) if you don't count the divisions that were used to find the smaller primes.

Arguably, what you're doing isn't a sieve (of Eratosthenes) at all, but a variant of trial division. The purpose of the sieve is to avoid division, and the various list-based implementations of it on the Haskell Wiki do no division. I think they are all still O(n²) up to log factors, but with a smaller constant factor.

If you are using nth only once to find the 10001st prime, then it doesn't contribute significantly to the run time. Using map nth [1..10000] to get a list of primes wouldn't make the algorithm any slower in the big-O sense, but might add noticeable overhead, and in any case should be avoided on general principle. Just use primes directly (take 10000 primes).

Instead of using head and tail, I would pattern match on the list: f (next:cands) = (next, filter (\n -> n `mod` next /= 0) cands).

The test x `mod` y /= 0 is equivalent to x `rem` y /= 0, and the latter is likely to be faster. mod (which computes the remainder from floor division, like Python's %) is better behaved mathematically, but rem (which computes the remainder from truncating division, like C's %) is the operation actually implemented in hardware on x86/x64 (and probably ARM, but I'm not sure). That said, the difference probably isn't noticeable for arbitrary-precision Integer, and you shouldn't be dividing anyway.

\$\endgroup\$

Your Answer

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

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