Performance of Haskell prime sieve

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.

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.