Sieve of Eratosthenes: Difference between revisions

Content added Content deleted

Inline

Revision as of 13:06, 18 February 2012

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. don't 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.

Note

It is important that the sieve algorithm be the actual algorithm used to find prime numbers for the task.

Cf

Ada

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.

<lang GLBasic>// Sieve of Eratosthenes (find primes) // GLBasic implementation

GLOBAL n%, k%, limit%, flags%[]

limit = 100 // search primes up to this number

DIM flags[limit+1] // GLBasic arrays start at 0

FOR n = 2 TO SQR(limit)

   IF flags[n] = 0
       FOR k = n*n TO limit STEP n
           flags[k] = 1
       NEXT
   ENDIF

   IF flags[n] = 0 THEN STDOUT n + ", "

ALGOL 68

<lang algol68>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))</lang>

Output:

FTTFTFTFFFTFTFFFTFTFFFTFFFFFTFTFFFFFTFFFTFTFFFTFFFFFTFFFFFTFTFFFFFTFFFTFTFFFFFTF

AutoHotkey

Search autohotkey.com: of Eratosthenes
Source: AutoHotkey forum by Laszlo <lang autohotkey> MsgBox % "12345678901234567890`n" Sieve(20)

Sieve(n) { ; Sieve of Eratosthenes => string of 0|1 chars, 1 at position k: k is prime

  Static zero := 48, one := 49 ; Asc("0"), Asc("1") 
  VarSetCapacity(S,n,one) 
  NumPut(zero,S,0,"char") 
  i := 2 
  Loop % sqrt(n)-1 { 
     If (NumGet(S,i-1,"char") = one) 
        Loop % n//i 
           If (A_Index > 1) 
              NumPut(zero,S,A_Index*i-1,"char") 
     i += 1+(i>2) 
  } 
  Return S

} </lang>

AWK

An initial array holds all numbers 2..max (which is entered on stdin); then all products of integers are deleted from it; the remaining are displayed in the unsorted appearance of a hash table.

$ awk '{for(i=2;i<=$1;i++) a[i]=1;
>       for(i=2;i<=sqrt($1);i++) for(j=2;j<=$1;j++) delete a[i*j];
>       for(i in a) printf i" "}'
100
71 53 17 5 73 37 19 83 47 29 7 67 59 11 97 79 89 31 13 41 23 2 61 43 3

The following variant does not unset non-primes, but set them to 0, to preserve order in output:

$ awk '{for(i=2;i<=$1;i++) a[i]=1;
>       for(i=2;i<=sqrt($1);i++) for(j=2;j<=$1;j++) a[i*j]=0;
>       for(i=2;i<=$1;i++) if(a[i])printf i" "}'
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

BASIC

Works with: FreeBASIC

Works with: RapidQ

<lang freebasic> 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 </lang>

Applesoft BASIC

<lang basic> 10 INPUT "ENTER NUMBER TO SEARCH TO: ";LIMIT 20 DIM FLAGS(LIMIT) 30 FOR N = 2 TO SQR (LIMIT) 40 IF FLAGS(N) < > 0 GOTO 80 50 FOR K = N * N TO LIMIT STEP N 60 FLAGS(K) = 1 70 NEXT K 80 NEXT N 90 REM DISPLAY THE PRIMES 100 FOR N = 2 TO LIMIT 110 IF FLAGS(N) = 0 THEN PRINT N;", "; 120 NEXT N </lang>

Locomotive Basic

<lang locobasic>10 DEFINT a-z 20 INPUT "Limit";limit 30 DIM f(limit) 40 FOR n=2 TO SQR(limit) 50 IF f(n)=1 THEN 90 60 FOR k=n*n TO limit STEP n 70 f(k)=1 80 NEXT k 90 NEXT n 100 FOR n=2 TO limit 110 IF f(n)=0 THEN PRINT n;","; 120 NEXT</lang>

BBC BASIC

<lang bbcbasic> limit% = 100000

     DIM sieve% limit%
     
     prime% = 2
     WHILE prime%^2 < limit%
       FOR I% = prime%*2 TO limit% STEP prime%
         sieve%?I% = 1
       NEXT
       REPEAT prime% += 1 : UNTIL sieve%?prime%=0
     ENDWHILE
     
     REM Display the primes:
     FOR I% = 1 TO limit%
       IF sieve%?I% = 0 PRINT I%;
     NEXT</lang>

bash

See solutions at UNIX Shell.

Befunge

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

C

Plain sieve, without any optimizations:

include <stdlib.h>
include <math.h>

char* eratosthenes (int n, int *c) {

 char* sieve;
 int i, j, m;

 *c = n-1;     /* primes count */

 /* 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) {
       if (sieve[j] == 0) {
         sieve[j] = 1; 
         --(*c);
       }
     }
   }
 }
 return sieve;

} </lang>

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

C++

<lang cpp> // yield all prime numbers less than 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); } </lang>

Full program:

Works with: Boost

<lang cpp>/**

  $ 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;

} </lang>

C#

Works with: C# version 3.0+

<lang csharp>using System; using System.Collections; using System.Collections.Generic; using System.Linq; using System.Text;

namespace sieve {

   class Program
   {
       static void Main(string[] args)
       {
           int maxprime = int.Parse(args[0]);
           ArrayList primelist = sieve(maxprime);

           foreach (int prime in primelist)
               Console.WriteLine(prime);

           Console.WriteLine("Count = " + primelist.Count);
           Console.ReadLine();
       }

       static ArrayList sieve(int arg_max_prime)
       {
           BitArray al = new BitArray(arg_max_prime, true);

           int lastprime = 1;
           int lastprimeSquare = 1;

           while (lastprimeSquare <= arg_max_prime)
           {
               lastprime++;

               while (!(bool)al[lastprime])
                   lastprime++;

               lastprimeSquare = lastprime * lastprime;

               for (int i = lastprimeSquare; i < arg_max_prime; i += lastprime)
                   if (i > 0)
                       al[i] = false;
           }

           ArrayList sieve_2_return = new ArrayList();

           for (int i = 2; i < arg_max_prime; i++)
               if (al[i])
                   sieve_2_return.Add(i);

           return sieve_2_return;
       }
   }

}</lang>

Clojure

<lang lisp>(defn primes-up-to [n]

 (let [n (int n)]
   "Returns a list of all primes from 2 to n"
   (let [root (int (Math/round (Math/floor (Math/sqrt n))))]
     (loop [i (int 3)
            a (int-array n)
            result (transient [2])]
       (if (>= i n)
         (persistent! result)
         (recur (+ i (int 2))
                (if (< i root)
                  (loop [arr a
                         inc (+ i i)
                         j (* i i)]
                    (if (>= j n)
                      arr
                      (recur (do (aset arr j (int 1)) arr)
                             inc
                             (+ j inc))))
                  a)
                (if (zero? (aget a i))
                  (conj! result i)
                  result)))))))</lang>

This code is wrong

(primes-up-to 29) => [2 3 5 7 11 13 17 19 23 25]

Or as a lazy sequence:

<lang lisp>(def prime-gen

    (let [primes (atom [])]
      (for [n (iterate inc 2)
            :when (not-any? #(zero? (rem n %))
                            (filter #(<= % (Math/sqrt n)) 
                                    @primes))]
        (do (swap! primes conj n)
            n))))</lang>

CMake

<lang cmake>function(eratosthenes var limit)

 # Check for integer overflow. With CMake using 32-bit signed integer,
 # this check fails when limit > 46340.
 if(NOT limit EQUAL 0)         # Avoid division by zero.
   math(EXPR i "(${limit} * ${limit}) / ${limit}")
   if(NOT limit EQUAL ${i})
     message(FATAL_ERROR "limit is too large, would cause integer overflow")
   endif()
 endif()

 # Use local variables prime_2, prime_3, ..., prime_${limit} as array.
 # Initialize array to y => yes it is prime.
 foreach(i RANGE 2 ${limit})
   set(prime_${i} y)
 endforeach(i)

 # Gather a list of prime numbers.
 set(list)
 foreach(i RANGE 2 ${limit})
   if(prime_${i})
     # Append this prime to list.
     list(APPEND list ${i})

     # For each multiple of i, set n => no it is not prime.
     # Optimization: start at i squared.
     math(EXPR square "${i} * ${i}")
     if(NOT square GREATER ${limit})   # Avoid fatal error.
       foreach(m RANGE ${square} ${limit} ${i})
         set(prime_${m} n)
       endforeach(m)
     endif()
   endif(prime_${i})
 endforeach(i)
 set(${var} ${list} PARENT_SCOPE)

endfunction(eratosthenes)</lang>

<lang cmake># Print all prime numbers through 100. eratosthenes(primes 100) message(STATUS "${primes}")</lang>

Common Lisp

<lang lisp>(defun sieve-of-eratosthenes (maximum)

 (let ((sieve (make-array (1+ maximum) :element-type 'bit
                                         :initial-element 0)))
   (loop for candidate from 2 to maximum
         when (zerop (bit sieve candidate))
           collect candidate
           and do (loop for composite from (expt candidate 2) 
                                        to maximum by candidate
                         do (setf (bit sieve composite) 1)))))</lang>

Working with odds only (above twice speedup):

<lang lisp>(defun sieve-odds (maximum)

(cons 2 
 (let ((maxi (floor (/ (- maximum 1) 2))))
  (let ((sieve (make-array (1+ maxi) :element-type 'bit 
                                       :initial-element 0)))
   (loop for i from 1 to maxi
         when (zerop (bit sieve i))
           collect (+ (* 2 i) 1)
           and do (loop for j from (* 2 i (+ i 1)) 
                                to maxi by (+ (* 2 i) 1)
                         do (setf (bit sieve j) 1)))))))</lang>

While formally a wheel, odds are uniformly spaced and do not require any special processing except for value translation. Wheels proper aren't uniformly spaced and are thus trickier.

D

<lang d>import std.stdio, std.algorithm, std.range;

uint[] primes(in uint n) /*pure nothrow*/ {

   if (n < 2) return [];

   auto F = new bool[n + 1];
   F[0 .. 2] = true;
   foreach (i; 2 .. cast(uint)(n ^^ 0.5))
       if (!F[i])
           for (uint j = i * i; j <= n; j += i)
               if (!F[j]) F[j] = true;

   return array(filter!(i => !F[i])(iota(n + 1)));

}

void main() {

   writeln(primes(50));

}</lang> Output:

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49]

