Sieve of Eratosthenes

From Rosetta Code
Jump to: navigation, search
Task
Sieve of Eratosthenes
You are encouraged to solve this task according to the task description, using any language you may know.
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

Contents

[edit] 68000 Assembly

Algorithm somewhat optimized: array omits 1, 2, all higher odd numbers. Optimized for storage: uses bit array for prime/composite flags.

Works with: [EASy68K v5.13.00]

Some of the macro code is derived from the examples included with EASy68K. See 68000 "100 Doors" listing for additional information.

*-----------------------------------------------------------
* Title  : BitSieve
* Written by : G. A. Tippery
* Date  : 2014-Feb-24, 2013-Dec-22
* Description: Prime number sieve
*-----------------------------------------------------------
ORG $1000
 
** ---- Generic macros ---- **
PUSH MACRO
MOVE.L \1,-(SP)
ENDM
 
POP MACRO
MOVE.L (SP)+,\1
ENDM
 
DROP MACRO
ADDQ #4,SP
ENDM
 
PUTS MACRO
** Print a null-terminated string w/o CRLF **
** Usage: PUTS stringaddress
** Returns with D0, A1 modified
MOVEQ #14,D0 ; task number 14 (display null string)
LEA \1,A1 ; address of string
TRAP #15 ; display it
ENDM
 
GETN MACRO
MOVEQ #4,D0 ; Read a number from the keyboard into D1.L.
TRAP #15
ENDM
 
** ---- Application-specific macros ---- **
 
val MACRO ; Used by bit sieve. Converts bit address to the number it represents.
ADD.L \1,\1 ; double it because odd numbers are omitted
ADDQ #3,\1 ; add offset because initial primes (1, 2) are omitted
ENDM
 
* ** ================================================================================ **
* ** Integer square root routine, bisection method **
* ** IN: D0, should be 0<D0<$1000 (65536) -- higher values MAY work, no guarantee
* ** OUT: D1
*
SquareRoot:
*
MOVEM.L D2-D4,-(SP) ; save registers needed for local variables
* DO == n
* D1 == a
* D2 == b
* D3 == guess
* D4 == temp
*
* a = 1;
* b = n;
MOVEQ #1,D1
MOVE.L D0,D2
* do {
REPEAT
* guess = (a+b)/2;
MOVE.L D1,D3
ADD.L D2,D3
LSR.L #1,D3
* if (guess*guess > n) { // inverse function of sqrt is square
MOVE.L D3,D4
MULU D4,D4 ; guess^2
CMP.L D0,D4
BLS .else
* b = guess;
MOVE.L D3,D2
BRA .endif
* } else {
.else:
* a = guess;
MOVE.L D3,D1
* } //if
.endif:
* } while ((b-a) > 1); ; Same as until (b-a)<=1 or until (a-b)>=1
MOVE.L D2,D4
SUB.L D1,D4 ; b-a
UNTIL.L D4 <LE> #1 DO.S
* return (a) ; Result is in D1
* } //LongSqrt()
MOVEM.L (SP)+,D2-D4 ; restore saved registers
RTS
*
* ** ================================================================================ **
 
 
** ======================================================================= **
*
** Prime-number Sieve of Eratosthenes routine using a big bit field for flags **
* Enter with D0 = size of sieve (bit array)
* Prints found primes 10 per line
* Returns # prime found in D6
*
* Register usage:
*
* D0 == n
* D1 == prime
* D2 == sqroot
* D3 == PIndex
* D4 == CIndex
* D5 == MaxIndex
* D6 == PCount
*
* A0 == PMtx[0]
*
* On return, all registers above except D0 are modified. Could add MOVEMs to save and restore D2-D6/A0.
*
 
** ------------------------ **
 
GetBit: ** sub-part of Sieve subroutine **
** Entry: bit # is on TOS
** Exit: A6 holds the byte number, D7 holds the bit number within the byte
** Note: Input param is still on TOS after return. Could have passed via a register, but
** wanted to practice with stack. :)
*
MOVE.L (4,SP),D7 ; get value from (pre-call) TOS
ASR.L #3,D7 ; /8
MOVEA D7,A6 ; byte #
MOVE.L (4,SP),D7 ; get value from (pre-call) TOS
AND.L #$7,D7 ; bit #
RTS
 
** ------------------------ **
 
Sieve:
MOVE D0,D5
SUBQ #1,D5
JSR SquareRoot ; sqrt D0 => D1
MOVE.L D1,D2
LEA PArray,A0
CLR.L D3
*
PrimeLoop:
MOVE.L D3,D1
val D1
MOVE.L D3,D4
ADD.L D1,D4
*
CxLoop: ; Goes through array marking multiples of d1 as composite numbers
CMP.L D5,D4
BHI ExitCx
PUSH D4 ; set D7 as bit # and A6 as byte pointer for D4'th bit of array
JSR GetBit
DROP
BSET D7,0(A0,A6.L) ; set bit to mark as composite number
ADD.L D1,D4 ; next number to mark
BRA CxLoop
ExitCx:
CLR.L D1 ; Clear new-prime-found flag
ADDQ #1,D3 ; Start just past last prime found
PxLoop: ; Searches for next unmarked (not composite) number
CMP.L D2,D3 ; no point searching past where first unmarked multiple would be past end of array
BHI ExitPx ; if past end of array
TST.L D1
BNE ExitPx ; if flag set, new prime found
PUSH D3 ; check D3'th bit flag
JSR GetBit ; sets D7 as bit # and A6 as byte pointer
DROP ; drop TOS
BTST D7,0(A0,A6.L) ; read bit flag
BNE IsSet ; If already tagged as composite
MOVEQ #-1,D1 ; Set flag that we've found a new prime
IsSet:
ADDQ #1,D3 ; next PIndex
BRA PxLoop
ExitPx:
SUBQ #1,D3 ; back up PIndex
TST.L D1 ; Did we find a new prime #?
BNE PrimeLoop ; If another prime # found, go process it
*
; fall through to print routine
 
** ------------------------ **
 
* Print primes found
*
* D4 == Column count
*
* Print header and assumed primes (#1, #2)
PUTS Header ; Print string @ Header, no CR/LF
MOVEQ #2,D6 ; Start counter at 2 because #1 and #2 are assumed primes
MOVEQ #2,D4
*
MOVEQ #0,D3
PrintLoop:
CMP.L D5,D3
BHS ExitPL
PUSH D3
JSR GetBit ; sets D7 as bit # and A6 as byte pointer
DROP ; drop TOS
BTST D7,0(A0,A6.L)
BNE NotPrime
* printf(" %6d", val(PIndex)
MOVE.L D3,D1
val D1
AND.L #$0000FFFF,D1
MOVEQ #6,D2
MOVEQ #20,D0 ; display signed RJ
TRAP #15
ADDQ #1,D4
ADDQ #1,D6
* *** Display formatting ***
* if((PCount % 10) == 0) printf("\n");
CMP #10,D4
BLO NoLF
PUTS CRLF
MOVEQ #0,D4
NoLF:
NotPrime:
ADDQ #1,D3
BRA PrintLoop
ExitPL:
RTS
 
** ======================================================================= **
 
N EQU 5000 ; *** Size of boolean (bit) array ***
SizeInBytes EQU (N+7)/8
*
START: ; first instruction of program
MOVE.L #N,D0 ; # to test
JSR Sieve
* printf("\n %d prime numbers found.\n", D6); ***
PUTS Summary1,A1
MOVE #3,D0 ; Display signed number in D1.L in decimal in smallest field.
MOVE.W D6,D1
TRAP #15
PUTS Summary2,A1
 
SIMHALT ; halt simulator
 
** ======================================================================= **
 
* Variables and constants here
 
ORG $2000
CR EQU 13
LF EQU 10
CRLF DC.B CR,LF,$00
 
PArray: DCB.B SizeInBytes,0
 
Header: DC.B CR,LF,LF,' Primes',CR,LF,' ======',CR,LF
DC.B ' 1 2',$00
 
Summary1: DC.B CR,LF,' ',$00
Summary2: DC.B ' prime numbers found.',CR,LF,$00
 
END START ; last line of source

[edit] ACL2

(defun nats-to-from (n i)
(declare (xargs :measure (nfix (- n i))))
(if (zp (- n i))
nil
(cons i (nats-to-from n (+ i 1)))))
 
(defun remove-multiples-up-to-r (factor limit xs i)
(declare (xargs :measure (nfix (- limit i))))
(if (or (> i limit)
(zp (- limit i))
(zp factor))
xs
(remove-multiples-up-to-r
factor
limit
(remove i xs)
(+ i factor))))
 
(defun remove-multiples-up-to (factor limit xs)
(remove-multiples-up-to-r factor limit xs (* factor 2)))
 
(defun sieve-r (factor limit)
(declare (xargs :measure (nfix (- limit factor))))
(if (zp (- limit factor))
(nats-to-from limit 2)
(remove-multiples-up-to factor (+ limit 1)
(sieve-r (1+ factor) limit))))
 
(defun sieve (limit)
(sieve-r 2 limit))

[edit] Ada

with Ada.Text_IO, Ada.Command_Line;
 
procedure Eratos is
 
Last: Positive := Positive'Value(Ada.Command_Line.Argument(1));
Prime: array(1 .. Last) of Boolean := (1 => False, others => True);
Base: Positive := 2;
Cnt: Positive;
begin
loop
exit when Base * Base > Last;
if Prime(Base) then
Cnt := Base + Base;
loop
exit when Cnt > Last;
Prime(Cnt) := False;
Cnt := Cnt + Base;
end loop;
end if;
Base := Base + 1;
end loop;
Ada.Text_IO.Put("Primes less or equal" & Positive'Image(Last) &" are:");
for Number in Prime'Range loop
if Prime(Number) then
Ada.Text_IO.Put(Positive'Image(Number));
end if;
end loop;
end Eratos;

Output:

> ./eratos 31
Primes less or equal 31 are : 2 3 5 7 11 13 17 19 23 29 31

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

Search autohotkey.com: of Eratosthenes
Source: AutoHotkey forum by Laszlo

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
}

[edit] AutoIt

#include <Array.au3>
$M = InputBox("Integer", "Enter biggest Integer")
Global $a[$M], $r[$M], $c = 1
For $i = 2 To $M -1
If Not $a[$i] Then
$r[$c] = $i
$c += 1
For $k = $i To $M -1 Step $i
$a[$k] = True
Next
EndIf
Next
$r[0] = $c - 1
ReDim $r[$c]
_ArrayDisplay($r)

[edit] 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. Here, the script is entered directly on the commandline, and input entered on stdin:

$ 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

Now with the script from a file, input from commandline as well as stdin, and input is checked for valid numbers:

 
# usage: gawk -v n=101 -f sieve.awk
 
function sieve(n) { # print n,":"
for(i=2; i<=n; i++) a[i]=1;
for(i=2; i<=sqrt(n);i++) for(j=2;j<=n;j++) a[i*j]=0;
for(i=2; i<=n; i++) if(a[i]) printf i" "
print ""
}
 
BEGIN { print "Sieve of Eratosthenes:"
if(n>1) sieve(n)
}
 
{ n=$1+0 }
n<2 { exit }
{ sieve(n) }
 
END { print "Bye!" }
 

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

[edit] Locomotive Basic

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

[edit] ZX Spectrum Basic

10 INPUT "Enter number to search to: ";l
20 DIM p(l)
30 FOR n=2 TO SQR l
40 IF p(n)<>0 THEN NEXT n
50 FOR k=n*n TO l STEP n
60 LET p(k)=1
70 NEXT k
80 NEXT n
90 REM Display the primes
100 FOR n=2 TO l
110 IF p(n)=0 THEN PRINT n;", ";
120 NEXT n

[edit] BBC BASIC

      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

[edit] bash

See solutions at UNIX Shell.

[edit] Befunge

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

[edit] Bracmat

This solution does not use an array. Instead, numbers themselves are used as variables. The numbers that are not prime are set (to the silly value "nonprime"). Finally all numbers up to the limit are tested for being initialised. The uninitialised (unset) ones must be the primes.

( ( eratosthenes
= n j i
.  !arg:?n
& 1:?i
& whl
' ( (1+!i:?i)^2:~>!n:?j
& ( !!i
| whl
' ( !j:~>!n
& nonprime:?!j
& !j+!i:?j
)
)
)
& 1:?i
& whl
' ( 1+!i:~>!n:?i
& (!!i|put$(!i " "))
)
)
& eratosthenes$100
)
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

[edit] C

Plain sieve, without any optimizations:
#include <stdlib.h>
#include <math.h>
 
char*
eratosthenes(int n, int *c)
{
char* sieve;
int i, j, m;
 
if(n < 2)
return NULL;
 
*c = n-1; /* primes count */
m = (int) sqrt((double) n);
 
/* calloc initializes to zero */
sieve = calloc(n+1,sizeof(char));
sieve[0] = 1;
sieve[1] = 1;
for(i = 2; i <= m; i++)
if(!sieve[i])
for (j = i*i; j <= n; j += i)
if(!sieve[j]){
sieve[j] = 1;
--(*c);
}
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.

Another example:

We first fill ones into an array and assume all numbers are prime. Then, in a loop, fill zeroes into those places where i * j is less than or equal to n (number of primes requested), which means they have multiples! To understand this better, look at the output of the following example. To print this back, we look for ones in the array and only print those spots.
#include <stdio.h>
#include <malloc.h>
void sieve(int *, int);
 
int main(int argc, char *argv)
{
int *array, n=10;
array =(int *)malloc(sizeof(int));
sieve(array,n);
return 0;
}
 
void sieve(int *a, int n)
{
int i=0, j=0;
 
for(i=2; i<=n; i++) {
a[i] = 1;
}
 
for(i=2; i<=n; i++) {
printf("\ni:%d", i);
if(a[i] == 1) {
for(j=i; (i*j)<=n; j++) {
printf ("\nj:%d", j);
printf("\nBefore a[%d*%d]: %d", i, j, a[i*j]);
a[(i*j)] = 0;
printf("\nAfter a[%d*%d]: %d", i, j, a[i*j]);
}
}
}
 
printf("\nPrimes numbers from 1 to %d are : ", n);
for(i=2; i<=n; i++) {
if(a[i] == 1)
printf("%d, ", i);
}
printf("\n\n");
}
Output:
i:2
j:2
Before a[2*2]: 1
After a[2*2]: 0
j:3
Before a[2*3]: 1
After a[2*3]: 0
j:4
Before a[2*4]: 1
After a[2*4]: 0
j:5
Before a[2*5]: 1
After a[2*5]: 0
i:3
j:3
Before a[3*3]: 1
After a[3*3]: 0
i:4
i:5
i:6
i:7
i:8
i:9
i:10
Primes numbers from 1 to 10 are : 2, 3, 5, 7,

[edit] C++

// 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);
}

Full program:

