Sieve of Eratosthenes
You are encouraged to solve this task according to the task description, using any language you may know.
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
68000 Assembly
Algorithm somewhat optimized: array omits 1, 2, all higher odd numbers. Optimized for storage: uses bit array for prime/composite flags.
Some of the macro code is derived from the examples included with EASy68K. See 68000 "100 Doors" listing for additional information. <lang 68000>*-----------------------------------------------------------
- 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</lang>
ACL2
<lang Lisp>(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))</lang>
Ada
<lang 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;</lang>
Output:
> ./eratos 31 Primes less or equal 31 are : 2 3 5 7 11 13 17 19 23 29 31
ALGOL 68
<lang algol68>BOOL prime = TRUE, non prime = FALSE; PROC eratosthenes = (INT n)[]BOOL: (
[n]BOOL sieve;
FOR i TO UPB sieve DO sieve[i] := prime OD;
INT m = ENTIER sqrt(n);
sieve[1] := non prime;
FOR i FROM 2 TO m DO
IF sieve[i] = prime THEN
FOR j FROM i*i BY i TO n DO
sieve[j] := non prime
OD
FI
OD;
sieve
);
print((eratosthenes(80),new line))</lang>
Output:
FTTFTFTFFFTFTFFFTFTFFFTFFFFFTFTFFFFFTFFFTFTFFFTFFFFFTFFFFFTFTFFFFFTFFFTFTFFFFFTF
AutoHotkey
Search autohotkey.com: of Eratosthenes
Source: AutoHotkey forum by Laszlo
<lang autohotkey>
MsgBox % "12345678901234567890`n" Sieve(20)
Sieve(n) { ; Sieve of Eratosthenes => string of 0|1 chars, 1 at position k: k is prime
Static zero := 48, one := 49 ; Asc("0"), Asc("1")
VarSetCapacity(S,n,one)
NumPut(zero,S,0,"char")
i := 2
Loop % sqrt(n)-1 {
If (NumGet(S,i-1,"char") = one)
Loop % n//i
If (A_Index > 1)
NumPut(zero,S,A_Index*i-1,"char")
i += 1+(i>2)
}
Return S
} </lang>
AutoIt
<lang 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) </lang>
AWK
An initial array holds all numbers 2..max (which is entered on stdin); then all products of integers are deleted from it; the remaining are displayed in the unsorted appearance of a hash table.
$ awk '{for(i=2;i<=$1;i++) a[i]=1;
> for(i=2;i<=sqrt($1);i++) for(j=2;j<=$1;j++) delete a[i*j];
> for(i in a) printf i" "}'
100
71 53 17 5 73 37 19 83 47 29 7 67 59 11 97 79 89 31 13 41 23 2 61 43 3
The following variant does not unset non-primes, but set them to 0, to preserve order in output:
$ awk '{for(i=2;i<=$1;i++) a[i]=1;
> for(i=2;i<=sqrt($1);i++) for(j=2;j<=$1;j++) a[i*j]=0;
> for(i=2;i<=$1;i++) if(a[i])printf i" "}'
100
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
BASIC
<lang freebasic> DIM n AS Integer, k AS Integer, limit AS Integer
INPUT "Enter number to search to: "; limit DIM flags(limit) AS Integer
FOR n = 2 TO SQR(limit)
IF flags(n) = 0 THEN
FOR k = n*n TO limit STEP n
flags(k) = 1
NEXT k
END IF
NEXT n
' Display the primes FOR n = 2 TO limit
IF flags(n) = 0 THEN PRINT n; ", ";
NEXT n </lang>
Applesoft BASIC
<lang basic> 10 INPUT "ENTER NUMBER TO SEARCH TO: ";LIMIT 20 DIM FLAGS(LIMIT) 30 FOR N = 2 TO SQR (LIMIT) 40 IF FLAGS(N) < > 0 GOTO 80 50 FOR K = N * N TO LIMIT STEP N 60 FLAGS(K) = 1 70 NEXT K 80 NEXT N 90 REM DISPLAY THE PRIMES 100 FOR N = 2 TO LIMIT 110 IF FLAGS(N) = 0 THEN PRINT N;", "; 120 NEXT N </lang>
Locomotive Basic
<lang locobasic> 10 DEFINT a-z 20 INPUT "Limit";limit 30 DIM f(limit) 40 FOR n=2 TO SQR(limit) 50 IF f(n)=1 THEN 90 60 FOR k=n*n TO limit STEP n 70 f(k)=1 80 NEXT k 90 NEXT n 100 FOR n=2 TO limit 110 IF f(n)=0 THEN PRINT n;","; 120 NEXT </lang>
ZX Spectrum Basic
<lang 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 </lang>
BBC BASIC
<lang bbcbasic> limit% = 100000
DIM sieve% limit%
prime% = 2
WHILE prime%^2 < limit%
FOR I% = prime%*2 TO limit% STEP prime%
sieve%?I% = 1
NEXT
REPEAT prime% += 1 : UNTIL sieve%?prime%=0
ENDWHILE
REM Display the primes:
FOR I% = 1 TO limit%
IF sieve%?I% = 0 PRINT I%;
NEXT</lang>
bash
See solutions at UNIX Shell.
Befunge
2>:3g" "-!v\ g30 < |!`"O":+1_:.:03p>03g+:"O"`| @ ^ p3\" ":< 2 234567890123456789012345678901234567890123456789012345678901234567890123456789
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. <lang bracmat>( ( 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 )</lang> Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
C
Plain sieve, without any optimizations:
<lang c>
- 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;
} </lang>
Possible optimizations include sieving only odd numbers (or more complex wheels), packing the sieve into bits to improve locality (and allow larger sieves), etc.
C++
<lang cpp> // yield all prime numbers less than limit. template<class UnaryFunction> void primesupto(int limit, UnaryFunction yield) {
std::vector<bool> is_prime(limit, true);
const int sqrt_limit = static_cast<int>(std::sqrt(limit));
for (int n = 2; n <= sqrt_limit; ++n)
if (is_prime[n]) {
yield(n);
for (unsigned k = n*n, ulim = static_cast<unsigned>(limit); k < ulim; k += n)
//NOTE: "unsigned" is used to avoid an overflow in `k+=n` for `limit` near INT_MAX
is_prime[k] = false;
}
for (int n = sqrt_limit + 1; n < limit; ++n) if (is_prime[n])
yield(n); } </lang>
Full program:
<lang cpp>/**
$ g++ -I/path/to/boost sieve.cpp -o sieve && sieve 10000000 */
- include <inttypes.h> // uintmax_t
- include <limits>
- include <cmath>
- include <iostream>
- include <sstream>
- include <vector>
- include <boost/lambda/lambda.hpp>
int main(int argc, char *argv[])
{
using namespace std; using namespace boost::lambda;
int limit = 10000;
if (argc == 2) {
stringstream ss(argv[--argc]);
ss >> limit;
if (limit < 1 or ss.fail()) {
cerr << "USAGE:\n sieve LIMIT\n\nwhere LIMIT in the range [1, "
<< numeric_limits<int>::max() << ")" << endl;
return 2; } }
// print primes less then 100 primesupto(100, cout << _1 << " "); cout << endl;
// find number of primes less then limit and their sum int count = 0; uintmax_t sum = 0; primesupto(limit, (var(sum) += _1, var(count) += 1));
cout << "limit sum pi(n)\n"
<< limit << " " << sum << " " << count << endl;
} </lang>
C#
<lang csharp>using System; using System.Collections; using System.Collections.Generic; using System.Linq; using System.Text;
namespace sieve {
class Program
{
static void Main(string[] args)
{
int maxprime = int.Parse(args[0]);
ArrayList primelist = sieve(maxprime);
foreach (int prime in primelist)
Console.WriteLine(prime);
Console.WriteLine("Count = " + primelist.Count);
Console.ReadLine();
}
static ArrayList sieve(int arg_max_prime)
{
BitArray al = new BitArray(arg_max_prime, true);
int lastprime = 1;
int lastprimeSquare = 1;
while (lastprimeSquare <= arg_max_prime)
{
lastprime++;
while (!(bool)al[lastprime])
lastprime++;
lastprimeSquare = lastprime * lastprime;
for (int i = lastprimeSquare; i < arg_max_prime; i += lastprime)
if (i > 0)
al[i] = false;
}
ArrayList sieve_2_return = new ArrayList();
for (int i = 2; i < arg_max_prime; i++)
if (al[i])
sieve_2_return.Add(i);
return sieve_2_return;
}
}
}</lang>
Chapel
This solution uses nested iterators to create new wheels at run time: <lang chapel>// 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;
}
}</lang>
The topmost sieve needs to be started with 2 (the smallest prime): <lang chapel>config const N = 30; for p in sieve(2) {
write(" ", p);
if p > N then {
writeln();
break;
}
}</lang>
Clojure
Calculates primes up to and including n. <lang clojure>(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)))))))</lang>
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: <lang clojure>(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)))))</lang>
Or as a lazy sequence (slightly faster): <lang clojure>(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)])))))))</lang>
CMake
<lang cmake>function(eratosthenes var limit)
# Check for integer overflow. With CMake using 32-bit signed integer,
# this check fails when limit > 46340.
if(NOT limit EQUAL 0) # Avoid division by zero.
math(EXPR i "(${limit} * ${limit}) / ${limit}")
if(NOT limit EQUAL ${i})
message(FATAL_ERROR "limit is too large, would cause integer overflow")
endif()
endif()
# Use local variables prime_2, prime_3, ..., prime_${limit} as array.
# Initialize array to y => yes it is prime.
foreach(i RANGE 2 ${limit})
set(prime_${i} y)
endforeach(i)
# Gather a list of prime numbers.
set(list)
foreach(i RANGE 2 ${limit})
if(prime_${i})
# Append this prime to list.
list(APPEND list ${i})
# For each multiple of i, set n => no it is not prime.
# Optimization: start at i squared.
math(EXPR square "${i} * ${i}")
if(NOT square GREATER ${limit}) # Avoid fatal error.
foreach(m RANGE ${square} ${limit} ${i})
set(prime_${m} n)
endforeach(m)
endif()
endif(prime_${i})
endforeach(i)
set(${var} ${list} PARENT_SCOPE)
endfunction(eratosthenes)</lang>
<lang cmake># Print all prime numbers through 100. eratosthenes(primes 100) message(STATUS "${primes}")</lang>
COBOL
<lang 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
.</lang>
Common Lisp
<lang lisp>(defun sieve-of-eratosthenes (maximum)
(let ((sieve (make-array (1+ maximum) :element-type 'bit
:initial-element 0)))
(loop for candidate from 2 to maximum
when (zerop (bit sieve candidate))
collect candidate
and do (loop for composite from (expt candidate 2)
to maximum by candidate
do (setf (bit sieve composite) 1)))))</lang>
Working with odds only (above twice speedup), and only test divide up to the square root of the maximum:
<lang lisp>(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)))))))</lang>
While formally a wheel, odds are uniformly spaced and do not require any special processing except for value translation. Wheels proper aren't uniformly spaced and are thus trickier.
D
Simpler Version
Prints all numbers less than the limit. <lang d>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;
}</lang>
- Output:
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
Faster Version
This version uses an array of bits (instead of booleans, that are represented with one byte), and skips even numbers. The output is the same. <lang d>import std.stdio, std.math, std.array;
size_t[] 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;
}</lang>
Extensible Version
(This version is used in the task Extensible prime generator.) <lang d>/// 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;
}
}</lang>
To see the output (that is the same), compile with -version=sieve_of_eratosthenes3_main.
Dart
<lang dart>// helper function to pretty print an Iterable String iterableToString(Iterable seq) {
String str = "[";
Iterator i = seq.iterator();
while(i.hasNext()) {
str += i.next();
if(i.hasNext()) {
str += ",";
}
}
return str + "]";
}
main() {
int limit = 1000; Set<int> sieve = new Set<int>();
for(int i = 2; i <= limit; i++) {
sieve.add(i);
}
for(int i = 2; i * i <= limit; i++) {
if(sieve.contains(i)) {
for(int j = i * i; j <= limit; j += i) {
sieve.remove(j);
}
}
}
var sortedValues = new List<int>.from(sieve);
sortedValues.sort((int arg1, int arg2) {
return arg1.compareTo(arg2);
});
print(iterableToString(sortedValues));
Expect.equals(168, sieve.length);
}</lang> Output:
[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997]
Delphi
<lang 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.</lang> Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
DWScript
<lang delphi>function Primes(limit : Integer) : array of Integer; var
n, k : Integer; sieve := new Boolean[limit+1];
begin
for n := 2 to Round(Sqrt(limit)) do begin
if not sieve[n] then begin
for k := n*n to limit step n do
sieve[k] := True;
end;
end;
for k:=2 to limit do
if not sieve[k] then
Result.Add(k);
end;
var r := Primes(50); var i : Integer; for i:=0 to r.High do
PrintLn(r[i]);</lang>
Dylan
With outer to sqrt and inner to p^2 optimizations: <lang dylan>define method primes(n)
let limit = floor(n ^ 0.5) + 1; let sieve = make(limited(<simple-vector>, of: <boolean>), size: n + 1, fill: #t); let last-prime = 2;
while (last-prime < limit)
for (x from last-prime ^ 2 to n by last-prime)
sieve[x] := #f;
end for;
block (found-prime)
for (n from last-prime + 1 below limit)
if (sieve[n] = #f)
last-prime := n;
found-prime()
end;
end;
last-prime := limit;
end block;
end while;
for (x from 2 to n)
if (sieve[x]) format-out("Prime: %d\n", x); end;
end;
end;</lang>
E
E's standard library doesn't have a step-by-N numeric range, so we'll define one, implementing the standard iteration protocol.
def rangeFromBelowBy(start, limit, step) {
return def stepper {
to iterate(f) {
var i := start
while (i < limit) {
f(null, i)
i += step
}
}
}
}
The sieve itself:
def eratosthenes(limit :(int > 2), output) {
def composite := [].asSet().diverge()
for i ? (!composite.contains(i)) in 2..!limit {
output(i)
composite.addAll(rangeFromBelowBy(i ** 2, limit, i))
}
}
Example usage:
? eratosthenes(12, println) # stdout: 2 # 3 # 5 # 7 # 11
eC
Note: this is not a Sieve of Eratosthenes; it is just trial division. <lang cpp> public class FindPrime {
Array<int> primeList { [ 2 ], minAllocSize = 64 };
int index;
index = 3;
bool HasPrimeFactor(int x)
{
int max = (int)floor(sqrt((double)x));
for(i : primeList)
{
if(i > max) break;
if(x % i == 0) return true;
}
return false;
}
public int GetPrime(int x)
{
if(x > primeList.count - 1)
{
for (; primeList.count != x; index += 2)
if(!HasPrimeFactor(index))
{
if(primeList.count >= primeList.minAllocSize) primeList.minAllocSize *= 2;
primeList.Add(index);
}
}
return primeList[x-1];
}
}
class PrimeApp : Application {
FindPrime fp { };
void Main()
{
int num = argc > 1 ? atoi(argv[1]) : 1;
PrintLn(fp.GetPrime(num));
}
} </lang>
Eiffel
<lang eiffel>class
APPLICATION
create
make
feature
make
-- Run application.
do
across primes_through (100) as ic loop print (ic.item.out + " ") end
end
primes_through (a_limit: INTEGER): LINKED_LIST [INTEGER]
-- Prime numbers through `a_limit'
require
valid_upper_limit: a_limit >= 2
local
l_tab: ARRAY [BOOLEAN]
do
create Result.make
create l_tab.make_filled (True, 2, a_limit)
across
l_tab as ic
loop
if ic.item then
Result.extend (ic.target_index)
across ((ic.target_index * ic.target_index) |..| l_tab.upper).new_cursor.with_step (ic.target_index) as id
loop
l_tab [id.item] := False
end
end
end
end
end</lang>
Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
Emacs Lisp
<lang 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)))))
</lang>
Straightforward implementation of sieve of Eratosthenes, 2 times faster:
<lang lisp> (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)))
</lang>
Erlang
<lang 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) )}. </lang>
- Output:
35> sieve_of_eratosthenes:primes_upto( 20 ). [2,3,5,7,11,13,17,19]
Euphoria
<lang euphoria>constant limit = 1000 sequence flags,primes flags = repeat(1, limit) for i = 2 to sqrt(limit) do
if flags[i] then
for k = i*i to limit by i do
flags[k] = 0
end for
end if
end for
primes = {} for i = 2 to limit do
if flags[i] = 1 then
primes &= i
end if
end for ? primes</lang>
Output:
{2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,
97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,
181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,
277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,
383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,
487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,
601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,
709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,
827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,
947,953,967,971,977,983,991,997}
F#
Functional
This is 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:
<lang fsharp>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)</lang>
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:
<lang fsharp>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) }</lang>
Imperative
The following code is written in functional style other than it uses a mutable bit array to sieve the composites:
<lang fsharp>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</lang>
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:
<lang fsharp>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 }</lang>
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:
<lang fsharp>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 }</lang>
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:
<lang fsharp> 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</lang>
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):
<lang fsharp> 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 }</lang>
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.
Forth
: prime? ( n -- ? ) here + c@ 0= ;
: composite! ( n -- ) here + 1 swap c! ;
: sieve ( n -- )
here over erase
2
begin
2dup dup * >
while
dup prime? if
2dup dup * do
i composite!
dup +loop
then
1+
repeat
drop
." Primes: " 2 do i prime? if i . then loop ;
100 sieve
Fortran
<lang fortran>program sieve
implicit none integer, parameter :: i_max = 100 integer :: i logical, dimension (i_max) :: is_prime
is_prime = .true. is_prime (1) = .false. do i = 2, int (sqrt (real (i_max))) if (is_prime (i)) is_prime (i * i : i_max : i) = .false. end do do i = 1, i_max if (is_prime (i)) write (*, '(i0, 1x)', advance = 'no') i end do write (*, *)
end program sieve</lang> Output: <lang>2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97</lang> Optimised using a pre-computed wheel based on 2: <lang fortran>program sieve_wheel_2
implicit none integer, parameter :: i_max = 100 integer :: i logical, dimension (i_max) :: is_prime
is_prime = .true. is_prime (1) = .false. is_prime (4 : i_max : 2) = .false. do i = 3, int (sqrt (real (i_max))), 2 if (is_prime (i)) is_prime (i * i : i_max : 2 * i) = .false. end do do i = 1, i_max if (is_prime (i)) write (*, '(i0, 1x)', advance = 'no') i end do write (*, *)
end program sieve_wheel_2</lang> Output: <lang>2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97</lang>
GAP
<lang gap>Eratosthenes := function(n)
local sieve, cur, mul;
sieve := [1 .. n]*0;
sieve[1] := 1;
cur := 1;
while cur <= n do
if sieve[cur] = 0 then
mul := cur*cur;
while mul <= n do
sieve[mul] := 1;
mul := mul + cur;
od;
fi;
cur := cur + 1;
od;
return Filtered([1 .. n], x -> sieve[x] = 0);
end;
Eratosthenes(100);
- [ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 ]</lang>
GLBasic
<lang GLBasic>// Sieve of Eratosthenes (find primes) // GLBasic implementation
GLOBAL n%, k%, limit%, flags%[]
limit = 100 // search primes up to this number
DIM flags[limit+1] // GLBasic arrays start at 0
FOR n = 2 TO SQR(limit)
IF flags[n] = 0
FOR k = n*n TO limit STEP n
flags[k] = 1
NEXT
ENDIF
NEXT
// Display the primes FOR n = 2 TO limit
IF flags[n] = 0 THEN STDOUT n + ", "
NEXT
KEYWAIT </lang>
Go
<lang go>package main
import (
"fmt"
)
func main() {
const limit = 201 // means sieve numbers < 201
// sieve
c := make([]bool, limit) // c for composite. false means prime candidate
c[1] = true // 1 not considered prime
p := 2
for {
// first allowed optimization: outer loop only goes to sqrt(limit)
p2 := p * p
if p2 >= limit {
break
}
// second allowed optimization: inner loop starts at sqr(p)
for i := p2; i < limit; i += p {
c[i] = true // it's a composite
}
// scan to get next prime for outer loop
for {
p++
if !c[p] {
break
}
}
}
// sieve complete. now print a representation.
for n := 1; n < limit; n++ {
if c[n] {
fmt.Print(" .")
} else {
fmt.Printf("%3d", n)
}
if n%20 == 0 {
fmt.Println("")
}
}
}</lang> Output:
. 2 3 . 5 . 7 . . . 11 . 13 . . . 17 . 19 . . . 23 . . . . . 29 . 31 . . . . . 37 . . . 41 . 43 . . . 47 . . . . . 53 . . . . . 59 . 61 . . . . . 67 . . . 71 . 73 . . . . . 79 . . . 83 . . . . . 89 . . . . . . . 97 . . . 101 .103 . . .107 .109 . . .113 . . . . . . . . . . . . .127 . . .131 . . . . .137 .139 . . . . . . . . .149 .151 . . . . .157 . . . . .163 . . .167 . . . . .173 . . . . .179 . 181 . . . . . . . . .191 .193 . . .197 .199 .
A fairly odd sieve tree method: <lang go>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())
}
}</lang>
Example with channels:
<lang go>
package main
import ( "fmt" )
func filter(n int, in chan int) chan int { out := make(chan int) go func() { for i := range in { if i%n != 0 { out <- i } } close(out) }()
return out }
func ints(max int) chan int { out := make(chan int) go func() { for i := 2; i <= max; i++ { out <- i } close(out) }()
return out }
func main() { ch := ints(201) for { i := <-ch if i == 0 { break } fmt.Println(i) ch = filter(i, ch) } } </lang>
Groovy
This solution uses a BitSet for compactness and speed, but in Groovy, BitSet has full List semantics. It also uses both the "square root of the boundary" shortcut and the "square of the prime" shortcut. <lang groovy>def sievePrimes = { bound ->
def isPrime = new BitSet(bound)
isPrime[0..1] = false
isPrime[2..bound] = true
(2..(Math.sqrt(bound))).each { pc ->
if (isPrime[pc]) {
((pc**2)..bound).step(pc) { isPrime[it] = false }
}
}
(0..bound).findAll { isPrime[it] }
}</lang>
Test: <lang groovy>println sievePrimes(100)</lang>
Output:
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
GW-BASIC
<lang qbasic>10 INPUT "ENTER NUMBER TO SEARCH TO: ";LIMIT 20 DIM FLAGS(LIMIT) 30 FOR N = 2 TO SQR (LIMIT) 40 IF FLAGS(N) < > 0 GOTO 80 50 FOR K = N * N TO LIMIT STEP N 60 FLAGS(K) = 1 70 NEXT K 80 NEXT N 90 REM DISPLAY THE PRIMES 100 FOR N = 2 TO LIMIT 110 IF FLAGS(N) = 0 THEN PRINT N;", "; 120 NEXT N</lang>
Haskell
Mutable unboxed arrays, odds only
Mutable array of unboxed Bools indexed by Ints, representing odds only:
<lang haskell>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 = []</lang>
This represents odds only in the array. Empirical orders of growth is 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.
Immutable arrays
This is the most straightforward: <lang haskell>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</lang>
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.
Immutable arrays, by segments
Works by segments between consecutive primes' squares. Should be the fastest of non-monadic code: <lang haskell>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]]</lang>
Basic list-based sieve
Straightforward implementation of the sieve of Ertosthenes in its original bounded form: <lang haskell>primesTo m = 2 : eratos [3,5..m] where
eratos (p : xs) | p*p>m = p : xs
| True = p : eratos (xs `minus` [p*p, p*p+2*p..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 </lang>
Its time complexity is similar to that of optimal trial division because of limitations of Haskell linked lists, where (minus a b) takes time proportional to length(union a b) and not (length b) as achieved in imperative setting with direct-access memory. Uses ordered list representation of sets.
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.
<lang haskell>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</lang>
List-based tree-merging incremental sieve
Linear merging structure can be further replaced with indefinitely deepening to the right tree-like structure (a-(b+((c+d)+( ((e+f)+(g+h)) + ... )))).
This finds primes in gaps between composites, and composites as union of the multiples of primes, merging in a tree-like fashion all the streams of primes' multiples, running with empirical orders of growth of 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): <lang haskell>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), when
| otherwise = gaps (k+2) xs -- k<=x && null(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</lang>
Here's the test entry on Ideone.com, a comparison with more versions, and a similar code with wheel optimization.
HicEst
<lang hicest>REAL :: N=100, sieve(N)
sieve = $ > 1 ! = 0 1 1 1 1 ... DO i = 1, N^0.5
IF( sieve(i) ) THEN
DO j = i^2, N, i
sieve(j) = 0
ENDDO
ENDIF
ENDDO
DO i = 1, N
IF( sieve(i) ) WRITE() i
ENDDO </lang>
Icon and Unicon
<lang Icon> procedure main()
sieve(100) end
procedure sieve(n)
local p,i,j
p:=list(n, 1)
every i:=2 to sqrt(n) & j:= i+i to n by i & p[i] == 1
do p[j] := 0
every write(i:=2 to n & p[i] == 1 & i)
end</lang>
Alternatively using sets <lang Icon> procedure main()
sieve(100) end
procedure sieve(n)
primes := set()
every insert(primes,1 to n)
every member(primes,i := 2 to n) do
every delete(primes,i + i to n by i)
delete(primes,1)
every write(!sort(primes))
end</lang>
J
sieve=: 3 : 0
m=. <.%:y
z=. $0
b=. y{.1
while. m>:j=. 1+b i. 0 do.
b=. b+.y$(-j){.1
z=. z,j
end.
z,1+I.-.b
)
"Wheels" may be implemented as follows:
sieve1=: 3 : 0
m=. <.%:y
z=. y (>:#]) 2 3 5 7
b=. 1,}.y$+./(*/z)$&>(-z){.&.>1
while. m>:j=. 1+b i. 0 do.
b=. b+.y$(-j){.1
z=. z,j
end.
z,1+I.-.b
)
The use of 2 3 5 7 as wheels provides a 20% time improvement for n=1000 and 2% for n=1e6 .
Java
<lang java5>import java.util.LinkedList;
public class Sieve{
public static LinkedList<Integer> sieve(int n){
if(n < 2) return new LinkedList<Integer>();
LinkedList<Integer> primes = new LinkedList<Integer>();
LinkedList<Integer> nums = new LinkedList<Integer>();
for(int i = 2;i <= n;i++){ //unoptimized
nums.add(i);
}
while(nums.size() > 0){
int nextPrime = nums.remove();
for(int i = nextPrime * nextPrime;i <= n;i += nextPrime){
nums.removeFirstOccurrence(i);
}
primes.add(nextPrime);
}
return primes;
}
}</lang>
To optimize by testing only odd numbers, replace the loop marked "unoptimized" with these lines: <lang java5>nums.add(2); for(int i = 3;i <= n;i += 2){
nums.add(i);
}</lang>
Version using a BitSet: <lang java5>import java.util.LinkedList; import java.util.BitSet;
public class Sieve{
public static LinkedList<Integer> sieve(int n){
LinkedList<Integer> primes = new LinkedList<Integer>();
BitSet nonPrimes = new BitSet(n+1);
for (int p = 2; p <= n ; p = nonPrimes.nextClearBit(p+1)) {
for (int i = p * p; i <= n; i += p)
nonPrimes.set(i);
primes.add(p);
}
return primes;
}
}</lang>
Infinite iterator
An iterator that will generate primes indefinitely (perhaps until it runs out of memory).
<lang java5>import java.util.Iterator; import java.util.PriorityQueue;
// generates all prime numbers public class InfiniteSieve implements Iterator<Integer> {
private static class NonPrimeSequence implements Comparable<NonPrimeSequence> {
int currentMultiple;
int prime;
public NonPrimeSequence(int p) {
prime = p;
currentMultiple = p*p; // start at square of prime
}
public int compareTo(NonPrimeSequence other) {
// sorted by value of current multiple
return ((Integer)currentMultiple).compareTo(other.currentMultiple);
}
}
private int i = 2; // priority queue of the sequences of non-primes // the priority queue allows us to get the "next" non-prime quickly private PriorityQueue<NonPrimeSequence> nonprimes = new PriorityQueue<NonPrimeSequence>();
public boolean hasNext() { return true; }
public Integer next() {
// skip non-prime numbers
for ( ; !nonprimes.isEmpty() && i == nonprimes.peek().currentMultiple; 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
while (nonprimes.peek().currentMultiple == i) {
NonPrimeSequence x = nonprimes.poll();
x.currentMultiple += x.prime;
nonprimes.offer(x);
}
}
// prime
// insert a NonPrimeSequence object into the priority queue
nonprimes.offer(new NonPrimeSequence(i));
return i++;
}
public static void main(String[] args) {
Iterator<Integer> sieve = new InfiniteSieve();
for (int i = 0; i < 8; i++) {
System.out.println(sieve.next());
}
}
}</lang>
- Output:
2 3 5 7 11 13 17 19
JavaScript
<lang 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);</lang> outputs:
2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97
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:
<lang javascript>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;
}</lang>
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:
<lang javascript>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;
})();</lang>
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:
<lang javascript>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);</lang>
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:
<lang javascript>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;
})();</lang>
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.
Liberty BASIC
<lang lb> 'Notice that arrays are globally visible to functions.
'The sieve() function uses the flags() array. 'This is a Sieve benchmark adapted from BYTE 1985 ' May, page 286
size = 7000
dim flags(7001)
start = time$("ms")
print sieve(size); " primes found."
print "End of iteration. Elapsed time in milliseconds: "; time$("ms")-start
end
function sieve(size)
for i = 0 to size
if flags(i) = 0 then
prime = i + i + 3
k = i + prime
while k <= size
flags(k) = 1
k = k + prime
wend
sieve = sieve + 1
end if
next i
end function</lang>
Logo
to sieve :limit
make "a (array :limit 2) ; initialized to empty lists
make "p []
for [i 2 :limit] [
if empty? item :i :a [
queue "p :i
for [j [:i * :i] :limit :i] [setitem :j :a :i]
]
]
output :p
end
print sieve 100 ; 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
Lua
<lang lua>function erato(n)
local t = {0, 2}
for i = 3, n, 2 do t[i], t[i+1] = i, 0 end
for i = 3, math.sqrt(n) do for j = i*i, n, 2*i do t[j] = 0 end end
return t
end</lang>
Lucid
prime
where
prime = 2 fby (n whenever isprime(n));
n = 3 fby n+2;
isprime(n) = not(divs) asa divs or prime*prime > N
where
N is current n;
divs = N mod prime eq 0;
end;
end
recursive
sieve( N )
where
N = 2 fby N + 1;
sieve( i ) =
i fby sieve ( i whenever i mod first i ne 0 ) ;
end
M4
<lang M4> define(`lim',100)dnl define(`for',
`ifelse($#,0,
``$0,
`ifelse(eval($2<=$3),1,
`pushdef(`$1',$2)$5`'popdef(`$1')$0(`$1',eval($2+$4),$3,$4,`$5')')')')dnl
for(`j',2,lim,1,
`ifdef(a[j],
`',
`j for(`k',eval(j*j),lim,j,
`define(a[k],1)')')')
</lang>
Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
Mathematica
<lang Mathematica> Eratosthenes[n_] := Module[{numbers = Range[n]},
Do[If[numbersi != 0, Do[numbersi j = 0, {j, 2, n/i}]], {i, 2, Sqrt[n]}]; Select[numbers, # > 1 &]]
Eratosthenes[100] </lang>
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: <lang Mathematica> Eratosthenes[n_] := Module[{numbers = Range[n]},
Do[If[numbersi != 0, Do[numbersj = 0, {j,i i,n,i}]],{i,2,Sqrt[n]}]; Select[numbers, # > 1 &]]
Eratosthenes[100] </lang>
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: <lang Mathematica> Eratosthenes2[n_] := Module[{numbers = Range[3, n, 2], limit = (n - 1)/2},
Do[c = numbersi; If[c != 0, Do[numbersj = 0, {j,(c c - 1)/2,limit,c}]], {i,1,(Sqrt[n] - 1)/2}]; Prepend[Select[numbers, # > 1 &], 2]]
Eratosthenes2[100] </lang>
MATLAB
Not a true Sieve of Eratosthenes but rather a Trial Division Sieve
<lang MATLAB>function primeList = sieveOfEratosthenes(lastNumber)
list = (2:lastNumber); %Construct list of numbers
primeList = []; %Preallocate prime list
while( list(1)^2 <lastNumber )
primeList = [primeList list(1)]; %add prime to the prime list
list( mod(list,list(1))==0 ) = []; %filter out all multiples of the current prime
end
primeList = [primeList list]; %The rest of the numbers in the list are primes
end</lang>
Sample Output: <lang MATLAB>sieveOfEratosthenes(30)
ans =
2 3 5 7 11 13 17 19 23 29</lang>
Somewhat optimized true Sieve of Eratosthenes
<lang MATLAB> function P = erato(x) % Sieve of Eratosthenes: returns all primes between 2 and x
P = [0 2:x] ; % Create vector with all ints between 2 and x where
% position 1 is hard-coded as 0 since 1 is not a prime.
for (n=2:sqrt(x)) % All primes factors lie between 2 and sqrt(x).
if P(n) % If the current value is not 0 (i.e. a prime),
P((2*n):n:x) = 0 ; % then replace all further multiples of it with 0.
end
end % At this point P is a vector with only primes and zeroes.
P = P(P ~= 0) ; % Remove all zeroes from P, leaving only the primes.
return </lang>
The optimization lies in fewer steps in the for loop, use of MATLAB's built-in array operations and no modulo calculation.
Limitation: your machine has to be able to allocate enough memory for an array of length x.
Maxima
<lang 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
)
)$ </lang>
MAXScript
fn eratosthenes n =
(
multiples = #()
print 2
for i in 3 to n do
(
if (findItem multiples i) == 0 then
(
print i
for j in (i * i) to n by i do
(
append multiples j
)
)
)
)
eratosthenes 100
Modula-3
Regular version
<lang modula3> MODULE Prime EXPORTS Main;
IMPORT IO;
CONST LastNum = 1000;
VAR a: ARRAY [2..LastNum] OF BOOLEAN;
BEGIN
FOR i := FIRST(a) TO LAST(a) DO a[i] := TRUE; END;
FOR i := FIRST(a) TO LAST(a) DO
IF a[i] THEN
IO.PutInt(i);
IO.Put(" ");
FOR j := FIRST(a) TO LAST(a) DO
IF j MOD i = 0 THEN
a[j] := FALSE;
END;
END;
END;
END;
IO.Put("\n");
END Prime. </lang>
Advanced version
This version uses more "advanced" types. <lang modula3>(* From the CM3 examples folder (comments removed). *)
MODULE Sieve EXPORTS Main;
IMPORT IO;
TYPE
Number = [2..1000]; Set = SET OF Number;
VAR
prime: Set := Set {FIRST(Number) .. LAST(Number)};
BEGIN
FOR i := FIRST(Number) TO LAST(Number) DO
IF i IN prime THEN
IO.PutInt(i);
IO.Put(" ");
FOR j := i TO LAST(Number) BY i DO
prime := prime - Set{j};
END;
END;
END;
IO.Put("\n");
END Sieve.</lang>
MUMPS
<lang MUMPS>ERATO1(HI)
;performs the Sieve of Erotosethenes up to the number passed in. ;This version sets an array containing the primes SET HI=HI\1 KILL ERATO1 ;Don't make it new - we want it to remain after we quit the function NEW I,J,P FOR I=2:1:(HI**.5)\1 FOR J=I*I:I:HI SET P(J)=1 FOR I=2:1:HI S:'$DATA(P(I)) ERATO1(I)=I KILL I,J,P QUIT</lang>
Example:
USER>SET MAX=100,C=0 DO ERATO1^ROSETTA(MAX) USER>WRITE !,"PRIMES BETWEEN 1 AND ",MAX,! FOR SET I=$ORDER(ERATO1(I)) Q:+I<1 WRITE I,", " PRIMES BETWEEN 1 AND 100 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73,79, 83, 89, 97,
NetRexx
Version 1 (slow)
<lang Rexx>/* NetRexx */
options replace format comments java crossref savelog symbols binary
parse arg loWatermark hiWatermark . if loWatermark = | loWatermark = '.' then loWatermark = 1 if hiWatermark = | hiWatermark = '.' then hiWatermark = 200
do
if \loWatermark.datatype('w') | \hiWatermark.datatype('w') then -
signal NumberFormatException('arguments must be whole numbers')
if loWatermark > hiWatermark then -
signal IllegalArgumentException('the start value must be less than the end value')
seive = sieveOfEratosthenes(hiWatermark) primes = getPrimes(seive, loWatermark, hiWatermark).strip
say 'List of prime numbers from' loWatermark 'to' hiWatermark 'via a "Sieve of Eratosthenes" algorithm:'
say ' 'primes.changestr(' ', ',')
say ' Count of primes:' primes.words
catch ex = Exception
ex.printStackTrace
end
return
method sieveOfEratosthenes(hn = long) public static binary returns Rexx
sv = Rexx(isTrue) sv[1] = isFalse ix = long jx = long
loop ix = 2 while ix * ix <= hn
if sv[ix] then loop jx = ix * ix by ix while jx <= hn
sv[jx] = isFalse
end jx
end ix
return sv
method getPrimes(seive = Rexx, lo = long, hi = long) private constant binary returns Rexx
primes = Rexx() loop p_ = lo to hi if \seive[p_] then iterate p_ primes = primes p_ end p_
return primes
method isTrue public constant binary returns boolean
return 1 == 1
method isFalse public constant binary returns boolean
return \isTrue
</lang>
- Output
List of prime numbers from 1 to 200 via a "Sieve of Eratosthenes" algorithm: 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199 Count of primes: 46
Version 2 (significantly, i.e. 10 times faster)
<lang NetRexx>/* 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</lang>
Nial
primes is sublist [ each (2 = sum eachright (0 = mod) [pass,count]), pass ] rest count
Using it
|primes 10 =2 3 5 7
Nimrod
Based on one of Python solutions: <lang nimrod>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,"")</lang>
- 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
Niue
<lang Niue>[ dup 2 < ] '<2 ; [ 1 + 'count ; [ <2 [ , ] when ] count times ] 'fill-stack ;
0 'n ; 0 'v ;
[ .clr 0 'n ; 0 'v ; ] 'reset ; [ len 1 - n - at 'v ; ] 'set-base ; [ n 1 + 'n ; ] 'incr-n ; [ mod 0 = ] 'is-factor ; [ dup * ] 'sqr ;
[ set-base
v sqr 2 at > not [ [ dup v = not swap v is-factor and ] remove-if incr-n run ] when ] 'run ;
[ fill-stack run ] 'sieve ;
( tests )
10 sieve .s ( => 2 3 5 7 9 ) reset newline 30 sieve .s ( => 2 3 5 7 11 13 17 19 23 29 ) </lang>
Oberon-2
<lang oberon2>MODULE Primes;
IMPORT Out, Math;
CONST N = 1000;
VAR a: ARRAY N OF BOOLEAN;
i, j, m: INTEGER;
BEGIN
(* Set all elements of a to TRUE. *)
FOR i := 1 TO N - 1 DO
a[i] := TRUE;
END;
(* Compute square root of N and convert back to INTEGER. *) m := ENTIER(Math.Sqrt(N));
FOR i := 2 TO m DO
IF a[i] THEN
FOR j := 2 TO (N - 1) DIV i DO
a[i*j] := FALSE;
END;
END;
END;
(* Print all the elements of a that are TRUE. *)
FOR i := 2 TO N - 1 DO
IF a[i] THEN
Out.Int(i, 4);
END;
END;
Out.Ln;
END Primes.</lang>
OCaml
Imperative
<lang ocaml>let sieve n =
let is_prime = Array.create n true in
let limit = truncate(sqrt (float n)) in
for i = 2 to pred limit do
if is_prime.(i) then
let j = ref (i*i) in
try while true do
is_prime.(!j) <- false;
j := !j + i;
done with _ -> ()
done;
(is_prime)</lang>
<lang ocaml>let primes n =
let _, primes =
Array.fold_left (fun (i,acc) -> function
| true -> (i+1, i::acc)
| false -> (i+1, acc))
(0, [])
(sieve n)
in
List.tl(List.tl(List.rev primes))</lang>
in the top-level:
# primes 100 ;; - : int list = [2; 3; 5; 7; 11; 13; 17; 19; 23; 29; 31; 37; 41; 43; 47; 53; 59; 61; 67; 71; 73; 79; 83; 89; 97]
Functional
<lang ocaml>(* first define some iterators *)
- let fold_iter f init a b =
let rec aux acc i =
if i > b
then (acc)
else aux (f acc i) (succ i)
in
aux init a ;;
val fold_iter : ('a -> int -> 'a) -> 'a -> int -> int -> 'a = <fun>
- let fold_step f init a b step =
let rec aux acc i =
if i > b
then (acc)
else aux (f acc i) (i + step)
in
aux init a ;;
val fold_step : ('a -> int -> 'a) -> 'a -> int -> int -> int -> 'a = <fun>
(* remove a given value from a list *)
- let remove li v =
let rec aux acc = function
| hd::tl when hd = v -> (List.rev_append acc tl)
| hd::tl -> aux (hd::acc) tl
| [] -> li
in
aux [] li ;;
val remove : 'a list -> 'a -> 'a list = <fun>
(* the main function *)
- let primes n =
let li =
(* create a list [from 2; ... until n] *)
List.rev(fold_iter (fun acc i -> (i::acc)) [] 2 n)
in
let limit = truncate(sqrt(float n)) in
fold_iter (fun li i ->
if List.mem i li (* test if (i) is prime *)
then (fold_step remove li (i*i) n i)
else li)
li 2 (pred limit)
;;
val primes : int -> int list = <fun>
- primes 200 ;;
- : int list = [2; 3; 5; 7; 11; 13; 17; 19; 23; 29; 31; 37; 41; 43; 47; 53; 59; 61; 67; 71;
73; 79; 83; 89; 97; 101; 103; 107; 109; 113; 127; 131; 137; 139; 149; 151; 157; 163; 167; 173; 179; 181; 191; 193; 197; 199]</lang>
Another functional version
This uses zero to denote struck-out numbers. It is slightly inefficient as it strikes-out multiples above p rather than p2
<lang ocaml># let rec strike_nth k n l = match l with
| [] -> [] | h :: t -> if k = 0 then 0 :: strike_nth (n-1) n t else h :: strike_nth (k-1) n t;;
val strike_nth : int -> int -> int list -> int list = <fun>
- let primes n =
let limit = truncate(sqrt(float n)) in
let rec range a b = if a > b then [] else a :: range (a+1) b in
let rec sieve_primes l = match l with
| [] -> []
| 0 :: t -> sieve_primes t
| h :: t -> if h > limit then List.filter ((<) 0) l else
h :: sieve_primes (strike_nth (h-1) h t) in
sieve_primes (range 2 n) ;;
val primes : int -> int list = <fun>
- primes 200;;
- : int list = [2; 3; 5; 7; 11; 13; 17; 19; 23; 29; 31; 37; 41; 43; 47; 53; 59; 61; 67; 71;
73; 79; 83; 89; 97; 101; 103; 107; 109; 113; 127; 131; 137; 139; 149; 151; 157; 163; 167; 173; 179; 181; 191; 193; 197; 199]</lang>
Oz
<lang oz>declare
fun {Sieve N}
S = {Array.new 2 N true}
M = {Float.toInt {Sqrt {Int.toFloat N}}}
in
for I in 2..M do
if S.I then for J in I*I..N;I do S.J := false end end
end
S
end
fun {Primes N}
S = {Sieve N}
in
for I in 2..N collect:C do
if S.I then {C I} end
end end
in
{Show {Primes 30}}</lang>
PARI/GP
<lang parigp>Eratosthenes(lim)={
my(v=Vectorsmall(lim\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", ")))
};</lang>
An alternate version:
<lang parigp>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))) };</lang>
Pascal
Note: Some Pascal implementations put quite low limits on the size of a set (e.g. Turbo Pascal doesn't allow more than 256 members). To compile on such an implementation, reduce the constant PrimeLimit accordingly. <lang pascal> program primes(output)
const
PrimeLimit = 1000;
var
primes: set of 1 .. PrimeLimit; n, k: integer; needcomma: boolean;
begin
{ calculate the primes }
primes := [2 .. PrimeLimit];
for n := 1 to trunc(sqrt(PrimeLimit)) do
begin
if n in primes
then
begin
k := n*n;
while k < PrimeLimit do
begin
primes := primes - [k];
k := k + n
end
end
end;
{ output the primes }
needcomma := false;
for n := 1 to PrimeLimit do
if n in primes
then
begin
if needcomma
then
write(', ');
write(n);
needcomma := true
end
end. </lang>
Perl
For highest performance and ease, typically a module would be used, such as Math::Prime::Util, Math::Prime::FastSieve, or Math::Prime::XS.
Classic Sieve
<lang perl>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;
}</lang>
Odds only (faster)
<lang perl>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;
}</lang>
Odds only, using vectors for lower memory use
<lang perl>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;
}</lang>
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. <lang perl>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;
}</lang>
Experimental
These are examples of golfing or unusual styles.
Golfing a bit, at the expense of speed: <lang perl>sub sieve{ my (@s, $i); grep { not $s[ $i = $_ ] and do { $s[ $i += $_ ]++ while $i <= $_[0]; 1 } } 2..$_[0] }
print join ", " => sieve 100;</lang>
Or with bit strings (much slower than the vector version above): <lang perl>sub sieve{ my ($s, $i); grep { not vec $s, $i = $_, 1 and do { (vec $s, $i += $_, 1) = 1 while $i <= $_[0]; 1 } } 2..$_[0] }
print join ", " => sieve 100;</lang>
A short recursive version: <lang perl>sub erat {
my $p = shift; return $p, $p**2 > $_[$#_] ? @_ : erat(grep $_%$p, @_)
}
print join ', ' => erat 2..100000;</lang>
Regexp (purely an example -- the regex engine limits it to only 32769):<lang perl>sub sieve { my ($s, $p) = "." . ("x" x shift);
1 while ++$p and $s =~ /^(.{$p,}?)x/g and $p = length($1) and $s =~ s/(.{$p})./${1}./g and substr($s, $p, 1) = "x"; $s }
print sieve(1000);</lang>
Extensible sieves
Here are two incremental versions, which allows one to create a tied array of primes: <lang perl>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;</lang> 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. <lang perl>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;</lang>
Perl 6
<lang perl6>sub sieve( Int $limit ) {
my @is-prime = False, False, True xx $limit - 1;
gather for @is-prime.kv -> $number, $is-prime {
if $is-prime {
take $number;
@is-prime[$_] = False if $_ %% $number for $number**2 .. $limit;
}
}
}
(sieve 100).join(",").say;</lang>
A recursive version: <lang perl6>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;</lang>
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):
<lang perl6>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 ... *;</lang>
PHP
<lang php> function iprimes_upto($limit) {
for ($i = 2; $i < $limit; $i++)
{
$primes[$i] = true;
}
for ($n = 2; $n < $limit; $n++)
{
if ($primes[$n]) { for ($i = $n*$n; $i < $limit; $i += $n) { $primes[$i] = false; } }
} return $primes;
} </lang>
PicoLisp
<lang PicoLisp>(de sieve (N)
(let Sieve (range 1 N)
(set Sieve)
(for I (cdr Sieve)
(when I
(for (S (nth Sieve (* I I)) S (nth (cdr S) I))
(set S) ) ) )
(filter bool Sieve) ) )</lang>
Output:
: (sieve 100) -> (2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)
PL/I
<lang PL/I> eratos: proc options (main) reorder;
dcl i fixed bin (31); dcl j fixed bin (31); dcl n fixed bin (31); dcl sn fixed bin (31);
dcl hbound builtin; dcl sqrt builtin;
dcl sysin file; dcl sysprint file;
get list (n); sn = sqrt(n);
begin;
dcl primes(n) bit (1) aligned init ((*)((1)'1'b));
i = 2;
do while(i <= sn);
do j = i ** 2 by i to hbound(primes, 1);
/* Adding a test would just slow down processing! */
primes(j) = '0'b;
end;
do i = i + 1 to sn until(primes(i)); end; end;
do i = 2 to hbound(primes, 1);
if primes(i) then
put data(i);
end;
end; end eratos; </lang>
Pop11
define eratostenes(n);
lvars bits = inits(n), i, j;
for i from 2 to n do
if bits(i) = 0 then
printf('' >< i, '%s\n');
for j from 2*i by i to n do
1 -> bits(j);
endfor;
endif;
endfor;
enddefine;
PowerShell
Basic procedure
It outputs immediately so that the number can be used by the pipeline. <lang PowerShell>function Sieve ( [int] $num ) {
$isprime = @{}
2..$num | Where-Object {
$isprime[$_] -eq $null } | ForEach-Object {
$_
$isprime[$_] = $true
$i=$_*$_
for ( ; $i -le $num; $i += $_ )
{ $isprime[$i] = $false }
}
}</lang>
Processing
Calculate the primes up to 1000000 with Processing, including a visualisation of the process. As an additional visual effect, the layout of the pixel could be changed from the line-by-line layout to a spiral-like layout starting in the middle of the screen. <lang java>int maxx,maxy; int max; boolean[] sieve;
void plot(int pos, boolean active) {
set(pos%maxx,pos/maxx, active?#000000:#ffffff);
}
void setup() {
size(1000, 1000, P2D);
frameRate(2);
maxx=width;
maxy=height;
max=width*height;
sieve=new boolean[max+1];
sieve[1]=false;
plot(0,false);
plot(1,false);
for(int i=2;i<=max;i++) {
sieve[i]=true;
plot(i,true);
}
}
int i=2;
void draw() {
if(!sieve[i]) {
while(i*i<max && !sieve[i]) {
i++;
}
}
if(sieve[i]) {
print(i+" ");
for(int j=i*i;j<=max;j+=i) {
if(sieve[j]) {
sieve[j]=false;
plot(j,false);
}
}
}
if(i*i<max) {
i++;
} else {
noLoop();
println("finished");
}
}</lang>
Prolog
Using facts to record composite numbers
The first two solutions use Prolog "facts" to record the composite (i.e. already-visited) numbers.
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.
<lang Prolog>% %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 ).
</lang> The above has been tested with SWI-Prolog and gprolog.
<lang Prolog>% 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. </lang>
Optimized approach
<lang 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.</lang>
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)
<lang 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, _).</lang>
PureBasic
Basic procedure
<lang PureBasic>For n=2 To Sqr(lim)
If Nums(n)=0
m=n*n
While m<=lim
Nums(m)=1
m+n
Wend
EndIf
Next n</lang>
Working example
<lang PureBasic>Dim Nums.i(0) Define l, n, m, lim
If OpenConsole()
; Ask for the limit to search, get that input and allocate a Array
Print("Enter limit for this search: ")
lim=Val(Input())
ReDim Nums(lim)
; Use a basic Sieve of Eratosthenes
For n=2 To Sqr(lim)
If Nums(n)=#False
m=n*n
While m<=lim
Nums(m)=#True
m+n
Wend
EndIf
Next n
;Present the result to our user
PrintN(#CRLF$+"The Prims up to "+Str(lim)+" are;")
m=0: l=Log10(lim)+1
For n=2 To lim
If Nums(n)=#False
Print(RSet(Str(n),l)+" ")
m+1
If m>72/(l+1)
m=0: PrintN("")
EndIf
EndIf
Next
Print(#CRLF$+#CRLF$+"Press ENTER to exit"): Input()
CloseConsole()
EndIf</lang>
Output may look like;
Enter limit for this search: 750 The Prims up to 750 are; 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691 701 709 719 727 733 739 743 Press ENTER to exit
Python
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.
Using set lookup
The version below uses a set to store the multiples. set objects are much faster (usually O(log n)) than lists (O(n)) for checking if a given object is a member. Using the set.update method avoids explicit iteration in the interpreter, giving a further speed improvement.
<lang python>def eratosthenes2(n):
multiples = set()
for i in range(2, n+1):
if i not in multiples:
print(i)
multiples.update(range(i*i, n+1, i))
eratosthenes2(100)</lang>
Using array lookup
The version below uses array lookup to test for primality. The function primes_upto() is a straightforward implementation of Sieve of Eratosthenesalgorithm. It returns prime numbers less than or equal to limit. <lang python>def primes_upto(limit):
is_prime = [False] * 2 + [True] * (limit - 1)
for n in 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]</lang>
Using generator
The following code may be slightly slower than using the array/list as above, but uses no memory for output: <lang python>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: <lang python>>>> list(iprimes_upto(15)) [2, 3, 5, 7, 11, 13]</lang>
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: <lang python>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]</lang>
Note that "range" needs to be changed to "xrange" for maximum speed with Python 2.
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:
<lang python>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)</lang>
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.
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:
<lang python>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])</lang>
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).
Using numpy
Below code adapted from literateprograms.org using numpy <lang python>import numpy def primes_upto2(limit):
is_prime = numpy.ones(limit + 1, dtype=numpy.bool)
for n in xrange(2, int(limit**0.5 + 1.5)):
if is_prime[n]:
is_prime[n*n::n] = 0
return numpy.nonzero(is_prime)[0][2:]</lang>
Performance note: there is no point to add wheels here, due to execution of p[n*n::n] = 0 and nonzero() takes us almost all time.
Also see Prime numbers and Numpy – Python.
Using wheels with numpy
Version with wheel based optimization: <lang python>from numpy import array, bool_, multiply, nonzero, ones, put, resize
def makepattern(smallprimes):
pattern = ones(multiply.reduce(smallprimes), dtype=bool_)
pattern[0] = 0
for p in smallprimes:
pattern[p::p] = 0
return pattern
def primes_upto3(limit, smallprimes=(2,3,5,7,11)):
sp = array(smallprimes)
if limit <= sp.max(): return sp[sp <= limit]
#
isprime = resize(makepattern(sp), limit + 1)
isprime[:2] = 0; put(isprime, sp, 1)
#
for n in range(sp.max() + 2, int(limit**0.5 + 1.5), 2):
if isprime[n]:
isprime[n*n::n] = 0
return nonzero(isprime)[0]</lang>
Examples: <lang python>>>> primes_upto3(10**6, smallprimes=(2,3)) # Wall time: 0.17 array([ 2, 3, 5, ..., 999961, 999979, 999983]) >>> primes_upto3(10**7, smallprimes=(2,3)) # Wall time: 2.13 array([ 2, 3, 5, ..., 9999971, 9999973, 9999991]) >>> primes_upto3(15) array([ 2, 3, 5, 7, 11, 13]) >>> primes_upto3(10**7, smallprimes=primes_upto3(15)) # Wall time: 1.31 array([ 2, 3, 5, ..., 9999971, 9999973, 9999991]) >>> primes_upto2(10**7) # Wall time: 1.39 <-- version without wheels array([ 2, 3, 5, ..., 9999971, 9999973, 9999991]) >>> primes_upto3(10**7) # Wall time: 1.30 array([ 2, 3, 5, ..., 9999971, 9999973, 9999991])</lang> The above-mentioned examples demonstrate that the given wheel based optimization does not show significant performance gain.
Infinite generator
A generator that will generate primes indefinitely (perhaps until it runs out of memory). Used as a library here.
<lang python>import heapq
- generates all prime numbers
def sieve():
# priority queue of the sequences of non-primes
# the priority queue allows us to get the "next" non-prime quickly
nonprimes = []
i = 2
while True:
if nonprimes and i == nonprimes[0][0]: # non-prime
while nonprimes[0][0] == i:
# for each sequence that generates this number,
# have it go to the next number (simply add the prime)
# and re-position it in the priority queue
x = nonprimes[0]
x[0] += x[1]
heapq.heapreplace(nonprimes, x)
else: # prime
# insert a 2-element list into the priority queue:
# [current multiple, prime]
# the first element allows sorting by value of current multiple
# we start with i^2
heapq.heappush(nonprimes, [i*i, i])
yield i
i += 1</lang>
Example:
>>> foo = sieve() >>> for i in range(8): ... print(next(foo)) ... 2 3 5 7 11 13 17 19
Note that the above implementation has a very low limit in use to the prime value of 46349 for Python 2.X due to the normal short integer "roll-over" overflows to negative numbers for the prime squares. Although the infinite precision by default integers of Python 3 fix this, it only masks a very poor algorithm that adds prime multiple cull sequences to the Priority Queue prematurely, resulting in a much larger queue than necessary and therefor much longer (log n) re-insertion times.
The below codes gain efficiency by postponing the adding of cull sequences until necessary and as well use a built-in dictionary which is hash table based and therefore amortized to O(1) access time rather than the log n access time for the Priority Queue.
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):
<lang python>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()))</lang>
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.
<lang python>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</lang>
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: <lang python>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
</lang>
R
This code is vectorised, so it runs quickly for n < one million. The allocation of the primes vector means memory usage is too high for n much larger than that. <lang R> sieve <- function(n) {
n <- as.integer(n)
if(n > 1e6) stop("n too large")
primes <- rep(TRUE, n)
primes[1] <- FALSE
last.prime <- 2L
for(i in last.prime:floor(sqrt(n)))
{
primes[seq.int(2L*last.prime, n, last.prime)] <- FALSE
last.prime <- last.prime + min(which(primes[(last.prime+1):n]))
}
which(primes)
}
sieve(1000) </lang>
Racket
Imperative versions
<lang Racket>
- lang racket
- ugly imperative version
(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) </lang>
<lang Racket>
- lang racket
- a little nicer, but still imperative
(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) </lang>
Infinite list using laziness
<lang Racket>
- lang lazy
- infinite list using lazy racket (see the language above)
(define nats (cons 1 (map add1 nats))) (define (sift n l) (filter (λ(x) (not (zero? (modulo x n)))) l)) (define (sieve l) (cons (first l) (sieve (sift (first l) (rest l))))) (define primes (sieve (rest nats)))
(define (take-upto n l)
(if (<= (car l) n) (cons (car l) (take-upto n (cdr l))) '()))
(!! (take-upto 100 primes)) </lang>
Infinite list using threads and channels
<lang Racket>
- lang racket
- infinite list using threads and channels (similar to newsqueak)
(define-syntax (bg-loop stx)
(syntax-case stx ()
[(_ expr ...)
(with-syntax ([out! (datum->syntax stx 'out!)])
#'(let ([out (make-channel)])
(define (out! x) (channel-put out x))
(thread (λ() (let loop () expr ... (loop))))
out))]))
(define nats (bg-loop (for ([i (in-naturals 1)]) (out! i)))) (define (filter pred? c)
(bg-loop (define x (channel-get c))
(when (pred? x) (out! x))))
(define (sift n c)
(filter (λ(x) (not (zero? (modulo x n)))) c))
(define (sieve c)
(bg-loop (define x (channel-get c))
(out! x)
(set! c (sift x c))))
(define primes (begin (channel-get nats) (sieve nats)))
(define (take-upto n c)
(let loop () (let ([x (channel-get c)]) (if (<= x n) (cons x (loop)) '()))))
(take-upto 100 primes) </lang>
Infinite list using generators
<lang Racket>
- lang racket
- yet another variation of the same algorithm, this time using generators
(require racket/generator)
(define nats (generator () (for ([i (in-naturals 1)]) (yield i)))) (define (filter pred g)
(generator () (for ([i (in-producer g #f)] #:when (pred i)) (yield i))))
(define (sift n g) (filter (λ(x) (not (zero? (modulo x n)))) g)) (define (sieve g)
(generator () (let loop ([g g]) (let ([x (g)]) (yield x) (loop (sift x g))))))
(define primes (begin (nats) (sieve nats)))
(define (take-upto n p)
(let loop () (let ([x (p)]) (if (<= x n) (cons x (loop)) '()))))
(take-upto 100 primes) </lang>
REXX
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).
<lang REXX>/*REXX program generates primes via the sieve of Eratosthenes algorithm.*/
parse arg H .; if H== then H=200 /*let the highest # be specified.*/
w = length(H); prime=right('prime',20) /*W is used for formatting output*/
@.=1 /*assume all numbers are prime. */
- =0 /*number of primes found so far. */
do j=2 to H /*all integers up to H inclusive.*/
if @.j then do; #=#+1 /*Prime? Then bump prime counter*/
say prime right(#,w) " ───► " right(j,w) /*show.*/
do m=j*j to H by j; @.m=0; end /*odd multiples*/
end /*plain*/ /*[↑] strike odd multiples ¬prime*/
end /*j*/
/*stick a fork in it, we're done.*/</lang>
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
wheel version with optional prime list suppression
This version skips striking the even numbers (as being not prime).
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.
<lang rexx>/*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. */
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 sieve striking. */
@.=1 /*assume all numbers are prime. */
- =1 /*number of primes found so far. */
do j=3 by 2 to H /*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=0; end /*odd multiples.*/
end /*plain*/ /*[↑] strike odd multiples ¬prime*/
end /*j*/
say; say right(#,w+20+2) 'prime's(#) "found." exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────S subroutine────────────────────────*/ s: if arg(1)==1 then return arg(3); return word(arg(2) 's',1) /*plural*/ /*──────────────────────────────────TH subroutine───────────────────────*/ th: procedure;parse arg x;x=abs(x);return word('th st nd rd',1+x//10*(x//100%10\==1)*(x//10<4))</lang> 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
13 +++ @.m=0 Error 5 running "C:\sieve_of_Eratosthenes.rex", line 13: System resources exhausted
The above 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 delayed by implementing further optimization to the wheel version.
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.
<lang rexx>/*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 sieve striking. */
@.=1 /*assume all numbers are prime. */
- =1 /*number of primes found so far. */
do j=3 by 2 to H /*odd integers up to H inclusive.*/
if @.j then do; #=#+1 /*Prime? Then bump prime counter*/
say prime right(#,w) " ───► " right(j,w) /*show.*/
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=0; end /*odd multiples.*/
end /*plain*/ /*[↑] strike odd multiples ¬prime*/
end /*j*/
/*stick a fork in it, we're done.*/</lang>
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.
Wheel Version restructured
<lang rexx>/*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</lang>
output is identical to the above versions.
Program execution time is slightly worse than that of the wheel version above
(tested with highest=2000000 and suppressed output: 2.4 seconds vs. 2.3 seconds).
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.
<lang ruby>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)</lang>
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
With a wheel
eratosthenes2 adds more optimizations, but the code is longer.
- The array
numsonly tracks odd numbers (skips multiples of 2). - The array
numsholds booleans instead of integers, and every multiple of 3 beginsfalse. - The outer loop skips multiples of 2 and 3.
- Both inner loops skip multiples of 2 and 3.
<lang ruby>def eratosthenes2(n)
# For odd i, if i is prime, nums[i >> 1] is true. # Set false for all multiples of 3. nums = [true, false, true] * ((n + 5) / 6) nums[0] = false # 1 is not prime. nums[1] = true # 3 is prime.
# Outer loop and both inner loops are skipping multiples of 2 and 3.
# Outer loop checks i * i > n, same as i > Math.sqrt(n).
i = 5
until (m = i * i) > 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)</lang>
This simple benchmark compares eratosthenes with eratosthenes2.
<lang ruby>require 'benchmark' Benchmark.bmbm {|x|
x.report("eratosthenes") { eratosthenes(1_000_000) }
x.report("eratosthenes2") { eratosthenes2(1_000_000) }
}</lang>
eratosthenes2 runs about 4 times faster than eratosthenes.
With the standard library
MRI 1.9.x implements the sieve of Eratosthenes at file prime.rb, class EratosthensesSeive (around line 421). This implementation optimizes for space, by packing the booleans into 16-bit integers. It also hardcodes all primes less than 256.
<lang ruby>require 'prime' p Prime::EratosthenesGenerator.new.take_while {|i| i <= 100}</lang>
Run BASIC
<lang runbasic>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</lang>
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,
Rust
Sieve of Eratosthenes - No optimization
<lang rust>use std::vec; use std::num; use std::iter;
fn int_sqrt(n: uint) -> uint {
num::sqrt(n as f64) as uint
}
fn simple_sieve(limit: uint) -> ~[uint] {
if limit < 2 {
return ~[];
}
let mut primes = 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[multiple] = false
}
}
}
iter::range_inclusive(2, limit).filter(|&n| primes[n]).collect()
}
fn main() {
println!("{:?}", simple_sieve(100))
}</lang>
- Output:
~[2u, 3u, 5u, 7u, 11u, 13u, 17u, 19u, 23u, 29u, 31u, 37u, 41u, 43u, 47u, 53u, 59u, 61u, 67u, 71u, 73u, 79u, 83u, 89u, 97u]
Scala
Genuine Eratosthenes sieve
<lang Scala>def sieveOfEratosthenes(nTo: Int) = {
val primes = collection.mutable.BitSet.empty.par ++ (2 to nTo)
for {
candidate <- 2 to Math.sqrt(nTo).toInt
if primes contains candidate
} primes --= candidate * candidate to nTo by candidate
primes
} // BitSet toList is shuffled when using the ParSet version. println(sieveOfEratosthenes(100).toList.sorted) println(sieveOfEratosthenes(100000000).size)</lang>
- Output:
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) 5761455
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.
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:
<lang Scala>import scala.annotation.tailrec
def sieveOfEratosthenes(nTo: Int) = {
val nonprimes: Array[Int] = new Array((nTo >>> 5) + 1)
def cullp(p: Int) = {
@tailrec def cull(c: Int): Unit = {
if (c <= nTo) nonprimes(c >>> 5) |= 1 << (c & 31)
cull(c + p)
}
cull(p * p)
}
def nextprime(i: Int): Int = {
if (i > nTo) i else if ((nonprimes(i >>> 5) & (1 << (i & 31))) != 0) nextprime(i + 1) else i
}
for {
p <- 2 to Math.sqrt(nTo).toInt
if (nonprimes(p >>> 5).&(1 << (p & 31))) == 0
} cullp(p)
Iterator.iterate(2)(c => nextprime(c + 1)).takeWhile(_ <= nTo)
}</lang>
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:
<lang Scala>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
}</lang>
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.
Odds-Only "infinite" generator sieve using Streams and Co-Inductive Streams
Using Streams, the "unfaithful sieve", i.e. sub-optimal trial division sieve. Not a sieve of Eratosthenes. <lang scala>def sieve(nums: Stream[Int]): Stream[Int] =
Stream.cons(nums.head, sieve((nums.tail).filter(_ % nums.head != 0))) val primes = 2 #:: sieve(Stream.from(3, 2))
println(primes take 10 toList) // //List(2, 3, 5, 7, 11, 13, 17, 19, 23, 29) println(primes takeWhile (_ < 30) toList) //List(2, 3, 5, 7, 11, 13, 17, 19, 23, 29)</lang>
- Output:
Both println statements give the same results
List(2, 3, 5, 7, 11, 13, 17, 19, 23, 29)
The above code is extremely inefficient for larger ranges, both because it tests for primality using computationally expensive divide (modulo) operations and because it sets up deferred tests for division by all of the primes up to each prime candidate, meaning that it has approximately a square law computational complexity with range.
The following code uses delayed recursion via Streams to implement the Richard Bird algorithm mentioned in the last part (the Epilogue) of the above linked article, 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 above "unfaithful" sieve.
<lang scala> def primesBird: Stream[Int] = {
def merge(xs: Stream[Int], ys: Stream[Int]): Stream[Int] = {
val (x, y) = (xs.head, ys.head)
if (y > x) Stream.cons(x, merge(xs.tail, ys))
else if (x > y) Stream.cons(y, merge(xs, ys.tail))
else Stream.cons(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]] =
Stream.cons(primeMltpls(ps.head), allMltpls(ps.tail))
def join(ams: Stream[Stream[Int]]): Stream[Int] =
Stream.cons(ams.head.head, merge(ams.head.tail, join(ams.tail)))
def oddPrimes(): Stream[Int] = { //implements "minua"
def oddPrms(n: Int, composites: Stream[Int]): Stream[Int] =
if (n >= composites.head) oddPrms(n + 2, composites.tail)
else Stream.cons(n, oddPrms(n + 2, composites))
//following uses a new recursive source of odd base primes
Stream.cons(3, oddPrms(5, join(allMltpls(oddPrimes()))))
}
Stream.cons(2, oddPrimes())
}</lang>
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: <lang scala> 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) }</lang>
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.
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:
<lang scala> 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
}</lang>
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.
Scheme
<lang scheme>; 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</lang>
<lang scheme>; 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)</lang> Output: <lang>(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)</lang> Optimised using a pre-computed wheel based on 2 (i.e. odds only): <lang scheme>(define (primes-wheel-2 limit)
(let ((stop (sqrt limit)))
(define (sieve lst)
(let ((p (car lst)))
(if (> p stop)
lst
(cons p (sieve (strike (cdr lst) (* p p) (* 2 p)))))))
(cons 2 (sieve (iota 3 limit 2)))))
(display (primes-wheel-2 100)) (newline)</lang> Output: <lang>(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)</lang>
Vector-based (faster), works with RRS: <lang scheme>; initialize v to vector of sequential integers (define (initialize! v)
(define (iter v n) (if (>= n (vector-length v))
(values)
(begin (vector-set! v n n) (iter v (+ n 1)))))
(iter v 0))
- set every nth element of vector v to 0,
- starting with element m
(define (strike! v m n)
(cond ((>= m (vector-length v)) (values))
(else (begin
(vector-set! v m 0)
(strike! v (+ m n) n)))))
- lowest non-zero index of vector v >= n
(define (nextprime v n)
(if (zero? (vector-ref v n))
(nextprime v (+ n 1))
(vector-ref v n)))
- remove elements satisfying pred? from list lst
(define (remove pred? lst)
(cond ((null? lst) '()) ((pred? (car lst))(remove pred? (cdr lst))) (else (cons (car lst) (remove pred? (cdr lst))))))
- the sieve itself
(define (sieve n)
(define stop (sqrt n))
(define (iter v p)
(cond
((> p stop) v)
(else
(begin
(strike! v (* p p) p)
(iter v (nextprime v (+ p 1)))))))
(let ((v (make-vector (+ n 1))))
(initialize! v)
(vector-set! v 1 0) ; 1 is not a prime
(remove zero? (vector->list (iter v 2)))))</lang>
Seed7
The program below computes the number of primes between 1 and 10000000: <lang seed7>$ 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;</lang>
Original source: [1]
SNOBOL4
Using strings instead of arrays, and the square/sqrt optimizations.
<lang SNOBOL4> define('sieve(n)i,j,k,p,str,res') :(sieve_end) sieve i = lt(i,n - 1) i + 1 :f(sv1)
str = str (i + 1) ' ' :(sieve)
sv1 str break(' ') . j span(' ') = :f(return)
sieve = sieve j ' '
sieve = gt(j ^ 2,n) sieve str :s(return) ;* Opt1
res =
str (arb ' ') @p ((j ^ 2) ' ') ;* Opt2
str len(p) . res = ;* Opt2
sv2 str break(' ') . k span(' ') = :f(sv3)
res = ne(remdr(k,j),0) res k ' ' :(sv2)
sv3 str = res :(sv1) sieve_end
- # Test and display
output = sieve(100)
end</lang>
Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
Tcl
<lang tcl>package require Tcl 8.5
proc sieve n {
if {$n < 2} {return {}}
# create a container to hold the sequence of numbers.
# use a dictionary for its speedy access (like an associative array)
# and for its insertion order preservation (like a list)
set nums [dict create]
for {set i 2} {$i <= $n} {incr i} {
# the actual value is never used
dict set nums $i ""
}
set primes [list]
while {[set nextPrime [lindex [dict keys $nums] 0]] <= sqrt($n)} {
dict unset nums $nextPrime
for {set i [expr {$nextPrime ** 2}]} {$i <= $n} {incr i $nextPrime} {
dict unset nums $i
}
lappend primes $nextPrime
}
return [concat $primes [dict keys $nums]]
}
puts [sieve 100] ;# 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97</lang>
UNIX Shell
With array
<lang bash>#!/usr/bin/zsh
function primes() { typeset -a a typeset i j
a[1]="" for (( i = 2; i <= $1; i++ )); do a[$i]=$i done
for (( i = 2; i * i <= $1; i++ )); do if [[ ! -z $a[$i] ]]; then for (( j = i * i; j <= $1; j += i )); do a[$j]="" done fi done print $a }
primes 1000</lang>
<lang bash>function primes { typeset a i=2 j m=$1 # No for (( ... )) loop in pdksh. Use while loop. while (( i <= m )); do a[$i]=$i (( i++ )) done
i=2 while (( j = i * i, j <= m )); do if [[ -n ${a[$i]} ]]; then while (( j <= m )); do unset a[$j] (( j += i )) done fi (( i++ )) done # No print command in bash. Use echo command. echo ${a[*]} }
primes 1000</lang>
Both scripts output a single long line.
2 3 5 7 11 13 17 19 23 ... 971 977 983 991 997
Using variables as fake array
Bourne Shell and Almquist Shell have no arrays. This script works with bash or dash (standard shell in Ubuntu), but uses no specifics of the shells, so it works with plain Bourne Shell as well.
<lang bash>#! /bin/sh
LIMIT=1000
- As a workaround for missing arrays, we use variables p2, p3, ...,
- p$LIMIT, to represent the primes. Values are true or false.
- eval p$i=true # Set value.
- eval \$p$i # Run command: true or false.
- A previous version of this script created a temporary directory and
- used files named 2, 3, ..., $LIMIT to represent the primes. We now use
- variables so that a killed script does not leave extra files. About
- performance, variables are about as slow as files.
i=2 while [ $i -le $LIMIT ] do
eval p$i=true # was touch $i i=`expr $i + 1`
done
i=2 while
j=`expr $i '*' $i` [ $j -le $LIMIT ]
do
if eval \$p$i # was if [ -f $i ]
then
while [ $j -le $LIMIT ]
do
eval p$j=false # was rm -f $j
j=`expr $j + $i`
done
fi
i=`expr $i + 1`
done
- was echo `ls|sort -n`
echo `i=2
while [ $i -le $LIMIT ]; do
eval \\$p$i && echo $i
i=\`expr $i + 1\`
done`</lang>
With piping
This version works by piping 1s and 0s through sed. The string of 1s and 0s represents the array of primes.
<lang bash># Fill $1 characters with $2. fill () { # This pipeline would begin # head -c $1 /dev/zero | ... # but some systems have no head -c. Use dd. dd if=/dev/zero bs=$1 count=1 2>/dev/null | tr '\0' $2 }
filter () { # Use sed to put an 'x' after each multiple of $1, remove # first 'x', and mark non-primes with '0'. sed -e s/$2/\&x/g -e s/x// -e s/.x/0/g | { if expr $1 '*' $1 '<' $3 > /dev/null; then filter `expr $1 + 1` .$2 $3 else cat fi } }
- Generate a sequence of 1s and 0s indicating primality.
oz () { fill $1 1 | sed s/1/0/ | filter 2 .. $1 }
- Echo prime numbers from 2 to $1.
prime () { # Escape backslash inside backquotes. sed sees one backslash. echo `oz $1 | sed 's/./&\\ /g' | grep -n 1 | sed s/:1//` }
prime 1000</lang>
C Shell
<lang csh># Sieve of Eratosthenes: Echoes all prime numbers through $limit. @ limit = 80
if ( ( $limit * $limit ) / $limit != $limit ) then echo limit is too large, would cause integer overflow. exit 1 endif
- Use $prime[2], $prime[3], ..., $prime[$limit] as array of booleans.
- Initialize values to 1 => yes it is prime.
set prime=( `repeat $limit echo 1` )
- Find and echo prime numbers.
@ i = 2 while ( $i <= $limit ) if ( $prime[$i] ) then echo $i
# For each multiple of i, set 0 => no it is not prime. # Optimization: start at i squared. @ m = $i * $i while ( $m <= $limit ) set prime[$m] = 0 @ m += $i end endif @ i += 1 end</lang>
Ursala
with no optimizations <lang Ursala>#import nat
sieve = ~<{0,1}&& iota; @NttPX ~&r->lx ^/~&rhPlC remainder@rlX~|@r</lang> test program: <lang Ursala>#cast %nL
example = sieve 50</lang> output:
<2,3,5,7,11,13,17,19,23,29,31,37,41,43,47>
Vala
Without any optimizations: <lang vala> using Gee;
ArrayList<int> primes(int limit){ var sieve = new ArrayList<bool>(); var prime_list = new ArrayList<int>();
for(int i = 0; i <= limit; i++){ sieve.add(true); }
sieve[0] = false; sieve[1] = false;
for (int i = 2; i <= limit/2; i++){ if (sieve[i] != false){ for (int j = 2; i*j <= limit; j++){ sieve[i*j] = false; } } }
for (int i = 0; i < sieve.size; i++){ if (sieve[i] != false){ prime_list.add(i); } }
return prime_list; } // end primes
public static void main(){ var prime_list = primes(50);
foreach(var prime in prime_list) stdout.printf("%s ", prime.to_string());
stdout.printf("\n"); } </lang>
Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47
Vedit macro language
This implementation uses an edit buffer as an array for flags. After the macro has been run, you can see how the primes are located in the array. Primes are marked with 'P' and non-primes with '-'. The first character position represents number 0.
#10 = Get_Num("Enter number to search to: ", STATLINE)
Buf_Switch(Buf_Free) // Use edit buffer as flags array
Ins_Text("--") // 0 and 1 are not primes
Ins_Char('P', COUNT, #10-1) // init rest of the flags to "prime"
for (#1 = 2; #1*#1 < #10; #1++) {
Goto_Pos(#1)
if (Cur_Char=='P') { // this is a prime number
for (#2 = #1*#1; #2 <= #10; #2 += #1) {
Goto_Pos(#2)
Ins_Char('-', OVERWRITE)
}
}
}
Sample output showing numbers in range 0 to 599.
--PP-P-P---P-P---P-P---P-----P-P-----P---P-P---P-----P-----P -P-----P---P-P-----P---P-----P-------P---P-P---P-P---P------ -------P---P-----P-P---------P-P-----P-----P---P-----P-----P -P---------P-P---P-P-----------P-----------P---P-P---P-----P -P---------P-----P-----P-----P-P-----P---P-P---------P------ -------P---P-P---P-------------P-----P---------P-P---P-----P -------P-----P-----P---P-----P-------P---P-------P---------P -P---------P-P-----P---P-----P-------P---P-P---P-----------P -------P---P-------P---P-----P-----------P-P---------------- -P-----P---------P-----P-----P-P-----P---------P-----P-----P
Visual Basic .NET
<lang vbnet>Dim n As Integer, k As Integer, limit As Integer Console.WriteLine("Enter number to search to: ") limit = Console.ReadLine Dim flags(limit) As Integer For n = 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</lang>
Vorpal
<lang vorpal> self.print_primes = method(m){
p = new() p.fill(0, m, 1, true)
count = 0
i = 2
while(i < m){
if(p[i] == true){
p.fill(i+i, m, i, false)
count = count + 1
}
i = i + 1
}
('primes: ' + count + ' in ' + m).print()
for(i = 2, i < m, i = i + 1){
if(p[i] == true){
( + i + ', ').put()
}
}
.print()
}
self.print_primes(100) </lang> Result:
primes: 25 in 100 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,
XPL0
<lang 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
];
];
]</lang>
Example output:
20 2 3 5 7 11 13 17 19
zkl
<lang zkl>fcn sieve(limit){
if (limit<2) return(T);
composite:=(0).pump(limit+1,Data,T(Void,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)) },
False)); // turn filter into a no-result loop
(2).filter(limit-1,composite.__sGet); // bytes still 1 are prime
}
sieve(53).println();</lang>
- Output:
L(2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53)
- Programming Tasks
- Prime Numbers
- Clarified and Needing Review
- 68000 Assembly
- ACL2
- Ada
- ALGOL 68
- AutoHotkey
- AutoIt
- AWK
- BASIC
- Applesoft BASIC
- Locomotive Basic
- ZX Spectrum Basic
- BBC BASIC
- Befunge
- Bracmat
- C
- C++
- C sharp
- Chapel
- Clojure
- CMake
- COBOL
- Common Lisp
- D
- Dart
- Delphi
- DWScript
- Dylan
- E
- EC
- EC examples needing attention
- Examples needing attention
- Eiffel
- Emacs Lisp
- Erlang
- Euphoria
- F Sharp
- Forth
- Fortran
- GAP
- GLBasic
- Go
- Groovy
- GW-BASIC
- Haskell
- HicEst
- Icon
- Unicon
- J
- Java
- JavaScript
- Liberty BASIC
- Logo
- Lua
- Lucid
- M4
- Mathematica
- MATLAB
- Maxima
- MAXScript
- Modula-3
- MUMPS
- NetRexx
- Nial
- Nial examples needing attention
- Nimrod
- Niue
- Niue examples needing attention
- Oberon-2
- OCaml
- Oz
- PARI/GP
- Pascal
- Perl
- Perl 6
- PHP
- PicoLisp
- PL/I
- Pop11
- PowerShell
- Processing
- Prolog
- PureBasic
- Python
- Numpy
- R
- Racket
- REXX
- Ruby
- Run BASIC
- Rust
- Scala
- Scheme
- Seed7
- SNOBOL4
- Tcl
- UNIX Shell
- C Shell
- Ursala
- Ursala examples needing attention
- Vala
- Gee
- Vedit macro language
- Visual Basic .NET
- Vorpal
- XPL0
- Zkl