Faster version

This version uses an array of bits (instead of booleans, that are represented with one byte), and skips even numbers. The output is the same: <lang d>import std.stdio, std.math, std.array;

size_t[] primes(in size_t n) /*pure nothrow*/ {

   enum size_t bpc = size_t.sizeof * 8;
   auto F = new size_t[((n + 2) / 2) / bpc + 1];
   F[] = size_t.max;

   size_t isSet(in size_t i) nothrow {
       immutable size_t offset = i / bpc;
       immutable size_t mask = 1 << (i % bpc);
       return F[offset] & mask;
   }

   void resetBit(in size_t i) nothrow {
       immutable size_t offset = i / bpc;
       immutable size_t mask = 1 << (i % bpc);
       if ((F[offset] & mask) != 0)
           F[offset] = F[offset] ^ mask;
   }

   for (size_t i = 3; i < cast(size_t)sqrt(n); i += 2)
       if (isSet((i - 3) / 2))
           for (size_t j = i * i; j <= n; j += 2 * i)
               resetBit((j - 3) / 2);

   Appender!(size_t[]) result;
   result.put(2);
   for(size_t i = 3; i <= n; i += 2)
       if (isSet((i - 3) / 2))
           result.put(i);
   return result.data;

}

void main() {

   writeln(primes(50));

}</lang>

Dart

<lang dart>// helper function to pretty print an Iterable String iterableToString(Iterable seq) {

 String str = "[";
 Iterator i = seq.iterator();
 while(i.hasNext()) {
   str += i.next();
   if(i.hasNext()) {
     str += ",";
   }
 }
 return str + "]";

}

main() {

 int limit = 1000;
 Set<int> sieve = new Set<int>();

 for(int i = 2; i <= limit; i++) {
   sieve.add(i);
 }
 for(int i = 2; i * i <= limit; i++) {
   if(sieve.contains(i)) {
     for(int j = i * i; j <= limit; j += i) {
       sieve.remove(j);
     }
   }
 }
 var sortedValues = new List<int>.from(sieve);
 sortedValues.sort((int arg1, int arg2) {
   return arg1.compareTo(arg2);
 });
 print(iterableToString(sortedValues));
 Expect.equals(168, sieve.length);

}</lang> Output:

[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,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997]

DWScript

<lang delphi>function Primes(limit : Integer) : array of Integer; var

  n, k : Integer;
  sieve := new Boolean[limit+1];

begin

  for n := 2 to Round(Sqrt(limit)) do begin
     if not sieve[n] then begin
        for k := n*n to limit step n do
           sieve[k] := True;
     end;
  end;
  
  for k:=2 to limit do
     if not sieve[k] then
        Result.Add(k);

end;

var r := Primes(50); var i : Integer; for i:=0 to r.High do

  PrintLn(r[i]);</lang>

Dylan