Works with: Boost
/**
$ 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] C#

Works with: C# version 3.0+
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;
}
}
}

[edit] Chapel

This solution uses nested iterators to create new wheels at run time:

// yield prime and remove all multiples of it from children sieves
iter sieve(prime):int {
 
yield prime;
 
var last = prime;
label candidates for candidate in sieve(prime+1) do {
for composite in last..candidate by prime do {
 
// candidate is a multiple of this prime
if composite == candidate then {
// remember size of last composite
last = composite;
// and try the next candidate
continue candidates;
}
}
 
// candidate cannot need to be removed by this sieve
// yield to parent sieve for checking
yield candidate;
}
}
The topmost sieve needs to be started with 2 (the smallest prime):
config const N = 30;
for p in sieve(2) {
write(" ", p);
if p > N then {
writeln();
break;
}
}

[edit] Clojure

Calculates primes up to and including n.
(defn primes-to
"Returns a list of all primes from 2 to n"
[n]
(let [n (int n)]
(let [root (-> n Math/sqrt int)]
(loop [i (int 2), a (boolean-array (inc n)), result (transient [])]
(if (> i n)
(persistent! result)
(recur (inc i)
(if (and (<= i root) (not (aget a i)))
(loop [arr a, j (* i i)]
(if (> j n)
arr
(recur (do (aset arr j true) arr)
(+ j i))))
a)
(if (not (aget a i))
(conj! result i)
result)))))))

This code is not very idiomatic, and it looks like something written in an imperative language.

The following codes don't work on the latest Clojure compiler and need fixing!

Here is a non-lazy version:
(defn primes [max-prime]
(let [sieve (fn [s n]
(if (<= (* n n) max-prime)
(recur (if (s n)
(reduce #(assoc %1 %2 false) s (range (* n n) (inc max-prime) n))
s)
(inc n))
s))]
(-> (vector-of :boolean) ; form the return vector
(reduce conj (range (inc max-prime)))
(assoc 0 false)
(assoc 1 false)
(sieve 2)
#(->> % count range (map vector %) (filter first) (map second)))))

Or as a lazy sequence (slightly faster):

(defn primes
"Computes all prime numbers up to a given number using sieve of Eratosthenes"
[n]
(loop [cs (range 2 number) ; candidates
ps [2]] ; results
(let [lp (last primes)
ncs (-> (range lp n lp) set (remove cs))]
(if (> lp (Math/sqrt n))
(concat ps ncs)
(recur ncs (concat ps [(first ncs)])))))))

[edit] 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)
# Print all prime numbers through 100.
eratosthenes(primes 100)
message(STATUS "${primes}")

[edit] COBOL

*> Please ignore the asterisks in the first column of the next comments,
*> which are kludges to get syntax highlighting to work.
IDENTIFICATION DIVISION.
PROGRAM-ID. Sieve-Of-Eratosthenes.
 
DATA DIVISION.
WORKING-STORAGE SECTION.
 
01 Max-Number USAGE UNSIGNED-INT.
01 Max-Prime USAGE UNSIGNED-INT.
 
01 Num-Group.
03 Num-Table PIC X VALUE "P"
OCCURS 1 TO 10000000 TIMES DEPENDING ON Max-Number
INDEXED BY Num-Index.
88 Is-Prime VALUE "P" FALSE "N".
 
01 Current-Prime USAGE UNSIGNED-INT.
 
01 I USAGE UNSIGNED-INT.
 
PROCEDURE DIVISION.
DISPLAY "Enter the limit: " WITH NO ADVANCING
ACCEPT Max-Number
DIVIDE Max-Number BY 2 GIVING Max-Prime
 
* *> Set Is-Prime of all non-prime numbers to false.
SET Is-Prime (1) TO FALSE
PERFORM UNTIL Max-Prime < Current-Prime
* *> Set current-prime to next prime.
ADD 1 TO Current-Prime
PERFORM VARYING Num-Index FROM Current-Prime BY 1
UNTIL Is-Prime (Num-Index)
END-PERFORM
MOVE Num-Index TO Current-Prime
 
* *> Set Is-Prime of all multiples of current-prime to
* *> false, starting from current-prime sqaured.
COMPUTE Num-Index = Current-Prime ** 2
PERFORM UNTIL Max-Number < Num-Index
SET Is-Prime (Num-Index) TO FALSE
SET Num-Index UP BY Current-Prime
END-PERFORM
END-PERFORM
 
* *> Display the prime numbers.
PERFORM VARYING Num-Index FROM 1 BY 1
UNTIL Max-Number < Num-Index
IF Is-Prime (Num-Index)
DISPLAY Num-Index
END-IF
END-PERFORM
 
GOBACK
.

[edit] Common 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)))))

Working with odds only (above twice speedup), and only test divide up to the square root of the maximum:

(defun sieve-odds (maximum) "sieve for odd numbers"
(cons 2
(let ((maxi (ash (1- maximum) -1)) (stop (ash (isqrt maximum) -1)))
(let ((sieve (make-array (1+ maxi) :element-type 'bit :initial-element 0)))
(loop for i from 1 to maxi
when (zerop (sbit sieve i))
collect (1+ (ash i 1))
and when (<= i stop) do
(loop for j from (ash (* i (1+ i)) 1) to maxi by (1+ (ash i 1))
do (setf (sbit sieve j) 1)))))))

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.

[edit] D

[edit] Simpler Version

Prints all numbers less than the limit.
import std.stdio, std.algorithm, std.range, std.functional;
 
uint[] sieve(in uint limit) nothrow @safe {
if (limit < 2)
return [];
auto composite = new bool[limit];
 
foreach (immutable uint n; 2 .. cast(uint)(limit ^^ 0.5) + 1)
if (!composite[n])
for (uint k = n * n; k < limit; k += n)
composite[k] = true;
 
//return iota(2, limit).filter!(not!composite).array;
return iota(2, limit).filter!(i => !composite[i]).array;
}
 
void main() {
50.sieve.writeln;
}
Output:
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]

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

import std.stdio, std.math, std.array;
 
size_t[] sieve(in size_t m) pure nothrow @safe {
if (m < 3)
return null;
immutable size_t n = m - 1;
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 @safe @nogc {
immutable size_t offset = i / bpc;
immutable size_t mask = 1 << (i % bpc);
return F[offset] & mask;
}
 
void resetBit(in size_t i) nothrow @safe @nogc {
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 <= sqrt(real(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 ~= 2;
for (size_t i = 3; i <= n; i += 2)
if (isSet((i - 3) / 2))
result ~= i;
return result.data;
}
 
void main() {
50.sieve.writeln;
}

[edit] Extensible Version

(This version is used in the task Extensible prime generator.)

/// Extensible Sieve of Eratosthenes.
struct Prime {
uint[] a = [2];
 
private void grow() pure nothrow @safe {
immutable p0 = a[$ - 1] + 1;
auto b = new bool[p0];
 
foreach (immutable di; a) {
immutable uint i0 = p0 / di * di;
uint i = (i0 < p0) ? i0 + di - p0 : i0 - p0;
for (; i < b.length; i += di)
b[i] = true;
}
 
foreach (immutable uint i, immutable bi; b)
if (!b[i])
a ~= p0 + i;
}
 
uint opCall(in uint n) pure nothrow @safe {
while (n >= a.length)
grow;
return a[n];
}
}
 
version (sieve_of_eratosthenes3_main) {
void main() {
import std.stdio, std.range, std.algorithm;
 
Prime prime;
uint.max.iota.map!prime.until!q{a > 50}.writeln;
}
}

To see the output (that is the same), compile with -version=sieve_of_eratosthenes3_main.

[edit] Dart

// helper function to pretty print an Iterable
String iterableToString(Iterable seq) {
String str = "[";
Iterator i = seq.iterator;
if (i.moveNext()) str += i.current.toString();
while(i.moveNext()) {
str += ", " + i.current.toString();
}
return str + "]";
}
 
main() {
int limit = 1000;
int strt = new DateTime.now().millisecondsSinceEpoch;
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);
int elpsd = new DateTime.now().millisecondsSinceEpoch - strt;
print("Found " + sieve.length.toString() + " primes up to " + limit.toString() +
" in " + elpsd.toString() + " milliseconds.");
print(iterableToString(sortedValues)); // expect sieve.length to be 168 up to 1000...
// Expect.equals(168, sieve.length);
}
Output:

Found 168 primes up to 1000 in 9 milliseconds. [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]

Although it has the characteristics of a true Sieve of Eratosthenes, the above code isn't very efficient due to the remove/modify operations on the Set. Due to these, the computational complexity isn't close to linear with increasing range and it is quite slow for larger sieve ranges compared to compiled languages, taking about four seconds to sieve to ten million.

[edit] faster bit-packed array odds-only solution

import 'dart:math';
 
List<int> SoEOdds(int limit) {
List<int> prms = new List();
if (limit < 2) return prms;
prms.add(2);
if (limit < 3) return prms;
int lmt = (limit - 3) >> 1;
int bfsz = (lmt >> 5) + 1;
int sqrtlmt = (sqrt(limit) - 3).floor() >> 1;
var buf = new List<int>();
for (int i = 0; i < bfsz; i++)
buf.add(0);
for (int i = 0; i <= sqrtlmt; i++)
if ((buf[i >> 5] & (1 << (i & 31))) == 0) {
int p = i + i + 3;
for (int j = (p * p - 3) >> 1; j <= lmt; j += p)
buf[j >> 5] |= 1 << (j & 31);
}
for (int i = 0; i <= lmt; i++)
if ((buf[i >> 5] & (1 << (i & 31))) == 0)
prms.add(i + i + 3);
return prms;
}
 
void main() {
int limit = 10000000;
int strt = new DateTime.now().millisecondsSinceEpoch;
List<int> primes = SoEOdds(limit);
int count = primes.length;
int elpsd = new DateTime.now().millisecondsSinceEpoch - strt;
print("Found " + count.toString() + " primes up to " + limit.toString() +
" in " + elpsd.toString() + " milliseconds.");
// print(iterableToString(primes)); // expect sieve.length to be 168 up to 1000...
}

The above code is somewhat faster at about ten seconds using the Dart VM to sieve to 100 million, although much faster at about 1.5 seconds run conventionally in Google Chrome using the JavaScript V8 engine, likely due to JavaScript using double floating point numbers for int's whereas the Dart VM uses arbitrary precision integers.

[edit] fast page segmented array infinite iterator (sieves odds-only)

Translation of: JavaScript
import 'dart:collection';
 
class _SoEPagedIterator implements Iterator<int> {
static const int _BFSZ = 1 << 16;
static const int _BFBTS = _BFSZ * 32;
static const int _BFRNG = _BFBTS * 2;
int _prime = null;
int _bi = -1;
int _lowi = 0;
List<int> _bpa = new List<int>();
Iterator<int> _bps;
List<int> _buf = new List<int>();
int get current => this._prime;
bool moveNext() {
// the following redundant local variable declaration is necessary to
// prevent the dart2js compiler from "tree-shaking" and eliminating some
// essential code from the below, which doesn't happen with the Dart VM compiler.
int lowi = this._lowi;
while (true) {
if (this._bi < 1) {
if (this._bi < 0) { this._bi++; this._prime = 2; break; }
int nxt = 3 + (this._lowi << 1) + _BFRNG;
this._buf.clear();
for (int i = 0; i < _BFSZ; i++) this._buf.add(0); // faster initialization:
if (lowi <= 0) { // special culling for first page as no base primes yet:
for (int i = 0, p = 3, sqr = 9; sqr < nxt; i++, p += 2, sqr = p * p)
if ((this._buf[i >> 5] & (1 << (i & 31))) == 0)
for (int j = (sqr - 3) >> 1; j < _BFBTS; j += p)
this._buf[j >> 5] |= 1 << (j & 31);
} else { // after the first page:
if (this._bpa.length == 0) { // if this is the first page after the zero one:
this._bps = new _SoEPagedIterator(); // initialize separate base primes stream:
this._bps.moveNext(); // advance to the only even prime of two
this._bps.moveNext(); // advance past 2 to the next prime of 3
}
// get enough base primes for the page range...
for (var lp = this._bps.current, sqr = lp * lp; sqr < nxt;
this._bps.moveNext(), lp = this._bps.current, sqr = lp * lp) this._bpa.add(lp);
for (var i = 0; i < this._bpa.length; i++) {
int p = this._bpa[i];
int s = (p * p - 3) >> 1;
if (s >= this._lowi) // adjust start index based on page lower limit...
s -= this._lowi;
else {
int r = (this._lowi - s) % p;
s = (r != 0) ? p - r : 0;
}
for (var j = s; j < _BFBTS; j += p)
this._buf[j >> 5] |= 1 << (j & 31);
}
}
}
while (this._bi < _BFBTS && ((this._buf[this._bi >> 5] & (1 << (this._bi & 31))) != 0))
this._bi++; // find next marker still with prime status
if (this._bi < _BFBTS) { // within buffer: output computed prime
this._prime = 3 + ((this._lowi + this._bi++) << 1); break; }
else { // beyond buffer range: advance buffer
this._bi = 0;
this._lowi += _BFBTS;
lowi = this._lowi;
}
} return true;
}
}
 
class SoEPagedOddsInfGen extends IterableBase<int> {
Iterator<int> get iterator { return new _SoEPagedIterator(); }
}
 
void main() {
int n = 1000000000;
int strt = new DateTime.now().millisecondsSinceEpoch;
int count = new SoEPagedOddsInfGen().takeWhile((p) => p <= n).length;
int elpsd = new DateTime.now().millisecondsSinceEpoch - strt;
print("For a range of " + n.toString() + ", " + count.toString() +
" primes found in " + elpsd.toString() + " milliseconds.");
}

This version calculates the 50,847,534 primes up to one billion in about 20 seconds under the Dart Virtual Machine (VM). Under the Google Chrome V8 JavaScript engine it should take the same time as the JavaScript from which it was translated of about five seconds, but takes about 14 seconds due to the dart2js compiler adding extra run time array buffer range checks to the innermost culling loops, even though the "check" compiler option was not selected.

Also note the comment at the beginning of the "moveNext()" method about the redundant local variable needed to be added in order for the code to run under JavaScript using Dart 1.5.1 (and possible other versions), which shouldn't happen when it runs fine under the Dart VM without that extra local variable (based only on the private class field _lowi).

[edit] Delphi

program erathostenes;
 
{$APPTYPE CONSOLE}
 
type
TSieve = class
private
fPrimes: TArray<boolean>;
procedure InitArray;
procedure Sieve;
function getNextPrime(aStart: integer): integer;
function getPrimeArray: TArray<integer>;
public
function getPrimes(aMax: integer): TArray<integer>;
end;
 
{ TSieve }
 
function TSieve.getNextPrime(aStart: integer): integer;
begin
result := aStart;
while not fPrimes[result] do
inc(result);
end;
 
function TSieve.getPrimeArray: TArray<integer>;
var
i, n: integer;
begin
n := 0;
setlength(result, length(fPrimes)); // init array with maximum elements
for i := 2 to high(fPrimes) do
begin
if fPrimes[i] then
begin
result[n] := i;
inc(n);
end;
end;
setlength(result, n); // reduce array to actual elements
end;
 
function TSieve.getPrimes(aMax: integer): TArray<integer>;
begin
setlength(fPrimes, aMax);
InitArray;
Sieve;
result := getPrimeArray;
end;
 
procedure TSieve.InitArray;
begin
for i := 2 to high(fPrimes) do
fPrimes[i] := true;
end;
 
procedure TSieve.Sieve;
var
i, n, max: integer;
begin
max := length(fPrimes);
i := 2;
while i < sqrt(max) do
begin
n := sqr(i);
while n < max do
begin
fPrimes[n] := false;
inc(n, i);
end;
i := getNextPrime(i + 1);
end;
end;
 
var
i: integer;
Sieve: TSieve;
 
begin
Sieve := TSieve.Create;
for i in Sieve.getPrimes(100) do
write(i, ' ');
Sieve.Free;
readln;
end.

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 


[edit] DWScript

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]);

[edit] Dylan

With outer to sqrt and inner to p^2 optimizations:

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;


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

This example is incorrect. It uses rem testing and so is a trial division algorithm, not a sieve of Eratosthenes. Please fix the code and remove this message.

Note: this is not a Sieve of Eratosthenes; it is just trial division.

 
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));
}
}
 

[edit] Eiffel

Works with: EiffelStudio version 6.6 beta (with provisional loop syntax)
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

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


[edit] Emacs Lisp

 
(defun sieve-set (limit)
(let ((xs (make-vector (1+ limit) 0)))
(loop for i from 2 to limit
when (zerop (aref xs i))
collect i
and do (loop for m from (* i i) to limit by i
do (aset xs m 1)))))
 

Straightforward implementation of sieve of Eratosthenes, 2 times faster:

 
(defun sieve (limit)
(let ((xs (vconcat [0 0] (number-sequence 2 limit))))
(loop for i from 2 to (sqrt limit)
when (aref xs i)
do (loop for m from (* i i) to limit by i
do (aset xs m 0)))
(remove 0 xs)))
 

[edit] Erlang

 
-module( sieve_of_eratosthenes ).
 
-export( [primes_upto/1] ).
 
primes_upto( N ) ->
Ns = lists:seq( 2, N ),
Dict = dict:from_list( [{X, potential_prime} || X <- Ns] ),
{Upto_sqrt_ns, _T} = lists:split( erlang:round(math:sqrt(N)), Ns ),
{N, Prime_dict} = lists:foldl( fun find_prime/2, {N, Dict}, Upto_sqrt_ns ),
lists:sort( dict:fetch_keys(Prime_dict) ).
 
 
 
find_prime( N, {Max, Dict} ) -> find_prime( dict:find(N, Dict), N, {Max, Dict} ).
 
find_prime( error, _N, Acc ) -> Acc;
find_prime( {ok, _Value}, N, {Max, Dict} ) -> {Max, lists:foldl( fun dict:erase/2, Dict, lists:seq(N*N, Max, N) )}.
 
Output:
35> sieve_of_eratosthenes:primes_upto( 20 ).
[2,3,5,7,11,13,17,19]

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

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}

[edit] F#

[edit] Functional

This is the idea behind Richard Bird's unbounded code presented in M. O'Neill's article in Haskell. It is about twice as much code as the Haskell code because F# does not have a built-in lazy list so that the effect must be constructed using a Co-Inductive Stream type and the use of recursive functions in combination with sequences:

type CIS<'T> = struct val v:'T val cont:unit->CIS<'T> //Co Inductive Stream for laziness
new (v,cont) = { v = v; cont = cont } end
let primes =
let rec pculls p cull = CIS(cull, fun() -> pculls p (cull + p))
let rec allculls (ps:CIS<_>) = //stream of streams of composite culls
CIS(pculls ps.v (ps.v * ps.v),fun() -> allculls (ps.cont()))
let rec (^^) (xs:CIS<uint32>) (ys:CIS<uint32>) = //union op for CIS<uint32>'
s
match compare xs.v ys.v with
| -1 -> CIS(xs.v, fun() -> xs.cont() ^^ ys) // <
| 0 -> CIS(xs.v, fun() -> xs.cont() ^^ ys.cont()) // ==
| _ -> CIS(ys.v, fun() -> xs ^^ ys.cont()) //must be > (= 1)
let rec join (cmpsts:CIS<CIS<_>>) =
CIS(cmpsts.v.v, fun() -> cmpsts.v.cont() ^^ join (cmpsts.cont()))
let rec mkPrms cnd (cmpsts:CIS<_>) =
let ncnd = cnd + 1u
if cnd >= cmpsts.v then mkPrms ncnd (cmpsts.cont()) //implements 'minus'
else CIS(cnd,fun()->mkPrms ncnd cmpsts) //found a prime
let rec basePrimes = CIS(2u, fun() -> mkPrms 3u initCmpsts)
and initCmpsts = join (allculls (basePrimes))
let genseq cis = Seq.unfold (fun (cs:CIS<_>) -> Some(cs.v, cs.cont())) cis
genseq (mkPrms 2u initCmpsts)

The above code sieves all numbers of two and up including all even numbers as per the page specification; the following code makes the very minor changes for an odds-only sieve, with a speedup of over a factor of two:

type CIS<'T> = struct val v:'T val cont:unit->CIS<'T> //Co Inductive Stream for laziness
new (v,cont) = { v = v; cont = cont } end
let primes =
let rec pculls p cull = CIS(cull, fun() -> pculls p (cull + 2u * p))
let rec allculls (ps:CIS<_>) = //stream of streams of composite culls
CIS(pculls ps.v (ps.v * ps.v),fun() -> allculls (ps.cont()))
let rec (^^) (xs:CIS<uint32>) (ys:CIS<uint32>) = //union op for CIS<uint32>'
s
match compare xs.v ys.v with
| -1 -> CIS(xs.v, fun() -> xs.cont() ^^ ys) // <
| 0 -> CIS(xs.v, fun() -> xs.cont() ^^ ys.cont()) // ==
| _ -> CIS(ys.v, fun() -> xs ^^ ys.cont()) //must be > (= 1)
let rec join (cmpsts:CIS<CIS<_>>) =
CIS(cmpsts.v.v, fun() -> cmpsts.v.cont() ^^ join (cmpsts.cont()))
let rec mkPrms cnd (cmpsts:CIS<_>) =
let ncnd = cnd + 2u
if cnd >= cmpsts.v then mkPrms ncnd (cmpsts.cont()) //implements 'minus'
else CIS(cnd,fun()->mkPrms ncnd cmpsts) //found a prime
let rec oddBasePrimes = CIS(3u, fun() -> mkPrms 5u initCmpsts)
and initCmpsts = join (allculls (oddBasePrimes))
let genseq cis = Seq.unfold (fun (cs:CIS<_>) -> Some(cs.v, cs.cont())) cis
seq { yield 2u; yield! genseq (mkPrms 3u initCmpsts) }

[edit] Imperative

The following code is written in functional style other than it uses a mutable bit array to sieve the composites:

let primes limit =
let buf = System.Collections.BitArray(int limit + 1, true)
let cull p = { p * p .. p .. limit } |> Seq.iter (fun c -> buf.[int c] <- false)
{ 2u .. uint32 (sqrt (double limit)) } |> Seq.iter (fun c -> if buf.[int c] then cull c)
{ 2u .. limit } |> Seq.map (fun i -> if buf.[int i] then i else 0u) |> Seq.filter ((<>) 0u)
 
[<EntryPoint>]
let main argv =
if argv = null || argv.Length = 0 then failwith "no command line argument for limit!!!"
printfn "%A" (primes (System.UInt32.Parse argv.[0]) |> Seq.length)
0 // return an integer exit code

Substituting the following minor changes to the code for the "primes" function will only deal with the odd prime candidates for a speed up of over a factor of two as well as a reduction of the buffer size by a factor of two:

let primes limit =
let lmtb,lmtbsqrt = (limit - 3u) / 2u, (uint32 (sqrt (double limit)) - 3u) / 2u
let buf = System.Collections.BitArray(int lmtb + 1, true)
let cull i = let p = i + i + 3u in let s = p * (i + 1u) + i in
{ s .. p .. lmtb } |> Seq.iter (fun c -> buf.[int c] <- false)
{ 0u .. lmtbsqrt } |> Seq.iter (fun i -> if buf.[int i] then cull i )
let oddprms = { 0u .. lmtb } |> Seq.map (fun i -> if buf.[int i] then i + i + 3u else 0u)
|> Seq.filter ((<>) 0u)
seq { yield 2u; yield! oddprms }

The following code uses other functional forms for the inner culling loops of the "primes function" to reduce the use of inefficient sequences so as to reduce the execution time by another factor of almost three:

let primes limit =
let lmtb,lmtbsqrt = (limit - 3u) / 2u, (uint32 (sqrt (double limit)) - 3u) / 2u
let buf = System.Collections.BitArray(int lmtb + 1, true)
let rec culltest i = if i <= lmtbsqrt then
let p = i + i + 3u in let s = p * (i + 1u) + i in
let rec cullp c = if c <= lmtb then buf.[int c] <- false; cullp (c + p)
(if buf.[int i] then cullp s); culltest (i + 1u) in culltest 0u
seq {yield 2u; for i = 0u to lmtb do if buf.[int i] then yield i + i + 3u }

Now much of the remaining execution time is just the time to enumerate the primes as can be seen by turning "primes" into a primes counting function by substituting the following for the last line in the above code doing the enumeration; this makes the code run about a further five times faster:

  let rec count i acc =
if i > int lmtb then acc else if buf.[i] then count (i + 1) (acc + 1) else count (i + 1) acc
count 0 1

Since the final enumeration of primes is the main remaining bottleneck, it is worth using a "roll-your-own" enumeration implemented as an object expression so as to save many inefficiencies in the use of the built-in seq computational expression by substituting the following code for the last line of the previous codes, which will decrease the execution time by a factor of over three (instead of almost five for the counting-only version, making it almost as fast):

  let nmrtr() =
let i = ref -2
let rec nxti() = i:=!i + 1;if !i <= int lmtb && not buf.[!i] then nxti() else !i <= int lmtb
let inline curr() = if !i < 0 then (if !i= -1 then 2u else failwith "Enumeration not started!!!")
else let v = uint32 !i in v + v + 3u
{ new System.Collections.Generic.IEnumerator<_> with
member this.Current = curr()
interface System.Collections.IEnumerator with
member this.Current = box (curr())
member this.MoveNext() = if !i< -1 then i:=!i+1;true else nxti()
member this.Reset() = failwith "IEnumerator.Reset() not implemented!!!"a
interface System.IDisposable with
member this.Dispose() = () }
{ new System.Collections.Generic.IEnumerable<_> with
member this.GetEnumerator() = nmrtr()
interface System.Collections.IEnumerable with
member this.GetEnumerator() = nmrtr() :> System.Collections.IEnumerator }

The various optimization techniques shown here can be used "jointly and severally" on any of the basic versions for various trade-offs between code complexity and performance. Not shown here are other techniques of making the sieve faster, including extending wheel factorization to much larger wheels such as 2/3/5/7, pre-culling the arrays, page segmentation, and multi-processing.

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

[edit] Fortran

Works with: Fortran version 90 and later
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

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

Optimised using a pre-computed wheel based on 2:

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

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

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

[edit] 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
NEXT
 
// Display the primes
FOR n = 2 TO limit
IF flags[n] = 0 THEN STDOUT n + ", "
NEXT
 
KEYWAIT
 

[edit] 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("")
}
}
}

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  .

A fairly odd sieve tree method:

package main
import "fmt"
 
type xint uint64
type xgen func()(xint)
 
func primes() func()(xint) {
pp, psq := make([]xint, 0), xint(25)
 
var sieve func(xint, xint)xgen
sieve = func(p, n xint) xgen {
m, next := xint(0), xgen(nil)
return func()(r xint) {
if next == nil {
r = n
if r <= psq {
n += p
return
}
 
next = sieve(pp[0] * 2, psq) // chain in
pp = pp[1:]
psq = pp[0] * pp[0]
 
m = next()
}
switch {
case n < m: r, n = n, n + p
case n > m: r, m = m, next()
default: r, n, m = n, n + p, next()
}
return
}
}
 
f := sieve(6, 9)
n, p := f(), xint(0)
 
return func()(xint) {
switch {
case p < 2: p = 2
case p < 3: p = 3
default:
for p += 2; p == n; {
p += 2
if p > n {
n = f()
}
}
pp = append(pp, p)
}
return p
}
}
 
 
func main() {
for i, p := 0, primes(); i < 100000; i++ {
fmt.Println(p())
}
}

See also the concurrent prime sieve example in the "Try Go" window at http://golang.org/

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

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

Test:

println sievePrimes(100)

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]

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

[edit] Haskell

[edit] Mutable unboxed arrays, odds only

Mutable array of unboxed Bools indexed by Ints, representing odds only:

import Control.Monad (forM_, when)
import Control.Monad.ST
import Data.Array.ST
import Data.Array.Unboxed
 
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 -- prime(i) = 2i+1
isPrime <- readArray sieve i -- ((2i+1)^2-1)`div`2 = 2i(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*i + 1 | (i,True) <- assocs $ sieveUO top]
| otherwise = []

This represents odds only in the array. Empirical orders of growth is O(n1.2) in n primes produced, and improving for bigger n‍ ‍s. Memory consumption is low (array seems to be packed) and growing about linearly with n.

[edit] Immutable arrays

This is the most straightforward:

import Data.Array.Unboxed
 
primesToA m = sieve 3 (array (3,m) [(i,odd i) | i<-[3..m]] :: UArray Int Bool)
where
sieve p a
| p*p > m = 2 : [i | (i,True) <- assocs a]
| a!p = sieve (p+2) $ a//[(i,False) | i <- [p*p, p*p+2*p..m]]
| otherwise = sieve (p+2) a

Its performance sharply depends on compiler optimizations. Compiled with -O2 flag in the presence of the explicit type signature, it is very fast in producing first few million primes. (//) is an array update operator.

[edit] Immutable arrays, by segments

Works by segments between consecutive primes' squares. Should be the fastest of non-monadic code:

import Data.Array.Unboxed
 
primesSA = 2 : prs ()
where
prs () = 3 : sieve 3 [] (prs ())
sieve x fs (p:ps) = [i*2 + x | (i,True) <- assocs a]
++ sieve (p*p) fs2 ps
where
q = (p*p-x)`div`2
fs2 = (p,0) : [(s, rem (y-q) s) | (s,y) <- fs]
a :: UArray Int Bool
a = accumArray (\ b c -> False) True (1,q-1)
[(i,()) | (s,y) <- fs, i <- [y+s, y+s+s..q]]

[edit] Basic list-based sieve

Straightforward implementation of the sieve of Eratosthenes in its original bounded form. This finds primes in gaps between the composites, and composites as an enumeration of each prime's multiples.

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..m])
 
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

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.

[edit] Unbounded list based 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.

primes = _Y $ ((2:) . minus [3..] 
. foldr (\x-> (x*x :) . union [x*x+x, x*x+2*x..]) [])
 
_Y g = g (_Y g) -- non-sharing multistage fixpoint combinator
-- = let x = g x in g x -- sharing two-stage fixpoint combinator
 
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

[edit] List-based tree-merging incremental sieve

Linear merging structure can further be replaced with an indefinitely deepening to the right tree-like structure, (a-(b+((c+d)+( ((e+f)+(g+h)) + ... )))).

This merges primes' multiples streams in a tree-like fashion, achieving theoretical time complexity which is only a log n factor above the optimal n log n log (log n), for n primes produced. Indeed, empirically it runs at about ~ n1.2 (for producing first few million primes), similarly to priority-queue–based version of M. O'Neill's, and with very low space complexity too (not counting the produced sequence of course):

primes :: [Int]   
primes = 2 : _Y ((3 :) . gaps 5 . _U . map(\p-> [p*p, p*p+2*p..]))
 
gaps k s@(x:xs) | k < x = k : gaps (k+2) s -- ~= ([k,k+2..] \\ s)
| otherwise = gaps (k+2) xs -- when s⊂[k,k+2..]
 
_U ((x:xs):t) = x : (union xs . _U . pairs) t -- ~= nub . sort . concat
where
pairs (xs:ys:t) = union xs ys : pairs t

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

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

[edit] Icon and Unicon

 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

Alternatively using sets

 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

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

This problem is a classic example of how J can be used to represent mathematical concepts.

J uses x|y (residue) to represent the operation of finding the remainder during integer division of y divided by x

   10|13
3

And x|/y (table) gives us a table with all possibilities from x and all possibilities from y.

   2 3 4 |/ 2 3 4
0 1 0
2 0 1
2 3 0

Meanwhile, |/~y (reflex) copies the right argument and uses it as the left argment.

   |/~ 0 1 2 3 4
0 1 2 3 4
0 0 0 0 0
0 1 0 1 0
0 1 2 0 1
0 1 2 3 0

(Bigger examples might make the patterns more obvious but they also take up more space.)

By the way, we can ask J to count out the first N integers for us using i. (integers):

   i. 5
0 1 2 3 4

Anyways, the 0s in that last table represent the Sieve of Eratosthenes (in a symbolic or mathematical sense), and we can use = (equal) to find them.

   0=|/~ i.5
1 0 0 0 0
1 1 1 1 1
1 0 1 0 1
1 0 0 1 0
1 0 0 0 1

Now all we need to do is add them up, using / (insert) in its single argument role to insert + between each row of that last table.

   +/0=|/~ i.5
5 1 2 2 3

The sieve wants the cases where we have two divisors:

   2=+/0=|/~ i.5
0 0 1 1 0

And we just want to know the positions of the 1s in that list, which we can find using I. (indices):

   I.2=+/0=|/~ i.5
2 3
I.2=+/0=|/~ 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

And we might want to express this sentence as a definition of a word that lets us use it with an arbitrary argument. There are a variety of ways of doing this. For example:

sieve0=: verb def 'I.2=+/0=|/~ i.y'

Of course, we can also express this in an even more elaborate fashion. The elaboration makes more efficient use of resources for large arguments, at the expense of less efficient use of resources for small arguments:

sieve1=: 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:

sieve2=: 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 but note that sieve2 is still 25 times slower than i.&.(p:inv) for n=1e6. Then again, the value of the sieve of eratosthenes was not efficiency but simplicity. So perhaps we should ignore resource consumption issues and instead focus on intermediate results for reasonably sized example problems?

   0=|/~ i.8
1 0 0 0 0 0 0 0
1 1 1 1 1 1 1 1
1 0 1 0 1 0 1 0
1 0 0 1 0 0 1 0
1 0 0 0 1 0 0 0
1 0 0 0 0 1 0 0
1 0 0 0 0 0 1 0
1 0 0 0 0 0 0 1

Columns with two "1" values correspond to prime numbers.

[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);
}

Version using a BitSet:

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

[edit] Infinite iterator

An iterator that will generate primes indefinitely (perhaps until it runs out of memory), but very slowly.

Translation of: Python
Works with: Java version 1.5+
import java.util.Iterator;
import java.util.PriorityQueue;
import java.math.BigInteger;
 
// generates all prime numbers
public class InfiniteSieve implements Iterator<BigInteger> {
 
private static class NonPrimeSequence implements Comparable<NonPrimeSequence> {
BigInteger currentMultiple;
BigInteger prime;
 
public NonPrimeSequence(BigInteger p) {
prime = p;
currentMultiple = p.multiply(p); // start at square of prime
}
@Override public int compareTo(NonPrimeSequence other) {
// sorted by value of current multiple
return currentMultiple.compareTo(other.currentMultiple);
}
}
 
private BigInteger i = BigInteger.valueOf(2);
// priority queue of the sequences of non-primes
// the priority queue allows us to get the "next" non-prime quickly
final PriorityQueue<NonPrimeSequence> nonprimes = new PriorityQueue<NonPrimeSequence>();
 
@Override public boolean hasNext() { return true; }
@Override public BigInteger next() {
// skip non-prime numbers
for ( ; !nonprimes.isEmpty() && i.equals(nonprimes.peek().currentMultiple); i = i.add(BigInteger.ONE)) {
// 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
while (nonprimes.peek().currentMultiple.equals(i)) {
NonPrimeSequence x = nonprimes.poll();
x.currentMultiple = x.currentMultiple.add(x.prime);
nonprimes.offer(x);
}
}
// prime
// insert a NonPrimeSequence object into the priority queue
nonprimes.offer(new NonPrimeSequence(i));
BigInteger result = i;
i = i.add(BigInteger.ONE);
return result;
}
 
public static void main(String[] args) {
Iterator<BigInteger> sieve = new InfiniteSieve();
for (int i = 0; i < 25; i++) {
System.out.println(sieve.next());
}
}
}
Output:
2
3
5
7
11
13
17
19

[edit] Infinite iterator with a faster algorithm (sieves odds-only)

The adding of each discovered prime's incremental step information to the mapping should be postponed until the candidate number reaches the primes square, as it is useless before that point. This drastically reduces the space complexity from O(n/log(n)) to O(sqrt(n/log(n))), in n primes produced, and also lowers the run time complexity due to the use of the hash table based HashMap, which is much more efficient for large ranges.

Translation of: Python
Works with: Java version 1.5+
import java.util.Iterator;
import java.util.HashMap;
 
// generates all prime numbers up to about 10 ^ 19 if one can wait 1000's of years or so...
public class SoEInfHashMap implements Iterator<Long> {
 
long candidate = 2;
Iterator<Long> baseprimes = null;
long basep = 3;
long basepsqr = 9;
// HashMap of the sequences of non-primes
// the hash map allows us to get the "next" non-prime reasonably quickly
// but further allows re-insertions to take amortized constant time
final HashMap<Long,Long> nonprimes = new HashMap<>();
 
@Override public boolean hasNext() { return true; }
@Override public Long next() {
// do the initial primes separately to initialize the base primes sequence
if (this.candidate <= 5L) if (this.candidate++ == 2L) return 2L; else {
this.candidate++; if (this.candidate == 5L) return 3L; else {
this.baseprimes = new SoEInfHashMap();
this.baseprimes.next(); this.baseprimes.next(); // throw away 2 and 3
return 5L;
} }
// skip non-prime numbers including squares of next base prime
for ( ; this.candidate >= this.basepsqr || //equals nextbase squared => not prime
nonprimes.containsKey(this.candidate); candidate += 2) {
// insert a square root prime sequence into hash map if to limit
if (candidate >= basepsqr) { // if square of base prime, always equal
long adv = this.basep << 1;
nonprimes.put(this.basep * this.basep + adv, adv);
this.basep = this.baseprimes.next();
this.basepsqr = this.basep * this.basep;
}
// else for each sequence that generates this number,
// have it go to the next number (simply add the advance)
// and re-position it in the hash map at an emply slot
else {
long adv = nonprimes.remove(this.candidate);
long nxt = this.candidate + adv;
while (this.nonprimes.containsKey(nxt)) nxt += adv; //unique keys
this.nonprimes.put(nxt, adv);
}
}
// prime
long tmp = candidate; this.candidate += 2; return tmp;
}
 
public static void main(String[] args) {
int n = 100000000;
long strt = System.currentTimeMillis();
SoEInfHashMap sieve = new SoEInfHashMap();
int count = 0;
while (sieve.next() <= n) count++;
long elpsd = System.currentTimeMillis() - strt;
System.out.println("Found " + count + " primes up to " + n + " in " + elpsd + " milliseconds.");
}
 
}
Output:
Found 5761455 primes up to 100000000 in 4297 milliseconds.

[edit] Infinite iterator with a very fast page segmentation algorithm (sieves odds-only)

Although somewhat faster than the previous infinite iterator version, the above code is still over 10 times slower than an infinite iterator based on array paged segmentation as in the following code, where the time to enumerate/iterate over the found primes (common to all the iterators) is now about half of the total execution time:

Translation of: JavaScript
Works with: Java version 1.5+
import java.util.Iterator;
import java.util.ArrayList;
 
// generates all prime numbers up to about 10 ^ 19 if one can wait 100's of years or so...
// practical range is about 10^14 in a week or so...
public class SoEPagedOdds implements Iterator<Long> {
private final int BFSZ = 1 << 16;
private final int BFBTS = BFSZ * 32;
private final int BFRNG = BFBTS * 2;
private long bi = -1;
private long lowi = 0;
private final ArrayList<Integer> bpa = new ArrayList<>();
private Iterator<Long> bps;
private final int[] buf = new int[BFSZ];
 
@Override public boolean hasNext() { return true; }
@Override public Long next() {
if (this.bi < 1) {
if (this.bi < 0) {
this.bi = 0;
return 2L;
}
//this.bi muxt be 0
long nxt = 3 + (this.lowi << 1) + BFRNG;
if (this.lowi <= 0) { // special culling for first page as no base primes yet:
for (int i = 0, p = 3, sqr = 9; sqr < nxt; i++, p += 2, sqr = p * p)
if ((this.buf[i >>> 5] & (1 << (i & 31))) == 0)
for (int j = (sqr - 3) >> 1; j < BFBTS; j += p)
this.buf[j >>> 5] |= 1 << (j & 31);
}
else { // after the first page:
for (int i = 0; i < this.buf.length; i++)
this.buf[i] = 0; // clear the sieve buffer
if (this.bpa.isEmpty()) { // if this is the first page after the zero one:
this.bps = new SoEPagedOdds(); // initialize separate base primes stream:
this.bps.next(); // advance past the only even prime of two
this.bpa.add(this.bps.next().intValue()); // get the next prime (3 in this case)
}
// get enough base primes for the page range...
for (long p = this.bpa.get(this.bpa.size() - 1), sqr = p * p; sqr < nxt;
p = this.bps.next(), this.bpa.add((int)p), sqr = p * p) ;
for (int i = 0; i < this.bpa.size() - 1; i++) {
long p = this.bpa.get(i);
long s = (p * p - 3) >>> 1;
if (s >= this.lowi) // adjust start index based on page lower limit...
s -= this.lowi;
else {
long r = (this.lowi - s) % p;
s = (r != 0) ? p - r : 0;
}
for (int j = (int)s; j < BFBTS; j += p)
this.buf[j >>> 5] |= 1 << (j & 31);
}
}
}
while ((this.bi < BFBTS) &&
((this.buf[(int)this.bi >>> 5] & (1 << ((int)this.bi & 31))) != 0))
this.bi++; // find next marker still with prime status
if (this.bi < BFBTS) // within buffer: output computed prime
return 3 + ((this.lowi + this.bi++) << 1);
else { // beyond buffer range: advance buffer
this.bi = 0;
this.lowi += BFBTS;
return this.next(); // and recursively loop
}
}
 
public static void main(String[] args) {
long n = 1000000000;
long strt = System.currentTimeMillis();
Iterator<Long> gen = new SoEPagedOdds();
int count = 0;
while (gen.next() <= n) count++;
long elpsd = System.currentTimeMillis() - strt;
System.out.println("Found " + count + " primes up to " + n + " in " + elpsd + " milliseconds.");
}
 
}
Output:
Found 50847534 primes up to 1000000000 in 3201 milliseconds.

[edit] JavaScript

function eratosthenes(limit) {
var primes = [];
if (limit >= 2) {
var sqrtlmt = Math.sqrt(limit) - 2;
var nums = new Array(); // start with an empty Array...
for (var i = 2; i <= limit; i++) // and
nums.push(i); // only initialize the Array once...
for (var i = 0; i <= sqrtlmt; i++) {
var p = nums[i]
if (p)
for (var j = p * p - 2; j < nums.length; j += p)
nums[j] = 0;
}
for (var i = 0; i < nums.length; i++) {
var p = nums[i];
if (p)
primes.push(p);
}
}
return primes;
}
 
var primes = eratosthenes(100);
 
if (typeof print == "undefined")
print = (typeof WScript != "undefined") ? WScript.Echo : alert;
print(primes);

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

Substituting the following code for the function for an odds-only algorithm using bit packing for the array produces code that is many times faster than the above:

function eratosthenes(limit) {
var prms = [];
if (limit >= 2) prms = [2];
if (limit >= 3) {
var sqrtlmt = (Math.sqrt(limit) - 3) >> 1;
var lmt = (limit - 3) >> 1;
var bfsz = (lmt >> 5) + 1
var buf = [];
for (var i = 0; i < bfsz; i++)
buf.push(0);
for (var i = 0; i <= sqrtlmt; i++)
if ((buf[i >> 5] & (1 << (i & 31))) == 0) {
var p = i + i + 3;
for (var j = (p * p - 3) >> 1; j <= lmt; j += p)
buf[j >> 5] |= 1 << (j & 31);
}
for (var i = 0; i <= lmt; i++)
if ((buf[i >> 5] & (1 << (i & 31))) == 0)
prms.push(i + i + 3);
}
return prms;
}

While the above code is quite fast especially using an efficient JavaScript engine such as Google Chrome's V8, it isn't as elegant as it could be using the features of the new EcmaScript6 specification when it comes out about the end of 2014 and when JavaScript engines including those of browsers implement that standard in that we might choose to implement an incremental algorithm iterators or generators similar to as implemented in Python or F# (yield). Meanwhile, we can emulate some of those features by using a simulation of an iterator class (which is easier than using a call-back function) for an "infinite" generator based on an Object dictionary as in the following odds-only code written as a JavaScript class:

var SoEIncClass = (function () {
function SoEIncClass() {
this.n = 0;
}
SoEIncClass.prototype.next = function () {
this.n += 2;
if (this.n < 7) { // initialization of sequence to avoid runaway:
if (this.n < 3) { // only even of two:
this.n = 1; // odds from here...
return 2;
}
if (this.n < 5)
return 3;
this.dict = {}; // n must be 5...
this.bps = new SoEIncClass(); // new source of base primes
this.bps.next(); // advance past the even prime of two...
this.p = this.bps.next(); // first odd prime (3 in this case)
this.q = this.p * this.p; // set guard
return 5;
} else { // past initialization:
var s = this.dict[this.n]; // may or may not be defined...
if (!s) { // not defined:
if (this.n < this.q) // haven't reached the guard:
return this.n; // found a prime
else { // n === q => not prime but at guard, so:
var p2 = this.p << 1; // the span odds-only is twice prime
this.dict[this.n + p2] = p2; // add next composite of prime to dict
this.p = this.bps.next();
this.q = this.p * this.p; // get next base prime guard
return this.next(); // not prime so advance...
}
} else { // is a found composite of previous base prime => not prime
delete this.dict[this.n]; // advance to next composite of this prime:
var nxt = this.n + s;
while (this.dict[nxt]) nxt += s; // find unique empty slot in dict
this.dict[nxt] = s; // to put the next composite for this base prime
return this.next(); // not prime so advance...
}
}
};
return SoEIncClass;
})();

The above code can be used to find the nth prime (which would require estimating the required range limit using the previous fixed range code) by using the following code:

var gen = new SoEIncClass(); 
for (var i = 1; i < 1000000; i++, gen.next());
var prime = gen.next();
 
if (typeof print == "undefined")
print = (typeof WScript != "undefined") ? WScript.Echo : alert;
print(prime);

to produce the following output (about five seconds using Google Chrome's V8 JavaScript engine):

15485863

The above code is considerably slower than the fixed range code due to the multiple method calls and the use of an object as a dictionary, which (using a hash table as its basis for most implementations) will have about a constant O(1) amortized time per operation but has quite a high constant overhead to convert the numeric indices to strings which are then hashed to be used as table keys for the look-up operations as compared to doing this more directly in implementations such as the Python dict with Python's built-in hashing functions for every supported type.

This can be implemented as an "infinite" odds-only generator using page segmentation for a considerable speed-up with the alternate JavaScript class code as follows:

var SoEPgClass = (function () {
function SoEPgClass() {
this.bi = -1; // constructor resets the enumeration to start...
}
SoEPgClass.prototype.next = function () {
if (this.bi < 1) {
if (this.bi < 0) {
this.bi++;
this.lowi = 0; // other initialization done here...
this.bpa = [];
return 2;
} else { // bi must be zero:
var nxt = 3 + (this.lowi << 1) + 262144;
this.buf = new Array();
for (var i = 0; i < 4096; i++) // faster initialization:
this.buf.push(0);
if (this.lowi <= 0) { // special culling for first page as no base primes yet:
for (var i = 0, p = 3, sqr = 9; sqr < nxt; i++, p += 2, sqr = p * p)
if ((this.buf[i >> 5] & (1 << (i & 31))) === 0)
for (var j = (sqr - 3) >> 1; j < 131072; j += p)
this.buf[j >> 5] |= 1 << (j & 31);
} else { // after the first page:
if (!this.bpa.length) { // if this is the first page after the zero one:
this.bps = new SoEPgClass(); // initialize separate base primes stream:
this.bps.next(); // advance past the only even prime of two
this.bpa.push(this.bps.next()); // get the next prime (3 in this case)
}
// get enough base primes for the page range...
for (var p = this.bpa[this.bpa.length - 1], sqr = p * p; sqr < nxt;
p = this.bps.next(), this.bpa.push(p), sqr = p * p) ;
for (var i = 0; i < this.bpa.length; i++) {
var p = this.bpa[i];
var s = (p * p - 3) >> 1;
if (s >= this.lowi) // adjust start index based on page lower limit...
s -= this.lowi;
else {
var r = (this.lowi - s) % p;
s = (r != 0) ? p - r : 0;
}
for (var j = s; j < 131072; j += p)
this.buf[j >> 5] |= 1 << (j & 31);
}
}
}
}
while (this.bi < 131072 && this.buf[this.bi >> 5] & (1 << (this.bi & 31)))
this.bi++; // find next marker still with prime status
if (this.bi < 131072) // within buffer: output computed prime
return 3 + ((this.lowi + this.bi++) << 1);
else { // beyond buffer range: advance buffer
this.bi = 0;
this.lowi += 131072;
return this.next(); // and recursively loop
}
};
return SoEPgClass;
})();

The above code is about fifty times faster (about five seconds to calculate 50 million primes to about a billion on the Google Chrome V8 JavaScript engine) than the above dictionary based code.

The speed for both of these "infinite" solutions will also respond to further wheel factorization techniques, especially for the dictionary based version where any added overhead to deal with the factorization wheel will be negligible compared to the dictionary overhead. The dictionary version would likely speed up about a factor of three or a little more with maximum wheel factorization applied; the page segmented version probably won't gain more than a factor of two and perhaps less due to the overheads of array look-up operations.

[edit] Liberty BASIC

    '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

[edit]

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

function erato(n)
if n < 2 then return {} end
local t = {0} -- clears '1'
local sqrtlmt = math.sqrt(n)
for i = 2, n do t[i] = 1 end
for i = 2, sqrtlmt do if t[i] ~= 0 then for j = i*i, n, i do t[j] = 0 end end end
local primes = {}
for i = 2, n do if t[i] ~= 0 then table.insert(primes, i) end end
return primes
end

The following changes the code to odds-only using the same large array-based algorithm:

function erato2(n)
if n < 2 then return {} end
if n < 3 then return {2} end
local t = {}
local lmt = (n - 3) / 2
local sqrtlmt = (math.sqrt(n) - 3) / 2
for i = 0, lmt do t[i] = 1 end
for i = 0, sqrtlmt do if t[i] ~= 0 then
local p = i + i + 3
for j = (p*p - 3) / 2, lmt, p do t[j] = 0 end end end
local primes = {2}
for i = 0, lmt do if t[i] ~= 0 then table.insert(primes, i + i + 3) end end
return primes
end

The following code implements an odds-only "infinite" generator style using a table as a hash table, including postponing adding base primes to the table:

function newEratoInf()
local _cand = 0; local _lstbp = 3; local _lstsqr = 9
local _composites = {}; local _bps = nil
local _self = {}
function _self.next()
if _cand < 9 then if _cand < 1 then _cand = 1; return 2
elseif _cand >= 7 then
--advance aux source base primes to 3...
_bps = newEratoInf()
_bps.next(); _bps.next() end end
_cand = _cand + 2
if _composites[_cand] == nil then -- may be prime
if _cand >= _lstsqr then -- if not the next base prime
local adv = _lstbp + _lstbp -- if next base prime
_composites[_lstbp * _lstbp + adv] = adv -- add cull seq
_lstbp = _bps.next(); _lstsqr = _lstbp * _lstbp -- adv next base prime
return _self.next()
else return _cand end -- is prime
else
local v = _composites[_cand]
_composites[_cand] = nil
local nv = _cand + v
while _composites[nv] ~= nil do nv = nv + v end
_composites[nv] = v
return _self.next() end
end
return _self
end
 
gen = newEratoInf()
count = 0
while gen.next() <= 10000000 do count = count + 1 end -- sieves to 10 million
print(count)
 

which outputs "664579" in about three seconds. As this code uses much less memory for a given range than the previous ones and retains efficiency better with range, it is likely more appropriate for larger sieve ranges.

[edit] Lucid

This example is incorrect. Not a true Sieve of Eratosthenes but rather a Trial Division Sieve Please fix the code and remove this message.
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

This example is incorrect. Not a true Sieve of Eratosthenes but rather a Trial Division Sieve Please fix the code and remove this message.
sieve( N )
   where
    N = 2 fby N + 1;
    sieve( i ) =
      i fby sieve ( i whenever i mod first i ne 0 ) ;
   end

[edit] 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)')')')
 

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

[edit] Mathematica

Eratosthenes[n_] := Module[{numbers = Range[n]},
Do[If[numbers[[i]] != 0, Do[numbers[[i j]] = 0, {j, 2, n/i}]], {i,
2, Sqrt[n]}];
Select[numbers, # > 1 &]]
 
Eratosthenes[100]

[edit] Slightly Optimized Version

The below has been improved to not require so many operations per composite number cull for about two thirds the execution time:

Eratosthenes[n_] := Module[{numbers = Range[n]},
Do[If[numbers[[i]] != 0, Do[numbers[[j]] = 0, {j,i i,n,i}]],{i,2,Sqrt[n]}];
Select[numbers, # > 1 &]]
 
Eratosthenes[100]

[edit] Sieving Odds-Only Version

The below has been further improved to only sieve odd numbers for a further reduction in execution time by a factor of over two:

Eratosthenes2[n_] := Module[{numbers = Range[3, n, 2], limit = (n - 1)/2}, 
Do[c = numbers[[i]]; If[c != 0,
Do[numbers[[j]] = 0, {j,(c c - 1)/2,limit,c}]], {i,1,(Sqrt[n] - 1)/2}];
Prepend[Select[numbers, # > 1 &], 2]]
 
Eratosthenes2[100]

[edit] MATLAB

This example is incorrect. Not a true Sieve of Eratosthenes but rather a Trial Division Sieve Please fix the code and remove this message.
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
Sample Output:
sieveOfEratosthenes(30)

ans =

    2     3     5     7    11    13    17    19    23    29

[edit] Somewhat optimized true Sieve of Eratosthenes

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

[edit] Maxima

sieve(n):=block(
[a:makelist(true,n),i:1,j],
a[1]:false,
do (
i:i+1,
unless a[i] do i:i+1,
if i*i>n then return(sublist_indices(a,identity)),
for j from i*i step i while j<=n do a[j]:false
)
)$

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

[edit] Regular version

This example is incorrect. Not a true Sieve of Eratosthenes but rather a Trial Division Sieve Please fix the code and remove this message.
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.

[edit] Advanced version

This version uses more "advanced" types.

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

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

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,

[edit] NetRexx

[edit] Version 1 (slow)

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

[edit] Version 2 (significantly, i.e. 10 times faster)

/* NetRexx ************************************************************
* Essential improvements:Use boolean instead of Rexx for sv
* and remove methods isTrue and isFalse
* 24.07.2012 Walter Pachl courtesy Kermit Kiser
**********************************************************************/

 
options replace format comments java crossref savelog symbols binary
 
