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.
Ada
This solution uses an array of bits representing boolean values to identify the prime numbers. The limit on the outer loop is set to the square root of the number received from the command line. This limit minimizes the number of unneeded iterations in the algorithm.
The lowest index of the array of boolean is set to 2 so that all calculations begin at 2.
<lang ada>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;</lang>
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.
<lang ada>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;</lang>
This version is the same as the previous version, but with a wheel of 2.
<lang ada>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;</lang>
ALGOL 68
<lang algol68>BOOL prime = TRUE, non prime = FALSE; PROC eratosthenes = (INT n)[]BOOL: (
[n]BOOL sieve;
FOR i TO UPB sieve DO sieve[i] := prime OD;
INT m = ENTIER sqrt(n);
sieve[1] := non prime;
FOR i FROM 2 TO m DO
IF sieve[i] = prime THEN
FOR j FROM i*i BY i TO n DO
sieve[j] := non prime
OD
FI
OD;
sieve
);
print((eratosthenes(80),new line))</lang>
Output:
FTTFTFTFFFTFTFFFTFTFFFTFFFFFTFTFFFFFTFFFTFTFFFTFFFFFTFFFFFTFTFFFFFTFFFTFTFFFFFTF
AutoHotkey
Search autohotkey.com: of Eratosthenes
Source: AutoHotkey forum by Laszlo
<lang autohotkey>
MsgBox % "12345678901234567890`n" Sieve(20)
Sieve(n) { ; Sieve of Eratosthenes => string of 0|1 chars, 1 at position k: k is prime
Static zero := 48, one := 49 ; Asc("0"), Asc("1")
VarSetCapacity(S,n,one)
NumPut(zero,S,0,"char")
i := 2
Loop % sqrt(n)-1 {
If (NumGet(S,i-1,"char") = one)
Loop % n//i
If (A_Index > 1)
NumPut(zero,S,A_Index*i-1,"char")
i += 1+(i>2)
}
Return S
} </lang>
AWK
An initial array holds all numbers 2..max (which is entered on stdin); then all products of integers are deleted from it; the remaining are displayed in the unsorted appearance of a hash table.
$ awk '{for(i=2;i<=$1;i++) a[i]=1;
> for(i=2;i<=sqrt($1);i++) for(j=2;j<=$1;j++) delete a[i*j];
> for(i in a) printf i" "}'
100
71 53 17 5 73 37 19 83 47 29 7 67 59 11 97 79 89 31 13 41 23 2 61 43 3
The following variant does not unset non-primes, but set them to 0, to preserve order in output:
$ awk '{for(i=2;i<=$1;i++) a[i]=1;
> for(i=2;i<=sqrt($1);i++) for(j=2;j<=$1;j++) a[i*j]=0;
> for(i=2;i<=$1;i++) if(a[i])printf i" "}'
100
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
BASIC
<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>
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>
Befunge
2>:3g" "-!v\ g30 < |!`"O":+1_:.:03p>03g+:"O"`| @ ^ p3\" ":< 2 234567890123456789012345678901234567890123456789012345678901234567890123456789
C
Plain sieve, without any optimizations:
<lang c>
- 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;
} </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>
Clojure
<lang lisp>(defn primes-up-to [n]
(let [n (int n)]
"Returns a list of all primes from 2 to n"
(let [root (int (Math/round (Math/floor (Math/sqrt n))))]
(loop [i (int 3)
a (int-array n)
result (transient [2])]
(if (>= i n)
(persistent! result)
(recur (+ i (int 2))
(if (< i root)
(loop [arr a
inc (+ i i)
j (* i i)]
(if (>= j n)
arr
(recur (do (aset arr j (int 1)) arr)
inc
(+ j inc))))
a)
(if (zero? (aget a i))
(conj! result i)
result)))))))</lang>
Or as a lazy sequence:
<lang lisp>(def prime-gen
(let [primes (atom [])]
(for [n (iterate inc 2)
:when (not-any? #(zero? (rem n %))
(filter #(<= % (Math/sqrt n))
@primes))]
(do (swap! primes conj n)
n))))</lang>
Common Lisp
<lang 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)))))</lang>
D
<lang d>import std.stdio, std.math, std.algorithm, std.range;
uint[] sieve(uint n) {
bool[] comp = new bool[n + 1];
comp[0 .. 2] = true;
foreach (i; 2 .. cast(uint)sqrt(n))
if (!comp[i])
for (uint j = i * i; j <= n; j += i)
comp[j] = true;
return array(filter!((i){ return !comp[i]; })(iota(n+1)));
}
void main() {
writeln(sieve(50));
}</lang> Output:
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49]
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
F#
Functional
This is not a Sieve of Eratosthenes, for the same reason as the naive Haskell version, because it simply filters the numbers.
#light
let not_divisible_by n =
List.filter (fun m -> m % n <> 0)
let sieve max =
let rec sieve_aux = function
| [] -> []
| x :: xs as l ->
if x * x > max then l else x :: sieve_aux (not_divisible_by x xs)
sieve_aux [2 .. max]
Forth
: prime? ( n -- ? ) here + c@ 0= ;
: composite! ( n -- ) here + 1 swap c! ;
: sieve ( n -- )
here over erase
2
begin
2dup dup * >
while
dup prime? if
2dup dup * do
i composite!
dup +loop
then
1+
repeat
drop
." Primes: " 2 do i prime? if i . then loop ;
100 sieve
Fortran
<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
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;
od;
return Filtered([1 .. n], x -> sieve[x] = 0);
end;
Eratosthenes(100);
- [ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 ]</lang>
Go
<lang go>
package main
import (
"fmt"
)
func main() {
const limit = 201 // means sieve numbers < 201
// sieve
c := make([]bool, limit) // c for composite. false means prime candidate
c[1] = true // 1 not considered prime
p := 2
for {
// first allowed optimization: outer loop only goes to sqrt(limit)
p2 := p * p
if p2 >= limit {
break
}
// second allowed optimization: inner loop starts at sqr(p)
for i := p2; i < limit; i += p {
c[i] = true // it's a composite
}
// scan to get next prime for outer loop
for {
p++
if !c[p] {
break
}
}
}
// sieve complete. now print a representation.
for n := 1; n < limit; n++ {
if c[n] {
fmt.Print(" .")
} else {
fmt.Printf("%3d", n)
}
if n%20 == 0 {
fmt.Println("")
}
}
} </lang> Output:
. 2 3 . 5 . 7 . . . 11 . 13 . . . 17 . 19 . . . 23 . . . . . 29 . 31 . . . . . 37 . . . 41 . 43 . . . 47 . . . . . 53 . . . . . 59 . 61 . . . . . 67 . . . 71 . 73 . . . . . 79 . . . 83 . . . . . 89 . . . . . . . 97 . . . 101 .103 . . .107 .109 . . .113 . . . . . . . . . . . . .127 . . .131 . . . . .137 .139 . . . . . . . . .149 .151 . . . . .157 . . . . .163 . . .167 . . . . .173 . . . . .179 . 181 . . . . . . . . .191 .193 . . .197 .199 .
Groovy
This solution uses a BitSet for compactness and speed, but in Groovy, BitSet has full List semantics. It also uses both the "square root of the boundary" shortcut and the "square of the prime" shortcut. <lang groovy>def sievePrimes = { bound ->
def isPrime = new BitSet(bound)
isPrime[0..1] = false
isPrime[2..<bound] = true
(2..(Math.sqrt(bound))).each { pc ->
if (isPrime[pc]) {
((pc**2)..<bound).step(pc) { isPrime[it] = false }
}
}
(0..<bound).findAll { isPrime[it] }
}</lang>
Test: <lang groovy>println sievePrimes(100)</lang>
Output:
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
Haskell
Imperative version, with unboxed arrays: <lang haskell> import Control.Monad (forM_, when) import Control.Monad.ST
import Data.Array.ST import Data.Array.Unboxed
primeSieve :: Integer -> UArray Integer Bool primeSieve 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
</lang>
Return primes from sieve as list: <lang haskell> primesTo :: Integer -> [Integer] primesTo n = [p | p <- [2..n], s ! p] where s = primeSieve n </lang> One often sees a purely functional version similar to <lang haskell> primes :: [Integer] primes = sieve [2..] where
sieve (p:xs) = p : sieve [x | x <- xs, x `mod` p > 0]
</lang> 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.
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>
JavaScript
<lang javascript>function erasthenes(limit) {
var primes = [];
if (limit >= 2) {
var nums = new Array(limit-1);
for (var i = 2; i <= limit; i++)
nums[i-2] = i;
var 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
}</lang> outputs:
2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97
Liberty BASIC
<lang lb> 'Notice that arrays are globally visible to functions.
'The sieve() function uses the flags() array. 'This is a Sieve benchmark adapted from BYTE 1985 ' May, page 286
size = 7000
dim flags(7001)
start = time$("ms")
print sieve(size); " primes found."
print "End of iteration. Elapsed time in milliseconds: "; time$("ms")-start
end
function sieve(size)
for i = 0 to size
if flags(i) = 0 then
prime = i + i + 3
k = i + prime
while k <= size
flags(k) = 1
k = k + prime
wend
sieve = sieve + 1
end if
next i
end function</lang>
Logo
to sieve :limit
make "a (array :limit 2) ; initialized to empty lists
make "p []
for [i 2 :limit] [
if empty? item :i :a [
queue "p :i
for [j [:i * :i] :limit :i] [setitem :j :a :i]
]
]
output :p
end
print sieve 100 ; 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
Lua
<lang lua>function erato(n)
local t = {0, 2}
for i = 3, n, 2 do t[i], t[i+1] = i, 0 end
for i = 3, math.sqrt(n) do for j = i*i, n, 2*i do t[j] = 0 end end
return t
end</lang>
Lucid
prime
where
prime = 2 fby (n whenever isprime(n));
n = 3 fby n+2;
isprime(n) = not(divs) asa divs or prime*prime > N
where
N is current n;
divs = N mod prime eq 0;
end;
end
recursive
sieve( N )
where
N = 2 fby N + 1;
sieve( i ) =
i fby sieve ( i whenever i mod first i ne 0 ) ;
end
M4
<lang M4> define(`lim',100)dnl define(`for',
`ifelse($#,0,
``$0,
`ifelse(eval($2<=$3),1,
`pushdef(`$1',$2)$5`'popdef(`$1')$0(`$1',eval($2+$4),$3,$4,`$5')')')')dnl
for(`j',2,lim,1,
`ifdef(a[j],
`',
`j for(`k',eval(j*j),lim,j,
`define(a[k],1)')')')
</lang>
Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
MATLAB
<lang MATLAB>function primeList = sieveOfEratosthenes(lastNumber)
list = (2:lastNumber); %Construct list of numbers
primeList = []; %Preallocate prime list
while( list(1)^2 <lastNumber )
primeList = [primeList list(1)]; %add prime to the prime list
list( mod(list,list(1))==0 ) = []; %filter out all multiples of the current prime
end
primeList = [primeList list]; %The rest of the numbers in the list are primes
end</lang>
Sample Output: <lang MATLAB>sieveOfEratosthenes(30)
ans =
2 3 5 7 11 13 17 19 23 29</lang>
Somewhat optimized
<lang MATLAB> function P = erato(x) % Sieve of Eratosthenes: returns all primes between 2 and x
P = [0 2:x] ; % Create vector with all ints between 2 and x where
% position 1 is hard-coded as 0 since 1 is not a prime.
for (n=2:sqrt(x)) % All primes factors lie between 2 and sqrt(x).
if P(n) % If the current value is not 0 (i.e. a prime),
P((2*n):n:x) = 0 ; % then replace all further multiples of it with 0.
end
end % At this point P is a vector with only primes and zeroes.
P = P(P ~= 0) ; % Remove all zeroes from P, leaving only the primes.
return </lang>
The optimization lies in fewer steps in the for loop, use of MATLAB's built-in array operations and no modulo calculation.
Limitation: your machine has to be able to allocate enough memory for an array of length x.
MAXScript
fn eratosthenes n =
(
multiples = #()
print 2
for i in 3 to n do
(
if (findItem multiples i) == 0 then
(
print i
for j in (i * i) to n by i do
(
append multiples j
)
)
)
)
eratosthenes 100
Modula-3
Regular version
<lang modula3> MODULE Prime EXPORTS Main;
IMPORT IO;
CONST LastNum = 1000;
VAR a: ARRAY [2..LastNum] OF BOOLEAN;
BEGIN
FOR i := FIRST(a) TO LAST(a) DO a[i] := TRUE; END;
FOR i := FIRST(a) TO LAST(a) DO
IF a[i] THEN
IO.PutInt(i);
IO.Put(" ");
FOR j := FIRST(a) TO LAST(a) DO
IF j MOD i = 0 THEN
a[j] := FALSE;
END;
END;
END;
END;
IO.Put("\n");
END Prime. </lang>
Advanced version
This version uses more "advanced" types. <lang modula3>(* From the CM3 examples folder (comments removed). *)
MODULE Sieve EXPORTS Main;
IMPORT IO;
TYPE
Number = [2..1000]; Set = SET OF Number;
VAR
prime: Set := Set {FIRST(Number) .. LAST(Number)};
BEGIN
FOR i := FIRST(Number) TO LAST(Number) DO
IF i IN prime THEN
IO.PutInt(i);
IO.Put(" ");
FOR j := i TO LAST(Number) BY i DO
prime := prime - Set{j};
END;
END;
END;
IO.Put("\n");
END Sieve.</lang>
MUMPS
<lang MUMPS>ERATO1(HI)
;performs the Sieve of Erotosethenes up to the number passed in. ;This version sets an array containing the primes SET HI=HI\1 KILL ERATO1 ;Don't make it new - we want it to remain after we quit the function NEW I,J,P FOR I=2:1:(HI**.5)\1 FOR J=I*I:I:HI SET P(J)=1 FOR I=2:1:HI S:'$DATA(P(I)) ERATO1(I)=I KILL I,J,P QUIT</lang>
Example:
USER>SET MAX=100,C=0 DO ERATO1^ROSETTA(MAX) USER>WRITE !,"PRIMES BETWEEN 1 AND ",MAX,! FOR SET I=$ORDER(ERATO1(I)) Q:+I<1 WRITE I,", " PRIMES BETWEEN 1 AND 100 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73,79, 83, 89, 97,
Nial
primes is sublist [ each (2 = sum eachright (0 = mod) [pass,count]), pass ] rest count
Using it
|primes 10 =2 3 5 7
Niue
<lang Niue>[ dup 2 < ] '<2 ; [ 1 + 'count ; [ <2 [ , ] when ] count times ] 'fill-stack ;
0 'n ; 0 'v ;
[ .clr 0 'n ; 0 'v ; ] 'reset ; [ len 1 - n - at 'v ; ] 'set-base ; [ n 1 + 'n ; ] 'incr-n ; [ mod 0 = ] 'is-factor ; [ dup * ] 'sqr ;
[ set-base
v sqr 2 at > not [ [ dup v = not swap v is-factor and ] remove-if incr-n run ] when ] 'run ;
[ fill-stack run ] 'sieve ;
( tests )
10 sieve .s ( => 2 3 5 7 9 ) reset newline 30 sieve .s ( => 2 3 5 7 11 13 17 19 23 29 ) </lang>
Oberon-2
<lang oberon2>MODULE Primes;
IMPORT Out, Math;
CONST N = 1000;
VAR a: ARRAY N OF BOOLEAN;
i, j, m: INTEGER;
BEGIN
(* Set all elements of a to TRUE. *)
FOR i := 1 TO N - 1 DO
a[i] := TRUE;
END;
(* Compute square root of N and convert back to INTEGER. *) m := ENTIER(Math.Sqrt(N));
FOR i := 2 TO m DO
IF a[i] THEN
FOR j := 2 TO (N - 1) DIV i DO
a[i*j] := FALSE;
END;
END;
END;
(* Print all the elements of a that are TRUE. *)
FOR i := 2 TO N - 1 DO
IF a[i] THEN
Out.Int(i, 4);
END;
END;
Out.Ln;
END Primes.</lang>
OCaml
Imperative
<lang ocaml>let sieve n =
let is_prime = Array.create n true in
let limit = truncate(sqrt (float n)) in
for i = 2 to pred limit do
if is_prime.(i) then
let j = ref (i*i) in
try while true do
is_prime.(!j) <- false;
j := !j + i;
done with _ -> ()
done;
(is_prime)
</lang>
<lang ocaml>let primes n =
let _, primes =
Array.fold_left (fun (i,acc) -> function
| true -> (i+1, i::acc)
| false -> (i+1, acc))
(0, [])
(sieve n)
in
List.tl(List.tl(List.rev primes))
</lang>
in the top-level:
# primes 100 ;; - : int list = [2; 3; 5; 7; 11; 13; 17; 19; 23; 29; 31; 37; 41; 43; 47; 53; 59; 61; 67; 71; 73; 79; 83; 89; 97]
Functional
<lang ocaml>(* first define some iterators *)
- let fold_iter f init a b =
let rec aux acc i =
if i > b
then (acc)
else aux (f acc i) (succ i)
in
aux init a ;;
val fold_iter : ('a -> int -> 'a) -> 'a -> int -> int -> 'a = <fun>
- let fold_step f init a b step =
let rec aux acc i =
if i > b
then (acc)
else aux (f acc i) (i + step)
in
aux init a ;;
val fold_step : ('a -> int -> 'a) -> 'a -> int -> int -> int -> 'a = <fun>
(* remove a given value from a list *)
- let remove li v =
let rec aux acc = function
| hd::tl when hd = v -> (List.rev_append acc tl)
| hd::tl -> aux (hd::acc) tl
| [] -> li
in
aux [] li ;;
val remove : 'a list -> 'a -> 'a list = <fun>
(* the main function *)
- let primes n =
let li =
(* create a list [from 2; ... until n] *)
List.rev(fold_iter (fun acc i -> (i::acc)) [] 2 n)
in
let limit = truncate(sqrt(float n)) in
fold_iter (fun li i ->
if List.mem i li (* test if (i) is prime *)
then (fold_step remove li (i*i) n i)
else li)
li 2 (pred limit)
;;
val primes : int -> int list = <fun>
- primes 200 ;;
- : int list = [2; 3; 5; 7; 11; 13; 17; 19; 23; 29; 31; 37; 41; 43; 47; 53; 59; 61; 67; 71;
73; 79; 83; 89; 97; 101; 103; 107; 109; 113; 127; 131; 137; 139; 149; 151; 157; 163; 167; 173; 179; 181; 191; 193; 197; 199]
</lang>
Another functional version
This uses zero to denote struck-out numbers. It is slightly inefficient as it strikes-out multiples above p rather than 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>Eratosthenes(lim)={
my(v=Vectorsmall(lim,unused,1));
forprime(p=2,sqrt(lim),
forstep(i=p^2,lim,p,
v[i]=0
)
);
for(i=1,lim,if(v[i],print1(i", ")))
};</lang>
Pascal
Note: Some Pascal implementations put quite low limits on the size of a set (e.g. Turbo Pascal doesn't allow more than 256 members). To compile on such an implementation, reduce the constant PrimeLimit accordingly. <lang pascal> program primes(output)
const
PrimeLimit = 1000;
var
primes: set of 1 .. PrimeLimit; n, k: integer; needcomma: boolean;
begin
{ calculate the primes }
primes := [2 .. PrimeLimit];
for n := 1 to trunc(sqrt(PrimeLimit)) do
begin
if n in primes
then
begin
k := n*n;
while k < PrimeLimit do
begin
primes := primes - [k];
k := k + n
end
end
end;
{ output the primes }
needcomma := false;
for n := 1 to PrimeLimit do
if n in primes
then
begin
if needcomma
then
write(', ');
write(n);
needcomma := true
end
end. </lang>
Perl
This is not a Sieve of Eratosthenes, because it simply filters the numbers. <lang perl>sub findprimes {
my $prime_limit=shift;
my @possible_primes;
my @primes;
my $prime;
if ($prime_limit < 2) {
return @primes;
}
@possible_primes=(2..$prime_limit);
do {
$prime=shift(@possible_primes);
push @primes, $prime;
@possible_primes = grep {$_%$prime!=0} @possible_primes;
} while ($possible_primes[0] < sqrt($prime_limit));
return (@primes,@possible_primes);
}
my $prime_limit; print "Input the number that you wish to test :"; chomp($prime_limit=<STDIN>); my @primes = findprimes($prime_limit); print "The primes are: @primes!!!\n"; </lang>
Version using bit strings: <lang perl>sub sieve {
my ($n) = @_;
my @primes; my $nonPrimes = ;
foreach my $p (2 .. $n) {
unless (vec($nonPrimes, $p, 1)) {
for (my $i = $p * $p; $i <= $n; $i += $p) {
vec($nonPrimes, $i, 1) = 1;
}
push @primes, $p;
}
}
@primes
}
my $prime_limit; print "Input the number that you wish to test :"; chomp($prime_limit = <STDIN>); my @primes = sieve($prime_limit); print "The primes are: @primes!!!\n";</lang>
Golfing a bit: <lang perl>sub sieve{ my (@s, $i); grep { not $s[ $i = $_ ] and do { $s[ $i += $_ ]++ while $i <= $_[0]; 1 } } 2..$_[0] }
print join ", " => sieve 100; </lang>
Or with bit strings: <lang perl>sub sieve{ my ($s, $i); grep { not vec $s, $i = $_, 1 and do { (vec $s, $i += $_, 1) = 1 while $i <= $_[0]; 1 } } 2..$_[0] }
print join ", " => sieve 100; </lang>
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>
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>
Prolog
Works with SWI-Prolog. <lang Prolog>sieve(N, Ls) :-
numlist(2, N, Ns),
sieve_(Ns, Ls).
divisor(D, N) :- N mod D =:= 0.
sieve_([], []). sieve_([D|Ls0], [D|Ls]) :-
exclude(divisor(D), Ls0, Ls1),
sieve_(Ls1, Ls).
</lang>
PureBasic
Basic procedure
<lang PureBasic>For n=2 To Sqr(lim)
If Nums(n)=0
m=n*n
While m<=lim
Nums(m)=1
m+n
Wend
EndIf
Next n</lang>
Working example
<lang PureBasic>Dim Nums.i(0) Define l, n, m, lim
If OpenConsole()
; Ask for the limit to search, get that input and allocate a Array
Print("Enter limit for this search: ")
lim=Val(Input())
ReDim Nums(lim)
; Use a basic Sieve of Eratosthenes
For n=2 To Sqr(lim)
If Nums(n)=#False
m=n*n
While m<=lim
Nums(m)=#True
m+n
Wend
EndIf
Next n
;Present the result to our user
PrintN(#CRLF$+"The Prims up to "+Str(lim)+" are;")
m=0: l=Log10(lim)+1
For n=2 To lim
If Nums(n)=#False
Print(RSet(Str(n),l)+" ")
m+1
If m>72/(l+1)
m=0: PrintN("")
EndIf
EndIf
Next
Print(#CRLF$+#CRLF$+"Press ENTER to exit"): Input()
CloseConsole()
EndIf</lang>
Output may look like;
Enter limit for this search: 750 The Prims up to 750 are; 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691 701 709 719 727 733 739 743 Press ENTER to exit
Python
Testing on list of primes up to that number
This is not a Sieve of Eratosthenes, it is just trial division. <lang python>def primes_trial(n):
primes = [2]
for x in range(2, n+1):
if all(x % p for p in primes):
primes.append(x)
return primes
print(primes_trial(100))</lang>
Using list of primes up to that number using list comprehension
This is not a Sieve of Eratosthenes, because it simply filters the numbers. <lang python>def GetPrimes(p): i = 1 while i**2 < p[-1]: p[i:] = [x for x in p[i:] if x % p[i-1] != 0] i += 1 return p
print GetPrimes(range(2,100))</lang>
Using list lookup
<lang python>def eratosthenes(n):
multiples = []
for i in xrange(2, n+1):
if i not in multiples:
print i
for j in xrange(i*i, n+1, i):
multiples.append(j)
eratosthenes(100)</lang>
Using set lookup
The version below uses a set to store the multiples. set objects are much faster (usually O(log n)) than lists (O(n)) for checking if a given object is a member. Using the set.update method avoids explicit iteration in the interpreter, giving a further speed improvement.
<lang python>def eratosthenes2(n):
multiples = set()
for i in xrange(2, n+1):
if i not in multiples:
print i
multiples.update(xrange(i*i, n+1, i))
eratosthenes2(100)</lang>
Using array lookup
The version below uses array lookup to test for primality. The function primes_upto() is a straightforward implementation of Sieve of Eratosthenesalgorithm. It returns prime numbers less than or equal to limit. <lang python>def primes_upto(limit):
is_prime = [False] * 2 + [True] * (limit - 1)
for n in xrange(int(limit**0.5 + 1.5)): # stop at ``sqrt(limit)``
if is_prime[n]:
for i in xrange(n*n, limit+1, n):
is_prime[i] = False
return [i for i, prime in enumerate(is_prime) if prime]</lang>
Using generator
<lang python> def iprimes_upto(limit):
is_prime = [False] * 2 + [True] * (limit - 1)
for n in range(limit + 1):
if is_prime[n]:
yield n
for i in range(n*n, limit+1, n): # start at ``n`` squared
is_prime[i] = False
</lang> Example: <lang python> >>> list(iprimes_upto(15)) [2, 3, 5, 7, 11, 13] </lang>
Using numpy
Belowed code adapted from literateprograms.org using numpy <lang python>import numpy def primes_upto2(limit):
is_prime = numpy.ones(limit + 1, dtype=numpy.bool)
for n in xrange(2, int(limit**0.5 + 1.5)):
if is_prime[n]:
is_prime[n*n::n] = 0
return numpy.nonzero(is_prime)[0][2:]
</lang> Performance note: there is no point to add wheels here, due to execution of p[n*n::n] = 0 and nonzero() takes us almost all time.
Using wheels
Version with wheel based optimization: <lang python>from numpy import array, bool_, multiply, nonzero, ones, put, resize
def makepattern(smallprimes):
pattern = ones(multiply.reduce(smallprimes), dtype=bool_)
pattern[0] = 0
for p in smallprimes:
pattern[p::p] = 0
return pattern
def primes_upto3(limit, smallprimes=(2,3,5,7,11)):
sp = array(smallprimes)
if limit <= sp.max(): return sp[sp <= limit]
#
isprime = resize(makepattern(sp), limit + 1)
isprime[:2] = 0; put(isprime, sp, 1)
#
for n in range(sp.max() + 2, int(limit**0.5 + 1.5), 2):
if isprime[n]:
isprime[n*n::n] = 0
return nonzero(isprime)[0]
</lang>
Examples: <lang python> >>> primes_upto3(10**6, smallprimes=(2,3)) # Wall time: 0.17 array([ 2, 3, 5, ..., 999961, 999979, 999983]) >>> primes_upto3(10**7, smallprimes=(2,3)) # Wall time: 2.13 array([ 2, 3, 5, ..., 9999971, 9999973, 9999991]) >>> primes_upto3(15) array([ 2, 3, 5, 7, 11, 13]) >>> primes_upto3(10**7, smallprimes=primes_upto3(15)) # Wall time: 1.31 array([ 2, 3, 5, ..., 9999971, 9999973, 9999991]) >>> primes_upto2(10**7) # Wall time: 1.39 <-- version without wheels array([ 2, 3, 5, ..., 9999971, 9999973, 9999991]) >>> primes_upto3(10**7) # Wall time: 1.30 array([ 2, 3, 5, ..., 9999971, 9999973, 9999991]) </lang> The above-mentioned examples demonstrate that the given wheel based optimization does not show significant performance gain.
Infinite generator
A generator that will generate primes indefinitely (perhaps until it runs out of memory).
<lang python>import heapq
- generates all prime numbers
def sieve():
# priority queue of the sequences of non-primes
# the priority queue allows us to get the "next" non-prime quickly
nonprimes = []
i = 2
while True:
if nonprimes and i == nonprimes[0][0]: # non-prime
while nonprimes[0][0] == i:
# for each sequence that generates this number,
# have it go to the next number (simply add the prime)
# and re-position it in the priority queue
x = nonprimes[0]
x[0] += x[1]
heapq.heapreplace(nonprimes, x)
else: # prime
# insert a 2-element list into the priority queue:
# [current multiple, prime]
# the first element allows sorting by value of current multiple
# we start with i^2
heapq.heappush(nonprimes, [i*i, i])
yield i
i += 1</lang>
Example:
>>> foo = sieve() >>> for i in range(20): ... print(next(foo)) ... 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71
Infinite generator with a faster algorithm
<lang python> def primes():
composites = {} # A mapping from composite numbers to the smallest prime
# that is a factor of it (its witness).
n = 2 # The current number being considered as a prime.
while True:
if n not in composites:
yield n # Not a composite, therefore prime.
composites[n**2] = n # The next unseen composite number is n squared.
else:
# n is composite. Find the next unseen composite number with the
# same witness as n.
witness = composites.pop(n)
next = n + witness
while next in composites:
next += witness
composites[next] = witness
n += 1
import itertools def primes_up_to(limit):
return list(itertools.takewhile(lambda p: p <= limit, primes()))
</lang>
R
This code is vectorised, so it runs quickly for n < one million. The allocation of the primes vector means memory usage is too high for n much larger than that. <lang R> sieve <- function(n) {
n <- as.integer(n)
if(n > 1e6) stop("n too large")
primes <- rep(TRUE, n)
primes[1] <- FALSE
last.prime <- 2L
for(i in last.prime:floor(sqrt(n)))
{
primes[seq.int(2L*last.prime, n, last.prime)] <- FALSE
last.prime <- last.prime + min(which(primes[(last.prime+1):n]))
}
which(primes)
}
sieve(1000) </lang>
REXX
<lang REXX> /*REXX program, generate primes via sieve of Eratosthenes algorithm*/
highItem=200 /*define ten score of #'s. */ widthH=length(highItem) /*width of biggest element#.*/ a.=1 /*assume all numbers prime. */ a.1=0 /*special case for unity, by*/
/*definition, 1 isn't prime.*/
do j=2 to highItem /*starting @ 2, strike mult.*/
if j*j>highItem then leave /*process up to û(highItem).*/
if \a.j then iterate /*this number isn't prime, */
/*it was marked as not prime*/
do k=j+j to highItem by j /*start with next multiple. */ a.k=0 /*mark J multiple not prime.*/ end /*k*/
end /*j*/
np=0 /*number of primes so far. */
do q=1 for highItem /*list all the primes found.*/
if a.q then do
np=np+1 /*bump the prime count. */
if np>1000 then iterate /*only list up to 1k primes.*/
say 'prime number' right(np,widthH) "is" right(q,widthH)
if np==1000 then say 'Prime number listing truncated.'
end
end /*q*/
say np 'prime numbers found up to' highItem </lang> Output:
prime number 1 is 2 prime number 2 is 3 prime number 3 is 5 prime number 4 is 7 prime number 5 is 11 prime number 6 is 13 prime number 7 is 17 prime number 8 is 19 prime number 9 is 23 prime number 10 is 29 prime number 11 is 31 prime number 12 is 37 prime number 13 is 41 prime number 14 is 43 prime number 15 is 47 prime number 16 is 53 prime number 17 is 59 prime number 18 is 61 prime number 19 is 67 prime number 20 is 71 prime number 21 is 73 prime number 22 is 79 prime number 23 is 83 prime number 24 is 89 prime number 25 is 97 prime number 26 is 101 prime number 27 is 103 prime number 28 is 107 prime number 29 is 109 prime number 30 is 113 prime number 31 is 127 prime number 32 is 131 prime number 33 is 137 prime number 34 is 139 prime number 35 is 149 prime number 36 is 151 prime number 37 is 157 prime number 38 is 163 prime number 39 is 167 prime number 40 is 173 prime number 41 is 179 prime number 42 is 181 prime number 43 is 191 prime number 44 is 193 prime number 45 is 197 prime number 46 is 199 46 prime numbers found up to 200
Ruby
<lang ruby>def eratosthenes(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 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]
Scala
<lang scala>def sieveOfEratosthenes(n: Int) = {
val primes = new scala.collection.mutable.BitSet(n)
primes ++= (2 to n)
val sqrt = Math.sqrt(n).toInt
for {
candidate <- 2 to sqrt
if primes contains candidate
} primes --= candidate * candidate to n by candidate
primes
}
println( sieveOfEratosthenes(30) )</lang>
BitSet(2, 3, 5, 7, 11, 13, 17, 19, 23, 29)
Using Streams (the "unfaithful sieve", see the link under the Haskell entry): <lang scala>def ints(n:Int):Stream[Int] = Stream.cons(n,ints(n+1)) def sieve(nums:Stream[Int]):Stream[Int] =
Stream.cons(nums.head, sieve((nums.tail) filter (_%nums.head != 0)))
def primes = sieve(ints(2))
println( primes take 10 toList ) println( primes takeWhile (_<30) toList ) </lang> Both println statements give the same results:
List(2, 3, 5, 7, 11, 13, 17, 19, 23, 29)
Scheme
<lang scheme>(define (iota start stop stride)
(if (> start stop)
(list)
(cons start (iota (+ start stride) stop stride))))
(define (strike lst start stride)
(cond ((null? lst) lst)
((= (car lst) start) (strike (cdr lst) (+ start stride) stride))
((> (car lst) start) (strike lst (+ start stride) stride))
(else (cons (car lst) (strike (cdr lst) start stride)))))
(define (primes limit)
(let ((stop (sqrt limit)))
(define (*primes lst)
(let ((prime (car lst)))
(if (> prime stop)
lst
(cons prime (*primes (strike (cdr lst) (* prime prime) prime))))))
(*primes (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: <lang scheme>(define (primes-wheel-2 limit)
(let ((stop (sqrt limit)))
(define (*primes lst)
(let ((prime (car lst)))
(if (> prime stop)
lst
(cons prime
(*primes (strike (cdr lst) (* prime prime) (* 2 prime)))))))
(cons 2 (*primes (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
<lang seed7>const func array boolean: eratosthenes (in integer: n) is func
result
var array boolean: sieve is 0 times FALSE;
local
var integer: i is 0;
var integer: j is 0;
begin
sieve := n times TRUE;
sieve[1] := FALSE;
for i range 2 to sqrt(n) do
if sieve[i] then
for j range i * i to n step i do
sieve[j] := FALSE;
end for;
end if;
end for;
end func;</lang>
Original source: [1]
SNOBOL4
Using strings instead of arrays, and the square/sqrt optimizations.
<lang SNOBOL4> define('sieve(n)i,j,k,p,str,res') :(sieve_end) sieve i = lt(i,n - 1) i + 1 :f(sv1)
str = str (i + 1) ' ' :(sieve)
sv1 str break(' ') . j span(' ') = :f(return)
sieve = sieve j ' '
sieve = gt(j ^ 2,n) sieve str :s(return) ;* Opt1
res =
str (arb ' ') @p ((j ^ 2) ' ') ;* Opt2
str len(p) . res = ;* Opt2
sv2 str break(' ') . k span(' ') = :f(sv3)
res = ne(remdr(k,j),0) res k ' ' :(sv2)
sv3 str = res :(sv1) sieve_end
- # Test and display
output = sieve(100)
end</lang>
Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
Tcl
<lang tcl>package require Tcl 8.5
proc sieve n {
if {$n < 2} {return {}}
# create a container to hold the sequence of numbers.
# use a dictionary for its speedy access (like an associative array)
# and for its insertion order preservation (like a list)
set nums [dict create]
for {set i 2} {$i <= $n} {incr i} {
# the actual value is never used
dict set nums $i ""
}
set primes [list]
while {[set nextPrime [lindex [dict keys $nums] 0]] <= sqrt($n)} {
dict unset nums $nextPrime
for {set i [expr {$nextPrime ** 2}]} {$i <= $n} {incr i $nextPrime} {
dict unset nums $i
}
lappend primes $nextPrime
}
return [concat $primes [dict keys $nums]]
}
puts [sieve 100] ;# 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97</lang>
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
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>
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
<lang vb>Dim n As Integer, k As Integer, limit As Integer Console.WriteLine("Enter number to search to: ") limit = Console.ReadLine Dim flags(limit) As Integer For n As Integer = 2 To Math.Sqrt(limit)
If flags(n) = 0 Then
For k As Integer = n * n To limit Step n
flags(k) = 1
Next k
End If
Next n
' Display the primes For n = 2 To limit
If flags(n) = 0 Then
Console.WriteLine(n)
End If
Next n</lang>
Vorpal
<lang vorpal> self.print_primes = method(m){
p = new() p.fill(0, m, 1, true)
count = 0
i = 2
while(i < m){
if(p[i] == true){
p.fill(i+i, m, i, false)
count = count + 1
}
i = i + 1
}
('primes: ' + count + ' in ' + m).print()
for(i = 2, i < m, i = i + 1){
if(p[i] == true){
( + i + ', ').put()
}
}
.print()
}
self.print_primes(100) </lang> Result:
primes: 25 in 100 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,
- Programming Tasks
- Prime Numbers
- Clarified and Needing Review
- Ada
- ALGOL 68
- AutoHotkey
- AWK
- BASIC
- Applesoft BASIC
- BBC BASIC
- Befunge
- C
- C++
- C sharp
- Clojure
- Common Lisp
- D
- Dylan
- E
- EC
- Eiffel
- F Sharp
- Forth
- Fortran
- GAP
- Go
- Groovy
- Haskell
- HicEst
- Icon
- Unicon
- J
- Java
- JavaScript
- Liberty BASIC
- Logo
- Lua
- Lucid
- M4
- MATLAB
- MAXScript
- Modula-3
- MUMPS
- Nial
- Niue
- Oberon-2
- OCaml
- Oz
- PARI/GP
- Pascal
- Perl
- Perl 6
- PHP
- PicoLisp
- PL/I
- Pop11
- PowerShell
- Prolog
- PureBasic
- Python
- Numpy
- R
- REXX
- Ruby
- Scala
- Scheme
- Seed7
- SNOBOL4
- Tcl
- UnixPipes
- Ursala
- Vedit macro language
- Visual Basic
- Vorpal