With outer to sqrt and inner to p^2 optimizations: <lang dylan>define method primes(n)

 let limit = floor(n ^ 0.5) + 1;
 let sieve = make(limited(<simple-vector>, of: <boolean>), size: n + 1, fill: #t);
 let last-prime = 2;

 while (last-prime < limit)
   for (x from last-prime ^ 2 to n by last-prime)
     sieve[x] := #f;
   end for;
   block (found-prime)
     for (n from last-prime + 1 below limit)
       if (sieve[n] = #f)
         last-prime := n;
         found-prime()
       end;
     end;
     last-prime := limit;
   end block;
 end while;

 for (x from 2 to n)
   if (sieve[x]) format-out("Prime: %d\n", x); end;
 end;

end;</lang>

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

eC

Note: this is not a Sieve of Eratosthenes; it is just trial division. <lang cpp> public class FindPrime {

  Array<int> primeList { [ 2 ], minAllocSize = 64 };
  int index;

  index = 3;

  bool HasPrimeFactor(int x)
  {
     int max = (int)floor(sqrt((double)x));
    
     for(i : primeList)
     {
        if(i > max) break;
        if(x % i == 0) return true;
     }
     return false;
  }

  public int GetPrime(int x)
  {
     if(x > primeList.count - 1)
     {
        for (; primeList.count != x; index += 2)
           if(!HasPrimeFactor(index))
           {
              if(primeList.count >= primeList.minAllocSize) primeList.minAllocSize *= 2;
              primeList.Add(index);
           }
     }
     return primeList[x-1];
  }

}

class PrimeApp : Application {

  FindPrime fp { };
  void Main()
  {
     int num = argc > 1 ? atoi(argv[1]) : 1;
     PrintLn(fp.GetPrime(num));
  }

} </lang>

Eiffel

Works with: EiffelStudio version 6.6 beta (with provisional loop syntax)

<lang eiffel>class

   APPLICATION

create

   make

feature

      make
           -- Run application.
       do
           across primes_through (100) as ic loop print (ic.item.out + " ") end
       end

   primes_through (a_limit: INTEGER): LINKED_LIST [INTEGER]
           -- Prime numbers through `a_limit'
       require
           valid_upper_limit: a_limit >= 2
       local
           l_tab: ARRAY [BOOLEAN]
       do
           create Result.make
           create l_tab.make_filled (True, 2, a_limit)
           across
               l_tab as ic
           loop
               if ic.item then
                   Result.extend (ic.target_index)
                   across ((ic.target_index * ic.target_index) |..| l_tab.upper).new_cursor.with_step (ic.target_index) as id
                   loop
                       l_tab [id.item] := False
                   end
               end
           end
       end

end</lang>

Output:

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

Euphoria

<lang euphoria>constant limit = 1000 sequence flags,primes flags = repeat(1, limit) for i = 2 to sqrt(limit) do

   if flags[i] then
       for k = i*i to limit by i do
           flags[k] = 0
       end for
   end if

end for

primes = {} for i = 2 to limit do

   if flags[i] = 1 then
       primes &= i
   end if

end for ? primes</lang>

Output:

{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,211,223,227,229,233,239,241,251,257,263,269,271,
277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,
383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,
487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,
601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,
709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,
827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,
947,953,967,971,977,983,991,997}

F#

Functional

This is not a Sieve of Eratosthenes, for the same reason as the naive Haskell version, because it simply filters the numbers.

#light

let not_divisible_by n =
    List.filter (fun m -> m % n <> 0)

let sieve max =
    let rec sieve_aux = function
        | [] -> []
        | x :: xs as l ->
        if x * x > max then l else x :: sieve_aux (not_divisible_by x xs)
    sieve_aux [2 .. max]

Forth

: prime? ( n -- ? ) here + c@ 0= ;
: composite! ( n -- ) here + 1 swap c! ;

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

100 sieve

Fortran

Works with: Fortran version 90 and later

<lang fortran>program sieve

 implicit none
 integer, parameter :: i_max = 100
 integer :: i
 logical, dimension (i_max) :: is_prime

 is_prime = .true.
 is_prime (1) = .false.
 do i = 2, int (sqrt (real (i_max)))
   if (is_prime (i)) is_prime (i * i : i_max : i) = .false.
 end do
 do i = 1, i_max
   if (is_prime (i)) write (*, '(i0, 1x)', advance = 'no') i
 end do
 write (*, *)

end program sieve</lang> Output: <lang>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</lang> Optimised using a pre-computed wheel based on 2: <lang fortran>program sieve_wheel_2

 implicit none
 integer, parameter :: i_max = 100
 integer :: i
 logical, dimension (i_max) :: is_prime

 is_prime = .true.
 is_prime (1) = .false.
 is_prime (4 : i_max : 2) = .false.
 do i = 3, int (sqrt (real (i_max))), 2
   if (is_prime (i)) is_prime (i * i : i_max : 2 * i) = .false.
 end do
 do i = 1, i_max
   if (is_prime (i)) write (*, '(i0, 1x)', advance = 'no') i
 end do
 write (*, *)

end program sieve_wheel_2</lang> Output: <lang>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</lang>

GAP

<lang gap>Eratosthenes := function(n)

  local sieve, cur, mul;
  sieve := [1 .. n]*0;
  sieve[1] := 1;
  cur := 1;
  while cur <= n do
     if sieve[cur] = 0 then
        mul := cur*cur;
        while mul <= n do
           sieve[mul] := 1;
           mul := mul + cur;
        od;
     fi;
     cur := cur + 1;
  od;
  return Filtered([1 .. n], x -> sieve[x] = 0);

end;

Eratosthenes(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 ]</lang>

Go

<lang go>package main

import (

   "fmt"

)

func main() {

   const limit = 201 // means sieve numbers < 201

   // sieve
   c := make([]bool, limit) // c for composite.  false means prime candidate
   c[1] = true              // 1 not considered prime
   p := 2
   for {
       // first allowed optimization:  outer loop only goes to sqrt(limit)
       p2 := p * p
       if p2 >= limit {
           break
       }
       // second allowed optimization:  inner loop starts at sqr(p)
       for i := p2; i < limit; i += p {
           c[i] = true // it's a composite

       }
       // scan to get next prime for outer loop
       for {
           p++
           if !c[p] {
               break
           }
       }
   }

   // sieve complete.  now print a representation.
   for n := 1; n < limit; n++ {
       if c[n] {
           fmt.Print("  .")
       } else {
           fmt.Printf("%3d", n)
       }
       if n%20 == 0 {
           fmt.Println("")
       }
   }

}</lang> Output:

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

Groovy

This solution uses a BitSet for compactness and speed, but in Groovy, BitSet has full List semantics. It also uses both the "square root of the boundary" shortcut and the "square of the prime" shortcut. <lang groovy>def sievePrimes = { bound ->

   def isPrime  = new BitSet(bound)
   isPrime[0..1] = false
   isPrime[2..<bound] = true
   (2..(Math.sqrt(bound))).each { pc ->
       if (isPrime[pc]) {
           ((pc**2)..<bound).step(pc) { isPrime[it] = false }
       }
   }
   (0..<bound).findAll { isPrime[it] }

}</lang>

Test: <lang groovy>println sievePrimes(100)</lang>

Output:

[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]

GW-BASIC

<lang qbasic>10 INPUT "ENTER NUMBER TO SEARCH TO: ";LIMIT 20 DIM FLAGS(LIMIT) 30 FOR N = 2 TO SQR (LIMIT) 40 IF FLAGS(N) < > 0 GOTO 80 50 FOR K = N * N TO LIMIT STEP N 60 FLAGS(K) = 1 70 NEXT K 80 NEXT N 90 REM DISPLAY THE PRIMES 100 FOR N = 2 TO LIMIT 110 IF FLAGS(N) = 0 THEN PRINT N;", "; 120 NEXT N</lang>

Haskell

Imperative version, with unboxed arrays

<lang haskell>import Control.Monad (forM_, when) import Control.Monad.ST import Data.Array.ST import Data.Array.Unboxed

primeSieve :: Integer -> UArray Integer Bool primeSieve top = runSTUArray $ do

   a <- newArray (2,top) True           -- :: ST s (STUArray s Integer Bool)
   let r = ceiling . sqrt $ fromInteger top
   forM_ [2..r] $ \i -> do
     ai <- readArray a i
     when ai $ do
       forM_ [i*i,i*i+i..top] $ \j -> do
         writeArray a j False
   return a

-- Return primes from sieve as list: primesTo :: Integer -> [Integer] primesTo top = [p | (p,True) <- assocs $ primeSieve top]</lang>

Unboxed arrays, on odds only

<lang haskell>sieveUO :: Int -> UArray Int Bool sieveUO top = runSTUArray $ do

   let m = (top-1) `div` 2
       r = floor . sqrt $ fromIntegral top + 1
   sieve <- newArray (1,m) True          -- :: ST s (STUArray s Int Bool)
   forM_ [1..r `div` 2] $ \i -> do
     isPrime <- readArray sieve i        -- ((2*i+1)^2-1)`div`2 = 2*i*(i+1)
     when isPrime $ do
       forM_ [2*i*(i+1), 2*i*(i+2)+1..m] $ \j -> do
         writeArray sieve j False
   return sieve

primesToUO :: Int -> [Int] primesToUO top | top > 1 = 2 : [2*p + 1 | (p,True) <- assocs $ sieveUO top]

              | otherwise = []</lang>

Using plain Int is faster but less general. This represents odds only in the array and is fastest. Empirical time complexity is $O(n^{1.2})$ in n primes produced, and improving for bigger n‍ ‍s. Memory consumption is low (array seems to be packed) but growing about linearly with n.

Basic list-based sieve

Straightforward implementation of the sieve of Ertosthenes in its original bounded form: <lang haskell>primesTo m = 2 : eratos [3,5..m] where

  eratos (p : xs) | p*p>m = p : xs
                  | True  = p : eratos (xs `minus` [p*p, p*p+2*p..])

minus a@(x:xs) b@(y:ys) = case compare x y of

        LT -> x : minus  xs b
        EQ ->     minus  xs ys
        GT ->     minus  a  ys

minus a b = a </lang> Its time complexity is similar to that of optimal trial division because of limitations of Haskell linked lists, where (minus a b) takes time proportional to length(union a b) and not (length b) as achieved in imperative setting with direct-access memory. Uses ordered list representation of sets.

List-based tree-merging incremental sieve

Basic version's (((a-b)-c)- ... ) workflow can be rearranged into (a-(b+(c+ ... ))). This is the idea behind Richard Bird's unbounded code presented in M. O'Neill's article. Linear merging structure can be further replaced with a tree-like one.

This finds primes in gaps between composites, and composites as union of the multiples of primes, merging in a tree-like fashion all the multiples streams, running with about same, as priority-queue–based version of M. O'Neill's, empirical time complexity of a steady $O(n^{1.2})$ , and in close to constant space too (not counting the produced sequence space): <lang haskell>primes :: [Integer] -- unbounded merging idea due to Richard Bird primes = 2 : g (fix g) -- double staged production idea due to Melissa O'Neill

where                      --  tree merging idea Dave Bayer / simplified formulation Will Ness 
 g xs = 3 : gaps 5 (unionAll [[p*p, p*p+2*p..] | p <- xs])

fix g = xs where xs = g xs -- global definition to prevent space leak

gaps k s@(x:xs) | k<x = k:gaps (k+2) s -- [k,k+2..]`minus`s, k<=x !

               | True =   gaps (k+2) xs

unionAll ((x:xs):t) = x : union xs (unionAll (pairs t))

 where 
   pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t

union a@(x:xs) b@(y:ys) = case compare x y of

        LT -> x : union  xs b
        EQ -> x : union  xs ys
        GT -> y : union  a  ys</lang>

Here's the test entry on Ideone.com, a comparison of more versions, and a similar code with wheel optimization.

HicEst

<lang hicest>REAL :: N=100, sieve(N)

sieve = $ > 1 ! = 0 1 1 1 1 ... DO i = 1, N^0.5

 IF( sieve(i) ) THEN
    DO j = i^2, N, i
      sieve(j) = 0
    ENDDO
 ENDIF

ENDDO

DO i = 1, N

 IF( sieve(i) ) WRITE() i

ENDDO </lang>

Icon and Unicon

<lang Icon> procedure main()

   sieve(100)
end

procedure sieve(n)
   local p,i,j
   p:=list(n, 1)
   every i:=2 to sqrt(n) & j:= i+i to n by i & p[i] == 1
     do p[j] := 0
   every write(i:=2 to n & p[i] == 1 & i)
end</lang>

Alternatively using sets <lang Icon> procedure main()

    sieve(100)
end

procedure sieve(n)
    primes := set()
    every insert(primes,1 to n)
    every member(primes,i := 2 to n) do
        every delete(primes,i + i to n by i)
    delete(primes,1)
    every write(!sort(primes))

end</lang>

J

Generally, this task should be accomplished in J using i.&.(p:inv) . Here we take an approach that's more comparable with the other examples on this page.

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" may be 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 .

Java

Works with: Java version 1.5+

<lang java5>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;
      }

}</lang>

To optimize by testing only odd numbers, replace the loop marked "unoptimized" with these lines: <lang java5>nums.add(2); for(int i = 3;i <= n;i += 2){

      nums.add(i);

}</lang>

Version using a BitSet: <lang java5>import java.util.LinkedList; import java.util.BitSet;

public class Sieve{

   public static LinkedList<Integer> sieve(int n){
       LinkedList<Integer> primes = new LinkedList<Integer>();
       BitSet nonPrimes = new BitSet(n+1);
       
       for (int p = 2; p <= n; p = nonPrimes.nextClearBit(p+1)) {
           for (int i = p * p; i <= n; i += p)
               nonPrimes.set(i);
           primes.add(p);
       }
       return primes;
   }

}</lang>

JavaScript

<lang javascript>function erasthenes(limit) {

   var primes = [];
   if (limit >= 2) {
       var nums = new Array(limit-1);
       for (var i = 2; i <= limit; i++)
           nums[i-2] = i;

       var last_prime;
       var idx = 0;
       while ((last_prime = nums[idx]) <= Math.sqrt(limit)) {
           if (last_prime != null)
               for (var i = idx + last_prime; i < limit - 1; i += last_prime) 
                   nums[i] = null;
           idx++;
       }
       for (var i = 0; i < nums.length; i++)
           if (nums[i] != null)
               primes.push(nums[i]);
   }
   return primes;

}

var primes = erasthenes(100);

if (typeof print == "undefined")

   print = (typeof WScript != "undefined") ? WScript.Echo : alert;

print(primes); </lang> outputs:

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

Liberty BASIC

<lang lb> 'Notice that arrays are globally visible to functions.

   'The sieve() function uses the flags() array.
   'This is a Sieve benchmark adapted from BYTE 1985
   ' May, page 286

   size = 7000
   dim flags(7001)
   start = time$("ms")
   print sieve(size); " primes found."
   print "End of iteration.  Elapsed time in milliseconds: "; time$("ms")-start
   end

   function sieve(size)
       for i = 0 to size
           if flags(i) = 0 then
               prime = i + i + 3
               k = i + prime
               while k <= size
                   flags(k) = 1
                   k = k + prime
               wend
               sieve = sieve + 1
           end if
       next i
   end function</lang>

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

Lua

<lang lua>function erato(n)

 local t = {0, 2}
 for i = 3, n, 2 do t[i], t[i+1] = i, 0 end
 for i = 3, math.sqrt(n) do for j = i*i, n, 2*i do t[j] = 0 end end
 return t

end</lang>

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

recursive

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

M4

<lang M4> define(`lim',100)dnl define(`for',

  `ifelse($#,0,
     ``$0,
     `ifelse(eval($2<=$3),1,
        `pushdef(`$1',$2)$5`'popdef(`$1')$0(`$1',eval($2+$4),$3,$4,`$5')')')')dnl

for(`j',2,lim,1,

  `ifdef(a[j],
     `',
     `j for(`k',eval(j*j),lim,j,
        `define(a[k],1)')')')

</lang>

Output:

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

Mathematica

<lang Mathematica> Eratosthenes[n_] :=

Module[{Numbers}, (* Make Numbers a local variable in the following (optional) *)
 Numbers = Table[i, {i, n}];

  For[i = 2, i <= Sqrt[n], i++, 
   If[Numbersi != 0, 
    For[j = 2, j <= n/i, j++,
     Numbersi j = 0;
    ]
   ]
  ];
 Select[Numbers, # > 1 &]
 ]

Eratosthenes[100] </lang>

MATLAB

<lang MATLAB>function primeList = sieveOfEratosthenes(lastNumber)

   list = (2:lastNumber); %Construct list of numbers
   primeList = []; %Preallocate prime list
   
   while( list(1)^2 <lastNumber )
       
       primeList = [primeList list(1)]; %add prime to the prime list
       list( mod(list,list(1))==0 ) = []; %filter out all multiples of the current prime
       
   end
   
   primeList = [primeList list]; %The rest of the numbers in the list are primes

end</lang>

Sample Output: <lang MATLAB>sieveOfEratosthenes(30)

ans =

    2     3     5     7    11    13    17    19    23    29</lang>

Somewhat optimized

<lang MATLAB> function P = erato(x) % Sieve of Eratosthenes: returns all primes between 2 and x

   P = [0 2:x] ;            % Create vector with all ints between 2 and x where
                            %   position 1 is hard-coded as 0 since 1 is not a prime.

   for (n=2:sqrt(x))        % All primes factors lie between 2 and sqrt(x).
      if P(n)               % If the current value is not 0 (i.e. a prime),
         P((2*n):n:x) = 0 ; % then replace all further multiples of it with 0.
      end
   end                      % At this point P is a vector with only primes and zeroes.

   P = P(P ~= 0) ;          % Remove all zeroes from P, leaving only the primes.

return </lang>

The optimization lies in fewer steps in the for loop, use of MATLAB's built-in array operations and no modulo calculation.

Limitation: your machine has to be able to allocate enough memory for an array of length x.

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

Modula-3

Regular version

<lang modula3> MODULE Prime EXPORTS Main;

IMPORT IO;

CONST LastNum = 1000;

VAR a: ARRAY [2..LastNum] OF BOOLEAN;

BEGIN

 FOR i := FIRST(a) TO LAST(a) DO
   a[i] := TRUE;
 END;

 FOR i := FIRST(a) TO LAST(a) DO
   IF a[i] THEN
     IO.PutInt(i);
     IO.Put(" ");
     FOR j := FIRST(a) TO LAST(a) DO
       IF j MOD i = 0 THEN
         a[j] := FALSE;
       END;
     END;
   END;
 END;
 IO.Put("\n");

END Prime. </lang>

Advanced version

This version uses more "advanced" types. <lang modula3>(* From the CM3 examples folder (comments removed). *)

MODULE Sieve EXPORTS Main;

IMPORT IO;

TYPE

 Number = [2..1000];
 Set = SET OF Number;

VAR

 prime: Set := Set {FIRST(Number) .. LAST(Number)};

BEGIN

 FOR i := FIRST(Number) TO LAST(Number) DO
   IF i IN prime THEN
     IO.PutInt(i);
     IO.Put(" ");

     FOR j := i TO LAST(Number) BY i DO
       prime := prime - Set{j};
     END;
   END;
 END;
 IO.Put("\n");

END Sieve.</lang>

MUMPS

<lang MUMPS>ERATO1(HI)

;performs the Sieve of Erotosethenes up to the number passed in.
;This version sets an array containing the primes
SET HI=HI\1
KILL ERATO1 ;Don't make it new - we want it to remain after we quit the function
NEW I,J,P
FOR I=2:1:(HI**.5)\1 FOR J=I*I:I:HI SET P(J)=1
FOR I=2:1:HI S:'$DATA(P(I)) ERATO1(I)=I
KILL I,J,P
QUIT</lang>

Example:

USER>SET MAX=100,C=0 DO ERATO1^ROSETTA(MAX) 
USER>WRITE !,"PRIMES BETWEEN 1 AND ",MAX,! FOR  SET I=$ORDER(ERATO1(I)) Q:+I<1  WRITE I,", "

PRIMES BETWEEN 1 AND 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,

NetRexx

<lang Rexx>/* NetRexx */

options replace format comments java crossref savelog symbols binary

parse arg loWatermark hiWatermark . if loWatermark = | loWatermark = '.' then loWatermark = 1 if hiWatermark = | hiWatermark = '.' then hiWatermark = 200

do

 if \loWatermark.datatype('w') | \hiWatermark.datatype('w') then -
   signal NumberFormatException('arguments must be whole numbers')
 if loWatermark > hiWatermark then -
   signal IllegalArgumentException('the start value must be less than the end value')

 seive = sieveOfEratosthenes(hiWatermark)
 primes = getPrimes(seive, loWatermark, hiWatermark).strip

 say 'List of prime numbers from' loWatermark 'to' hiWatermark 'via a "Sieve of Eratosthenes" algorithm:'
 say '  'primes.changestr(' ', ',')
 say '  Count of primes:' primes.words

catch ex = Exception

 ex.printStackTrace

end

return

method sieveOfEratosthenes(hn = long) public static binary returns Rexx

 sv = Rexx(isTrue)
 sv[1] = isFalse
 ix = long
 jx = long

 loop ix = 2 while ix * ix <= hn
   if sv[ix] then loop jx = ix * ix by ix while jx <= hn
     sv[jx] = isFalse
     end jx
   end ix

 return sv

method getPrimes(seive = Rexx, lo = long, hi = long) private constant binary returns Rexx

 primes = Rexx()
 loop p_ = lo to hi
   if \seive[p_] then iterate p_
   primes = primes p_
   end p_

 return primes

method isTrue public constant binary returns boolean

 return 1 == 1

method isFalse public constant binary returns boolean

 return \isTrue

</lang>

Output

List of prime numbers from 1 to 200 via a "Sieve of Eratosthenes" algorithm:
  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
  Count of primes: 46

Nial

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

Using it

|primes 10
=2 3 5 7

Niue

<lang Niue>[ dup 2 < ] '<2 ; [ 1 + 'count ; [ <2 [ , ] when ] count times ] 'fill-stack ;

0 'n ; 0 'v ;

[ .clr 0 'n ; 0 'v ; ] 'reset ; [ len 1 - n - at 'v ; ] 'set-base ; [ n 1 + 'n ; ] 'incr-n ; [ mod 0 = ] 'is-factor ; [ dup * ] 'sqr ;

[ set-base

 v sqr 2 at > not 
 [ [ dup v = not swap v is-factor and ] remove-if incr-n run ] when ] 'run ;

[ fill-stack run ] 'sieve ;

( tests )

10 sieve .s ( => 2 3 5 7 9 ) reset newline 30 sieve .s ( => 2 3 5 7 11 13 17 19 23 29 ) </lang>

Oberon-2

<lang oberon2>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.</lang>

OCaml

Imperative

<lang ocaml>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)</lang>

<lang ocaml>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))</lang>

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]