parse arg loWatermark hiWatermark .
if loWatermark = '' | loWatermark = '.' then loWatermark = 1
if hiWatermark = '' | hiWatermark = '.' then hiWatermark = 200000
 
startdate=Date Date()
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')
sieve = sieveOfEratosthenes(hiWatermark)
primes = getPrimes(sieve, loWatermark, hiWatermark).strip
if hiWatermark = 200 Then do
say 'List of prime numbers from' loWatermark 'to' hiWatermark
say ' 'primes.changestr(' ', ',')
end
catch ex = Exception
ex.printStackTrace
end
enddate=Date Date()
Numeric Digits 20
say (enddate.getTime-startdate.getTime)/1000 'seconds elapsed'
say ' Count of primes:' primes.words
 
return
 
method sieveOfEratosthenes(hn = int) -
public static binary returns boolean[]
true = boolean 1
false = boolean 0
sv = boolean[hn+1]
sv[1] = false
 
ix = int
jx = int
 
loop ix=2 to hn
sv[ix]=true
end ix
 
loop ix = 2 while ix * ix <= hn
if sv[ix] then loop jx = ix * ix by ix while jx <= hn
sv[jx] = false
end jx
end ix
 
return sv
 
method getPrimes(sieve = boolean[], lo = int, hi = int) -
private constant binary Returns Rexx
p_ = int
primes = Rexx('')
loop p_ = lo to hi
if \sieve[p_] then iterate p_
primes = primes p_
end p_
 
