Sorry, I'm a little late to this thread, but the topic of
sieve [] = []
sieve (x:xs) = x : sieve [y | y <- xs, y `mod` x /= 0]
(as posted by Rafael Almeida) is somewhat dear to my heart, as I
wrote a paper about it last summer and sent it in to JFP. Still
waiting for a reply though.
Let's go back to the beginning with the classic complaint and the
excuses...
Creighton Hogg wrote:
So a friend and I were thinking about making code faster in
Haskell, and I was wondering if there was a way to improve the
following method of generating the list of all prime numbers. It
takes about 13 seconds to run, meanwhile my friend's C version took
0.1.
Here come the excuses, like this one from apfelmus,
While Haskell makes it possible to express very complicated
algorithms in simple and elegant ways, you have to expect to pay a
constant factor (roughly 2x-10x) when competing against the same
algorithm in low-level C.
and this one from Nicolas Frisby,
I have yet to need my Haskell to perform well
Matthew Brecknell came up with something much faster, namely
primes :: [Int]
primes = 2 : filter isPrime [3,5..] where
f x p r = x < p*p || mod x p /= 0 && r
isPrime x = foldr (f x) True primes
FWIW, another way to write this code (without needing to think about
how "fold bails early") is
primes = 2: [x | x <−[3,5..], all (\p −> x `mod` p > 0)
(factorsToTry x)]
where
factorsToTry x = takeWhile (\p −> p*p <= x) primes
Both of these algorithms best the "sieve" we began with and run
quickly, but as you can see (possibly more clearly from my
rephrasing), this algorithm is not actually the Sieve of Eratosthenes
-- it's actually a classic naive primes algorithm which checks a
number for primality by trying to divide it by every prime up to its
square root.
But that's okay, because our initial algorithm ISN'T THE REAL SIEVE
EITHER. Markus Fischmann hits the nail on the head when he says
The characteristics of a sieve is, that it uses the already found
primes to generate a list of non-primes that is then removed from a
list of candiates for primeness.
But then we get distracted by a discussion about avoiding division.
It's true that the real sieve avoids division, but it is NOT true
that every algorithm that avoids division is the sieve. The thread
ends with this algorithm from Yitzchak Gale:
-- Delete the elements of the first list from the second list,
-- where both lists are assumed sorted and without repetition.
deleteOrd :: Ord a => [a] -> [a] -> [a]
deleteOrd xs@(x:xs') ys@(y:ys')
| x > y = y : deleteOrd xs ys'
| x < y = deleteOrd xs' ys
| otherwise = deleteOrd xs' ys'
deleteOrd _ ys = ys
sieve (x:xs) = x : sieve (deleteOrd [x+x,x+x+x..] xs)
sieve _ = []
primes = sieve [2..]
Which seems reasonable, until you realize that every prime p we come
up with is still "considered" by k deleteOrd filters, where k is the
number of primes that preceeded it.
So, let's recap: the original algorithm is beautiful and simple, but
it is NOT the actual Sieve of Eratosthenes, NOT because it uses
modulus, but because fundamentally, at the highest level, it is a
different algorithm. At the risk of beating a dead horse, let's see
why it's not the real sieve.
What makes the sieve an efficient algorithm are the details of *what*
gets "crossed off", *when*, and *how many times*. Suppose that we are
finding the first 100 primes (i.e., 2 through 541), and have just
discovered that 17 is prime. We will begin crossing off at 289 (i.e.,
17 * 17) and cross off the multiples 289, 306, 323, ... , 510, 527,
making fifteen crossings off in total. Notice that we cross off 306
(17 * 18), even though it is a multiple of both 2 and 3 and has thus
already been crossed off twice. (The starting point of 17 * 17 is a
pleasing, but actually *minor*, optimization for the *genuine* sieve.)
The algorithm is efficient because each composite number, c, gets
crossed off f times, where f is the number of unique factors of c
less than sqrt(c). The average value for f increases slowly, being
less than 3 for the first 10^12 composites, and less than 4 for the
first 10^34.
None of the algorithms we've discussed correspond to the above
algorithm. It is not merely that they don't perform "optimizations",
such as starting at the square of the prime, or that some of then use
a divisibility check rather than using a simple increment. Rather, at
a fundamental level, they all work differently than the real sieve.
Following our earlier example, after finding that 17 is prime, the
"phony" algorithm will check to see if 19 is divisible by 17 (in the
case of Yitzchak's algorithm, divisibility is checked by comparing
against 17*2), followed by 23, 29, 31, ... , 523, 527, checking a
total of ninety-nine numbers for divisibility by 17. In fact, even if
it did (somehow) begin at 289, it would examine all forty-five
numbers that are not multiples of the primes prior to 17 (289, 293,
307, ... , 523, 527). It is also worth realizing that the set of
numbers examined by the phony algorithm is not a superset of the ones
examined by the real sieve (as we saw, the real algorithm crossed off
306, which the phony algorithm skips). In general, the speed of the
phony algorithm depends on the number of primes that are *not*
factors of each composite, whereas the speed of the real algorithm
depends on the number of (unique) primes that *are*.
Hopefully at this point I've piqued your interest. Maybe you're
wondering what our original sieve algorithm is, if it isn't the real
sieve and the mod issue is just a red herring...? Maybe you're
wondering what the the real sieve looks like, coded in Haskell. Or
maybe you're still not convinced about any of it. If so, I invite
you to read a draft of my JFP paper, available at:
http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf
Best Regards,
Melissa.
_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