Functional

<lang ocaml>(* 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 remove 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]</lang>

Another functional version

This uses zero to denote struck-out numbers. It is slightly inefficient as it strikes-out multiples above p rather than p²

<lang ocaml># let rec strike_nth k n l = match l with

 | [] -> []
 | h :: t ->
   if k = 0 then 0 :: strike_nth (n-1) n t
   else h :: strike_nth (k-1) n t;;

val strike_nth : int -> int -> int list -> int list = <fun>

let primes n =

 let limit = truncate(sqrt(float n)) in
 let rec range a b = if a > b then [] else a :: range (a+1) b in
 let rec sieve_primes l = match l with
   | [] -> []
   | 0 :: t -> sieve_primes t
   | h :: t -> if h > limit then List.filter ((<) 0) l else
       h :: sieve_primes (strike_nth (h-1) h t) in
 sieve_primes (range 2 n) ;;

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]</lang>

Oz

Translation of: Haskell

<lang oz>declare

 fun {Sieve N}
    S = {Array.new 2 N true}
    M = {Float.toInt {Sqrt {Int.toFloat N}}}
 in
    for I in 2..M do

if S.I then for J in I*I..N;I do S.J := false end end

    end
    S
 end

 fun {Primes N}
    S = {Sieve N}
 in
    for I in 2..N collect:C do