return primes

[edit] Nial

This example is incorrect. It uses rem testing and so is a trial division algorithm, not a sieve of Eratosthenes. Please fix the code and remove this message.
primes is sublist [ each (2 = sum eachright (0 = mod) [pass,count]), pass ] rest count

Using it

|primes 10
=2 3 5 7

[edit] Nimrod

Based on one of Python solutions:

import math, strutils
 
var is_prime: seq[Bool] = @[]
is_prime.add(False)
is_prime.add(False)
 
iterator iprimes_upto(limit: int): int =
for n in high(is_prime) .. limit+2: is_prime.add(True)
for n in 2 .. limit + 1:
if is_prime[n]:
yield n
for i in countup((n *% n), limit+1, n): # start at ``n`` squared
try:
is_prime[i] = False
except EInvalidIndex: break
 
 
echo("Primes are:")
for x in iprimes_upto(200):
write(stdout, x, " ")
writeln(stdout,"")
Output:
Primes 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

[edit] Niue

This example is incorrect. It uses rem testing and so is a trial division algorithm, not a sieve of Eratosthenes. Please fix the code and remove this message.
[ 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 )

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

[edit] Another functional version

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

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

[edit] Oz

Translation of: Haskell
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}}

[edit] PARI/GP

Eratosthenes(lim)={
my(v=Vectorsmall(lim\1,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", ")))
};

