Sieve of Eratosthenes

From Rosetta Code

Programming Task
This is a programming task. It lays out a problem which Rosetta Code users are encouraged to solve, using languages they know.

Code examples should be formatted along the lines of one of the existing prototypes.

This task has been clarified. Its programming examples are in need of review to ensure that they still fit the requirements of the task.

The Sieve of Eratosthenes is a simple algorithm that finds the prime numbers up to a given integer. Implement this algorithm, with the only allowed optimization that the outer loop can stop at the square root of the limit, and the inner loop may start at the square of the prime just found. That means especially that you shouldn't optimize by using pre-computed wheels, i.e. assume you need only to cross out odd numbers (wheel based on 2), numbers equal to 1 or 5 modulo 6 (wheel based on 2 and 3), or similar wheels based on low primes.

If there's an easy way to add such a wheel based optimization, implement this as an alternative version.

[edit] Ada

This example may be incorrect due to a recent change in the task requirements. Please verify it and remove this message. If the example does not match the requirements, replace this message with Template:incorrect or fix the code yourself.

This solution uses an array of bits representing boolean values to identify the prime numbers. The limit on the outer loop is set to the square root of the number received from the command line. This limit minimizes the number of unneeded iterations in the algorithm.

The lowest index of the array of boolean is set to 2 so that all calculations begin at 2.

with Ada.Text_Io; use Ada.Text_Io;
with Ada.Command_Line; use Ada.Command_Line;
with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;
 