if S.I then {C I} end

    end
 end

in

 {Show {Primes 30}}</lang>

PARI/GP

<lang parigp>Eratosthenes(lim)={

 my(v=Vectorsmall(lim,unused,1));
 forprime(p=2,sqrt(lim),
   forstep(i=p^2,lim,p,
     v[i]=0
   )
 );
 for(i=1,lim,if(v[i],print1(i", ")))

};</lang>

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. <lang pascal> 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. </lang>

Perl

<lang perl>$max= $ARGV[0]; @primes= (); @tested= (1); $j= 1; while ($j < $max) {

  next if $tested[$j++];
  push @primes, $j;
  for ($k= $j; $k <= $max; $k+=$j) {
     $tested[$k-1]= 1;
  }

} print "@primes\n";</lang>

Version using bit strings: <lang perl>sub sieve {

   my ($n) = @_;

   my @primes;
   my $nonPrimes = ;

   foreach my $p (2 .. $n) {
       unless (vec($nonPrimes, $p, 1)) {
           for (my $i = $p * $p; $i <= $n; $i += $p) {
               vec($nonPrimes, $i, 1) = 1;
           }
           push @primes, $p;
       }
   }
   @primes

}

my $prime_limit; print "Input the number that you wish to test :"; chomp($prime_limit = <STDIN>); my @primes = sieve($prime_limit); print "The primes are: @primes!!!\n";</lang>

Golfing a bit: <lang perl>sub sieve{ my (@s, $i); grep { not $s[ $i = $_ ] and do { $s[ $i += $_ ]++ while $i <= $_[0]; 1 } } 2..$_[0] }

print join ", " => sieve 100; </lang>

Or with bit strings: <lang perl>sub sieve{ my ($s, $i); grep { not vec $s, $i = $_, 1 and do { (vec $s, $i += $_, 1) = 1 while $i <= $_[0]; 1 } } 2..$_[0] }

print join ", " => sieve 100; </lang>