An alternate version:

Sieve(n)=
{
v=vector(n,unused,1);
for(i=2,sqrt(n),
if(v[i],
forstep(j=i^2,n,i,v[j]=0)));
for(i=2,n,if(v[i],print1(i)))
};

[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

For highest performance and ease, typically a module would be used, such as Math::Prime::Util, Math::Prime::FastSieve, or Math::Prime::XS.

[edit] Classic Sieve

sub sieve {
my $n = shift;
my @composite;
for my $i (2 .. int(sqrt($n))) {
if (!$composite[$i]) {
for (my $j = $i*$i; $j <= $n; $j += $i) {
$composite[$j] = 1;
}
}
}
my @primes;
for my $i (2 .. $n) {
$composite[$i] || push @primes, $i;
}
@primes;
}

[edit] Odds only (faster)

sub sieve2 {
my($n) = @_;
return @{([],[],[2],[2,3],[2,3])[$n]} if $n <= 4;
 
my @composite;
for (my $t = 3; $t*$t <= $n; $t += 2) {
if (!$composite[$t]) {
for (my $s = $t*$t; $s <= $n; $s += $t*2)
{ $composite[$s]++ }
}
}
my @primes = (2);
for (my $t = 3; $t <= $n; $t += 2) {
$composite[$t] || push @primes, $t;
}
@primes;
}

[edit] Odds only, using vectors for lower memory use

sub dj_vector {
my($end) = @_;
return @{([],[],[2],[2,3],[2,3])[$end]} if $end <= 4;
$end-- if ($end & 1) == 0; # Ensure end is odd
 
my ($sieve, $n, $limit, $s_end) = ( '', 3, int(sqrt($end)), $end >> 1 );
while ( $n <= $limit ) {
for (my $s = ($n*$n) >> 1; $s <= $s_end; $s += $n) {
vec($sieve, $s, 1) = 1;
}
do { $n += 2 } while vec($sieve, $n >> 1, 1) != 0;
}
my @primes = (2);
do { push @primes, 2*$_+1 if !vec($sieve,$_,1) } for (1..int(($end-1)/2));
@primes;
}

[edit] Odds only, using strings for best performance

Compared to array versions, about 2x faster (with 5.16.0 or later) and lower memory. Much faster than the experimental versions below. It's possible a mod-6 or mod-30 wheel could give more improvement, though possibly with obfuscation. The best next step for performance and functionality would be segmenting.

sub dj_string {
my($end) = @_;
return @{([],[],[2],[2,3],[2,3])[$end]} if $end <= 4;
$end-- if ($end & 1) == 0;
my $s_end = $end >> 1;
 
my $whole = int( ($end>>1) / 15); # prefill with 3 and 5 marked
my $sieve = '100010010010110' . '011010010010110' x $whole;
substr($sieve, ($end>>1)+1) = '';
my ($n, $limit, $s) = ( 7, int(sqrt($end)), 0 );
while ( $n <= $limit ) {
for ($s = ($n*$n) >> 1; $s <= $s_end; $s += $n) {
substr($sieve, $s, 1) = '1';
}
do { $n += 2 } while substr($sieve, $n>>1, 1);
}
# If you just want the count, it's very fast:
# my $count = 1 + $sieve =~ tr/0//;
my @primes = (2, 3, 5);
($s, $n) = (3, 7-2);
while ( (my $nexts = 1 + index($sieve, "0", $s)) > 0 ) {
$n += 2 * ($nexts - $s);
$s = $nexts;
push @primes, $n;
}
@primes;
}

[edit] Experimental

These are examples of golfing or unusual styles.

Golfing a bit, at the expense of speed:

sub sieve{ my (@s, $i);
grep { not $s[ $i = $_ ] and do
{ $s[ $i += $_ ]++ while $i <= $_[0]; 1 }
} 2..$_[0]
}
 
print join ", " => sieve 100;

Or with bit strings (much slower than the vector version above):

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;

A short recursive version:

sub erat {
my $p = shift;
return $p, $p**2 > $_[$#_] ? @_ : erat(grep $_%$p, @_)
}
 
print join ', ' => erat 2..100000;
Regexp (purely an example -- the regex engine limits it to only 32769):
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);

[edit] Extensible sieves

Here are two incremental versions, which allows one to create a tied array of primes:

use strict;
use warnings;
package Tie::SieveOfEratosthenes;
 
sub TIEARRAY {
my $class = shift;
bless \$class, $class;
}
 
# If set to true, produces copious output. Observing this output
# is an excellent way to gain insight into how the algorithm works.
use constant DEBUG => 0;
 
# If set to true, causes the code to skip over even numbers,
# improving runtime. It does not alter the output content, only the speed.
use constant WHEEL2 => 0;
 
BEGIN {
 
# This is loosely based on the Python implementation of this task,
# specifically the "Infinite generator with a faster algorithm"
 
my @primes = (2, 3);
my $ps = WHEEL2 ? 1 : 0;
my $p = $primes[$ps];
my $q = $p*$p;
my $incr = WHEEL2 ? 2 : 1;
my $candidate = $primes[-1] + $incr;
my %sieve;
 
print "Initial: p = $p, q = $q, candidate = $candidate\n" if DEBUG;
 
sub FETCH {
my $n = pop;
return if $n < 0;
return $primes[$n] if $n <= $#primes;
OUTER: while( 1 ) {
 
# each key in %sieve is a composite number between
# p and p-squared. Each value in %sieve is $incr x the prime
# which acted as a 'seed' for that key. We use the value
# to step through multiples of the seed-prime, until we find
# an empty slot in %sieve.
while( my $s = delete $sieve{$candidate} ) {
print "$candidate a multiple of ".($s/$incr).";\t\t" if DEBUG;
my $composite = $candidate + $s;
$composite += $s while exists $sieve{$composite};
print "The next stored multiple of ".($s/$incr)." is $composite\n" if DEBUG;
$sieve{$composite} = $s;
$candidate += $incr;
}
 
print "Candidate $candidate is not in sieve\n" if DEBUG;
 
while( $candidate < $q ) {
print "$candidate is prime\n" if DEBUG;
push @primes, $candidate;
$candidate += $incr;
next OUTER if exists $sieve{$candidate};
}
 
die "Candidate = $candidate, p = $p, q = $q" if $candidate > $q;
print "Candidate $candidate is equal to $p squared;\t" if DEBUG;
 
# Thus, it is now time to add p to the sieve,
my $step = $incr * $p;
my $composite = $q + $step;
$composite += $step while exists $sieve{$composite};
print "The next multiple of $p is $composite\n" if DEBUG;
$sieve{$composite} = $step;
 
# and fetch out a new value for p from our primes array.
$p = $primes[++$ps];
$q = $p * $p;
 
# And since $candidate was equal to some prime squared,
# it's obviously composite, and we need to increment it.
$candidate += $incr;
print "p is $p, q is $q, candidate is $candidate\n" if DEBUG;
} continue {
return $primes[$n] if $n <= $#primes;
}
}
 
}
 
if( !caller ) {
tie my (@prime_list), 'Tie::SieveOfEratosthenes';
my $limit = $ARGV[0] || 100;
my $line = "";
for( my $count = 0; $prime_list[$count] < $limit; ++$count ) {
$line .= $prime_list[$count]. ", ";
next if length($line) <= 70;
if( $line =~ tr/,// > 1 ) {
$line =~ s/^(.*,) (.*, )/$2/;
print $1, "\n";
} else {
print $line, "\n";
$line = "";
}
}
$line =~ s/, \z//;
print $line, "\n" if $line;
}
 
1;

This one is based on the vector sieve shown earlier, but adds to a list as needed, just sieving in the segment. Slightly faster and half the memory vs. the previous incremental sieve. It uses the same API -- arguably we should be offset by one so $primes[$n] returns the $n'th prime.

use strict;
use warnings;
package Tie::SieveOfEratosthenes;
 
sub TIEARRAY {
my $class = shift;
my @primes = (2,3,5,7);
return bless \@primes, $class;
}
 
sub prextend { # Extend the given list of primes using a segment sieve
my($primes, $to) = @_;
$to-- unless $to & 1; # Ensure end is odd
return if $to < $primes->[-1];
my $sqrtn = int(sqrt($to)+0.001);
prextend($primes, $sqrtn) if $primes->[-1] < $sqrtn;
my($segment, $startp) = ('', $primes->[-1]+1);
my($s_beg, $s_len) = ($startp >> 1, ($to>>1) - ($startp>>1));
for my $p (@$primes) {
last if $p > $sqrtn;
if ($p >= 3) {
my $p2 = $p*$p;
if ($p2 < $startp) { # Bump up to next odd multiple of p >= startp
my $f = 1+int(($startp-1)/$p);
$p2 = $p * ($f | 1);
}
for (my $s = ($p2>>1)-$s_beg; $s <= $s_len; $s += $p) {
vec($segment, $s, 1) = 1; # Mark composites in segment
}
}
}
# Now add all the primes found in the segment to the list
do { push @$primes, 1+2*($_+$s_beg) if !vec($segment,$_,1) } for 0 .. $s_len;
}
 
sub FETCHSIZE { 0x7FFF_FFFF } # Allows foreach to work
sub FETCH {
my($primes, $n) = @_;
return if $n < 0;
# Keep expanding the list as necessary, 5% larger each time.
prextend($primes, 1000+int(1.05*$primes->[-1])) while $n > $#$primes;
return $primes->[$n];
}
 
if( !caller ) {
tie my @prime_list, 'Tie::SieveOfEratosthenes';
my $limit = $ARGV[0] || 100;
print $prime_list[0];
my $i = 1;
while ($prime_list[$i] < $limit) { print " ", $prime_list[$i++]; }
print "\n";
}
 
1;

[edit] Perl 6

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;

A recursive version:

multi erat(Int $N) { erat 2 .. $N }
multi erat(@a where @a[0] > sqrt @a[*-1]) { @a }
multi erat(@a) { @a[0], erat(@a.grep: * % @a[0]) }
 
say erat 100;

Of course, upper limits are for wusses. Here's a version using infinite streams, that just keeps going until you ^C it (works under Niecza):

role Filter[Int $factor] {
method next { repeat until $.value % $factor { callsame } }
}
 
class Stream {
has Int $.value is rw = 1;
 
method next { ++$.value }
method filter { self but Filter[$.value] }
}
 
.next, say .value for Stream.new, *.filter ... *;

[edit] 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;
}
 

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

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)

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

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

[edit] Basic procedure

It outputs immediately so that the number can be used by the pipeline.

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

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

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");
}
}

[edit] Prolog

[edit] Using facts to record composite numbers

The first two solutions use Prolog "facts" to record the composite (i.e. already-visited) numbers.