procedure Sieve is 
   package Nat_Io is new Ada.Text_Io.Integer_Io(Natural);
   use Nat_Io;
   type Bool_Array is array(Positive range <>) of Boolean;
   pragma pack (Bool_Array);
 
   procedure Get_Sieve (A : out Bool_Array) is     
      Test_Num : Positive;  
      Limit : Positive := Positive(sqrt(float(A'Last)));
   begin
      A := (Others => True); -- set all values to true 
      for Num in A'First..Limit loop
         if A(Num) then
            Test_Num := Num * Num;
            while Test_Num <= A'Last loop
               A(Test_Num) := False;
               Test_Num := Test_Num + Num;
            end loop;
         end if;
      end loop;
   end Get_Sieve;
 
   pragma Inline(Get_Sieve);
 
   N : Positive := 2;  
begin
   if Argument_Count > 0 then
      N := Positive'Value(Argument(1));
   end if;
   declare
      Vals : Bool_Array(2..N);
   begin
      Get_Sieve(Vals);
      for I in Vals'range loop
         if Vals(I) then
            Put_Line(Positive'Image(I));
         end if;
      end loop;
   end;
end Sieve;

This version uses tasking to try to improve the performance on multiprocessor systems which are becoming common. It is still an implementation of the Sieve of Eratosthenes: the list of numbers is in the main task, and the list of primes is in the list of Siever tasks.

with Ada.Command_Line;
with Ada.Text_IO;
 
procedure Sieve is
   task type Siever is
      entry Check (Value : in Positive);
      entry Print;
   end Siever;
 
   subtype Siever_Task is Siever;
 
   type Sieve_Ptr is access Siever;
 
   task body Siever is
      Number    : Positive;
      Candidate : Positive;
      Next      : Sieve_Ptr;
   begin -- Siever
      accept Check (Value : in Positive) do
         Number := Value;
      end Check;
 
      All_Values : loop
         select
            accept Check (Value : in Positive) do
               Candidate := Value;
            end Check;
 
            if Candidate rem Number /= 0 then
               if Next = null then
                  Next := new Siever_Task;
               end if;
 
               Next.Check (Value => Candidate);
            end if;
         or
            accept Print do
               Ada.Text_IO.Put (Item => Integer'Image (Number) );
 
               if Next /= null then
                  Next.Print;
               end if;
            end Print;
 
            exit All_Values;
         end select;
      end loop All_Values;
   end Siever;
 
   Head : Sieve_Ptr := new Siever;
   Max  : Positive  := 20;
begin -- Sieve
   if Ada.Command_Line.Argument_Count > 0 then
      Max := Integer'Value (Ada.Command_Line.Argument (1) );
   end if;
 
   All_Values : for I in 2 .. Max loop
      Head.Check (Value => I);
   end loop All_Values;
 
   Head.Print;
end Sieve;

This version is the same as the previous version, but with a wheel of 2.

with Ada.Command_Line;
with Ada.Text_IO;
 
procedure Sieve_With_Wheel is
   task type Siever is
      entry Check (Value : in Positive);
      entry Print;
   end Siever;
 
   subtype Siever_Task is Siever;
 
   type Sieve_Ptr is access Siever;
 
   task body Siever is
      Number    : Positive;
      Candidate : Positive;
      Next      : Sieve_Ptr;
   begin -- Siever
      accept Check (Value : in Positive) do
         Number := Value;
      end Check;
 
      All_Values : loop
         select
            accept Check (Value : in Positive) do
               Candidate := Value;
            end Check;
 
            if Candidate rem Number /= 0 then
               if Next = null then
                  Next := new Siever_Task;
               end if;
 
               Next.Check (Value => Candidate);
            end if;
         or
            accept Print do
               Ada.Text_IO.Put (Item => Integer'Image (Number) );
 
               if Next /= null then
                  Next.Print;
               end if;
            end Print;
 
            exit All_Values;
         end select;
      end loop All_Values;
   end Siever;
 
   Head  : Sieve_Ptr := new Siever;
   Max   : Positive  := 20;
   Value : Positive  := 3;
begin -- Sieve_With_Wheel
   if Ada.Command_Line.Argument_Count > 0 then
      Max := Integer'Value (Ada.Command_Line.Argument (1) );
   end if;
 
   All_Values : loop
      exit All_Values when Value > Max;
 
      Head.Check (Value => Value);
      Value := Value + 2;
   end loop All_Values;
 
   Ada.Text_IO.Put (Item => " 2");
   Head.Print;
end Sieve_With_Wheel;

[edit] ALGOL 68

BOOL prime = TRUE, non prime = FALSE;
PROC eratosthenes = (INT n)[]BOOL:
(
  [n]BOOL sieve;
  FOR i TO UPB sieve DO sieve[i] := prime OD;
  INT m = ENTIER sqrt(n);
  sieve[1] := non prime;
  FOR i FROM 2 TO m DO
    IF sieve[i] = prime THEN
      FOR j FROM i*i BY i TO n DO
        sieve[j] := non prime
      OD
    FI
  OD;
  sieve
);

print((eratosthenes(80),new line))

Output:

FTTFTFTFFFTFTFFFTFTFFFTFFFFFTFTFFFFFTFFFTFTFFFTFFFFFTFFFFFTFTFFFFFTFFFTFTFFFFFTF

[edit] BASIC

Works with: FreeBASIC

Works with: RapidQ

 
DIM n AS Integer, k AS Integer, limit AS Integer
 
INPUT "Enter number to search to: "; limit
DIM flags(limit) AS Integer
 
FOR n = 2 TO SQR(limit)
    IF flags(n) = 0 THEN
        FOR k = n*n TO limit STEP n
            flags(k) = 1
        NEXT k
    END IF
NEXT n
 
' Display the primes
FOR n = 2 TO limit
    IF flags(n) = 0 THEN PRINT n; ", ";
NEXT n

[edit] Befunge

2>:3g" "-!v\  g30          <
 |!`"O":+1_:.:03p>03g+:"O"`|
 @               ^  p3\" ":<
2 234567890123456789012345678901234567890123456789012345678901234567890123456789

[edit] C

Plain sieve, without any optimizations:

#include <stdlib.h>
#include <math.h>
 
char* eratosthenes (int n)
{
  char* sieve;
  int i, j, m;
 
  /* calloc initializes to zero */
  sieve = calloc (n+1, sizeof(char));
  m = (int) sqrt ((double) n);
  sieve[0] = 1;
  sieve[1] = 1;
  for (i = 2; i <= m; i++) {
    if (sieve[i] == 0) {
      for (j = i*i; j <= n; j += i) {
        sieve[j] = 1;
      }
    }
  }
  return sieve;
}

Possible optimizations include sieving only odd numbers (or more complex wheels), packing the sieve into bits to improve locality (and allow larger sieves), etc.

[edit] C++

Works with: Boost

// yield all prime numbers less then limit. 
template<class UnaryFunction>
void primesupto(int limit, UnaryFunction yield)
{
  std::vector<bool> is_prime(limit, true);
 
  const int sqrt_limit = static_cast<int>(std::sqrt(limit));
  for (int n = 2; n <= sqrt_limit; ++n)
    if (is_prime[n]) {
	yield(n);
 
	for (unsigned k = n*n, ulim = static_cast<unsigned>(limit); k < ulim; k += n) 
      //NOTE: "unsigned" is used to avoid an overflow in `k+=n` for `limit` near INT_MAX
	  is_prime[k] = false;
    }
 
  for (int n = sqrt_limit + 1; n < limit; ++n)
    if (is_prime[n])
	yield(n);
}

Full program:

/**
   $ g++ -I/path/to/boost sieve.cpp -o sieve && sieve 10000000
 */
#include <inttypes.h> // uintmax_t
#include <limits>
#include <cmath>
#include <iostream>
#include <sstream>
#include <vector>
 
#include <boost/lambda/lambda.hpp>
 
 
int main(int argc, char *argv[])
{
  using namespace std;
  using namespace boost::lambda;
 
  int limit = 10000;
  if (argc == 2) {
    stringstream ss(argv[--argc]);
    ss >> limit;
 
    if (limit < 1 or ss.fail()) {
      cerr << "USAGE:\n  sieve LIMIT\n\nwhere LIMIT in the range [1, " 
	   << numeric_limits<int>::max() << ")" << endl;
      return 2;
    }
  }
 
  // print primes less then 100
  primesupto(100, cout << _1 << " ");
  cout << endl;  
 
  // find number of primes less then limit and their sum
  int count = 0;
  uintmax_t sum = 0;
  primesupto(limit, (var(sum) += _1, var(count) += 1));
 
  cout << "limit sum pi(n)\n" 
       << limit << " " << sum << " " << count << endl;
}

[edit] Common Lisp

(defun">defun primer ()
 ((lambda (&rest q)
  (lambda ()
  ((lambda (u n) (funcall u u n))
   (lambda (u n) (if (some (lambda (x y) (and (/= 1 y) (= (mod x y) 0)))
                           ((lambda (x &aux (y (list">list x))) (setf (cdr y) y)) n) q)
                     (funcall u u (1+ n)) (setf q (cons n q) n n))) (1+ (car q))))) 1))
 
;; Example use:
(defparameter *png* (primer))
(loop repeat 50 collect (funcall *png*))

[edit] E

E's standard library doesn't have a step-by-N numeric range, so we'll define one, implementing the standard iteration protocol.

def rangeFromBelowBy(start, limit, step) {
  return def stepper {
    to iterate(f) {
      var i := start
      while (i < limit) {
        f(null, i)
        i += step
      }
    }
  }
}

The sieve itself:

def eratosthenes(limit :(int > 2), output) {
  def composite := [].asSet().diverge()
  for i ? (!composite.contains(i)) in 2..!limit {
    output(i)
    composite.addAll(rangeFromBelowBy(i ** 2, limit, i))
  }
}

Example usage:

? eratosthenes(12, println)
# stdout: 2
#         3
#         5
#         7
#         11

[edit] Forth

: sieve ( n -- )
  here over erase
  2
  begin
    2dup dup * >
  while
    here over + c@ 0= if
      2dup dup * do
        1 here i + c!
      dup +loop
    then
    1+
  repeat
  drop
  ." Primes: " 2 do here i + c@ 0= if i . then loop ;

100 sieve

[edit] Fortran

Fortran 90 and later

PROGRAM SIEVE
  IMPLICIT NONE
  INTEGER :: i, j, number
  LOGICAL, ALLOCATABLE :: pflags(:)

  WRITE(*,*) "Enter number to search to"
  READ (*,*) number

  ALLOCATE (pflags(2:number))
  pflags = .TRUE.    

  DO i = 2, SQRT(REAL(number))
    IF (pflags(i)) THEN
      DO j = i+i, number, i
        pflags(j) = .FALSE.
      END DO
    END IF
  END DO

  DO i = 2, number
    IF (pflags(i)) WRITE(*,*) i
  END DO

  DEALLOCATE (pflags)

END PROGRAM SIEVE

Since 2 is trivially the only even prime number this version just sieves odd numbers

ALLOCATE (pflags((number-1)/2))
pflags = .TRUE.    
 
DO i = 1, SQRT(REAL(number-1)/2)
  IF (pflags(i)) THEN
    DO j = 3*i+1, (number-1)/2, 2*i+1
      pflags(j) = .FALSE.
    END DO
  END IF
END DO

DO i = 1, (number-1)/2
  IF (pflags(i)) WRITE(*,*) 2*i + 1
END DO

[edit] Haskell

Imperative version, with unboxed arrays:

import Control.Monad
import Control.Monad.ST

import Data.Array.ST
import Data.Array.Unboxed

primeSieve n = runST sieve where
  newSieve :: ST s (STUArray s Integer Bool)
  newSieve = newArray (2,n) True
  sieve :: ST s (UArray Integer Bool)
  sieve = do
    s <- newSieve
    let m = ceiling $ sqrt $ fromInteger n
    forM_ [2..m] $ \i -> do
      si <- readArray s i
      when si $ do
        forM_ [i*i,i*i+i..n] $ \j -> do
        writeArray s j False
    unsafeFreeze s

Return primes from sieve as list:

primesTo n = [p | p <- [2..n], s ! p] where s = primeSieve n

One often sees a purely functional version similar to

primes :: [Integer]
primes = sieve [2..] where
  sieve (p:xs) = p : sieve [x | x <- xs, x `mod` p > 0]

As Melissa O'Neill demonstrates in The Genuine Sieve of Eratosthenes, this version is not better than trial division. She also presents an actual purely functional version that uses priority queues, and is much easier to combine with wheels (pre-calculated sieves for small primes) than the imperative version. Even combined with a simple wheel, the speed of this version is similar to the imperative version.

[edit] J

sieve=: 3 : 0
 m=. <.%:y
 z=. $0
 b=. y{.1
 while. m>j=. 1+b i. 0 do.
  b=. b+.y$(-j){.1
  z=. z,j
 end.
 z,1+I.-.b
)

"Wheels" are readily implemented as follows:

sieve1=: 3 : 0
 m=. <.%:y
 z=. y (>:#]) 2 3 5 7
 b=. 1,}.y$+./(*/z)$&>(-z){.&.>1
 while. m>j=. 1+b i. 0 do.
  b=. b+.y$(-j){.1
  z=. z,j
 end.
 z,1+I.-.b
)

The use of 2 3 5 7 as wheels provides a 20% time improvement for n=1000 and 2% for n=1e6 .

[edit] Java

Works with: Java version 1.5+

import java.util.LinkedList;
 
public class Sieve{
       public static LinkedList<Integer> sieve(int n){
               if(n < 2) return new LinkedList<Integer>();
               LinkedList<Integer> primes = new LinkedList<Integer>();
               LinkedList<Integer> nums = new LinkedList<Integer>();
 
               for(int i = 2;i <= n;i++){ //unoptimized
                       nums.add(i);
               }
 
               while(nums.size() > 0){
                       int nextPrime = nums.remove();
                       for(int i = nextPrime * nextPrime;i <= n;i += nextPrime){
                               nums.removeFirstOccurrence(i);
                       }
                       primes.add(nextPrime);
               }
               return primes;
       }
}

To optimize by testing only odd numbers, replace the loop marked "unoptimized" with these lines:

nums.add(2);
for(int i = 3;i <= n;i += 2){
       nums.add(i);
}

[edit] Logo

to sieve :limit
  make "a (array :limit 2)     ; initialized to empty lists
  make "p []
  for [i 2 :limit] [
    if empty? item :i :a [
      queue "p :i
      for [j [:i * :i] :limit :i] [setitem :j :a :i]
    ]
  ]
  output :p
end
print sieve 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

[edit] Lucid

prime
 where
    prime = 2 fby (n whenever isprime(n));
    n = 3 fby n+2;
    isprime(n) = not(divs) asa divs or prime*prime > N
                    where
                      N is current n;
                      divs = N mod prime eq 0;
                    end;
 end

[edit] recursive

sieve( N )
   where
    N = 2 fby N + 1;
    sieve( i ) =
      i fby sieve ( i whenever i mod first i ne 0 ) ;
   end

[edit] MAXScript

fn eratosthenes n =
(
    multiples = #()
    print 2
    for i in 3 to n do
    (
        if (findItem multiples i) == 0 then
        (
            print i
            for j in (i * i) to n by i do
            (
                append multiples j
            )
        )
    )
)

eratosthenes 100

[edit] Nial

primes is sublist [ each (2 = sum eachright (0 = mod) [pass,count]), pass ] rest count

Using it

|primes 10
=2 3 5 7

[edit] Oberon-2

MODULE Primes;

   IMPORT Out, Math;

   CONST N = 1000;

   VAR a: ARRAY N OF BOOLEAN;
      i, j, m: INTEGER;

BEGIN
   (* Set all elements of a to TRUE. *)
   FOR i := 1 TO N - 1 DO
      a[i] := TRUE;
   END;

   (* Compute square root of N and convert back to INTEGER. *)
   m := ENTIER(Math.Sqrt(N));

   FOR i := 2 TO m DO
      IF a[i] THEN
         FOR j := 2 TO (N - 1) DIV i DO 
            a[i*j] := FALSE;
         END;
      END;
   END;

   (* Print all the elements of a that are TRUE. *)
   FOR i := 2 TO N - 1 DO
      IF a[i] THEN
         Out.Int(i, 4);
      END;
   END;
   Out.Ln;
END Primes.

[edit] OCaml

[edit] Imperative

let sieve n =
  let is_prime = Array.create n true in
  let limit = truncate(sqrt (float n)) in
  for i = 2 to pred limit do
    if is_prime.(i) then
      let j = ref (i*i) in
      try while true do
        is_prime.(!j) <- false;
        j := !j + i;
      done with _ -> ()
  done;
  (is_prime)

let primes n =
  let _, primes =
    Array.fold_left (fun (i,acc) -> function
        | true -> (i+1, i::acc)
        | false -> (i+1, acc))
      (0, [])
      (sieve n)
  in
  List.tl(List.tl(List.rev primes))

in the top-level:

# primes 100 ;;
- : int list =
[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]

[edit] Functional

(* first define some iterators *)
# let fold_iter f init a b =
    let rec aux acc i =
      if i > b
      then (acc)
      else aux (f acc i) (succ i)
    in
    aux init a ;;
val fold_iter : ('a -> int -> 'a) -> 'a -> int -> int -> 'a = <fun>
 
# let fold_step f init a b step =
    let rec aux acc i =
      if i > b
      then (acc)
      else aux (f acc i) (i + step)
    in
    aux init a ;;
val fold_step : ('a -> int -> 'a) -> 'a -> int -> int -> int -> 'a = <fun>
 
(* remove a given value from a list *)
# let remove li v =
    let rec aux acc = function
      | hd::tl when hd = v -> (List.rev_append acc tl)
      | hd::tl -> aux (hd::acc) tl
      | [] -> li
    in
    aux [] li ;;
val remove : 'a list -> 'a -> 'a list = <fun>
 
(* the main function *)
# let primes n =
    let li =
      (* create a list [from 2; ... until n] *)
      List.rev(fold_iter (fun acc i -> (i::acc)) [] 2 n)
    in
    let limit = truncate(sqrt(float n)) in
    fold_iter (fun li i ->
        if List.mem i li  (* test if (i) is prime *)
        then (fold_step (fun acc j -> remove acc j) li (i*i) n i)
        else li)
      li 2 (pred limit)
  ;;
val primes : int -> int list = <fun>
 
# primes 200 ;;
- : int list =
[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; 101; 103; 107; 109; 113; 127; 131; 137; 139; 149; 151;
 157; 163; 167; 173; 179; 181; 191; 193; 197; 199]

[edit] Pascal

Note: Some Pascal implementations put quite low limits on the size of a set (e.g. Turbo Pascal doesn't allow more than 256 members). To compile on such an implementation, reduce the constant PrimeLimit accordingly.

 
program primes(output)
 
const
 PrimeLimit = 1000;
 
var
 primes: set of 1 .. PrimeLimit;
 n, k: integer;
 needcomma: boolean;
 
begin
 { calculate the primes }
 primes := [2 .. PrimeLimit];
 for n := 1 to trunc(sqrt(PrimeLimit)) do
  begin
   if n in primes
    then
     begin
      k := n*n;
      while k < PrimeLimit do
       begin
        primes := primes - [k];
        k := k + n
       end
     end
  end;
 
  { output the primes }
  needcomma := false;
  for n := 1 to PrimeLimit do
   if n in primes
    then
     begin
      if needcomma
       then
        write(', ');
      write(n);
      needcomma := true
     end
end.

[edit] Perl

sub findprimes {
 my $prime_limit=shift;
 my @possible_primes;
 my @primes;
 my $prime;
 if ($prime_limit < 2) {
  return \@primes;
 }
 @possible_primes=(2..$prime_limit);
 do {
  $prime=shift(@possible_primes);
  push @primes, $prime;
  @possible_primes = grep {$_%$prime!=0} @possible_primes;
 } while ($possible_primes[0] < sqrt($prime_limit));	
 return [@primes,@possible_primes];	
}
 
my $prime_limit;
print "Input the number that you wish to test :";
chomp($prime_limit=<STDIN>);
my $primes=&findprimes($prime_limit);
print "The primes are: @{$primes}!!!\n";

[edit] Pop11

define eratostenes(n);
lvars bits = inits(n), i, j;
for i from 2 to n do
   if bits(i) = 0 then
      printf('' >< i, '%s\n');
      for j from 2*i by i to n do
         1 -> bits(j);
      endfor;
   endif;
endfor;
enddefine;

[edit] Python

[edit] Testing on list of primes up to that number

for x in range(1,100000)
    z=0
    for y in p:
        if x%y==0:
            z=1
    if z==0:
        p.append(x)

[edit] Using list lookup

def eratosthenes(n):
    multiples = []
    print 2
    for i in xrange(3, n+1):
        if not i in multiples:
            print i
            for j in xrange(i*i, n+1, i):
                multiples.append(j)
 
eratosthenes(100)

[edit] Using set lookup

The version below uses a set to store the multiples. set objects are much faster (usually O(log n)) than lists (O(n)) for checking if a given object is a member. Using the set.update method avoids explicit iteration in the interpreter, giving a further speed improvement.

def eratosthenes2(n):
    multiples = set()
    print 2
    for i in xrange(3, n+1):
        if not i in multiples:
            print i
            multiples.update(xrange(i*i, n+1, i))
 
eratosthenes2(100)

[edit] Using array lookup

The version below uses array lookup to test for primality. The function primes_upto() is a straightforward implementation of Sieve of Eratosthenesalgorithm. It returns prime numbers less than or equal to limit.

def primes_upto(limit):
    is_prime = [False] * 2 + [True] * (limit - 1) 
    for n in xrange(int(limit**0.5 + 1.5)): # stop at ``sqrt(limit)``
        if is_prime[n]:
            for i in xrange(n*n, limit+1, n):
                is_prime[i] = False
    return [i for i, prime in enumerate(is_prime) if prime]

[edit] Using generator

 
def iprimes_upto(limit):
    is_prime = [False] * 2 + [True] * (limit - 1)
    for n in range(limit + 1):
        if is_prime[n]:
            yield n
            for i in range(n*n, limit+1, n): # start at ``n`` squared
                is_prime[i] = False

Example:

 
>>> list(iprimes_upto(15))
[2, 3, 5, 7, 11, 13]

[edit] Using numpy

Library: numpy

Belowed code adapted from literateprograms.org using numpy

import numpy
def primes_upto2(limit):
    is_prime = numpy.ones(limit + 1, dtype=numpy.bool)
    for n in xrange(2, int(limit**0.5 + 1.5)): 
        if is_prime[n]:
            is_prime[n*n::n] = 0
    return numpy.nonzero(is_prime)[0][2:]

Performance note: there is no point to add wheels here, due to execution of p[n*n::n] = 0 and nonzero() takes us almost all time.

[edit] Using wheels

Version with wheel based optimization:

from numpy import array, bool_, multiply, nonzero, ones, put, resize
#
def makepattern(smallprimes):
    pattern = ones(multiply.reduce(smallprimes), dtype=bool_)
    pattern[0] = 0
    for p in smallprimes:
        pattern[p::p] = 0
    return pattern
#
def primes_upto3(limit, smallprimes=(2,3,5,7,11)):    
    sp = array(smallprimes)
    if limit <= sp.max(): return sp[sp <= limit]
    #
    isprime = resize(makepattern(sp), limit + 1) 
    isprime[:2] = 0; put(isprime, sp, 1) 
    #
    for n in range(sp.max() + 2, int(limit**0.5 + 1.5), 2): 
        if isprime[n]:
            isprime[n*n::n] = 0 
    return nonzero(isprime)[0]

Examples:

>>> primes_upto3(10**6, smallprimes=(2,3)) # Wall time: 0.17
array([     2,      3,      5, ..., 999961, 999979, 999983])
>>> primes_upto3(10**7, smallprimes=(2,3))            # Wall time: '''2.13'''
array([      2,       3,       5, ..., 9999971, 9999973, 9999991])
>>> primes_upto3(15)
array([ 2,  3,  5,  7, 11, 13])
>>> primes_upto3(10**7, smallprimes=primes_upto3(15)) # Wall time: '''1.31'''
array([      2,       3,       5, ..., 9999971, 9999973, 9999991])
>>> primes_upto2(10**7)                               # Wall time: '''1.39''' <-- version ''without'' wheels
array([      2,       3,       5, ..., 9999971, 9999973, 9999991])
>>> primes_upto3(10**7)                               # Wall time: '''1.30'''
array([      2,       3,       5, ..., 9999971, 9999973, 9999991])

The above-mentioned examples demonstrate that the given wheel based optimization does not show significant performance gain.

[edit] Infinite generator

A generator that will generate primes indefinitely (until it runs out of memory).

import heapq
 
# generates p*p, p*(p+1), p*(p+2), ...
#  i.e. generates a sequence of non-primes that are multiples of p
# not a Python generator though
# it supports "peeking"
class MultiplesGen:
    def __init__(self, p):
        self.p = p
        self.q = p * p
 
    def get(self):
        return self.q
 
    def step(self):
        self.q += self.p
 
    def __cmp__(self, other):
        # order these sequences by their first not-yet-used element
        return cmp(self.q, other.q)
 
# generates all prime numbers
def sieve():
    # priority queue of above generators for all primes generated so far
    # i.e. priority queue of sequences of non-primes
    # the priority queue allows us to get the "next" non-prime quickly
    nonprimes = []
 
    i = 2
    while True:
        if nonprimes and i == nonprimes[0].get(): # non-prime
            while nonprimes[0].get() == i:
                # for each generator that generates this number
                # have it go to the next number
                # and re-position it in the priority queue
                x = heapq.heappop(nonprimes)
                x.step()
                heapq.heappush(nonprimes, x)
 
        else: # prime
            heapq.heappush(nonprimes, MultiplesGen(i))
            yield i
 
        i += 1