A sort, recursive version: <lang perl>sub erat {

   my $p = shift;
   return $p, $p**2 > $_[$#_] ? @_ : erat(grep $_%$p, @_)

}

print join ', ' => erat 2..100000; </lang>

Regexp:<lang perl>sub sieve { my ($s, $p) = "." . ("x" x shift);

1 while ++$p and $s =~ /^(.{$p,}?)x/g and $p = length($1) and $s =~ s/(.{$p})./${1}./g and substr($s, $p, 1) = "x"; $s }

print sieve(1000);</lang>

Perl 6

A very short and slow solution: <lang perl6> sub erat { grep none(@^a[^sqrt @a] X* @a), @a }

say erat 2..100 </lang>

A recursive version: <lang perl6> sub erat(@n is copy) {

   ($_ = @n.shift),
   $_**2 > @n[*-1] ?? @n !! &?ROUTINE(grep * % $_, @n)

}

say erat 2..1000 </lang>

An other version:

<lang perl6>sub sieve( Int $limit ) {

   my @is-prime = False, False, True xx $limit - 1;

   gather for @is-prime.kv -> $number, $is-prime {
       if $is-prime {
           take $number;
           @is-prime[$_] = False if $_ %% $number for $number**2 .. $limit; 
       }
   }

}

(sieve 100).join(",").say;</lang>

PHP

<lang php> function iprimes_upto($limit) {

   for ($i = 2; $i < $limit; $i++)
   {

$primes[$i] = true;

   }
   
   for ($n = 2; $n < $limit; $n++)
   {

if ($primes[$n]) { for ($i = $n*$n; $i < $limit; $i += $n) { $primes[$i] = false; } }

   }
   
   return $primes;

} </lang>

PicoLisp

<lang PicoLisp>(de sieve (N)

  (let Sieve (range 1 N)
     (set Sieve)
     (for I (cdr Sieve)
        (when I
           (for (S (nth Sieve (* I I)) S (nth (cdr S) I))
              (set S) ) ) )
     (filter bool Sieve) ) )</lang>

Output:

: (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)

PL/I

<lang PL/I> eratos: proc options (main) reorder;

dcl i fixed bin (31); dcl j fixed bin (31); dcl n fixed bin (31); dcl sn fixed bin (31);

dcl hbound builtin; dcl sqrt builtin;

dcl sysin file; dcl sysprint file;

get list (n); sn = sqrt(n);

begin;

 dcl primes(n) bit (1) aligned init ((*)((1)'1'b));

 i = 2;

 do while(i <= sn);
   do j = i ** 2 by i to hbound(primes, 1);
     /* Adding a test would just slow down processing! */
     primes(j) = '0'b; 
    end;

   do i = i + 1 to sn until(primes(i));
   end;
 end;

 do i = 2 to hbound(primes, 1);
   if primes(i) then
     put data(i);
 end;

end; end eratos; </lang>

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;

PowerShell

Basic procedure

It outputs immediately so that the number can be used by the pipeline. <lang PowerShell>function Sieve ( [int] $num ) {

   $isprime = @{}
   2..$num | Where-Object {
       $isprime[$_] -eq $null } | ForEach-Object {
       $_
       $isprime[$_] = $true
       $i=$_*$_
       for ( ; $i -le $num; $i += $_ )
       { $isprime[$i] = $false }
   }

}</lang>

Processing

Calculate the primes up to 1000000 with Processing, including a visualisation of the process. As an additional visual effect, the layout of the pixel could be changed from the line-by-line layout to a spiral-like layout starting in the middle of the screen. <lang java>int maxx,maxy; int max; boolean[] sieve;

void plot(int pos, boolean active) {

 set(pos%maxx,pos/maxx, active?#000000:#ffffff);

}

void setup() {

 size(1000, 1000, P2D);
 frameRate(2);
 maxx=width;
 maxy=height;
 max=width*height;
 sieve=new boolean[max+1];

 sieve[1]=false;
 plot(0,false);
 plot(1,false);
 for(int i=2;i<=max;i++) {
   sieve[i]=true;
   plot(i,true);
 }

}

int i=2;

void draw() {

 if(!sieve[i]) {
   while(i*i<max && !sieve[i]) {
     i++;
   }
 }
 if(sieve[i]) {
   print(i+" ");
   for(int j=i*i;j<=max;j+=i) {
     if(sieve[j]) {
       sieve[j]=false;
       plot(j,false);
     }
   }
 }
 if(i*i<max) {
   i++;
 } else {
   noLoop();
   println("finished");
 }

}</lang>

Prolog

using facts as a set of non-prime numbers

Builds a set of non-prime numbers (multiples of already visited numbers) pretty much the same way as imperative languages, by counting up in increments of p for each prime p (or 2p for odd primes). Works with SWI-Prolog.

<lang Prolog>sieve(N, S) :-

   retractall(mult(_)),
   sieve_(2,N,S).

sieve_(I,N,Xs) :-

   I =< N, !,
   I1 is I+1,
   (   mult(I)
   ->  sieve_(I1,N,Xs)
   ;   ISq is I*I,
       ( I=:=2 -> add_mults(I,ISq,N) 
         ; DI is 2*I, add_mults(DI,ISq,N) ),
       Xs=[ I|Ys ],
       sieve_(I1,N,Ys)
   ).

sieve_(I,N,[]) :- I > N.

add_mults(DI,I,N) :-

   I =< N, !,
   ( mult(I) -> true ; assert(mult(I)) ),
   I1 is I+DI,
   add_mults(DI,I1,N).

add_mults(_,I,N) :- I > N.</lang>

PureBasic

Basic procedure

<lang PureBasic>For n=2 To Sqr(lim)

 If Nums(n)=0
   m=n*n
   While m<=lim
     Nums(m)=1
     m+n
   Wend
 EndIf

Next n</lang>

Working example

<lang PureBasic>Dim Nums.i(0) Define l, n, m, lim

If OpenConsole()

 ; Ask for the limit to search, get that input and allocate a Array
 Print("Enter limit for this search: ")
 lim=Val(Input())
 ReDim Nums(lim)
 
 ; Use a basic Sieve of Eratosthenes
 For n=2 To Sqr(lim)
   If Nums(n)=#False
     m=n*n
     While m<=lim
       Nums(m)=#True
       m+n
     Wend
   EndIf
 Next n
 
 ;Present the result to our user
 PrintN(#CRLF$+"The Prims up to "+Str(lim)+" are;")
 m=0: l=Log10(lim)+1
 For n=2 To lim
   If Nums(n)=#False
     Print(RSet(Str(n),l)+" ")
     m+1
     If m>72/(l+1)
       m=0: PrintN("")
     EndIf
   EndIf
 Next
 
 Print(#CRLF$+#CRLF$+"Press ENTER to exit"): Input()
 CloseConsole()

EndIf</lang>

Output may look like;

Enter limit for this search: 750

The Prims up to 750 are;
   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  211  223  227  229  233  239  241  251  257  263  269  271  277  281
 283  293  307  311  313  317  331  337  347  349  353  359  367  373  379
 383  389  397  401  409  419  421  431  433  439  443  449  457  461  463
 467  479  487  491  499  503  509  521  523  541  547  557  563  569  571
 577  587  593  599  601  607  613  617  619  631  641  643  647  653  659
 661  673  677  683  691  701  709  719  727  733  739  743

Press ENTER to exit

Python

Testing on list of primes up to that number

This is not a Sieve of Eratosthenes, it is just trial division. <lang python>def primes_trial(n):

   primes = [2]
   
   for x in range(2, n+1):
       if all(x % p for p in primes):
           primes.append(x)
   
   return primes

print(primes_trial(100))</lang>

Using list of primes up to that number using list comprehension

This is not a Sieve of Eratosthenes, because it simply filters the numbers. <lang python>def GetPrimes(p): i = 1 while i**2 < p[-1]: p[i:] = [x for x in p[i:] if x % p[i-1] != 0] i += 1 return p

print GetPrimes(range(2,100))</lang>

Using list lookup

<lang python>def eratosthenes(n):

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

eratosthenes(100)</lang>

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.

<lang python>def eratosthenes2(n):

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

eratosthenes2(100)</lang>

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. <lang python>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]</lang>

Using generator

<lang python> 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

</lang> Example: <lang python> >>> list(iprimes_upto(15)) [2, 3, 5, 7, 11, 13] </lang>

Using numpy

Library: numpy

Belowed code adapted from literateprograms.org using numpy <lang python>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:]

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

Using wheels

Version with wheel based optimization: <lang python>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]

</lang>

Examples: <lang python> >>> 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]) </lang> The above-mentioned examples demonstrate that the given wheel based optimization does not show significant performance gain.

Infinite generator

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

Works with: Python version 2.6, 3.x

<lang python>import heapq

generates all prime numbers

def sieve():

   # priority queue of the 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][0]: # non-prime
           while nonprimes[0][0] == i:
               # for each sequence that generates this number,
               # have it go to the next number (simply add the prime)
               # and re-position it in the priority queue
               x = nonprimes[0]
               x[0] += x[1]
               heapq.heapreplace(nonprimes, x)
       
       else: # prime
           # insert a 2-element list into the priority queue:
           # [current multiple, prime]
           # the first element allows sorting by value of current multiple
           # we start with i^2
           heapq.heappush(nonprimes, [i*i, i])
           yield i
       
       i += 1</lang>

Example:

>>> foo = sieve()
>>> for i in range(20):
...     print(next(foo))
... 
2
3
5
7
11
13
17
19
23
29
31
37
41
43
47
53
59
61
67
71

Infinite generator with a faster algorithm

Works with: Python version 2.4 through 3.1

<lang python> def primes():

   composites = {} # A mapping from composite numbers to the smallest prime
                   # that is a factor of it (its witness).
   n = 2 # The current number being considered as a prime.
   while True:
       if n not in composites:
           yield n # Not a composite, therefore prime.
           composites[n**2] = n # The next unseen composite number is n squared.
       else:
           # n is composite. Find the next unseen composite number with the
           # same witness as n.
           witness = composites.pop(n)
           next = n + witness
           while next in composites:
               next += witness
           composites[next] = witness
       n += 1

import itertools def primes_up_to(limit):

   return list(itertools.takewhile(lambda p: p <= limit, primes()))

</lang>

R

This code is vectorised, so it runs quickly for n < one million. The allocation of the primes vector means memory usage is too high for n much larger than that. <lang R> sieve <- function(n) {

  n <- as.integer(n)
  if(n > 1e6) stop("n too large")
  primes <- rep(TRUE, n)
  primes[1] <- FALSE
  last.prime <- 2L
  for(i in last.prime:floor(sqrt(n)))
  {
     primes[seq.int(2L*last.prime, n, last.prime)] <- FALSE
     last.prime <- last.prime + min(which(primes[(last.prime+1):n]))
  }
  which(primes)

}

sieve(1000) </lang>

REXX

<lang REXX> /*REXX program, generate primes via sieve of Eratosthenes algorithm*/

highItem=200 /*define ten score of #'s. */ widthH=length(highItem) /*width of biggest element#.*/ a.=1 /*assume all numbers prime. */ a.1=0 /*special case for unity, by*/

                                      /*definition, 1 isn't prime.*/

 do j=2 to highItem                   /*starting @ 2, strike mult.*/
 if j*j>highItem then leave           /*process up to û(highItem).*/
 if \a.j then iterate                 /*this number isn't prime,  */
                                      /*it was marked as not prime*/

   do k=j+j to highItem by j          /*start with next multiple. */
   a.k=0                              /*mark J multiple not prime.*/
   end   /*k*/

 end     /*j*/

np=0 /*number of primes so far. */

 do q=1 for highItem                  /*list all the primes found.*/
 if a.q then do
             np=np+1                  /*bump the prime count.     */
             if np>1000 then iterate  /*only list up to 1k primes.*/
             say 'prime number' right(np,widthH) "is" right(q,widthH)
             if np==1000 then say 'Prime number listing truncated.'
             end
 end   /*q*/

say np 'prime numbers found up to' highItem </lang> Output:

prime number   1 is   2
prime number   2 is   3
prime number   3 is   5
prime number   4 is   7
prime number   5 is  11
prime number   6 is  13
prime number   7 is  17
prime number   8 is  19
prime number   9 is  23
prime number  10 is  29
prime number  11 is  31
prime number  12 is  37
prime number  13 is  41
prime number  14 is  43
prime number  15 is  47
prime number  16 is  53
prime number  17 is  59
prime number  18 is  61
prime number  19 is  67
prime number  20 is  71
prime number  21 is  73
prime number  22 is  79
prime number  23 is  83
prime number  24 is  89
prime number  25 is  97
prime number  26 is 101
prime number  27 is 103
prime number  28 is 107
prime number  29 is 109
prime number  30 is 113
prime number  31 is 127
prime number  32 is 131
prime number  33 is 137
prime number  34 is 139
prime number  35 is 149
prime number  36 is 151
prime number  37 is 157
prime number  38 is 163
prime number  39 is 167
prime number  40 is 173
prime number  41 is 179
prime number  42 is 181
prime number  43 is 191
prime number  44 is 193
prime number  45 is 197
prime number  46 is 199
46 prime numbers found up to 200