[edit] Elementary approach: multiplication-free, division-free, mod-free, and cut-free

The basic Eratosthenes sieve depends on nothing more complex than counting. In celebration of this simplicity, the first approach to the problem taken here is free of multiplication and division, as well as Prolog's non-logical "cut".

It uses the predicate between/4 to avoid division, and composite/1 to record integers that are found to be composite.

% %sieve( +N, -Primes ) is true if Primes is the list of consecutive primes
% that are less than or equal to N
sieve( N, [2|Rest]) :-
retractall( composite(_) ),
sieve( N, 2, Rest ) -> true. % only one solution
 
% sieve P, find the next non-prime, and then recurse:
sieve( N, P, [I|Rest] ) :-
sieve_once(P, N),
(P = 2 -> P2 is P+1; P2 is P+2),
between(P2, N, I),
(composite(I) -> fail; sieve( N, I, Rest )).
 
% It is OK if there are no more primes less than or equal to N:
sieve( N, P, [] ).
 
sieve_once(P, N) :-
forall( between(P, N, P, IP),
(composite(IP) -> true ; assertz( composite(IP) )) ).
 
 
% To avoid division, we use the iterator
% between(+Min, +Max, +By, -I)
% where we assume that By > 0
% This is like "for(I=Min; I <= Max; I+=By)" in C.
between(Min, Max, By, I) :-
Min =< Max,
A is Min + By,
(I = Min; between(A, Max, By, I) ).
 
 
% Some Prolog implementations require the dynamic predicates be
% declared:
 
:- dynamic( composite/1 ).
 

The above has been tested with SWI-Prolog and gprolog.

% SWI-Prolog:
 
?- time( (sieve(100000,P), length(P,N), writeln(N), last(P, LP), writeln(LP) )).
% 1,323,159 inferences, 0.862 CPU in 0.921 seconds (94% CPU, 1534724 Lips)
P = [2, 3, 5, 7, 11, 13, 17, 19, 23|...],
N = 9592,
LP = 99991.
 

[edit] Optimized approach

Works with SWI-Prolog.

sieve(N, [2|PS]) :-       % PS is list of odd primes up to N
retractall(mult(_)),
sieve_O(3,N,PS).
 
