Sieve Of Euler

February 25, 2011

Here’s our solution:

(define (euler n)
    (let loop ((xs (range 3 n 2)) (ps '(2)))
      (let ((p (car xs)))
        (if (< n (* p p))
            (append (reverse ps) xs)
            (loop (minus (cdr xs) (sub-list n p xs))
                  (cons p ps))))))

Here xs is the first list, (sub-list n p xs) is the new list, and the subtraction is done by the minus function:

(define (minus xs ys)
  (let loop ((xs xs) (ys ys) (zs '()))
    (cond ((null? ys) (append (reverse zs) xs))
          ((equal? (car xs) (car ys))
            (loop (cdr xs) (cdr ys) zs))
          (else (loop (cdr xs) ys (cons (car xs) zs))))))

We use auxiliary functions to compute the new list to be subtracted and to curry the multiply function:

(define (sub-list n p xs)
  (let loop ((xs xs) (zs '()))
    (let ((px (* p (car xs))))
      (if (< n px) (reverse zs)
        (loop (cdr xs) (cons px zs))))))

(define (times p) (lambda (x) (* p x)))

The three optimizations are included: We consider only odd numbers by specifying (range 3 n 2) and starting with (2) in the ps result list. We stop early when (< n (* p p)), appending the remaining list, which must necessarily contain only primes, to the primes already found. And we stop computing the new list when the product (* p x) exceeds n.

Here’s an example:

> (euler 100)
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)

Though Euler’s sieve touches each number only once, it is slower than Eratosthenes’ sieve because it uses more operations and because so many intermediate lists must be created, then destroyed. The codepad solution shows a timing test in which Euler’s sieve was about twenty times slower than Eratosthenes’ sieve.

We used range from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/dcEnoKkS.

Pages: 1 2

16 Responses to “Sieve Of Euler”

  1. 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
    
  2. […] today’s Programming Praxis exercise, our goal is to implement the Sieve of Euler to generate primes. […]

  3. Graham said

    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 ps
    

    Strictly 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.

  4. Hannes Tydén said

    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]
    
  5. slabounty said

    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)
    
  6. Rainer said

    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.
    */

  7. programmingpraxis said

    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?

  8. Rainer said

    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

  9. […] crible d’Euler, à la différence du crible d’Ératosthène, n’élimine les nombres composites […]

  10. 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.

  11. Mike said

    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)
    
  12. Khanh Nguyen said

    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 30
    
  13. C++ 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

  14. Rohan said

    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) ?

Leave a comment