Sieve of Eratosthenes
From Rosetta Code
Sieve of Eratosthenes is a programming task. Visitors like you are encouraged to solve it according to the task description, using any language they may happen to 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. assume you need only to cross out odd numbers (wheel based on 2), numbers equal to 1 or 5 modulo 6 (wheel based on 2 and 3), or similar wheels based on low primes.
If there's an easy way to add such a wheel based optimization, implement this as an alternative version.
[edit] Ada
This solution uses an array of bits representing boolean values to identify the prime numbers. The limit on the outer loop is set to the square root of the number received from the command line. This limit minimizes the number of unneeded iterations in the algorithm.
The lowest index of the array of boolean is set to 2 so that all calculations begin at 2.
with Ada.Text_Io; use Ada.Text_Io;
with Ada.Command_Line; use Ada.Command_Line;
with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;
procedure Sieve is
package Nat_Io is new Ada.Text_Io.Integer_Io(Natural);
use Nat_Io;
type Bool_Array is array(Positive range <>) of Boolean;
pragma pack (Bool_Array);
procedure Get_Sieve (A : out Bool_Array) is
Test_Num : Positive;
Limit : Positive := Positive(sqrt(float(A'Last)));
begin
A := (Others => True); -- set all values to true
for Num in A'First..Limit loop
if A(Num) then
Test_Num := Num * Num;
while Test_Num <= A'Last loop
A(Test_Num) := False;
Test_Num := Test_Num + Num;
end loop;
end if;
end loop;
end Get_Sieve;
pragma Inline(Get_Sieve);
N : Positive := 2;
begin
if Argument_Count > 0 then
N := Positive'Value(Argument(1));
end if;
declare
Vals : Bool_Array(2..N);
begin
Get_Sieve(Vals);
for I in Vals'range loop
if Vals(I) then
Put_Line(Positive'Image(I));
end if;
end loop;
end;
end Sieve;
This version uses tasking to try to improve the performance on multiprocessor systems which are becoming common. It is still an implementation of the Sieve of Eratosthenes: the list of numbers is in the main task, and the list of primes is in the list of Siever tasks.
with Ada.Command_Line;
with Ada.Text_IO;
procedure Sieve is
task type Siever is
entry Check (Value : in Positive);
entry Print;
end Siever;
subtype Siever_Task is Siever;
type Sieve_Ptr is access Siever;
task body Siever is
Number : Positive;
Candidate : Positive;
Next : Sieve_Ptr;
begin -- Siever
accept Check (Value : in Positive) do
Number := Value;
end Check;
All_Values : loop
select
accept Check (Value : in Positive) do
Candidate := Value;
end Check;
if Candidate rem Number /= 0 then
if Next = null then
Next := new Siever_Task;
end if;
Next.Check (Value => Candidate);
end if;
or
accept Print do
Ada.Text_IO.Put (Item => Integer'Image (Number) );
if Next /= null then
Next.Print;
end if;
end Print;
exit All_Values;
end select;
end loop All_Values;
end Siever;
Head : Sieve_Ptr := new Siever;
Max : Positive := 20;
begin -- Sieve
if Ada.Command_Line.Argument_Count > 0 then
Max := Integer'Value (Ada.Command_Line.Argument (1) );
end if;
All_Values : for I in 2 .. Max loop
Head.Check (Value => I);
end loop All_Values;
Head.Print;
end Sieve;
This version is the same as the previous version, but with a wheel of 2.
with Ada.Command_Line;
with Ada.Text_IO;
procedure Sieve_With_Wheel is
task type Siever is
entry Check (Value : in Positive);
entry Print;
end Siever;
subtype Siever_Task is Siever;
type Sieve_Ptr is access Siever;
task body Siever is
Number : Positive;
Candidate : Positive;
Next : Sieve_Ptr;
begin -- Siever
accept Check (Value : in Positive) do
Number := Value;
end Check;
All_Values : loop
select
accept Check (Value : in Positive) do
Candidate := Value;
end Check;
if Candidate rem Number /= 0 then
if Next = null then
Next := new Siever_Task;
end if;
Next.Check (Value => Candidate);
end if;
or
accept Print do
Ada.Text_IO.Put (Item => Integer'Image (Number) );
if Next /= null then
Next.Print;
end if;
end Print;
exit All_Values;
end select;
end loop All_Values;
end Siever;
Head : Sieve_Ptr := new Siever;
Max : Positive := 20;
Value : Positive := 3;
begin -- Sieve_With_Wheel
if Ada.Command_Line.Argument_Count > 0 then
Max := Integer'Value (Ada.Command_Line.Argument (1) );
end if;
All_Values : loop
exit All_Values when Value > Max;
Head.Check (Value => Value);
Value := Value + 2;
end loop All_Values;
Ada.Text_IO.Put (Item => " 2");
Head.Print;
end Sieve_With_Wheel;
[edit] ALGOL 68
BOOL prime = TRUE, non prime = FALSE;
PROC eratosthenes = (INT n)[]BOOL:
(
[n]BOOL sieve;
FOR i TO UPB sieve DO sieve[i] := prime OD;
INT m = ENTIER sqrt(n);
sieve[1] := non prime;
FOR i FROM 2 TO m DO
IF sieve[i] = prime THEN
FOR j FROM i*i BY i TO n DO
sieve[j] := non prime
OD
FI
OD;
sieve
);
print((eratosthenes(80),new line))
Output:
FTTFTFTFFFTFTFFFTFTFFFTFFFFFTFTFFFFFTFFFTFTFFFTFFFFFTFFFFFTFTFFFFFTFFFTFTFFFFFTF
[edit] AutoHotkey
Search autohotkey.com: of Eratosthenes
Source: AutoHotkey forum by Laszlo
MsgBox % "12345678901234567890`n" Sieve(20)
Sieve(n) { ; Sieve of Eratosthenes => string of 0|1 chars, 1 at position k: k is prime
Static zero := 48, one := 49 ; Asc("0"), Asc("1")
VarSetCapacity(S,n,one)
NumPut(zero,S,0,"char")
i := 2
Loop % sqrt(n)-1 {
If (NumGet(S,i-1,"char") = one)
Loop % n//i
If (A_Index > 1)
NumPut(zero,S,A_Index*i-1,"char")
i += 1+(i>2)
}
Return S
}
[edit] 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
[edit] BASIC
Works with: FreeBASIC Works with: RapidQ
DIM n AS Integer, k AS Integer, limit AS Integer
INPUT "Enter number to search to: "; limit
DIM flags(limit) AS Integer
FOR n = 2 TO SQR(limit)
IF flags(n) = 0 THEN
FOR k = n*n TO limit STEP n
flags(k) = 1
NEXT k
END IF
NEXT n
' Display the primes
FOR n = 2 TO limit
IF flags(n) = 0 THEN PRINT n; ", ";
NEXT n
[edit] Befunge
2>:3g" "-!v\ g30 < |!`"O":+1_:.:03p>03g+:"O"`| @ ^ p3\" ":< 2 234567890123456789012345678901234567890123456789012345678901234567890123456789
[edit] C
Plain sieve, without any optimizations:
#include <stdlib.h>
#include <math.h>
char* eratosthenes (int n)
{
char* sieve;
int i, j, m;
/* calloc initializes to zero */
sieve = calloc (n+1, sizeof(char));
m = (int) sqrt ((double) n);
sieve[0] = 1;
sieve[1] = 1;
for (i = 2; i <= m; i++) {
if (sieve[i] == 0) {
for (j = i*i; j <= n; j += i) {
sieve[j] = 1;
}
}
}
return sieve;
}
Possible optimizations include sieving only odd numbers (or more complex wheels), packing the sieve into bits to improve locality (and allow larger sieves), etc.
[edit] C++
// yield all prime numbers less then limit.
template<class UnaryFunction>
void primesupto(int limit, UnaryFunction yield)
{
std::vector<bool> is_prime(limit, true);
const int sqrt_limit = static_cast<int>(std::sqrt(limit));
for (int n = 2; n <= sqrt_limit; ++n)
if (is_prime[n]) {
yield(n);
for (unsigned k = n*n, ulim = static_cast<unsigned>(limit); k < ulim; k += n)
//NOTE: "unsigned" is used to avoid an overflow in `k+=n` for `limit` near INT_MAX
is_prime[k] = false;
}
for (int n = sqrt_limit + 1; n < limit; ++n)
if (is_prime[n])
yield(n);
}
Full program:
Works with: Boost
/**
$ g++ -I/path/to/boost sieve.cpp -o sieve && sieve 10000000
*/
#include <inttypes.h> // uintmax_t
#include <limits>
#include <cmath>
#include <iostream>
#include <sstream>
#include <vector>
#include <boost/lambda/lambda.hpp>
int main(int argc, char *argv[])
{
using namespace std;
using namespace boost::lambda;
int limit = 10000;
if (argc == 2) {
stringstream ss(argv[--argc]);
ss >> limit;
if (limit < 1 or ss.fail()) {
cerr << "USAGE:\n sieve LIMIT\n\nwhere LIMIT in the range [1, "
<< numeric_limits<int>::max() << ")" << endl;
return 2;
}
}
// print primes less then 100
primesupto(100, cout << _1 << " ");
cout << endl;
// find number of primes less then limit and their sum
int count = 0;
uintmax_t sum = 0;
primesupto(limit, (var(sum) += _1, var(count) += 1));
cout << "limit sum pi(n)\n"
<< limit << " " << sum << " " << count << endl;
}
[edit] Common Lisp
(defun sieve-of-eratosthenes (maximum)
(let ((composites (make-array (1+ maximum) :element-type 'bit
:initial-element 0)))
(loop for candidate from 2 to maximum
when (zerop (bit composites candidate))
collect candidate
and do (loop for composite from (expt candidate 2) to maximum by candidate
do (setf (bit composites composite) 1)))))
[edit] D
D V.2 code:
import std.stdio: writeln;
import std.conv: to;
int[] sieve(int m) {
bool[] isPrime = new bool[m];
isPrime[] = true;
for (int i = 2; i < isPrime.length; i++)
if (isPrime[i])
for (int k = i * 2; k < isPrime.length; k += i)
isPrime[k] = false;
int[] result;
foreach (i, p; isPrime[2 .. $])
if (p)
result ~= i + 2;
return result;
}
void main(string[] args) {
int n = args.length > 1 ? to!int(args[1]) : 100;
writeln("Primes up to ", n, ": ", sieve(n));
}
[edit] E
E's standard library doesn't have a step-by-N numeric range, so we'll define one, implementing the standard iteration protocol.
def rangeFromBelowBy(start, limit, step) {
return def stepper {
to iterate(f) {
var i := start
while (i < limit) {
f(null, i)
i += step
}
}
}
}
The sieve itself:
def eratosthenes(limit :(int > 2), output) {
def composite := [].asSet().diverge()
for i ? (!composite.contains(i)) in 2..!limit {
output(i)
composite.addAll(rangeFromBelowBy(i ** 2, limit, i))
}
}
Example usage:
? eratosthenes(12, println) # stdout: 2 # 3 # 5 # 7 # 11
[edit] eC
Note: this is not a Sieve of Eratosthenes; it is just trial division.
public class FindPrime
{
Array<int> primeList { [ 2 ], minAllocSize = 64 };
int index;
index = 3;
bool HasPrimeFactor(int x)
{
int max = (int)floor(sqrt((double)x));
for(i : primeList)
{
if(i > max) break;
if(x % i == 0) return true;
}
return false;
}
public int GetPrime(int x)
{
if(x > primeList.count - 1)
{
for (; primeList.count != x; index += 2)
if(!HasPrimeFactor(index))
{
if(primeList.count >= primeList.minAllocSize) primeList.minAllocSize *= 2;
primeList.Add(index);
}
}
return primeList[x-1];
}
}
class PrimeApp : Application
{
FindPrime fp { };
void Main()
{
int num = argc > 1 ? atoi(argv[1]) : 1;
PrintLn(fp.GetPrime(num));
}
}
[edit] Eiffel
Works with: EiffelStudio version 6.6 beta (with provisional loop syntax)
class
APPLICATION
create
make
feature
make
-- Run application.
do
across primes_through (100) as ic loop print (ic.item.out + " ") end
end
primes_through (a_limit: INTEGER): LINKED_LIST [INTEGER]
-- Prime numbers through `a_limit'
require
valid_upper_limit: a_limit >= 2
local
l_tab: ARRAY [BOOLEAN]
do
create Result.make
create l_tab.make_filled (True, 2, a_limit)
across
l_tab as ic
loop
if ic.item then
Result.extend (ic.target_index)
across ((ic.target_index * ic.target_index) |..| l_tab.upper).new_cursor.with_step (ic.target_index) as id
loop
l_tab [id.item] := False
end
end
end
end
end
Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
[edit] F Sharp
[edit] Functional
This is not a Sieve of Eratosthenes, for the same reason as the naive Haskell version, because it simply filters the numbers.
#light
let not_divisible_by n =
List.filter (fun m -> m % n <> 0)
let sieve max =
let rec sieve_aux = function
| [] -> []
| x :: xs as l ->
if x * x > max then l else x :: sieve_aux (not_divisible_by x xs)
sieve_aux [2 .. max]
[edit] Forth
: prime? ( n -- ? ) here + c@ 0= ;
: composite! ( n -- ) here + 1 swap c! ;
: sieve ( n -- )
here over erase
2
begin
2dup dup * >
while
dup prime? if
2dup dup * do
i composite!
dup +loop
then
1+
repeat
drop
." Primes: " 2 do i prime? if i . then loop ;
100 sieve
[edit] Fortran
Fortran 90 and later
PROGRAM SIEVE
IMPLICIT NONE
INTEGER :: i, j, number
LOGICAL, ALLOCATABLE :: pflags(:)
WRITE(*,*) "Enter number to search to"
READ (*,*) number
ALLOCATE (pflags(2:number))
pflags = .TRUE.
DO i = 2, SQRT(REAL(number))
IF (pflags(i)) THEN
DO j = i+i, number, i
pflags(j) = .FALSE.
END DO
END IF
END DO
DO i = 2, number
IF (pflags(i)) WRITE(*,*) i
END DO
DEALLOCATE (pflags)
END PROGRAM SIEVE
Since 2 is trivially the only even prime number this version just sieves odd numbers
ALLOCATE (pflags((number-1)/2))
pflags = .TRUE.
DO i = 1, SQRT(REAL(number-1)/2)
IF (pflags(i)) THEN
DO j = 3*i+1, (number-1)/2, 2*i+1
pflags(j) = .FALSE.
END DO
END IF
END DO
DO i = 1, (number-1)/2
IF (pflags(i)) WRITE(*,*) 2*i + 1
END DO
[edit] Haskell
Imperative version, with unboxed arrays:
import Control.Monad (forM_, when)
import Control.Monad.ST
import Data.Array.ST
import Data.Array.Unboxed
primeSieve :: Integer -> UArray Integer Bool
primeSieve n = runSTUArray $ do
s <- newArray (2,n) True :: ST s (STUArray s Integer Bool)
let m = ceiling $ sqrt $ fromInteger n
forM_ [2..m] $ \i -> do
si <- readArray s i
when si $ do
forM_ [i*i,i*i+i..n] $ \j -> do
writeArray s j False
return s
Return primes from sieve as list:
primesTo :: Integer -> [Integer]
primesTo n = [p | p <- [2..n], s ! p] where s = primeSieve n
One often sees a purely functional version similar to
primes :: [Integer]
primes = sieve [2..] where
sieve (p:xs) = p : sieve [x | x <- xs, x `mod` p > 0]
As Melissa O'Neill demonstrates in The Genuine Sieve of Eratosthenes, this version is not better than trial division. She also presents an actual purely functional version that uses priority queues, and is much easier to combine with wheels (pre-calculated sieves for small primes) than the imperative version. Even combined with a simple wheel, the speed of this version is similar to the imperative version.
[edit] 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
[edit] J
sieve=: 3 : 0
m=. <.%:y
z=. $0
b=. y{.1
while. m>j=. 1+b i. 0 do.
b=. b+.y$(-j){.1
z=. z,j
end.
z,1+I.-.b
)
"Wheels" are readily implemented as follows:
sieve1=: 3 : 0
m=. <.%:y
z=. y (>:#]) 2 3 5 7
b=. 1,}.y$+./(*/z)$&>(-z){.&.>1
while. m>j=. 1+b i. 0 do.
b=. b+.y$(-j){.1
z=. z,j
end.
z,1+I.-.b
)
The use of 2 3 5 7 as wheels provides a 20% time improvement for n=1000 and 2% for n=1e6 .
[edit] Java
Works with: Java version 1.5+
import java.util.LinkedList;
public class Sieve{
public static LinkedList<Integer> sieve(int n){
if(n < 2) return new LinkedList<Integer>();
LinkedList<Integer> primes = new LinkedList<Integer>();
LinkedList<Integer> nums = new LinkedList<Integer>();
for(int i = 2;i <= n;i++){ //unoptimized
nums.add(i);
}
while(nums.size() > 0){
int nextPrime = nums.remove();
for(int i = nextPrime * nextPrime;i <= n;i += nextPrime){
nums.removeFirstOccurrence(i);
}
primes.add(nextPrime);
}
return primes;
}
}
To optimize by testing only odd numbers, replace the loop marked "unoptimized" with these lines:
nums.add(2);
for(int i = 3;i <= n;i += 2){
nums.add(i);
}
Version using a BitSet:
import java.util.LinkedList;
import java.util.BitSet;
public class Sieve{
public static LinkedList<Integer> sieve(int n){
LinkedList<Integer> primes = new LinkedList<Integer>();
BitSet nonPrimes = new BitSet(n+1);
for (int p = 2; p <= n; p = nonPrimes.nextClearBit(p+1)) {
for (int i = p * p; i <= n; i += p)
nonPrimes.set(i);
primes.add(p);
}
return primes;
}
}
[edit] JavaScript
function erasthenes(limit) {
var primes = new Array();
if (limit >= 2) {
var nums = new Array(limit-1);
for (var i = 2; i <= limit; i++)
nums[i-2] = i;
var idx = 0;
while ((last_prime = nums[idx]) <= Math.sqrt(limit)) {
if (last_prime != null) {
for (var i = idx + last_prime; i < limit - 1; i += last_prime)
nums[i] = null;
}
idx++;
}
for (var i = 0; i < nums.length; i++)
if (nums[i] != null)
primes.push(nums[i]);
}
return primes;
}
var primes = erasthenes(100);
try {
WScript.Echo(primes); // testing in WSH
} catch(err) {
print(primes); // testing in Rhino
}
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
[edit] Logo
to sieve :limit
make "a (array :limit 2) ; initialized to empty lists
make "p []
for [i 2 :limit] [
if empty? item :i :a [
queue "p :i
for [j [:i * :i] :limit :i] [setitem :j :a :i]
]
]
output :p
end
print sieve 100 ; 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
[edit] 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
[edit] Lucid
prime
where
prime = 2 fby (n whenever isprime(n));
n = 3 fby n+2;
isprime(n) = not(divs) asa divs or prime*prime > N
where
N is current n;
divs = N mod prime eq 0;
end;
end
[edit] recursive
sieve( N )
where
N = 2 fby N + 1;
sieve( i ) =
i fby sieve ( i whenever i mod first i ne 0 ) ;
end
[edit] M4
define(`lim',100)dnl
define(`for',
`ifelse($#,0,
``$0'',
`ifelse(eval($2<=$3),1,
`pushdef(`$1',$2)$5`'popdef(`$1')$0(`$1',eval($2+$4),$3,$4,`$5')')')')dnl
for(`j',2,lim,1,
`ifdef(a[j],
`',
`j for(`k',eval(j*j),lim,j,
`define(a[k],1)')')')
Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
[edit] MAXScript
fn eratosthenes n =
(
multiples = #()
print 2
for i in 3 to n do
(
if (findItem multiples i) == 0 then
(
print i
for j in (i * i) to n by i do
(
append multiples j
)
)
)
)
eratosthenes 100
[edit] Modula-3
[edit] Regular version
MODULE Prime EXPORTS Main;
IMPORT IO;
CONST LastNum = 1000;
VAR a: ARRAY [2..LastNum] OF BOOLEAN;
BEGIN
FOR i := FIRST(a) TO LAST(a) DO
a[i] := TRUE;
END;
FOR i := FIRST(a) TO LAST(a) DO
IF a[i] THEN
IO.PutInt(i);
IO.Put(" ");
FOR j := FIRST(a) TO LAST(a) DO
IF j MOD i = 0 THEN
a[j] := FALSE;
END;
END;
END;
END;
IO.Put("\n");
END Prime.
[edit] Advanced version
This version uses more "advanced" types.
(* From the CM3 examples folder (comments removed). *)
MODULE Sieve EXPORTS Main;
IMPORT IO;
TYPE
Number = [2..1000];
Set = SET OF Number;
VAR
prime: Set := Set {FIRST(Number) .. LAST(Number)};
BEGIN
FOR i := FIRST(Number) TO LAST(Number) DO
IF i IN prime THEN
IO.PutInt(i);
IO.Put(" ");
FOR j := i TO LAST(Number) BY i DO
prime := prime - Set{j};
END;
END;
END;
IO.Put("\n");
END Sieve.
[edit] Nial
primes is sublist [ each (2 = sum eachright (0 = mod) [pass,count]), pass ] rest count
Using it
|primes 10 =2 3 5 7
[edit] Oberon-2
MODULE Primes;
IMPORT Out, Math;
CONST N = 1000;
VAR a: ARRAY N OF BOOLEAN;
i, j, m: INTEGER;
BEGIN
(* Set all elements of a to TRUE. *)
FOR i := 1 TO N - 1 DO
a[i] := TRUE;
END;
(* Compute square root of N and convert back to INTEGER. *)
m := ENTIER(Math.Sqrt(N));
FOR i := 2 TO m DO
IF a[i] THEN
FOR j := 2 TO (N - 1) DIV i DO
a[i*j] := FALSE;
END;
END;
END;
(* Print all the elements of a that are TRUE. *)
FOR i := 2 TO N - 1 DO
IF a[i] THEN
Out.Int(i, 4);
END;
END;
Out.Ln;
END Primes.
[edit] OCaml
[edit] Imperative
let sieve n =
let is_prime = Array.create n true in
let limit = truncate(sqrt (float n)) in
for i = 2 to pred limit do
if is_prime.(i) then
let j = ref (i*i) in
try while true do
is_prime.(!j) <- false;
j := !j + i;
done with _ -> ()
done;
(is_prime)
let primes n =
let _, primes =
Array.fold_left (fun (i,acc) -> function
| true -> (i+1, i::acc)
| false -> (i+1, acc))
(0, [])
(sieve n)
in
List.tl(List.tl(List.rev primes))
in the top-level:
# primes 100 ;; - : int list = [2; 3; 5; 7; 11; 13; 17; 19; 23; 29; 31; 37; 41; 43; 47; 53; 59; 61; 67; 71; 73; 79; 83; 89; 97]
[edit] Functional
(* first define some iterators *)
# let fold_iter f init a b =
let rec aux acc i =
if i > b
then (acc)
else aux (f acc i) (succ i)
in
aux init a ;;
val fold_iter : ('a -> int -> 'a) -> 'a -> int -> int -> 'a = <fun>
# let fold_step f init a b step =
let rec aux acc i =
if i > b
then (acc)
else aux (f acc i) (i + step)
in
aux init a ;;
val fold_step : ('a -> int -> 'a) -> 'a -> int -> int -> int -> 'a = <fun>
(* remove a given value from a list *)
# let remove li v =
let rec aux acc = function
| hd::tl when hd = v -> (List.rev_append acc tl)
| hd::tl -> aux (hd::acc) tl
| [] -> li
in
aux [] li ;;
val remove : 'a list -> 'a -> 'a list = <fun>
(* the main function *)
# let primes n =
let li =
(* create a list [from 2; ... until n] *)
List.rev(fold_iter (fun acc i -> (i::acc)) [] 2 n)
in
let limit = truncate(sqrt(float n)) in
fold_iter (fun li i ->
if List.mem i li (* test if (i) is prime *)
then (fold_step remove li (i*i) n i)
else li)
li 2 (pred limit)
;;
val primes : int -> int list = <fun>
# primes 200 ;;
- : int list =
[2; 3; 5; 7; 11; 13; 17; 19; 23; 29; 31; 37; 41; 43; 47; 53; 59; 61; 67; 71;
73; 79; 83; 89; 97; 101; 103; 107; 109; 113; 127; 131; 137; 139; 149; 151;
157; 163; 167; 173; 179; 181; 191; 193; 197; 199]
[edit] Another functional version
This uses zero to denote struck-out numbers. It is slightly inefficient as it strikes-out multiples above p rather than p2
# let rec strike_nth k n l = match l with
| [] -> []
| h :: t ->
if k = 0 then 0 :: strike_nth (n-1) n t
else h :: strike_nth (k-1) n t;;
val strike_nth : int -> int -> int list -> int list = <fun>
# let primes n =
let limit = truncate(sqrt(float n)) in
let rec range a b = if a > b then [] else a :: range (a+1) b in
let rec sieve_primes l = match l with
| [] -> []
| 0 :: t -> sieve_primes t
| h :: t -> if h > limit then List.filter ((<) 0) l else
h :: sieve_primes (strike_nth (h-1) h t) in
sieve_primes (range 2 n) ;;
val primes : int -> int list = <fun>
# primes 200;;
- : int list =
[2; 3; 5; 7; 11; 13; 17; 19; 23; 29; 31; 37; 41; 43; 47; 53; 59; 61; 67; 71;
73; 79; 83; 89; 97; 101; 103; 107; 109; 113; 127; 131; 137; 139; 149; 151;
157; 163; 167; 173; 179; 181; 191; 193; 197; 199]
[edit] Oz
Translation of: Haskell
declare
fun {Sieve N}
S = {Array.new 2 N true}
M = {Float.toInt {Sqrt {Int.toFloat N}}}
in
for I in 2..M do
if S.I then
for J in I*I..N;I do
S.J := false
end
end
end
S
end
fun {Primes N}
S = {Sieve N}
in
for I in 2..N collect:C do
if S.I then {C I} end
end
end
in
{Show {Primes 30}}
[edit] Pascal
Note: Some Pascal implementations put quite low limits on the size of a set (e.g. Turbo Pascal doesn't allow more than 256 members). To compile on such an implementation, reduce the constant PrimeLimit accordingly.
program primes(output)
const
PrimeLimit = 1000;
var
primes: set of 1 .. PrimeLimit;
n, k: integer;
needcomma: boolean;
begin
{ calculate the primes }
primes := [2 .. PrimeLimit];
for n := 1 to trunc(sqrt(PrimeLimit)) do
begin
if n in primes
then
begin
k := n*n;
while k < PrimeLimit do
begin
primes := primes - [k];
k := k + n
end
end
end;
{ output the primes }
needcomma := false;
for n := 1 to PrimeLimit do
if n in primes
then
begin
if needcomma
then
write(', ');
write(n);
needcomma := true
end
end.
[edit] Perl
This is not a Sieve of Eratosthenes, because it simply filters the numbers.
sub findprimes {
my $prime_limit=shift;
my @possible_primes;
my @primes;
my $prime;
if ($prime_limit < 2) {
return @primes;
}
@possible_primes=(2..$prime_limit);
do {
$prime=shift(@possible_primes);
push @primes, $prime;
@possible_primes = grep {$_%$prime!=0} @possible_primes;
} while ($possible_primes[0] < sqrt($prime_limit));
return (@primes,@possible_primes);
}
my $prime_limit;
print "Input the number that you wish to test :";
chomp($prime_limit=<STDIN>);
my @primes = findprimes($prime_limit);
print "The primes are: @primes!!!\n";
Version using bit strings:
sub sieve {
my ($n) = @_;
my @primes;
my $nonPrimes = '';
foreach my $p (2 .. $n) {
unless (vec($nonPrimes, $p, 1)) {
for (my $i = $p * $p; $i <= $n; $i += $p) {
vec($nonPrimes, $i, 1) = 1;
}
push @primes, $p;
}
}
@primes
}
my $prime_limit;
print "Input the number that you wish to test :";
chomp($prime_limit = <STDIN>);
my @primes = sieve($prime_limit);
print "The primes are: @primes!!!\n";
Golfing a bit:
sub sieve{ my (@s, $i);
grep { not $s[ $i = $_ ] and do
{ $s[ $i += $_ ]++ while $i <= $_[0]; 1 }
} 2..$_[0]
}
print join ", " => sieve 100;
Or with bit strings:
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;
[edit] PHP
function iprimes_upto($limit)
{
for ($i = 2; $i < $limit; $i++)
{
$primes[$i] = true;
}
for ($n = 2; $n < $limit; $n++)
{
if ($primes[$n])
{
for ($i = $n*$n; $i < $limit; $i += $n)
{
$primes[$i] = false;
}
}
}
return $primes;
}
[edit] PicoLisp
(de sieve (N)
(let Sieve (range 1 N)
(set Sieve)
(for I (cdr Sieve)
(when I
(for (S (nth Sieve (* I I)) S (nth (cdr S) I))
(set S) ) ) )
(filter bool Sieve) ) )
Output:
: (sieve 100) -> (2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)
[edit] PL/I
eratos: procedure options (main, reorder);
declare n fixed binary (31);
get list (n);
begin;
declare s(n) bit (1) unaligned;
declare (false initial ('0'b), true initial ('1'b)) bit (1);
declare (f, i, j) fixed binary;
s = true;
f = 2;
loop:
do j = 2 to sqrt(n);
do i = f to n by f;
if s(i) then s(i) = false;
end;
do i = f+1 to sqrt(n) by 1;
if s(i) then do; f = i; iterate loop; end;
end;
leave;
end;
put skip list ('The primes are:');
do i = 2 to n;
if s(i) then put skip list (i);
end;
end;
end eratos;
[edit] Pop11
define eratostenes(n);
lvars bits = inits(n), i, j;
for i from 2 to n do
if bits(i) = 0 then
printf('' >< i, '%s\n');
for j from 2*i by i to n do
1 -> bits(j);
endfor;
endif;
endfor;
enddefine;
[edit] PureBasic
[edit] Basic procedure
For n=2 To Sqr(lim)
If Nums(n)=0
m=n*n
While m<=lim
Nums(m)=1
m+n
Wend
EndIf
Next n
[edit] Working example
When compiled for Windows x86 using PureBasic 4.41 , this program is only 7 kB.
Dim Nums.i(0)
Define l, n, m, lim
If OpenConsole()
; Ask for the limit to search, get that input and allocate a Array
Print("Enter limit for this search: ")
lim=Val(Input())
ReDim Nums(lim)
; Use a basic Sieve of Eratosthenes
For n=2 To Sqr(lim)
If Nums(n)=#False
m=n*n
While m<=lim
Nums(m)=#True
m+n
Wend
EndIf
Next n
;Present the result to our user
PrintN(#CRLF$+"The Prims up to "+Str(lim)+" are;")
m=0: l=Log10(lim)+1
For n=2 To lim
If Nums(n)=#False
Print(RSet(Str(n),l)+" ")
m+1
If m>72/(l+1)
m=0: PrintN("")
EndIf
EndIf
Next
Print(#CRLF$+#CRLF$+"Press ENTER to exit"): Input()
CloseConsole()
EndIf
Output may look like;
Enter limit for this search: 750 The Prims up to 750 are; 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691 701 709 719 727 733 739 743 Press ENTER to exit
[edit] Python
[edit] Testing on list of primes up to that number
This is not a Sieve of Eratosthenes, it is just trial division.
for x in range(1,100000)
z=0
for y in p:
if x%y==0:
z=1
if z==0:
p.append(x)
[edit] Using list of primes up to that number using list comprehension
This is not a Sieve of Eratosthenes, because it simply filters the numbers.
def GetPrimes(p):
i = 1
while i**2 < p[-1]:
p[i:] = [x for x in p[i:] if x % p[i-1] != 0]
i += 1
return p
print GetPrimes(range(2,100))
[edit] Using list lookup
def eratosthenes(n):
multiples = []
print 2
for i in xrange(3, n+1):
if i not in multiples:
print i
for j in xrange(i*i, n+1, i):
multiples.append(j)
eratosthenes(100)
[edit] Using set lookup
The version below uses a set to store the multiples. set objects are much faster (usually O(log n)) than lists (O(n)) for checking if a given object is a member. Using the set.update method avoids explicit iteration in the interpreter, giving a further speed improvement.
def eratosthenes2(n):
multiples = set()
print 2
for i in xrange(3, n+1):
if i not in multiples:
print i
multiples.update(xrange(i*i, n+1, i))
eratosthenes2(100)
[edit] Using array lookup
The version below uses array lookup to test for primality. The function primes_upto() is a straightforward implementation of Sieve of Eratosthenesalgorithm. It returns prime numbers less than or equal to limit.
def primes_upto(limit):
is_prime = [False] * 2 + [True] * (limit - 1)
for n in xrange(int(limit**0.5 + 1.5)): # stop at ``sqrt(limit)``
if is_prime[n]:
for i in xrange(n*n, limit+1, n):
is_prime[i] = False
return [i for i, prime in enumerate(is_prime) if prime]
[edit] Using generator
def iprimes_upto(limit):
is_prime = [False] * 2 + [True] * (limit - 1)
for n in range(limit + 1):
if is_prime[n]:
yield n
for i in range(n*n, limit+1, n): # start at ``n`` squared
is_prime[i] = False
Example:
>>> list(iprimes_upto(15))
[2, 3, 5, 7, 11, 13]
[edit] Using numpy
Library: numpy Belowed code adapted from literateprograms.org using numpy
import numpy
def primes_upto2(limit):
is_prime = numpy.ones(limit + 1, dtype=numpy.bool)
for n in xrange(2, int(limit**0.5 + 1.5)):
if is_prime[n]:
is_prime[n*n::n] = 0
return numpy.nonzero(is_prime)[0][2:]
Performance note: there is no point to add wheels here, due to execution of p[n*n::n] = 0 and nonzero() takes us almost all time.
[edit] Using wheels
Version with wheel based optimization:
from numpy import array, bool_, multiply, nonzero, ones, put, resize
#
def makepattern(smallprimes):
pattern = ones(multiply.reduce(smallprimes), dtype=bool_)
pattern[0] = 0
for p in smallprimes:
pattern[p::p] = 0
return pattern
#
def primes_upto3(limit, smallprimes=(2,3,5,7,11)):
sp = array(smallprimes)
if limit <= sp.max(): return sp[sp <= limit]
#
isprime = resize(makepattern(sp), limit + 1)
isprime[:2] = 0; put(isprime, sp, 1)
#
for n in range(sp.max() + 2, int(limit**0.5 + 1.5), 2):
if isprime[n]:
isprime[n*n::n] = 0
return nonzero(isprime)[0]
Examples:
>>> primes_upto3(10**6, smallprimes=(2,3)) # Wall time: 0.17
array([ 2, 3, 5, ..., 999961, 999979, 999983])
>>> primes_upto3(10**7, smallprimes=(2,3)) # Wall time: '''2.13'''
array([ 2, 3, 5, ..., 9999971, 9999973, 9999991])
>>> primes_upto3(15)
array([ 2, 3, 5, 7, 11, 13])
>>> primes_upto3(10**7, smallprimes=primes_upto3(15)) # Wall time: '''1.31'''
array([ 2, 3, 5, ..., 9999971, 9999973, 9999991])
>>> primes_upto2(10**7) # Wall time: '''1.39''' <-- version ''without'' wheels
array([ 2, 3, 5, ..., 9999971, 9999973, 9999991])
>>> primes_upto3(10**7) # Wall time: '''1.30'''
array([ 2, 3, 5, ..., 9999971, 9999973, 9999991])
The above-mentioned examples demonstrate that the given wheel based optimization does not show significant performance gain.
[edit] Infinite generator
A generator that will generate primes indefinitely (perhaps until it runs out of memory).
Works with: Python version 2.6, 3.x
import heapq
# generates all prime numbers
def sieve():
# priority queue of the sequences of non-primes
# the priority queue allows us to get the "next" non-prime quickly
nonprimes = []
i = 2
while True:
if nonprimes and i == nonprimes[0][0]: # non-prime
while nonprimes[0][0] == i:
# for each sequence that generates this number,
# have it go to the next number (simply add the prime)
# and re-position it in the priority queue
x = nonprimes[0]
x[0] += x[1]
heapq.heapreplace(nonprimes, x)
else: # prime
# insert a 2-element list into the priority queue:
# [current multiple, prime]
# the first element allows sorting by value of current multiple
# we start with i^2
heapq.heappush(nonprimes, [i*i, i])
yield i
i += 1
Example:
>>> foo = sieve() >>> for i in range(20): ... print(next(foo)) ... 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71
[edit] Infinite generator with a faster algorithm
Works with: Python version 2.4 through 3.1
def primes():
composites = {} # A mapping from composite numbers to the smallest prime
# that is a factor of it (its witness).
n = 2 # The current number being considered as a prime.
while True:
if n not in composites:
yield n # Not a composite, therefore prime.
composites[n**2] = n # The next unseen composite number is n squared.
else:
# n is composite. Find the next unseen composite number with the
# same witness as n.
witness = composites.pop(n)
next = n + witness
while next in composites:
next += witness
composites[next] = witness
n += 1
import itertools
def primes_up_to(limit):
return list(itertools.takewhile(lambda p: p <= limit, primes()))
[edit] R
This code is vectorised, so it runs quickly for n < one million. The allocation of the primes vector means memory usage is too high for n much larger than that.
sieve <- function(n)
{
n <- as.integer(n)
if(n > 1e6) stop("n too large")
primes <- rep(TRUE, n)
primes[1] <- FALSE
last.prime <- 2L
for(i in last.prime:floor(sqrt(n)))
{
primes[seq.int(2L*last.prime, n, last.prime)] <- FALSE
last.prime <- last.prime + min(which(primes[(last.prime+1):n]))
}
which(primes)
}
sieve(1000)
[edit] Ruby
def erasthenes(n)
nums = [0, 0] + (2..n).to_a
(2**2..n).step(2) {|m| nums[m] = 0}
(3..Math.sqrt(n).to_i).step(2) do |i|
if nums[i].nonzero?
(i**2..n).step(i) {|m| nums[m] = 0}
end
end
nums.find_all {|m| m.nonzero?}
end
p erasthenes(100)
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
[edit] Scala
def sieveOfEratosthenes(n: Int) = {
val primes = new scala.collection.mutable.BitSet(n)
primes ++= (2 to n)
val sqrt = Math.sqrt(n).toInt
for {
candidate <- 2 to sqrt
if primes contains candidate
} primes --= candidate * candidate to n by candidate
primes
}
[edit] Scheme
Works with: Chicken Scheme
(use srfi-1) ; for FILTER
(define (range a b)
(if (> a b)
'()
(cons a (range (+ a 1) b))))
(define (divisible? x y)
(= (remainder x y) 0))
(define (sieve list)
(if (null? list)
'()
(cons
(car list)
(sieve (filter (lambda (x) (not (divisible? x (car list)))) (cdr list))))))
(define (primes num)
(sieve (range 2 num)))
Output:
#;1> (primes 100) (2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)
[edit] Seed7
const func array boolean: eratosthenes (in integer: n) is func
result
var array boolean: sieve is 0 times FALSE;
local
var integer: i is 0;
var integer: j is 0;
begin
sieve := n times TRUE;
sieve[1] := FALSE;
for i range 2 to sqrt(n) do
if sieve[i] then
for j range i * i to n step i do
sieve[j] := FALSE;
end for;
end if;
end for;
end func;
Original source: [1]
[edit] Tcl
package require Tcl 8.5
proc sieve n {
if {$n < 2} {return {}}
# create a container to hold the sequence of numbers.
# use a dictionary for its speedy access (like an associative array)
# and for its insertion order preservation (like a list)
set nums [dict create]
for {set i 2} {$i <= $n} {incr i} {
# the actual value is never used
dict set nums $i ""
}
set primes [list]
while {[set nextPrime [lindex [dict keys $nums] 0]] <= sqrt($n)} {
dict unset nums $nextPrime
for {set i [expr {$nextPrime ** 2}]} {$i <= $n} {incr i $nextPrime} {
dict unset nums $i
}
lappend primes $nextPrime
}
return [concat $primes [dict keys $nums]]
}
puts [sieve 100] ;# 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
[edit] UnixPipes
# Uses file a as a buffer for sequences.
echo > a && yes \ | cat -n | tail +2|
while read x ; do
( (echo 1 ; cat -s a) | xargs -n 1 expr $x % | grep "^0$" | wc -l | grep "^ *1 *$" > /dev/null && echo $x |tee -a a )
done
[edit] Ursala
with no optimizations
#import nat
sieve = ~<{0,1}&& iota; @NttPX ~&r->lx ^/~&rhPlC remainder@rlX~|@r
test program:
#cast %nL
example = sieve 50
output:
<2,3,5,7,11,13,17,19,23,29,31,37,41,43,47>
[edit] Vedit macro language
This implementation uses an edit buffer as an array for flags. After the macro has been run, you can see how the primes are located in the array. Primes are marked with 'P' and non-primes with '-'. The first character position represents number 0.
#10 = Get_Num("Enter number to search to: ", STATLINE)
Buf_Switch(Buf_Free) // Use edit buffer as flags array
Ins_Text("--") // 0 and 1 are not primes
Ins_Char('P', COUNT, #10-1) // init rest of the flags to "prime"
for (#1 = 2; #1*#1 < #10; #1++) {
Goto_Pos(#1)
if (Cur_Char=='P') { // this is a prime number
for (#2 = #1*#1; #2 <= #10; #2 += #1) {
Goto_Pos(#2)
Ins_Char('-', OVERWRITE)
}
}
}
Sample output showing numbers in range 0 to 599.
--PP-P-P---P-P---P-P---P-----P-P-----P---P-P---P-----P-----P -P-----P---P-P-----P---P-----P-------P---P-P---P-P---P------ -------P---P-----P-P---------P-P-----P-----P---P-----P-----P -P---------P-P---P-P-----------P-----------P---P-P---P-----P -P---------P-----P-----P-----P-P-----P---P-P---------P------ -------P---P-P---P-------------P-----P---------P-P---P-----P -------P-----P-----P---P-----P-------P---P-------P---------P -P---------P-P-----P---P-----P-------P---P-P---P-----------P -------P---P-------P---P-----P-----------P-P---------------- -P-----P---------P-----P-----P-P-----P---------P-----P-----P
[edit] Visual Basic
Dim n As Integer, k As Integer, limit As Integer
Console.WriteLine("Enter number to search to: ")
limit = Console.ReadLine
Dim flags(limit) As Integer
For n As Integer = 2 To Math.Sqrt(limit)
If flags(n) = 0 Then
For k As Integer = n * n To limit Step n
flags(k) = 1
Next k
End If
Next n
' Display the primes
For n = 2 To limit
If flags(n) = 0 Then
Console.WriteLine(n)
End If
Next n
[edit] Vorpal
self.print_primes = method(m){
p = new()
p.fill(0, m, 1, true)
count = 0
i = 2
while(i < m){
if(p[i] == true){
p.fill(i+i, m, i, false)
count = count + 1
}
i = i + 1
}
('primes: ' + count + ' in ' + m).print()
for(i = 2, i < m, i = i + 1){
if(p[i] == true){
('' + i + ', ').put()
}
}
''.print()
}
self.print_primes(100)
Result:
primes: 25 in 100 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,
Categories: Programming Tasks | Prime Numbers | Clarified and Needing Review | Ada | ALGOL 68 | AutoHotkey | AWK | BASIC | Befunge | C | C++ | Common Lisp | D | E | EC | Eiffel | F Sharp | Forth | Fortran | Haskell | Icon | J | Java | JavaScript | Logo | Lua | Lucid | M4 | MAXScript | Modula-3 | Nial | Oberon-2 | OCaml | Oz | Pascal | Perl | PHP | PicoLisp | PL/I | Pop11 | PureBasic | Python | Numpy | R | Ruby | Scala | Scheme | Seed7 | Tcl | UnixPipes | Ursala | Vedit macro language | Visual Basic | Vorpal