Ruby

eratosthenes starts with nums = [0, 0, 2, 3, 4, 5, ..., n], then marks (zeroes out) multiples of 2, 3, 5, 7, ... there, then returns all non-zero numbers which are the primes.

<lang ruby>def eratosthenes(n)

 nums = [0, 0] + (2..n).to_a
 (2..Math.sqrt(n).to_i).each do |i|
   if nums[i].nonzero?
     (i**2..n).step(i) {|m| nums[m] = 0}
   end
 end
 nums.find_all {|m| m.nonzero?}

end

p eratosthenes(100)</lang>

[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]

With a wheel

eratosthenes2 adds more optimizations, but the code is longer.

The array nums only tracks odd numbers (skips multiples of 2).
The array nums holds booleans instead of integers, and every multiple of 3 begins false.
The outer loop skips multiples of 2 and 3.
Both inner loops skip multiples of 2 and 3.

<lang ruby>def eratosthenes2(n)

 # For odd i, if i is prime, nums[i >> 1] is true.
 # Set false for all multiples of 3.
 nums = [true, false, true] * ((n + 5) / 6)
 nums[0] = false  # 1 is not prime.
 nums[1] = true   # 3 is prime.

 # Outer loop and both inner loops are skipping multiples of 2 and 3.
 # Outer loop checks i * i > n, same as i > Math.sqrt(n).
 i = 5
 until(m = i * i
       m > n) do
   if nums[i >> 1]
     i_times_2 = i << 1
     i_times_4 = i << 2
     while m <= n
       nums[m >> 1] = false
       m += i_times_2
       nums[m >> 1] = false
       m += i_times_4  # When i = 5, skip 45, 75, 105, ...
     end
   end
   i += 2
   if nums[i >> 1]
     m = i * i
     i_times_2 = i << 1
     i_times_4 = i << 2
     while m <= n
       nums[m >> 1] = false
       m += i_times_4  # When i = 7, skip 63, 105, 147, ...
       nums[m >> 1] = false
       m += i_times_2
     end
   end
   i += 4
 end

 primes = [2]
 nums.each_index {|i| primes << (i * 2 + 1) if nums[i]}
 primes.pop while primes.last > n
 return primes

end

p eratosthenes2(100)</lang>

This simple benchmark compares eratosthenes with eratosthenes2.

<lang ruby>require 'benchmark' Benchmark.bmbm {|x|

 x.report("eratosthenes") { eratosthenes(1_000_000) }
 x.report("eratosthenes2") { eratosthenes2(1_000_000) }

}</lang>

eratosthenes2 runs about 4 times faster than eratosthenes.

With the standard library

MRI 1.9.x implements the sieve of Eratosthenes at file prime.rb, class EratosthensesSeive (around line 421). This implementation optimizes for space, by packing the booleans into 16-bit integers. It also hardcodes all primes less than 256.

<lang ruby>require 'prime' p Prime::EratosthenesGenerator.new.take_while {|i| i <= 100}</lang>

Scala

<lang scala>def sieveOfEratosthenes(n: Int) = {

 val primes = new scala.collection.mutable.BitSet(n)
 primes ++= (2 to n)
 val sqrt = Math.sqrt(n).toInt
 for {
   candidate <- 2 to sqrt
   if primes contains candidate
 } primes --= candidate * candidate to n by candidate
 primes

}

println( sieveOfEratosthenes(30) )</lang>

BitSet(2, 3, 5, 7, 11, 13, 17, 19, 23, 29)

Using Streams (the "unfaithful sieve", see the link under the Haskell entry): <lang scala>def ints(n:Int):Stream[Int] = Stream.cons(n,ints(n+1)) def sieve(nums:Stream[Int]):Stream[Int] =

   Stream.cons(nums.head, sieve((nums.tail) filter (_%nums.head != 0)))

def primes = sieve(ints(2))

println( primes take 10 toList ) println( primes takeWhile (_<30) toList ) </lang> Both println statements give the same results:

List(2, 3, 5, 7, 11, 13, 17, 19, 23, 29)

Scheme

Works with: Scheme version R

^{5}

RS

<lang scheme>(define (iota start stop stride)

 (if (> start stop)
     (list)
     (cons start (iota (+ start stride) stop stride))))

(define (strike lst start stride)

 (cond ((null? lst) lst)
       ((= (car lst) start) (strike (cdr lst) (+ start stride) stride))
       ((> (car lst) start) (strike lst (+ start stride) stride))
       (else (cons (car lst) (strike (cdr lst) start stride)))))

(define (primes limit)

 (let ((stop (sqrt limit)))
   (define (sieve lst)
     (let ((p (car lst)))
       (if (> p stop)
           lst
           (cons p (sieve (strike (cdr lst) (* p p) p))))))
   (sieve (iota 2 limit 1))))

(display (primes 100)) (newline)</lang> Output: <lang>(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)</lang> Optimised using a pre-computed wheel based on 2 (i.e. odds only): <lang scheme>(define (primes-wheel-2 limit)

 (let ((stop (sqrt limit)))
   (define (sieve lst)
     (let ((p (car lst)))
       (if (> p stop)
           lst
           (cons p (sieve (strike (cdr lst) (* p p) (* 2 p)))))))
   (cons 2 (sieve (iota 3 limit 2)))))

(display (primes-wheel-2 100)) (newline)</lang> Output: <lang>(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)</lang> Vector-based (faster), works with R $^{5}$ RS: <lang scheme>; initialize v to vector of sequential integers (define (initialize! v)

 (define (iter v n) (if (>= n (vector-length v)) 
                        (values) 
                        (begin (vector-set! v n n) (iter v (+ n 1)))))
 (iter v 0))

set every nth element of vector v to 0,
starting with element m

(define (strike! v m n)

 (cond ((>= m (vector-length v)) (values))
       (else (begin
               (vector-set! v m 0)
               (strike! v (+ m n) n)))))

lowest non-zero index of vector v >= n

(define (nextprime v n)

 (if (zero? (vector-ref v n))
     (nextprime v (+ n 1))
     (vector-ref v n)))

remove elements satisfying pred? from list lst