sieve_O(I,N,PS) :- % sieve odds from I up to N to get PS
I =< N, !, I1 is I+2,
( mult(I) -> sieve_O(I1,N,PS)
; ( I =< N / I ->
ISq is I*I, DI is 2*I, add_mults(DI,ISq,N)
; true
),
PS = [I|T],
sieve_O(I1,N,T)
).
sieve_O(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.

[edit] Using a priority queue

Uses a ariority queue, from the paper "The Genuine Sieve of Eratosthenes" by Melissa O'Neill. Works with YAP (Yet Another Prolog)

?- use_module(library(heaps)).
 
prime(2).
prime(N) :- prime_heap(N, _).
 
prime_heap(3, H) :- list_to_heap([9-6], H).
prime_heap(N, H) :-
prime_heap(M, H0), N0 is M + 2,
next_prime(N0, H0, N, H).
 
next_prime(N0, H0, N, H) :-
\+ min_of_heap(H0, N0, _),
N = N0, Composite is N*N, Skip is N+N,
add_to_heap(H0, Composite, Skip, H).
next_prime(N0, H0, N, H) :-
min_of_heap(H0, N0, _),
adjust_heap(H0, N0, H1), N1 is N0 + 2,
next_prime(N1, H1, N, H).
 
adjust_heap(H0, N, H) :-
min_of_heap(H0, N, _),
get_from_heap(H0, N, Skip, H1),
Composite is N + Skip, add_to_heap(H1, Composite, Skip, H2),
adjust_heap(H2, N, H).
adjust_heap(H, N, H) :-
\+ min_of_heap(H, N, _).

[edit] PureBasic

[edit] Basic procedure

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

[edit] Working example

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

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

[edit] Python

Note that the examples use range instead of xrange for Python 3 and Python 2 compatability, but when using Python 2 xrange is the nearest equivalent to Python 3's implementation of range and should be substituted for ranges with a very large number of items.

[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()
for i in range(2, n+1):
if i not in multiples:
print(i)
multiples.update(range(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 range(int(limit**0.5 + 1.5)): # stop at ``sqrt(limit)``
if is_prime[n]:
for i in range(n*n, limit+1, n):
is_prime[i] = False
return [i for i, prime in enumerate(is_prime) if prime]

[edit] Using generator

The following code may be slightly slower than using the array/list as above, but uses no memory for output:

def iprimes_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 range(n * n, limit + 1, n): # start at ``n`` squared
is_prime[i] = False
for i in xrange(limit + 1):
if is_prime[i]: yield i
Example:
>>> list(iprimes_upto(15))
[2, 3, 5, 7, 11, 13]

[edit] Odds-only version of the array sieve above

The following code is faster than the above array version using only odd composite operations (for a factor of over two) and because it has been optimized to use slice operations for composite number culling to avoid extra work by the interpreter:

def primes2(limit):
if limit < 2: return []
if limit < 3: return [2]
lmtbf = (limit - 3) // 2
buf = [True] * (lmtbf + 1)
for i in range((int(limit ** 0.5) - 3) // 2 + 1):
if buf[i]:
p = i + i + 3
s = p * (i + 1) + i
buf[s::p] = [False] * ((lmtbf - s) // p + 1)
return [2] + [i + i + 3 for i, v in enumerate(buf) if v]

Note that "range" needs to be changed to "xrange" for maximum speed with Python 2.

[edit] Odds-only version of the generator version above

The following code is faster than the above generator version using only odd composite operations (for a factor of over two) and because it has been optimized to use slice operations for composite number culling to avoid extra work by the interpreter:

def iprimes2(limit):
yield 2
if limit < 3: return
lmtbf = (limit - 3) // 2
buf = [True] * (lmtbf + 1)
for i in range((int(limit ** 0.5) - 3) // 2 + 1):
if buf[i]:
p = i + i + 3
s = p * (i + 1) + i
buf[s::p] = [False] * ((lmtbf - s) // p + 1)
for i in range(lmtbf + 1):
if buf[i]: yield (i + i + 3)

Note that this version may actually run slightly faster than the equivalent array version with the advantage that the output doesn't require any memory.

Also note that "range" needs to be changed to "xrange" for maximum speed with Python 2.

[edit] Factorization wheel235 version of the generator version

This uses a 235 factorial wheel for further reductions in operations; the same techniques can be applied to the array version as well; it runs slightly faster and uses slightly less memory as compared to the odds-only algorithms:

def primes235(limit):
yield 2; yield 3; yield 5
if limit < 7: return
modPrms = [7,11,13,17,19,23,29,31]
gaps = [4,2,4,2,4,6,2,6,4,2,4,2,4,6,2,6] # 2 loops for overflow
ndxs = [0,0,0,0,1,1,2,2,2,2,3,3,4,4,4,4,5,5,5,5,5,5,6,6,7,7,7,7,7,7]
lmtbf = (limit + 23) // 30 * 8 - 1 # integral number of wheels rounded up
lmtsqrt = (int(limit ** 0.5) - 7)
lmtsqrt = lmtsqrt // 30 * 8 + ndxs[lmtsqrt % 30] # round down on the wheel
buf = [True] * (lmtbf + 1)
for i in range(lmtsqrt + 1):
if buf[i]:
ci = i & 7; p = 30 * (i >> 3) + modPrms[ci]
s = p * p - 7; p8 = p << 3
for j in range(8):
c = s // 30 * 8 + ndxs[s % 30]
buf[c::p8] = [False] * ((lmtbf - c) // p8 + 1)
s += p * gaps[ci]; ci += 1
for i in range(lmtbf - 6 + (ndxs[(limit - 7) % 30])): # adjust for extras
if buf[i]: yield (30 * (i >> 3) + modPrms[i & 7])

Note: Much of the time (almost two thirds for this last case for Python 2.7.6) for any of these array/list or generator algorithms is used in the computation and enumeration of the final output in the last line(s), so any slight changes to those lines can greatly affect execution time. For Python 3 this enumeration is about twice as slow as Python 2 (Python 3.3 slow and 3.4 slower) for an even bigger percentage of time spent just outputting the results. This slow enumeration means that there is little advantage to versions that use even further wheel factorization, as the composite number culling is a small part of the time to enumerate the results.

If just the count of the number of primes over a range is desired, then converting the functions to prime counting functions by changing the final enumeration lines to "return buf.count(True)" will save a lot of time.

Note that "range" needs to be changed to "xrange" for maximum speed with Python 2 where Python 2's "xrange" is a better choice for very large sieve ranges.
Timings were done primarily in Python 2 although source is Python 2/3 compatible (shows range and not xrange).

[edit] Using numpy

Library: numpy

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

Also see Prime numbers and Numpy – Python.

[edit] Using wheels with numpy

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 (perhaps until it runs out of memory). Used as a library here.

Works with: Python version 2.6+, 3.x
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

Example:

>>> foo = sieve()
>>> for i in range(8):
...     print(next(foo))
... 
2
3
5
7
11
13
17
19

[edit] Infinite generator with a faster algorithm

The adding of each discovered prime's incremental step info to the mapping should be postponed until the prime's square is seen amongst the candidate numbers, as it is useless before that point. This drastically reduces the space complexity from O(n) to O(sqrt(n/log(n))), in n primes produced, and also lowers the run time complexity quite low (this test entry in Python 2.7 and this test entry in Python 3.x shows about O(n^1.08) empirical order of growth which is very close to the theoretical value of O(n log(n) log(log(n))), in n primes produced):

Works with: Python version 2.6+, 3.x
def primes():
yield 2; yield 3; yield 5; yield 7;
bps = (p for p in primes()) # additional primes supply
p = next(bps) and next(bps) # discard 2, then get 3
q = p * p # 9 - square of next prime to be put
sieve = {} # into sieve dict
n = 9 # the initial candidate number
while True:
if n not in sieve: # is not a multiple of previously recorded primes
if n < q:
yield n # n is prime
else:
p2 = p + p # n == p * p: for prime p, add p * p + 2 * p,
sieve[q + p2] = p2 # with 2 * p as incremental step
p = next(bps); q = p * p # advance base prime and next prime to be put
else:
s = sieve.pop(n); nxt = n + s # n is composite, advance
while nxt in sieve: nxt += s # ensure each entry is unique
sieve[nxt] = s # next non-marked multiple of this prime
n += 2 # work on odds only
 
import itertools
def primes_up_to(limit):
return list(itertools.takewhile(lambda p: p <= limit, primes()))

[edit] Fast infinite generator using a wheel

Although theoretically over three times faster than odds-only, the following code using a 2/3/5/7 wheel is only about 1.5 times faster than the above odds-only code due to the extra overheads in code complexity. The test link for Python 2.7 and test link for Python 3.x show about the same empirical order of growth as the odds-only implementation above once the range grows enough so the dict operations become amortized to a constant factor.

Works with: Python version 2.6+, 3.x
def primes():
for p in [2,3,5,7]: yield p # base wheel primes
gaps1 = [ 2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,8 ]
gaps = gaps1 + [ 6,4,6,2,4,6,2,6,6,4,2,4,6,2,6,4,2,4,2,10,2,10 ] # wheel2357
def wheel_prime_pairs():
yield (11,0); bps = wheel_prime_pairs() # additional primes supply
p, pi = next(bps); q = p * p # adv to get 11 sqr'd is 121 as next square to put
sieve = {}; n = 13; ni = 1 # into sieve dict; init cndidate, wheel ndx
while True:
if n not in sieve: # is not a multiple of previously recorded primes
if n < q: yield (n, ni) # n is prime with wheel modulo index
else:
npi = pi + 1 # advance wheel index
if npi > 47: npi = 0
sieve[q + p * gaps[pi]] = (p, npi) # n == p * p: put next cull position on wheel
p, pi = next(bps); q = p * p # advance next prime and prime square to put
else:
s, si = sieve.pop(n)
nxt = n + s * gaps[si] # move current cull position up the wheel
si = si + 1 # advance wheel index
if si > 47: si = 0
while nxt in sieve: # ensure each entry is unique by wheel
nxt += s * gaps[si]
si = si + 1 # advance wheel index
if si > 47: si = 0
sieve[nxt] = (s, si) # next non-marked multiple of a prime
nni = ni + 1 # advance wheel index
if nni > 47: nni = 0
n += gaps[ni]; ni = nni # advance on the wheel
for p, pi in wheel_prime_pairs(): yield p # strip out indexes

Further gains of about 1.5 times in speed can be made using the same code by only changing the tables and a few constants for a further constant factor gain of about 1.5 times in speed by using a 2/3/5/7/11/13/17 wheel (with the gaps list 92160 elements long) computed for a slight constant overhead time as per the test link for Python 2.7 and test link for Python 3.x. Further wheel factorization will not really be worth it as the gains will be small (if any and not losses) and the gaps table huge - it is already too big for efficient use by 32-bit Python 3 and the wheel should likely be stopped at 13:

def primes():
whlPrms = [2,3,5,7,11,13,17] # base wheel primes
for p in whlPrms: yield p
def makeGaps():
buf = [True] * (3 * 5 * 7 * 11 * 13 * 17 + 1) # all odds plus extra for o/f
for p in whlPrms:
if p < 3:
continue # no need to handle evens
strt = (p * p - 19) >> 1 # start position (divided by 2 using shift)
while strt < 0: strt += p
buf[strt::p] = [False] * ((len(buf) - strt - 1) // p + 1) # cull for p
whlPsns = [i + i for i,v in enumerate(buf) if v]
return [whlPsns[i + 1] - whlPsns[i] for i in range(len(whlPsns) - 1)]
gaps = makeGaps() # big wheel gaps
def wheel_prime_pairs():
yield (19,0); bps = wheel_prime_pairs() # additional primes supply
p, pi = next(bps); q = p * p # adv to get 11 sqr'd is 121 as next square to put
sieve = {}; n = 23; ni = 1 # into sieve dict; init cndidate, wheel ndx
while True:
if n not in sieve: # is not a multiple of previously recorded primes
if n < q: yield (n, ni) # n is prime with wheel modulo index
else:
npi = pi + 1 # advance wheel index
if npi > 92159: npi = 0
sieve[q + p * gaps[pi]] = (p, npi) # n == p * p: put next cull position on wheel
p, pi = next(bps); q = p * p # advance next prime and prime square to put
else:
s, si = sieve.pop(n)
nxt = n + s * gaps[si] # move current cull position up the wheel
si = si + 1 # advance wheel index
if si > 92159: si = 0
while nxt in sieve: # ensure each entry is unique by wheel
nxt += s * gaps[si]
si = si + 1 # advance wheel index
if si > 92159: si = 0
sieve[nxt] = (s, si) # next non-marked multiple of a prime
nni = ni + 1 # advance wheel index
if nni > 92159: nni = 0
n += gaps[ni]; ni = nni # advance on the wheel
for p, pi in wheel_prime_pairs(): yield p # strip out indexes
 

[edit] 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.
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)

[edit] Racket

[edit] Imperative versions

Ugly imperative version:

#lang racket
 
(define (sieve n)
(define non-primes '())
(define primes '())
(for ([i (in-range 2 (add1 n))])
(unless (member i non-primes)
(set! primes (cons i primes))
(for ([j (in-range (* i i) (add1 n) i)])
(set! non-primes (cons j non-primes)))))
(reverse primes))
 
(sieve 100)

A little nicer, but still imperative:

#lang racket
 
(define (sieve n)
(define primes (make-vector (add1 n) #t))
(for* ([i (in-range 2 (add1 n))]
#:when (vector-ref primes i)
[j (in-range (* i i) (add1 n) i)])
(vector-set! primes j #f))
(for/list ([n (in-range 2 (add1 n))]
#:when (vector-ref primes n))
n))
 
(sieve 100)

[edit] Infinite list of primes

[edit] Using laziness

These examples use infinite lists (streams) to implement the sieve of Eratosthenes in a functional way, and producing all prime numbers. The following functions are used as a prefix for pieces of code that follow:

#lang lazy
(define (ints-from i d) (cons i (ints-from (+ i d) d)))
(define (when-bigger n l f)
(if (< (car l) n) (cons (car l) (when-bigger n (cdr l) f)) (f l)))
(define (diff l1 l2)
(let ([x1 (car l1)] [x2 (car l2)])
(cond [(< x1 x2) (cons x1 (diff (cdr l1) l2 ))]
[(> x1 x2) (diff l1 (cdr l2)) ]
[else (diff (cdr l1) (cdr l2)) ])))
(define (merge l1 l2)
(let ([x1 (car l1)] [x2 (car l2)])
(cond [(< x1 x2) (cons x1 (merge (cdr l1) l2 ))]
[(> x1 x2) (cons x2 (merge l1 (cdr l2)))]
[else (cons x1 (merge (cdr l1) (cdr l2)))])))
[edit] Basic sieve
(define (sieve l)
(define x (car l))
(cons x (sieve (diff (cdr l) (ints-from (+ x x) x)))))
(define primes (sieve (ints-from 2 1)))
(!! (take 25 primes))

Runs at ~ n^2.1 empirically, for n <= 1500 primes produced.

[edit] With merged composites

Note that the first number, 2, is done manually outside of the sieve loop -- that's not for optimization, it's done to ensure that the non-primes list is never empty, which simplifies the code since all functions assume non-empty infinite lists.

(define (sieve l non-primes)
(let ([x (car l)] [np (car non-primes)])
(cond [(< np x) (sieve l (cdr non-primes))]
[(= np x) (sieve (cdr l) (cdr non-primes))]
[else (cons x (sieve (cdr l) (merge (ints-from (* x x) x)
non-primes)))])))
(define primes (cons 2 (sieve (ints-from 3 2) (ints-from 2 2))))

Runs at ~ n^2.4, empirically, producing up to n=700 primes.

[edit] Basic sieve Optimized with postponed processing

Since a prime's multiples that count start from its square, we should only start removing them when we reach that square.

(define (sieve l prs)
(define p (car prs))
(define q (* p p))
(when-bigger q l (λ(t) (sieve (diff t (ints-from q p)) (cdr prs)))))
(define primes (cons 2 (sieve (ints-from 3 1) primes)))

Runs at ~ n^1.4 up to n=10,000. The initial 2 in the self-referential primes definition is needed to prevent a "black hole".

[edit] Merged composites Optimized with postponed processing

Since a prime's multiples that count start from its square, we should only add them when we reach that square. Improves the time complexity to ~ n^1.5, tested producing up to n=10,000 primes.

(define (sieve l q prs non-primes)
(let ([x (car l)])
(if (= x q)  ; also p*p == q
(let ([p (car prs)] [p2 (cadr prs)])
(sieve (cdr l) (* p2 p2) (cdr prs)
(if (not non-primes)
(ints-from q (+ p p))
(merge non-primes (ints-from q (+ p p))))))
(let ([np (car non-primes)])  ; here x < q
(cond [(< np x) (sieve l q prs (cdr non-primes)) ]
[(= np x) (sieve (cdr l) q prs (cdr non-primes)) ]
[else (cons x (sieve (cdr l) q prs non-primes ))])))))
(define primes (append (list 2 3 5 7)
(sieve (ints-from 9 2) 9 (cdr primes) #f)))

(list 2 3 5 7) in the primes definition is not part of a bigger wheel optimization, but just for a simpler coding of the sieve function. Though the processing starts from odds only, ignoring all evens above 2, so is effectively employing the simplest kind of a wheel.

[edit] Implementation of Richard Bird's algorithm

Appears in M.O'Neill's paper. Achieves on its own what is done laboriously in the postponed version above, and should be yet more efficient, because it folds to the right and so builds the right-leaning structure of merges, at run time, where the more frequently-producing streams of multiples appear higher in that structure, so have less merge nodes to percolate through.

(define (minus l1 l2)         ; assumes l2 is a subset of l1
(if (= (car l1) (car l2))
(minus (cdr l1) (cdr l2))
(cons (car l1) (minus (cdr l1) l2))))
(define primes
(cons 2 (minus (ints-from 3 1)
(foldr (λ(p r) (define q (* p p))
(cons q (merge (ints-from (+ q p) p) r)))
'() primes))))

[edit] Using threads and channels

Same algorithm as "merged composites" above (without the optimization), but now using threads and channels to produce a channel of all prime numbers (similar to newsqueak). The macro at the top is a convenient wrapper around definitions of channels using a thread that feeds them.

#lang racket
(define-syntax (define-thread-loop stx)
(syntax-case stx ()
[(_ (name . args) expr ...)
(with-syntax ([out! (datum->syntax stx 'out!)])
#'(define (name . args)
(define out (make-channel))
(define (out! x) (channel-put out x))
(thread (λ() (let loop () expr ... (loop))))
out))]))
(define-thread-loop (ints-from i d) (out! i) (set! i (+ i d)))
(define-thread-loop (merge c1 c2)
(let loop ([x1 (channel-get c1)] [x2 (channel-get c2)])
(cond [(< x1 x2) (out! x1) (loop (channel-get c1) x2)]
[(> x1 x2) (out! x2) (loop x1 (channel-get c2))]
[else (out! x1) (loop (channel-get c1) (channel-get c2))])))
(define-thread-loop (sieve l non-primes)
(let loop ([x (channel-get l)] [np (channel-get non-primes)])
(cond [(< np x) (loop x (channel-get non-primes))]
[(= np x) (loop (channel-get l) (channel-get non-primes))]
[else (out! x)
(set! non-primes (merge (ints-from (* x x) x) non-primes))])))
(define-thread-loop (cons x l)
(out! x) (let loop () (out! (channel-get l)) (loop)))
(define primes (cons 2 (sieve (ints-from 3 2) (ints-from 2 2))))
(for/list ([i 25] [x (in-producer channel-get eof primes)]) x)

[edit] Using generators

Yet another variation of the same algorithm as above, this time using generators.

#lang racket
(require racket/generator)
(define (ints-from i d)
(generator () (let loop ([i i]) (yield i) (loop (+ i d)))))
(define (merge g1 g2)
(generator ()
(let loop ([x1 (g1)] [x2 (g2)])
(cond [(< x1 x2) (yield x1) (loop (g1) x2)]
[(> x1 x2) (yield x2) (loop x1 (g2))]
[else (yield x1) (loop (g1) (g2))]))))
(define (sieve l non-primes)
(generator ()
(let loop ([x (l)] [np (non-primes)])
(cond [(< np x) (loop x (non-primes))]
[(= np x) (loop (l) (non-primes))]
[else (yield x)
(set! non-primes (merge (ints-from (* x x) x) non-primes))
(loop (l) (non-primes))]))))
(define (cons x l) (generator () (yield x) (let loop () (yield (l)) (loop))))
(define primes (cons 2 (sieve (ints-from 3 2) (ints-from 2 2))))
(for/list ([i 25] [x (in-producer primes)]) x)

[edit] REXX

[edit] no wheel version

The first three REXX versions make use of a sparse stemmed array:   [@.].
As the stemmed array gets heavily populated, the number of entries may slow down the REXX interpreter
substantially, depending upon the efficacy of the hashing technique being used for REXX variables (setting/retrieving).

/*REXX program generates primes via the sieve of Eratosthenes algorithm.*/
parse arg H .; if H=='' then H=200 /*was the high limit specified? */
w=length(H); @prime=right('prime',20) /*W is used for formatting output*/
@.=. /*assume all numbers are prime. */
#=0 /*number of primes found so far. */
do j=2 for H-1 /*all integers up to H inclusive.*/
if @.j==. then do; #=#+1 /*Prime? Then bump prime counter*/
say @prime right(#,w) " ───► " right(j,w)
do m=j*j to H by j; @.m=; end /*m*/
end /*[↑] strike odd multiples ¬prime*/
end /*j*/ /* ─── */
/*stick a fork in it, we're done.*/
say; say right(#,w+length(@prime)+1) 'primes found.'

output when using the input default of:   200

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

                      46 primes found.

[edit] wheel version, optional prime list suppression

This version skips striking the even numbers   (as being not prime),   2   is handled as a special case.

Also supported is the suppression of listing the primes if the H (high limit) is negative.
Also added is a final message indicating the number of primes found.

/*REXX pgm gens primes via a  wheeled  sieve of Eratosthenes  algorithm.*/
parse arg H .; if H=='' then H=200 /*let the highest # be specified.*/
tell=h>0; H=abs(H); w=length(H) /*neg H suppresses prime listing.*/
if 2<=H & tell then say right(1,w+20)'st prime ───► ' right(2,w)
skip=0 /*skips top part of sieve marking*/
#= w<=H /*number of primes found so far. */
@.=. /*assume all numbers are prime. */
 
do j=3 by 2 for (H-2)%2 /*odd integers up to H inclusive.*/
if @.j==. then do; #=#+1 /*Prime? Then bump prime counter*/
if tell then say right(#,w+20)th(#) 'prime ───► ' right(j,w)
if skip then iterate /*should the top part be skipped?*/
jj=j*j /*compute the square of J. __ */
if jj>H then skip=1 /*indicate skipping if j > √ H. */
do m=jj to H by j+j; @.m=; end /*m*/
end /*[↑] strike odd multiples ¬prime*/
end /*j*/ /* ─── */
 
say; say right(#,w+20) 'prime's(#) "found."
exit /*stick a fork in it, we're done.*/
/*──────────────────────────────────one─liner subroutines────────────────────────────────────────*/
s: if arg(1)==1 then return arg(3); return word(arg(2) 's',1) /*pluralizer.*/
th: procedure; parse arg x; x=abs(x); return word('th st nd rd',1+x//10*(x//100%10\==1)*(x//10<4))

output when using the input default of:   200

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

                       46 primes found.

output when using the input of:   -1000

                       168 primes found.

output when using the input of:   -10000

                       1229 primes found.

output when using the input of:   -100000

                        9592 primes found.

output when using the input of:   -1000000

                        78498 primes found.

output when using the input of:   -10000000

                        664579 primes found.

output when using the input of:   -100000000

    16 +++       @.m=
Error 5 running "C:\sieve_of_Eratosthenes.rex", line 16: System resources exhausted

The above (using Regina 3.82 under Windows/XP) shows one of the weaknesses of this implementation of the Sieve of Eratosthenes:   it must keep an array of all (if not most) values which is used to strike out composite numbers.

The   System resources exhausted   error can be postponed by implementing further optimizations (expanding the wheel with low primes).

[edit] wheel version

This version skips striking the even numbers   (as being not prime).
It also uses a short-circuit test for striking out composites ≤ √target.

/*REXX pgm gens primes via a  wheeled  sieve of Eratosthenes  algorithm.*/
parse arg H .; if H=='' then H=200 /*high# can be specified on C.L. */
w=length(H); @prime=right('prime',20) /*W is used for formatting output*/
if 2<=H then say @prime right(1,w) " ───► " right(2,w)
skip=0 /*skips top part of sieve marking*/
#= 2<=H /*number of primes found so far. */
@.=. /*assume all numbers are prime. */
do j=3 by 2 for (H-2)%2 /*odd integers up to H inclusive.*/
if @.j==. then do; #=#+1 /*Prime? Then bump prime counter*/
say @prime right(#,w) " ───► " right(j,w)
if skip then iterate /*should the top part be skipped?*/
jj=j*j /*compute the square of J. __ */
if jj>H then skip=1 /*indicate skipping if j > √ H.*/
do m=jj to H by j+j; @.m=; end /*m*/
end /*[↑] strike odd multiples ¬prime*/
end /*j*/ /* ─── */
/*stick a fork in it, we're done.*/
say; say right(#,w+length(@prime)+1) 'primes found.'

output is identical to the first (non-wheel) version;   program execution is twice as fast as the non-wheel version.
The addition of the short-circuit test makes it about another 20% faster.

[edit] Wheel Version restructured

/*REXX program generates primes via sieve of Eratosthenes algorithm.
* 21.07.2012 Walter Pachl derived from above Rexx version
* avoid symbols @ and # (not supported by ooRexx)
* avoid negations (think positive)
**********************************************************************/

highest=200 /*define highest number to use. */
is_prime.=1 /*assume all numbers are prime. */
w=length(highest) /*width of the biggest number, */
/* it's used for aligned output.*/
Do j=3 To highest By 2, /*strike multiples of odd ints. */
While j*j<=highest /* up to sqrt(highest) */
If is_prime.j Then Do
Do jm=j*3 To highest By j+j /*start with next odd mult. of J */
is_prime.jm=0 /*mark odd mult. of J not prime. */
End
End
End
np=0 /*number of primes shown */
Call tell 2
Do n=3 To highest By 2 /*list all the primes found. */
If is_prime.n Then Do
Call tell n
End
End
Exit
tell: Parse Arg prime
np=np+1
Say ' prime number' right(np,w) " --> " right(prime,w)
Return

output is identical to the above versions.

[edit] Ruby

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

def eratosthenes(n)
nums = [nil, nil, *2..n]
(2..Math.sqrt(n)).each do |i|
(i**2..n).step(i){|m| nums[m] = nil} if nums[i]
end
nums.compact
end
 
p 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]

[edit] 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.
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) > n
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
primes
end
 
p eratosthenes2(100)

This simple benchmark compares eratosthenes with eratosthenes2.

require 'benchmark'
Benchmark.bmbm {|x|
x.report("eratosthenes") { eratosthenes(1_000_000) }
x.report("eratosthenes2") { eratosthenes2(1_000_000) }
}

eratosthenes2 runs about 4 times faster than eratosthenes.

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

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

[edit] Run BASIC

input "Gimme the limit:"; limit
dim flags(limit)
for i = 2 to limit
for k = i*i to limit step i
flags(k) = 1
next k
if flags(i) = 0 then print i;", ";
next i
Gimme the limit:?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] Rust

Library: Rust

[edit] Sieve of Eratosthenes - No optimization

use std::iter;
 
fn int_sqrt(n: uint) -> uint {
(n as f64).sqrt() as uint
}
 
fn simple_sieve(limit: uint) -> Vec<uint> {
 
if limit < 2 {
return Vec::new();
}
 
let mut primes: Vec<bool> = Vec::from_elem(limit + 1, true);
 
for prime in iter::range_inclusive(2, int_sqrt(limit) + 1) {
if primes[prime] {
for multiple in iter::range_step(prime * prime, limit + 1, prime) {
*primes.get_mut(multiple) = false;
}
}
}
 
iter::range_inclusive(2, limit).filter(|&n| primes[n]).collect()
}
 
fn main() {
println!("{}", simple_sieve(100))
}
 
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]

[edit] Scala

[edit] Genuine Eratosthenes sieve

import scala.annotation.tailrec
import scala.collection.parallel.mutable
import scala.compat.Platform
 
object GenuineEratosthenesSieve extends App {
def sieveOfEratosthenes(limit: Int) = {
val (primes: mutable.ParSet[Int], sqrtLimit) = (mutable.ParSet.empty ++ (2 to limit), math.sqrt(limit).toInt)
@tailrec
def prim(candidate: Int): Unit = {
if (candidate <= sqrtLimit) {
if (primes contains candidate) primes --= candidate * candidate to limit by candidate
prim(candidate + 1)
}
}
prim(2)
primes
}
// BitSet toList is shuffled when using the ParSet version. So it has to be sorted before using it as a sequence.
 
assert(sieveOfEratosthenes(15099480).size == 976729)
println(s"Successfully completed without errors. [total ${Platform.currentTime - executionStart} ms]")
}
Output:
Successfully completed without errors. [total 39807 ms]

Process finished with exit code 0

While concise, the above code is quite slow but a little faster not using the ParSet (take out the '.par' for speed), in which case the sorting ('sorted') is not necessary for an additional small gain in speed; the above code is slow because of all the overhead in processing the bit packed "BitSet" bib-by-bit using complex "higher-order" method calls.

[edit] Using tail recursion

The below code using a primitive array (bit packed) and tail recursion to avoid some of the enumeration delays is almost three times faster:

import scala.compat.Platform
import scala.annotation.tailrec
 
object SoEwithOutBitSet extends App {
val expected = 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)
 
def sieveOfEratosthenes(limit: Int) = {
require(limit <= 2147437383, s"Limit exceeded with calling $limit")
val (composites, sqrtLimit) = (new Array[Int]((limit >>> 5) + 1), Math.sqrt(limit).toInt)
 
def getBit(c: Int): Boolean = (composites(c >>> 5) & (1 << (c & 31))) == 0
@tailrec
def loop(c: Int): Unit = {
def cullp(p: Int) = {
@tailrec
def cull(c: Int): Unit = {
if (c <= limit) {
composites(c >>> 5) |= 1 << (c & 31)
cull(c + p)
}
}
cull(p * p)
}
 
if (c <= sqrtLimit) {
if (getBit(c)) cullp(c)
loop(c + 1)
}
}
loop(2)
 
@tailrec
def nextPrime(i: Int): Int = if ((i > limit) || getBit(i)) i else nextPrime(i + 1)
 
Iterator.iterate(2)(p => nextPrime(p + 1)).takeWhile(_ <= limit)
}
 
{
val soe101 = sieveOfEratosthenes(101).toList
assert(soe101 == expected)
assert(soe101.indexOf(31) == 10)
}
println("Generating more than 100 million primes, moment please..")
// sieveOfEratosthenes(2147437291).dropWhile(_ < 2147437283).foreach(println)
assert(sieveOfEratosthenes(2147437291).size == 105095435)
println(s"Successfully completed without errors. [total ${Platform.currentTime - executionStart} ms]") // ~~ 1 minute
}
Output:
Generating more than 100 million primes, moment please..
Successfully completed without errors. [total 63270 ms]

Process finished with exit code 0

[edit] Odds-only page-segmented "infinite" generator version using tail recursion

The above code still uses an amount of memory proportional to the range of the sieve (although bit packed). As well as only sieving odd candidates, the following code uses a fixed range buffer that is about the size of the CPU L1 cache plus only storage for the base primes up to the square root of the range for a large potential saving in RAM memory used. The use of innermost tail recursive loops rather than "higher order" functions from iterators also greatly reduces execution time, with much of the remaining time used just to enumerate the primes output:

import scala.collection.mutable.ArrayBuffer
import scala.annotation.tailrec
 
def SoEPg(): Iterator[Int] = {
val basePrimesArr: ArrayBuffer[Int] = ArrayBuffer.empty
 
def makePg(lowp: Int, basePrimes: Iterator[Int]) = {
val nxtlowp = lowp + 262144
val low = (lowp - 3) >>> 1
val buf: Array[Int] = new Array(131072)
 
def nextPrime(c: Int) = {
@tailrec def nexti(i: Int): Int = {
if ((i < 131072) && ((buf(i >>> 5) & (1 << (i & 31)))) != 0) nexti(i + 1)
else i
}
(nexti((c - lowp) >>> 1) << 1) + lowp
}
 
def cullp(s: Int, p: Int) = {
@tailrec def cull(c: Int): Unit = {
if (c < 131072) {
buf(c >>> 5) |= 1 << (c & 31)
cull(c + p)
}
}
cull(s)
}
 
val bps = if (low <= 0) { //page zero has no base primes...
Iterator.iterate(3)(p => nextPrime(p + 2))
.takeWhile(_ <= Math.sqrt(262145).toInt).foreach(p => cullp((p * p - 3) >>> 1, p))
basePrimes
} else {
val nbps = if (basePrimesArr.length == 0) {
val nbasePrimes = SoEPg()
nbasePrimes.next() // new source primes > 2
basePrimesArr += nbasePrimes.next()
nbasePrimes
} else basePrimes
 
def fillBasePrimes(bp: Int): Int = {
if (bp * bp < nxtlowp) {
val np = nbps.next(); basePrimesArr += np
fillBasePrimes(np)
} else bp
}
 
fillBasePrimes(basePrimesArr(basePrimesArr.length - 1))
for (p <- basePrimesArr) {
val b = (p * p - 3) >>> 1
val s = if (b >= low) b - low else { val r = (low - b) % p; if (r != 0) p - r else 0 }
cullp(s, p)
}
nbps
}
(Iterator.iterate(nextPrime(lowp))(p => nextPrime(p + 2)).takeWhile(_ < nxtlowp), nxtlowp, bps)
}
 
Iterator.single(2) ++
Iterator.iterate(makePg(3, Iterator.empty)) { case (_, nlp, b) => makePg(nlp, b) }.
map{ case (itr, _, _) => itr }.flatten
}

As the above and all following sieves are "infinite", they all require an extra range limiting condition to produce a finite output, such as the addition of ".takeWhile(_ <= topLimit)" where "topLimit" is the specified range.

While the above code is reasonably fast, much of the execution time is consumed by the use of the built-in functions and iterators for concise code, especially in the use of iterators for primes output. A longer version using "roll-your-own" iteration or using a class with built in combined "higher-order" functions would be faster for many tasks using primes such as summing over a range, counting, et cetera.

[edit] Odds-Only "infinite" generator sieve using Streams and Co-Inductive Streams

The following code uses delayed recursion via Streams to implement the Richard Bird algorithm mentioned in the last part (the Epilogue) of M.O'Neill's paper, which is a true incremental Sieve of Eratosthenes. It is nowhere near as fast as the array based solutions due to the overhead of functionally chasing the merging of the prime multiple streams; this also means that the empirical performance is not according to the usual Sieve of Eratosthenes approximations due to this overhead increasing as the log of the sieved range, but it is much better than the "unfaithful" sieve.
  def oddPrimes: Stream[Int] = {
def merge(xs: Stream[Int], ys: Stream[Int]): Stream[Int] = {
val (x, y) = (xs.head, ys.head)
 
if (y > x) x #:: merge(xs.tail, ys) else if (x > y) y #:: merge(xs, ys.tail) else x #:: merge(xs.tail, ys.tail)
}
 
def primeMltpls(p: Int): Stream[Int] = Stream.iterate(p * p)(_ + p + p)
 
def allMltpls(ps: Stream[Int]): Stream[Stream[Int]] = primeMltpls(ps.head) #:: allMltpls(ps.tail)
 
def join(ams: Stream[Stream[Int]]): Stream[Int] = ams.head.head #:: merge(ams.head.tail, join(ams.tail))
 
def oddPrms(n: Int, composites: Stream[Int]): Stream[Int] =
if (n >= composites.head) oddPrms(n + 2, composites.tail) else n #:: oddPrms(n + 2, composites)
 
//following uses a new recursive source of odd base primes
3 #:: oddPrms(5, join(allMltpls(oddPrimes)))
}
 
val birdPrimes = 2 #:: oddPrimes
Now this algorithm doesn't really need the memoization and full laziness as offered by Streams, so an implementation and use of a Co-Inductive Stream (CIS) class is sufficient and reduces execution time by almost a factor of two:
  class CIS[A](val start: A, val continue: () => CIS[A])
 
def primesBirdCIS: Iterator[Int] = {
def merge(xs: CIS[Int], ys: CIS[Int]): CIS[Int] = {
val (x, y) = (xs.start, ys.start)
 
if (y > x) new CIS(x, () => merge(xs.continue(), ys))
else if (x > y) new CIS(y, () => merge(xs, ys.continue()))
else new CIS(x, () => merge(xs.continue(), ys.continue()))
}
 
def primeMltpls(p: Int): CIS[Int] = {
def nextCull(cull: Int): CIS[Int] = new CIS[Int](cull, () => nextCull(cull + 2 * p))
nextCull(p * p)
}
 
def allMltpls(ps: CIS[Int]): CIS[CIS[Int]] =
new CIS[CIS[Int]](primeMltpls(ps.start), () => allMltpls(ps.continue()))
def join(ams: CIS[CIS[Int]]): CIS[Int] = {
new CIS[Int](ams.start.start, () => merge(ams.start.continue(), join(ams.continue())))
}
 
def oddPrimes(): CIS[Int] = {
def oddPrms(n: Int, composites: CIS[Int]): CIS[Int] = { //"minua"
if (n >= composites.start) oddPrms(n + 2, composites.continue())
else new CIS[Int](n, () => oddPrms(n + 2, composites))
}
//following uses a new recursive source of odd base primes
new CIS(3, () => oddPrms(5, join(allMltpls(oddPrimes()))))
}
 
Iterator.single(2) ++ Iterator.iterate(oddPrimes())(_.continue()).map(_.start)
}
Further gains in performance for these last two implementations can be had by using further wheel factorization and "tree folding/merging" as per this Haskell implementation.

[edit] Odds-Only "infinite" generator sieve using a hash table (HashMap)

As per the "unfaithful sieve" article linked above, the incremental "infinite" Sieve of Eratosthenes can be implemented using a hash table instead of a Priority Queue or Map (Binary Heap) as were used in that article. The following implementation postpones the adding of base prime representations to the hash table until necessary to keep the hash table small:

  def SoEInc: Iterator[Int] = {
val nextComposites = scala.collection.mutable.HashMap[Int, Int]()
def oddPrimes: Iterator[Int] = {
val basePrimes = SoEInc
basePrimes.next()
basePrimes.next() // skip the two and three prime factors
@tailrec def makePrime(state: (Int, Int, Int)): (Int, Int, Int) = {
val (candidate, nextBasePrime, nextSquare) = state
if (candidate >= nextSquare) {
val adv = nextBasePrime << 1
nextComposites += ((nextSquare + adv) -> adv)
val np = basePrimes.next()
makePrime((candidate + 2, np, np * np))
} else if (nextComposites.contains(candidate)) {
val adv = nextComposites(candidate)
nextComposites -= (candidate) += (Iterator.iterate(candidate + adv)(_ + adv).
dropWhile(nextComposites.contains(_)).next() -> adv)
makePrime((candidate + 2, nextBasePrime, nextSquare))
} else (candidate, nextBasePrime, nextSquare)
}
Iterator.iterate((5, 3, 9)) { case (c, p, q) => makePrime((c + 2, p, q)) }
.map { case (p, _, _) => p }
}
List(2, 3).toIterator ++ oddPrimes
}

The above could be implemented using Streams or Co-Inductive Streams to pass the continuation parameters as passed here in a tuple but there would be no real difference in speed and there is no need to use the implied laziness. As compared to the versions of the Bird (or tree folding) Sieve of Eratosthenes, this has the expected same computational complexity as the array based versions, but is about 20 times slower due to the constant overhead of processing the key value hashing. Memory use is quite low, only being the hash table entries for each of the base prime values less than the square root of the last prime enumerated multiplied by the size of each hash entry (about 12 bytes in this case) plus a "load factor" percentage overhead in hash table size to minimize hash collisions (about twice as large as entries actually used by default on average).

The Scala implementable of a mutable HashMap is slower than the java.util.HashMap one by a factor of almost two, but the Scala version is used here to keep the code more portable (as to CLR). One can also quite easily convert this code to use the immutable Scala HashMap, but the code runs about four times slower due to the required "copy on update" operations for immutable objects.

This algorithm is very responsive to further application of wheel factorization, which can make it run up to about four times faster for the composite number culling operations; however, that is not enough to allow it to catch up to the array based sieves.

[edit] Scheme

Works with: Scheme version R5RS
; Tail-recursive solution :
(define (sieve n)
(define (aux u v)
(let ((p (car v)))
(if (> (* p p) n)
(let rev-append ((u u) (v v))
(if (null? u) v (rev-append (cdr u) (cons (car u) v))))
(aux (cons p u)
(let wheel ((u '()) (v (cdr v)) (a (* p p)))
(cond ((null? v) (reverse u))
((= (car v) a) (wheel u (cdr v) (+ a p)))
((> (car v) a) (wheel u v (+ a p)))
(else (wheel (cons (car v) u) (cdr v) a))))))))
(aux '(2)
(let range ((v '()) (k (if (odd? n) n (- n 1))))
(if (< k 3) v (range (cons k v) (- k 2))))))
 
; > (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)
; > (length (sieve 10000000))
; 664579
; Simpler solution, with the penalty that none of 'iota, 'strike or 'sieve is tail-recursive :
(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)

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)

Optimised using a pre-computed wheel based on 2 (i.e. odds only):

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

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)

Vector-based (faster), works with R5RS:

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

Using SICP-style head-forced streams. Works with MIT-Scheme.

 ;;;; Stream Implementation
(define (head s) (car s))
(define (tail s) ((cdr s)))
(define-syntax s-cons
(syntax-rules () ((s-cons h t) (cons h (lambda () t)))))
 
;;;; Stream Utility Functions
(define (from-By x s)
(s-cons x (from-By (+ x s) s)))
(define (take n s)
(cond
((> n 1) (cons (head s) (take (- n 1) (tail s))))
((= n 1) (list (head s))) ;; don't force it too soon!!
(else ()))) ;; so (take 4 (s-map / (from-By 4 -1))) works
(define (drop n s)
(cond
((> n 0) (drop (- n 1) (tail s)))
(else s)))
(define (s-map f s)
(s-cons (f (head s)) (s-map f (tail s))))
(define (s-diff s1 s2)
(let ((h1 (head s1)) (h2 (head s2)))
(cond
((< h1 h2) (s-cons h1 (s-diff (tail s1) s2 )))
((< h2 h1) (s-diff s1 (tail s2)))
(else (s-diff (tail s1) (tail s2))))))
(define (s-union s1 s2)
(let ((h1 (head s1)) (h2 (head s2)))
(cond
((< h1 h2) (s-cons h1 (s-union (tail s1) s2 )))
((< h2 h1) (s-cons h2 (s-union s1 (tail s2))))
(else (s-cons h1 (s-union (tail s1) (tail s2)))))))
 
;;;; all primes' multiples are removed, merged through a tree of unions
;;;; runs in ~ n^1.15 run time in producing n = 100K .. 1M primes
(define (primes-stream)
(define (mults p) (from-By (* p p) (* 2 p)))
(define (no-mults-From from)
(s-diff (from-By from 2)
(s-tree-join (s-map mults odd-primes))))
(define odd-primes
(s-cons 3 (no-mults-From 5)))
(s-cons 2 (no-mults-From 3)))
 
;;;; join an ordered stream of streams (here, of primes' multiples)
;;;; into one ordered stream, via an infinite right-deepening tree
(define (s-tree-join sts) ;; sts -> s
(define (join-With of-Tail sts) ;; sts -> s
(s-cons (head (head sts))
(s-union (tail (head sts)) (of-Tail (tail sts)))))
(define (pairs sts) ;; sts -> sts
(s-cons (join-With head sts) (pairs (tail (tail sts)))))
(join-With (lambda (t) (s-tree-join (pairs t))) sts))

Print 10 last primes of the first thousand primes:

(display (take 10 (drop 990 (primes-stream)))) 
;
(7841 7853 7867 7873 7877 7879 7883 7901 7907 7919)

[edit] Seed7

The program below computes the number of primes between 1 and 10000000:

$ include "seed7_05.s7i";
 
const func set of integer: eratosthenes (in integer: n) is func
result
var set of integer: sieve is EMPTY_SET;
local
var integer: i is 0;
var integer: j is 0;
begin
sieve := {2 .. n};
for i range 2 to sqrt(n) do
if i in sieve then
for j range i ** 2 to n step i do
excl(sieve, j);
end for;
end if;
end for;
end func;
 
const proc: main is func
begin
writeln(card(eratosthenes(10000000)));
end func;

Original source: [1]

[edit] SNOBOL4

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

        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

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

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

[edit] UNIX Shell

[edit] With array

Works with: zsh
#!/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
Works with: bash
Works with: ksh93
Works with: pdksh
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

Both scripts output a single long line.

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

[edit] 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
#! /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`

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

[edit] C Shell

Translation of: CMake
# 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

[edit] Ursala

This example is incorrect. It probably (remainder) uses rem testing and so is a trial division algorithm, not a sieve of Eratosthenes. Please fix the code and remove this message.

with no optimizations

#import nat
 
sieve = ~<{0,1}&& iota; @NttPX ~&r->lx ^/~&rhPlC remainder@rlX~|@r

test program:

#cast %nL
 
example = sieve 50
Output:
<2,3,5,7,11,13,17,19,23,29,31,37,41,43,47>

[edit] Vala

Library: Gee
Without any optimizations:
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");
}
{{out}
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47

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

[edit] Visual Basic .NET

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 = 2 To Math.Sqrt(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
Console.WriteLine(n)
End If
Next n

[edit] 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)
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, 

[edit] XPL0

include c:\cxpl\codes;                  \intrinsic 'code' declarations
int Size, Prime, I, Kill;
char Flag;
[Size:= IntIn(0);
Flag:= Reserve(Size+1);
for I:= 2 to Size do Flag(I):= true;
for I:= 2 to Size do
if Flag(I) then \found a prime
[Prime:= I;
IntOut(0, Prime); CrLf(0);
Kill:= Prime + Prime; \first multiple to kill
while Kill <= Size do
[Flag(Kill):= false; \zero a non-prime
Kill:= Kill + Prime; \next multiple
];
];
]
Example output:
20

2 3 5 7 11 13 17

19

[edit] zkl

fcn sieve(limit){
if (limit<2) return(T);
composite:=(0).pump(limit+1,Data,1); // bucket of bytes set to 1 (prime)
(2).filter(limit.toFloat().sqrt()+1,T(composite.__sGet, // if prime, zero multiples
'wrap(n){ [n*n .. limit,n].pump(Void,composite.__sSet.fp(0)) }, //composite[n*p]=0
False)); // turn filter into a no-result loop
(2).filter(limit-1,composite.__sGet); // bytes still 1 are prime
}
sieve(53).println();

The filter method, when given multiple filters, acts like a conditional and. Here, the first filter checks the table for a prime, if so, the second filter does some side effects and the third filter ensures that no items make it through the filter (False(<anything>)-->False) so that the filter returns an empty list, minimizing garbage.

Output:
L(2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53)
Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox