Fermat vs Eratosthenes
Benchmarking the generation of prime numbers using Fermat's Little Theorem and the sieve of Eratosthenes in Haskell.
I was thinking about Fermat’s Little Theorem and I wondered whether this primality test could be used to generate prime numbers faster than the Sieve of Eratosthenes which is typically used to market Haskell to imperative programmers.
Fermat’s Little Theorem says that “if \(p\) is a prime number, then for any integer \(a\), the number \(a^p−a\) is an integer multiple of \(p\).”
This is equivalent to \(a^p \equiv a \ \ (\mod p).\) which can be simplified as: \(a^{p-1} \equiv 1 \ \ (\mod p).\)
Assuming \(a\) to be \(2\), this can be expressed as follows:
To generate a sequence up to an integer max
, we include the above function in fermat
:
fermat :: Integer -> [Integer]
fermat max = 2:filter isPrime [3..max]
where isPrime x = (2^(x-1)) `mod` x == 1
Let us now turn our attention to the sieve of Eratosthenes. Melissa E. O’Neill argues that the algorithm typically cited in introductions to lazy functional programming is not the genuine sieve of Eratosthenes. This applies to Haskell.org’s implementation which we will call falseEratosthenes
:
falseEratosthenes :: Integer -> [Integer]
falseEratosthenes max = filterPrime $ [2..max]
where filterPrime [] = []
filterPrime (p:xs) =
p : filterPrime [x | x <- xs, x `mod` p /= 0]
The filterPrime
function may be understood, using a pseudo-strict evaluation strategy, for integers \(\{2..10\}\), as follows:
filterPrime [2,3,4,5,6,7,8,9,10]
= 2 : filterPrime [ x | x <- [3,4,5,6,7,8,9,10], x `mod` 2 /= 0]
= 2 : filterPrime [3,5,7,9]
= 2 : 3 : filterPrime [ x | x <- [5,7,9], x `mod` 3 p/= 0]
= 2 : 3 : filterPrime [5,7]
= 2 : 3 : 5 : filterPrime [x | x <- [7], x `mod` 5 p /= 0]
= 2 : 3 : 5 : filterPrime [7]
= 2 : 3 : 5 : 7 : filterPrime []
= 2 : 3 : 5 : 7 : []
= [2,3,5,7]
The genuine version, which is algorithmically closer to Eratosthenes’ method, is presented here:
import qualified Data.Set as PQ
genuineEratosthenes :: Integer -> [Integer]
genuineEratosthenes max = 2:sieve [3,5..max]
where
sieve (x:xs) = x : sieve' xs (insertprime x xs PQ.empty)
sieve' [] _ = []
sieve' (x:xs) table
| nextComposite == x = sieve' xs (adjust x table)
| otherwise = x : sieve' xs (insertprime x xs table)
where
(nextComposite,_) = PQ.findMin table
adjust x table
| n == x = adjust x (PQ.insert (n', ns) newPQ)
| otherwise = table
where
Just ((n, n':ns), newPQ) = PQ.minView table
insertprime p xs = PQ.insert (p*p, map (*p) xs)
Finally, we use Criterion to measure which of the three functions is faster at finding primes between 2 and 100,000:
main = defaultMain [
bgroup "primes"
[ bench "Fermat"
$ whnf (length . fermat) 100000
, bench "False Eratosthenes"
$ whnf (length . falseEratosthenes) 100000
, bench "Genuine Eratosthenes"
$ whnf (length . genuineEratosthenes) 100000
]
]
Intuitively, I had expected Fermat’s approach to be faster given than I reasoned that the trade-off was accuracy for speed given that Fermat’s method eventually picks up false positives—that is to say, some of the numbers that satisfy the theorem’s property are composite. The results, though, showed a different picture, especially when considering the genuine version of the sieve.
In tabular form:
Function | Time Taken |
---|---|
Fermat | 10.6s |
False Eratosthenes | 4.96s |
Genuine Eratosthenes | 112ms |
There are, arguably, many optimisations that may be applicable to all versions of these implementations; for instance, the power in Fermat’s function may be implemented using bit shifting which might be, CPU cycle-wise, cheaper. My intention though, was to observe the speed of classic, well-understood algorithms. I haven’t analysed the behaviour of these functions in larger integer ranges in which primes become more scarce either.
The conclusion is that, for a relatively small number of primes, there is no need to shoulder with the pseudoprimes that sometimes the Fermat’s method produces; the genuine sieve of Eratosthenes is faster and produces actual primes.