(define (remove pred? lst)

 (cond 
   ((null? lst) '())
   ((pred? (car lst))(remove pred? (cdr lst)))
   (else (cons (car lst) (remove pred? (cdr lst))))))

the sieve itself

(define (sieve n)

 (define stop (sqrt n))
 (define (iter v p)
   (cond 
     ((> p stop) v)
     (else 
      (begin
        (strike! v (* p p) p)
        (iter v (nextprime v (+ p 1)))))))
 
 (let ((v (make-vector (+ n 1))))
   (initialize! v)
   (vector-set! v 1 0) ; 1 is not a prime
   (remove zero? (vector->list (iter v 2)))))</lang>

Seed7

<lang seed7>const func array boolean: eratosthenes (in integer: n) is func

 result
   var array boolean: sieve is 0 times FALSE;
 local
   var integer: i is 0;
   var integer: j is 0;
 begin 
   sieve := n times TRUE;
   sieve[1] := FALSE;
   for i range 2 to sqrt(n) do
     if sieve[i] then
       for j range i * i to n step i do
         sieve[j] := FALSE;
       end for;
     end if;
   end for;
 end func;</lang>

Original source: [1]

SNOBOL4

Using strings instead of arrays, and the square/sqrt optimizations.

<lang SNOBOL4> define('sieve(n)i,j,k,p,str,res') :(sieve_end) sieve i = lt(i,n - 1) i + 1 :f(sv1)

       str = str (i + 1) ' ' :(sieve)

sv1 str break(' ') . j span(' ') = :f(return)

       sieve = sieve j ' '
       sieve = gt(j ^ 2,n) sieve str :s(return) ;* Opt1
       res = 
       str (arb ' ') @p ((j ^ 2) ' ') ;* Opt2
       str len(p) . res = ;* Opt2

sv2 str break(' ') . k span(' ') = :f(sv3)

       res = ne(remdr(k,j),0) res k ' ' :(sv2)

sv3 str = res :(sv1) sieve_end

# Test and display

       output = sieve(100)

end</lang>

Output:

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

Tcl

<lang tcl>package require Tcl 8.5

proc sieve n {

   if {$n < 2} {return {}}
   
   # create a container to hold the sequence of numbers.
   # use a dictionary for its speedy access (like an associative array) 
   # and for its insertion order preservation (like a list)
   set nums [dict create]
   for {set i 2} {$i <= $n} {incr i} {
       # the actual value is never used
       dict set nums $i ""
   }
   
   set primes [list]
   while {[set nextPrime [lindex [dict keys $nums] 0]] <= sqrt($n)} {
       dict unset nums $nextPrime
       for {set i [expr {$nextPrime ** 2}]} {$i <= $n} {incr i $nextPrime} {
           dict unset nums $i
       }
       lappend primes $nextPrime
   }
   return [concat $primes [dict keys $nums]]

}

puts [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</lang>

UNIX Shell

With array

Works with: zsh

<lang bash>#!/usr/bin/zsh

function primes() { typeset -a a typeset i j

a[1]="" for (( i = 2; i <= $1; i++ )); do a[$i]=$i done

for (( i = 2; i * i <= $1; i++ )); do if [[ ! -z $a[$i] ]]; then for (( j = i * i; j <= $1; j += i )); do a[$j]="" done fi done print $a }

primes 1000</lang>

Works with: bash

Works with: ksh93

Works with: pdksh

<lang bash>function primes { typeset a i=2 j m=$1 # No for (( ... )) loop in pdksh. Use while loop. while (( i <= m )); do a[$i]=$i (( i++ )) done

i=2 while (( j = i * i, j <= m )); do if [[ -n ${a[$i]} ]]; then while (( j <= m )); do unset a[$j] (( j += i )) done fi (( i++ )) done # No print command in bash. Use echo command. echo ${a[*]} }

primes 1000</lang>

Both scripts output a single long line.

2 3 5 7 11 13 17 19 23 ... 971 977 983 991 997

Using variables as fake array

Bourne Shell and Almquist Shell have no arrays. This script works with bash or dash (standard shell in Ubuntu), but uses no specifics of the shells, so it works with plain Bourne Shell as well.

Works with: Bourne Shell

<lang bash>#! /bin/sh

LIMIT=1000

As a workaround for missing arrays, we use variables p2, p3, ...,
p$LIMIT, to represent the primes. Values are true or false.
eval p$i=true # Set value.
eval \$p$i # Run command: true or false.
A previous version of this script created a temporary directory and
used files named 2, 3, ..., $LIMIT to represent the primes. We now use
variables so that a killed script does not leave extra files. About
performance, variables are about as slow as files.

i=2 while [ $i -le $LIMIT ] do

   eval p$i=true               # was touch $i
   i=`expr $i + 1`

done

i=2 while

   j=`expr $i '*' $i`
   [ $j -le $LIMIT ]

do

   if eval \$p$i               # was if [ -f $i ]
   then
       while [ $j -le $LIMIT ]
       do
           eval p$j=false      # was rm -f $j
           j=`expr $j + $i`
       done
   fi
   i=`expr $i + 1`

done

was echo `ls|sort -n`

echo `i=2

     while [ $i -le $LIMIT ]; do
         eval \\$p$i && echo $i
         i=\`expr $i + 1\`
     done`</lang>

With piping

This version works by piping 1s and 0s through sed. The string of 1s and 0s represents the array of primes.

Works with: Bourne Shell

<lang bash># Fill $1 characters with $2. fill () { # This pipeline would begin # head -c $1 /dev/zero | ... # but some systems have no head -c. Use dd. dd if=/dev/zero bs=$1 count=1 2>/dev/null | tr '\0' $2 }

filter () { # Use sed to put an 'x' after each multiple of $1, remove # first 'x', and mark non-primes with '0'. sed -e s/$2/\&x/g -e s/x// -e s/.x/0/g | { if expr $1 '*' $1 '<' $3 > /dev/null; then filter `expr $1 + 1` .$2 $3 else cat fi } }

Generate a sequence of 1s and 0s indicating primality.

oz () { fill $1 1 | sed s/1/0/ | filter 2 .. $1 }

Echo prime numbers from 2 to $1.

prime () { # Escape backslash inside backquotes. sed sees one backslash. echo `oz $1 | sed 's/./&\\ /g' | grep -n 1 | sed s/:1//` }

prime 1000</lang>

C Shell

Translation of: CMake

<lang csh># Sieve of Eratosthenes: Echoes all prime numbers through $limit. @ limit = 80

if ( ( $limit * $limit ) / $limit != $limit ) then echo limit is too large, would cause integer overflow. exit 1 endif

Use $prime[2], $prime[3], ..., $prime[$limit] as array of booleans.
Initialize values to 1 => yes it is prime.

set prime=( `repeat $limit echo 1` )

Find and echo prime numbers.

@ i = 2 while ( $i <= $limit ) if ( $prime[$i] ) then echo $i

# For each multiple of i, set 0 => no it is not prime. # Optimization: start at i squared. @ m = $i * $i while ( $m <= $limit ) set prime[$m] = 0 @ m += $i end endif @ i += 1 end</lang>

Ursala

with no optimizations <lang Ursala>#import nat

sieve = ~<{0,1}&& iota; @NttPX ~&r->lx ^/~&rhPlC remainder@rlX~|@r</lang> test program: <lang Ursala>#cast %nL

example = sieve 50</lang> output:

<2,3,5,7,11,13,17,19,23,29,31,37,41,43,47>

Vala

Library: Gee

Without any optimizations: <lang vala> using Gee;

ArrayList<int> primes(int limit){ var sieve = new ArrayList<bool>(); var prime_list = new ArrayList<int>();

for(int i = 0; i <= limit; i++){ sieve.add(true); }

sieve[0] = false; sieve[1] = false;

for (int i = 2; i <= limit/2; i++){ if (sieve[i] != false){ for (int j = 2; i*j <= limit; j++){ sieve[i*j] = false; } } }

for (int i = 0; i < sieve.size; i++){ if (sieve[i] != false){ prime_list.add(i); } }

return prime_list; } // end primes

public static void main(){ var prime_list = primes(50);

foreach(var prime in prime_list) stdout.printf("%s ", prime.to_string());

stdout.printf("\n"); } </lang>

Output:

2 3 5 7 11 13 17 19 23 29 31 37 41 43 47

Vedit macro language

This implementation uses an edit buffer as an array for flags. After the macro has been run, you can see how the primes are located in the array. Primes are marked with 'P' and non-primes with '-'. The first character position represents number 0.

#10 = Get_Num("Enter number to search to: ", STATLINE)
Buf_Switch(Buf_Free)                    // Use edit buffer as flags array
Ins_Text("--")                          // 0 and 1 are not primes
Ins_Char('P', COUNT, #10-1)             // init rest of the flags to "prime"
for (#1 = 2; #1*#1 < #10; #1++) {
    Goto_Pos(#1)
    if (Cur_Char=='P') {                // this is a prime number
        for (#2 = #1*#1; #2 <= #10; #2 += #1) {
            Goto_Pos(#2)
            Ins_Char('-', OVERWRITE)
        }
    }
}

Sample output showing numbers in range 0 to 599.

--PP-P-P---P-P---P-P---P-----P-P-----P---P-P---P-----P-----P
-P-----P---P-P-----P---P-----P-------P---P-P---P-P---P------
-------P---P-----P-P---------P-P-----P-----P---P-----P-----P
-P---------P-P---P-P-----------P-----------P---P-P---P-----P
-P---------P-----P-----P-----P-P-----P---P-P---------P------
-------P---P-P---P-------------P-----P---------P-P---P-----P
-------P-----P-----P---P-----P-------P---P-------P---------P
-P---------P-P-----P---P-----P-------P---P-P---P-----------P
-------P---P-------P---P-----P-----------P-P----------------
-P-----P---------P-----P-----P-P-----P---------P-----P-----P

Visual Basic .NET

<lang vbnet>Dim n As Integer, k As Integer, limit As Integer Console.WriteLine("Enter number to search to: ") limit = Console.ReadLine Dim flags(limit) As Integer For n As Integer = 2 To Math.Sqrt(limit)

   If flags(n) = 0 Then
       For k As Integer = 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
       Console.WriteLine(n)
   End If

Next n</lang>

Vorpal

<lang vorpal> self.print_primes = method(m){

  p = new()
  p.fill(0, m, 1, true)

  count = 0
  i = 2
  while(i < m){
     if(p[i] == true){
        p.fill(i+i, m, i, false)
        count = count + 1
     }
     i = i + 1
  }
  ('primes: ' + count + ' in ' + m).print()
  for(i = 2, i < m, i = i + 1){
     if(p[i] == true){
        ( + i + ', ').put()
     }
  }
  .print()

}

self.print_primes(100) </lang> Result:

primes: 25 in 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,