Sieve Of Euler
February 25, 2011
The ancient Sieve of Eratosthenes that computes the list of prime numbers is inefficient in the sense that some composite numbers are struck out more than once; for instance, 21 is struck out by both 3 and 7. The great Swiss mathematician Leonhard Euler invented a sieve that strikes out each composite number exactly once, at the cost of some additional bookkeeping. It works like this:
First, make a list of numbers from 2, as large as you wish; call the maximum number n.
Second, extract the first number from the list, make a new list in which each element of the original list, including the first, is multiplied by the extracted first number.
Third, “subtract” the new list from the original, keeping in an output list only those numbers in the original list that do not appear in the new list.
Fourth, output the first number from the list, which is prime, and repeat the second, third and fourth steps on the reduced list excluding its first element, continuing until the input list is exhausted.
For example, start with the list 2 3 4 5 … 30. Then the new list is 4 6 8 10 … 60. Subtracting gives the list 2 3 5 7 9 … 29. Now 2 is prime and the process repeats on the list 3 5 7 9 … 29. At the next step, the new list is 9 15 21 27 … 87, subtracting gives the list 3 5 7 11 13 … 29, now 2 and 3 are prime and the process repeats on the list 5 7 11 13 … 29. Likewise for the primes 5 and 7, and since 7 · 7 > 30, the process stops, with the remaining list 11 13 17 19 23 29, so the complete list of primes less than 30 is 2 3 5 7 11 13 17 19 23 29.
Just as in the Sieve of Eratosthenes, you can speed up the Sieve of Euler by considering only odd numbers, by stopping once the first item in the list is greater than the square root of n, and by computing the new list in the second step only as far as n.
Your task is to write a program that implements the Sieve of Euler, then compare its performance to the Sieve of Eratosthenes. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
My Haskell solution (see http://bonsaicode.wordpress.com/2011/02/25/programming-praxis-sieve-of-euler/ for a version with comments):
import Data.List.Ordered import Data.Time.Clock primes_euler :: [Integer] primes_euler = 2 : euler [3,5..] where euler ~(x:xs) = x : euler (minus xs $ map (* x) (x:xs)) primes_eratosthenes :: [Integer] primes_eratosthenes = 2 : 3 : sieve (tail primes_eratosthenes) 3 [] where sieve ~(p:ps) x fs = let q=p*p in foldr (flip minus) [x+2,x+4..q-2] [[y+s,y+2*s..q] | (s,y) <- fs] ++ sieve ps q ((2*p,q):[(s,q-rem (q-y) s) | (s,y) <- fs]) benchPrimes :: Num a => [a] -> IO () benchPrimes f = do start <- getCurrentTime print . sum $ take 10000 f print . flip diffUTCTime start =<< getCurrentTime main :: IO () main = do benchPrimes primes_eratosthenes benchPrimes primes_euler[…] today’s Programming Praxis exercise, our goal is to implement the Sieve of Euler to generate primes. […]
Here’s my implementation of Euler’s Sieve:
def _drop_mults(xs): for i in xrange(0, len(xs), xs[0]): xs[i] = 0 return filter(None, xs) def euler(n): xs = range(3, n, 2) ps = [2] while xs: p = xs[0] ps.append(p) xs = _drop_mults(xs) return psStrictly speaking, it doesn’t build the second list and subtract—the best way I could find to do this was
with sets, but the time complexity was still larger than it neede to be—but it does filter out all multiples
of the first item once and only once. I tried to get a faster version going using iterators (kind of like lazy
lists), but it turned out to be slower. Eratosthenes still wins; see codepad for timing, where I
pit this against the Python Cookbook’s
erat2()implementation.Unoptimized:
# Also available as gist: https://gist.github.com/844494 def sieve_of_euler(list) if list.empty? list else out = list - (list.map { |e| e * list.first }) [ out.first ] + sieve_of_euler(out[1..-1]) end end sieve_of_euler((2..30).to_a) # => [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]In ruby without any speedups …
def euler_sieve(list) prime = [] while list != [] do new = list.map { |i| list[0] * i } sub = list - new prime << sub.shift list = sub end prime end p euler_sieve((2..30).to_a)My try in REXX
n = 100 start = time('E') x = euler_primes(n) z = time('R')-start say 'Euler needed' z 'seconds for', x 'primes in' n 'numbers.' start = time('E') x = eratosthenes_primes(n) z = time('R')-start say 'Eratosthenes needed' z 'seconds for', x 'primes in' n 'numbers.' exit euler_primes: procedure parse arg n ur_liste = '' wk_liste = '' aus = '' do x = 2 to n ur_liste = ur_liste x end do while length(ur_liste) > 0 fak = word(ur_liste, 1) wk_liste = multipliziere(ur_liste, fak) ur_liste = entferne(wk_liste, ur_liste) parse value ur_liste with first ur_liste aus = aus first end return words(aus) multipliziere: procedure parse arg liste,faktor neu_liste = '' do while length(liste) > 0 parse value liste with wort liste neu_liste = neu_liste (wort * faktor) end return strip(neu_liste) entferne: procedure parse arg liste2,liste1 do x = 1 to words(liste1) wort = word(liste1, x) if wordpos(wort, liste2) > 0 then liste1 = delword(liste1, x, 1) end return liste1 eratosthenes_primes: procedure parse arg n ur_liste = '' do x = 2 to n ur_liste = ur_liste x end ur_liste = delete_multiples(ur_liste, n) return words(ur_liste) delete_multiples: procedure parse arg liste,n z = 2 do while (z*z) <= n do y = (z+1) to n if y // z == 0 then do p = wordpos(y,liste) if p > 0 then liste = delword(liste, p, 1) end end z = z + 1 end return liste/*
Euler needed 0.010000 seconds for 25 primes in 100 numbers.
Eratosthenes needed 0.006000 seconds for 25 primes in 100 numbers.
*/
Here’s my implementation in C++ (see repository at https://bitbucket.org/mekajfire/prime-number-generators).
Rainer: I used REXX back in my OS/2 days, probably twenty years ago. I wrote a critique of REXX at http://groups.google.com/group/comp.lang.rexx/msg/7584cb0c6d4d2981. Is the language still active?
Hello Phil,
I used it professionally on the AS/400 and like to think of it as a great language for teaching. I will write you a private mail.
Best regards
Rainer
My Python Solution
[…] crible d’Euler, à la différence du crible d’Ératosthène, n’élimine les nombres composites […]
Without any optimizations, the Sieve of Eratosthenes is incredibly superior. Here’s the output of my program when obtaining the list of primes below 15,000:
Sieve of Eratosthenes ran in 0 s
Sieve of Euler ran in 4 s
Here’s the Java code for my implementation of the Sieve of Euler:
public List<Integer> getEulerSieve(int maximum){ List<Integer> sieve = generateListFor(maximum); List<Integer> primes = new ArrayList<Integer>(); while(sieve.size() > 0){ List<Integer> temp = new ArrayList<Integer>(); for(int i = 0; i < sieve.size(); i++) temp.add(sieve.get(i) * sieve.get(0)); for(Integer composite : temp) sieve.remove(composite); if(sieve.size() > 0){ primes.add(sieve.get(0)); sieve.remove(0); } } return primes; }And here’s my crude implementation of the Sieve of Eratosthenes:
public List<Integer> getEratosthenesSieve(int maximum) { List<Integer> primes = new ArrayList<Integer>(); int[] sieve = generateArrayFor(maximum); int[] tempSieve = Arrays.copyOf(sieve, maximum - 1); for(int i = 0; i < sieve.length - 1; i++){ for(int j = i + tempSieve[i]; j < maximum - 1; j += tempSieve[i]){ if(j > maximum - 1 || sieve[j] == Integer.MIN_VALUE) continue; sieve[j] = Integer.MIN_VALUE; } } for (Integer prime : sieve) { if(prime != Integer.MIN_VALUE) primes.add(prime); } return primes; }The complete implementation can be found here.
After spending way too much time on this task, my fastest version of the Sieve of Eratsthenes is consistently 2 to 3 times faster than my fastest Sieve of Euler.
On my HP 2GHz core2 duo laptop:
erat14(1e7) took about 1.2 seconds,
euler18(1e7) took about 3.4 seconds.
For reference, erat2(1e7) from the Python Cookbook took about 7.4 seconds.
def euler18(n): ps = range(1,n,2) ps[0] = 0 root_n = int(sqrt(n)) limit = (root_n - 1)/2 + 1 for p in ifilter(None, islice(ps,limit)): for q in ifilter(None, ps[((n-1)/p - 1)/2: (p-1)/2 - 1: -1]): ps[(p*q - 1)/2] = 0 return [2] + filter(None, ps) def erat14(n): primes = range(1,n,2) primes[0] = 0 zeros = [0]*(len(primes)) limit = int(sqrt(n)) for p in takewhile(lambda x: x <= limit, ifilter(None,primes)): start = (p*p-1)/2 primes[start::p] = zeros[start::p] primes[0] = 2 return filter(None,primes)Mine in F#
let euler_sieves (n:int) = let L = 2 :: [3 .. 2 .. n] let rec sieves L = match L with | [] -> [] | h::t when h * h > n -> L | h::t -> h :: (t |> List.filter (fun x -> x % h <> 0) |> sieves)\ sieves L euler_sieves 30C++ implementation. Made sure to iterate over each element exactly once, which is the crux of this algorithm.
To achieve this I had to store the ‘non-deleted’ previous and next number of each number, at the end of each iteration. So, that in next time i can go the correct next in constant time.
https://github.com/puneetmnit/playground/blob/master/Algorithms/sieve_of_euler.cpp
Regarding the related problem of the number of primes less than N;
I wonder whether the fact that the Euler sieve only excludes numbers once;
can provide a basis for a general expression for the number of compound numbers less than N ?;
eg f(N) so that the number of primes less than N can be expressed as = N- f(N) ?