Sieve of Eratosthenes
From Rosetta Code
You are encouraged to solve this task according to the task description, using any language you may know.
The Sieve of Eratosthenes is a simple algorithm that finds the prime numbers up to a given integer. Implement this algorithm, with the only allowed optimization that the outer loop can stop at the square root of the limit, and the inner loop may start at the square of the prime just found. That means especially that you shouldn't optimize by using pre-computed wheels, i.e. don't assume you need only to cross out odd numbers (wheel based on 2), numbers equal to 1 or 5 modulo 6 (wheel based on 2 and 3), or similar wheels based on low primes.
If there's an easy way to add such a wheel based optimization, implement this as an alternative version.
- Note
- It is important that the sieve algorithm be the actual algorithm used to find prime numbers for the task.
- Cf
[edit] ACL2
(defun nats-to-from (n i)
(declare (xargs :measure (nfix (- n i))))
(if (zp (- n i))
nil
(cons i (nats-to-from n (+ i 1)))))
(defun remove-multiples-up-to-r (factor limit xs i)
(declare (xargs :measure (nfix (- limit i))))
(if (or (> i limit)
(zp (- limit i))
(zp factor))
xs
(remove-multiples-up-to-r
factor
limit
(remove i xs)
(+ i factor))))
(defun remove-multiples-up-to (factor limit xs)
(remove-multiples-up-to-r factor limit xs (* factor 2)))
(defun sieve-r (factor limit)
(declare (xargs :measure (nfix (- limit factor))))
(if (zp (- limit factor))
(nats-to-from limit 2)
(remove-multiples-up-to factor (+ limit 1)
(sieve-r (1+ factor) limit))))
(defun sieve (limit)
(sieve-r 2 limit))
[edit] Bracmat
This solution does not use an array. Instead, numbers themselves are used as variables. The numbers that are not prime are set (to the silly value "nonprime"). Finally all numbers up to the limit are tested for being initialised. The uninitialised (unset) ones must be the primes.
( ( eratosthenes
= n j i
. !arg:?n
& 1:?i
& whl
' ( (1+!i:?i)^2:~>!n:?j
& ( !!i
| whl
' ( !j:~>!n
& nonprime:?!j
& !j+!i:?j
)
)
)
& 1:?i
& whl
' ( 1+!i:~>!n:?i
& (!!i|put$(!i " "))
)
)
& eratosthenes$100
)
Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
[edit] GLBasic
// Sieve of Eratosthenes (find primes)
// GLBasic implementation
GLOBAL n%, k%, limit%, flags%[]
limit = 100 // search primes up to this number
DIM flags[limit+1] // GLBasic arrays start at 0
FOR n = 2 TO SQR(limit)
IF flags[n] = 0
FOR k = n*n TO limit STEP n
flags[k] = 1
NEXT
ENDIF
NEXT
// Display the primes
FOR n = 2 TO limit
IF flags[n] = 0 THEN STDOUT n + ", "
NEXT
KEYWAIT
[edit] Ada
with Ada.Text_IO, Ada.Command_Line;
procedure Eratos is
Last: Positive := Positive'Value(Ada.Command_Line.Argument(1));
Prime: array(1 .. Last) of Boolean := (1 => False, others => True);
Base: Positive := 2;
Cnt: Positive;
begin
loop
exit when Base * Base > Last;
if Prime(Base) then
Cnt := Base + Base;
loop
exit when Cnt > Last;
Prime(Cnt) := False;
Cnt := Cnt + Base;
end loop;
end if;
Base := Base + 1;
end loop;
Ada.Text_IO.Put("Primes less or equal" & Positive'Image(Last) &" are:");
for Number in Prime'Range loop
if Prime(Number) then
Ada.Text_IO.Put(Positive'Image(Number));
end if;
end loop;
end Eratos;
Output:
> ./eratos 31 Primes less or equal 31 are : 2 3 5 7 11 13 17 19 23 29 31
[edit] ALGOL 68
BOOL prime = TRUE, non prime = FALSE;
PROC eratosthenes = (INT n)[]BOOL:
(
[n]BOOL sieve;
FOR i TO UPB sieve DO sieve[i] := prime OD;
INT m = ENTIER sqrt(n);
sieve[1] := non prime;
FOR i FROM 2 TO m DO
IF sieve[i] = prime THEN
FOR j FROM i*i BY i TO n DO
sieve[j] := non prime
OD
FI
OD;
sieve
);
print((eratosthenes(80),new line))
Output:
FTTFTFTFFFTFTFFFTFTFFFTFFFFFTFTFFFFFTFFFTFTFFFTFFFFFTFFFFFTFTFFFFFTFFFTFTFFFFFTF
[edit] AutoHotkey
Search autohotkey.com: of Eratosthenes
Source: AutoHotkey forum by Laszlo
MsgBox % "12345678901234567890`n" Sieve(20)
Sieve(n) { ; Sieve of Eratosthenes => string of 0|1 chars, 1 at position k: k is prime
Static zero := 48, one := 49 ; Asc("0"), Asc("1")
VarSetCapacity(S,n,one)
NumPut(zero,S,0,"char")
i := 2
Loop % sqrt(n)-1 {
If (NumGet(S,i-1,"char") = one)
Loop % n//i
If (A_Index > 1)
NumPut(zero,S,A_Index*i-1,"char")
i += 1+(i>2)
}
Return S
}
[edit] AutoIt
#include <Array.au3>
$M = InputBox("Integer", "Enter biggest Integer")
Global $a[$M], $r[$M], $c = 1
For $i = 2 To $M -1
If Not $a[$i] Then
$r[$c] = $i
$c += 1
For $k = $i To $M -1 Step $i
$a[$k] = True
Next
EndIf
Next
$r[0] = $c - 1
ReDim $r[$c]
_ArrayDisplay($r)
[edit] AWK
An initial array holds all numbers 2..max (which is entered on stdin); then all products of integers are deleted from it; the remaining are displayed in the unsorted appearance of a hash table.
$ 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
DIM n AS Integer, k AS Integer, limit AS Integer
INPUT "Enter number to search to: "; limit
DIM flags(limit) AS Integer
FOR n = 2 TO SQR(limit)
IF flags(n) = 0 THEN
FOR k = n*n TO limit STEP n
flags(k) = 1
NEXT k
END IF
NEXT n
' Display the primes
FOR n = 2 TO limit
IF flags(n) = 0 THEN PRINT n; ", ";
NEXT n
[edit] Applesoft BASIC
10 INPUT "ENTER NUMBER TO SEARCH TO: ";LIMIT
20 DIM FLAGS(LIMIT)
30 FOR N = 2 TO SQR (LIMIT)
40 IF FLAGS(N) < > 0 GOTO 80
50 FOR K = N * N TO LIMIT STEP N
60 FLAGS(K) = 1
70 NEXT K
80 NEXT N
90 REM DISPLAY THE PRIMES
100 FOR N = 2 TO LIMIT
110 IF FLAGS(N) = 0 THEN PRINT N;", ";
120 NEXT N
[edit] Locomotive Basic
10 DEFINT a-z
20 INPUT "Limit";limit
30 DIM f(limit)
40 FOR n=2 TO SQR(limit)
50 IF f(n)=1 THEN 90
60 FOR k=n*n TO limit STEP n
70 f(k)=1
80 NEXT k
90 NEXT n
100 FOR n=2 TO limit
110 IF f(n)=0 THEN PRINT n;",";
120 NEXT
[edit] ZX Spectrum Basic
10 INPUT "Enter number to search to: ";l
20 DIM p(l)
30 FOR n=2 TO SQR l
40 IF p(n)<>0 THEN NEXT n
50 FOR k=n*n TO l STEP n
60 LET p(k)=1
70 NEXT k
80 NEXT n
90 REM Display the primes
100 FOR n=2 TO l
110 IF p(n)=0 THEN PRINT n;", ";
120 NEXT n
[edit] BBC BASIC
limit% = 100000
DIM sieve% limit%
prime% = 2
WHILE prime%^2 < limit%
FOR I% = prime%*2 TO limit% STEP prime%
sieve%?I% = 1
NEXT
REPEAT prime% += 1 : UNTIL sieve%?prime%=0
ENDWHILE
REM Display the primes:
FOR I% = 1 TO limit%
IF sieve%?I% = 0 PRINT I%;
NEXT
[edit] bash
See solutions at UNIX Shell.
[edit] Befunge
2>:3g" "-!v\ g30 < |!`"O":+1_:.:03p>03g+:"O"`| @ ^ p3\" ":< 2 234567890123456789012345678901234567890123456789012345678901234567890123456789
[edit] C
Plain sieve, without any optimizations:
#include <stdlib.h>
#include <math.h>
char* eratosthenes (int n, int *c)
{
char* sieve;
int i, j, m;
if (n < 2) { *c = 0; return NULL; }
*c = n-1; /* primes count */
m = (int) sqrt ((double) n);
/* calloc initializes to zero */
sieve = calloc (n+1, sizeof(char));
sieve[0] = 1;
sieve[1] = 1;
for (i = 2; i <= m; i++) {
if (sieve[i] == 0) {
for (j = i*i; j <= n; j += i) {
if (sieve[j] == 0) {
sieve[j] = 1;
--(*c);
}
}
}
}
return sieve;
}
Possible optimizations include sieving only odd numbers (or more complex wheels), packing the sieve into bits to improve locality (and allow larger sieves), etc.
[edit] C++
// yield all prime numbers less than limit.
template<class UnaryFunction>
void primesupto(int limit, UnaryFunction yield)
{
std::vector<bool> is_prime(limit, true);
const int sqrt_limit = static_cast<int>(std::sqrt(limit));
for (int n = 2; n <= sqrt_limit; ++n)
if (is_prime[n]) {
yield(n);
for (unsigned k = n*n, ulim = static_cast<unsigned>(limit); k < ulim; k += n)
//NOTE: "unsigned" is used to avoid an overflow in `k+=n` for `limit` near INT_MAX
is_prime[k] = false;
}
for (int n = sqrt_limit + 1; n < limit; ++n)
if (is_prime[n])
yield(n);
}
Full program:
/**
$ g++ -I/path/to/boost sieve.cpp -o sieve && sieve 10000000
*/
#include <inttypes.h> // uintmax_t
#include <limits>
#include <cmath>
#include <iostream>
#include <sstream>
#include <vector>
#include <boost/lambda/lambda.hpp>
int main(int argc, char *argv[])
{
using namespace std;
using namespace boost::lambda;
int limit = 10000;
if (argc == 2) {
stringstream ss(argv[--argc]);
ss >> limit;
if (limit < 1 or ss.fail()) {
cerr << "USAGE:\n sieve LIMIT\n\nwhere LIMIT in the range [1, "
<< numeric_limits<int>::max() << ")" << endl;
return 2;
}
}
// print primes less then 100
primesupto(100, cout << _1 << " ");
cout << endl;
// find number of primes less then limit and their sum
int count = 0;
uintmax_t sum = 0;
primesupto(limit, (var(sum) += _1, var(count) += 1));
cout << "limit sum pi(n)\n"
<< limit << " " << sum << " " << count << endl;
}
[edit] C#
using System;
using System.Collections;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace sieve
{
class Program
{
static void Main(string[] args)
{
int maxprime = int.Parse(args[0]);
ArrayList primelist = sieve(maxprime);
foreach (int prime in primelist)
Console.WriteLine(prime);
Console.WriteLine("Count = " + primelist.Count);
Console.ReadLine();
}
static ArrayList sieve(int arg_max_prime)
{
BitArray al = new BitArray(arg_max_prime, true);
int lastprime = 1;
int lastprimeSquare = 1;
while (lastprimeSquare <= arg_max_prime)
{
lastprime++;
while (!(bool)al[lastprime])
lastprime++;
lastprimeSquare = lastprime * lastprime;
for (int i = lastprimeSquare; i < arg_max_prime; i += lastprime)
if (i > 0)
al[i] = false;
}
ArrayList sieve_2_return = new ArrayList();
for (int i = 2; i < arg_max_prime; i++)
if (al[i])
sieve_2_return.Add(i);
return sieve_2_return;
}
}
}
[edit] Clojure
Calculates primes up to but not including n.
(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)))))))
This code is not very idiomatic, looks very much like a procedural language. Here is a non-lazy version:
(defn primes [max-prime]
(let [new-sieve (fn []
(let [v (reduce conj (vector-of :boolean) (range (inc max-prime)))]
(assoc (assoc v 0 false) 1 false))),
sieve (fn [S]
(loop [s S, n 2]
(let [n**2 (* n n)]
(if (<= n**2 max-prime)
(recur (if (s n)
(reduce #(assoc %1 %2 false) s (range n**2 (inc max-prime) n))
s)
(inc n))
s)))),
bitset-contents (fn [S]
(map second (filter first (map vector S (range (count S))))))]
(bitset-contents (sieve (new-sieve)))))
Or as a lazy sequence:
(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))))
[edit] CMake
function(eratosthenes var limit)
# Check for integer overflow. With CMake using 32-bit signed integer,
# this check fails when limit > 46340.
if(NOT limit EQUAL 0) # Avoid division by zero.
math(EXPR i "(${limit} * ${limit}) / ${limit}")
if(NOT limit EQUAL ${i})
message(FATAL_ERROR "limit is too large, would cause integer overflow")
endif()
endif()
# Use local variables prime_2, prime_3, ..., prime_${limit} as array.
# Initialize array to y => yes it is prime.
foreach(i RANGE 2 ${limit})
set(prime_${i} y)
endforeach(i)
# Gather a list of prime numbers.
set(list)
foreach(i RANGE 2 ${limit})
if(prime_${i})
# Append this prime to list.
list(APPEND list ${i})
# For each multiple of i, set n => no it is not prime.
# Optimization: start at i squared.
math(EXPR square "${i} * ${i}")
if(NOT square GREATER ${limit}) # Avoid fatal error.
foreach(m RANGE ${square} ${limit} ${i})
set(prime_${m} n)
endforeach(m)
endif()
endif(prime_${i})
endforeach(i)
set(${var} ${list} PARENT_SCOPE)
endfunction(eratosthenes)
# Print all prime numbers through 100.
eratosthenes(primes 100)
message(STATUS "${primes}")
[edit] Common Lisp
(defun sieve-of-eratosthenes (maximum)
(let ((sieve (make-array (1+ maximum) :element-type 'bit
:initial-element 0)))
(loop for candidate from 2 to maximum
when (zerop (bit sieve candidate))
collect candidate
and do (loop for composite from (expt candidate 2)
to maximum by candidate
do (setf (bit sieve composite) 1)))))
Working with odds only (above twice speedup), and only test divide up to the square root of the maximum:
(defun sieve-odds (maximum) "sieve for odd numbers"
(cons 2
(let ((maxi (ash (1- maximum) -1)) (stop (ash (isqrt maximum) -1)))
(let ((sieve (make-array (1+ maxi) :element-type 'bit :initial-element 0)))
(loop for i from 1 to maxi
when (zerop (sbit sieve i))
collect (1+ (ash i 1))
and when (<= i stop) do
(loop for j from (ash (* i (1+ i)) 1) to maxi by (1+ (ash i 1))
do (setf (sbit sieve j) 1)))))))
While formally a wheel, odds are uniformly spaced and do not require any special processing except for value translation. Wheels proper aren't uniformly spaced and are thus trickier.
[edit] D
[edit] Simpler Version
Prints all numbers less than the limit.
import std.stdio, std.algorithm, std.range, std.functional;
uint[] sieve(in uint limit) {
if (limit < 2)
return [];
auto composite = new bool[limit];
foreach (immutable uint n; 2 .. cast(uint)(limit ^^ 0.5) + 1)
if (!composite[n])
for (uint k = n * n; k < limit; k += n)
composite[k] = true;
//return iota(2, limit).filter!(not!composite)().array();
return iota(2, limit).filter!(i => !composite[i])().array();
}
void main() {
writeln(sieve(50));
}
- Output:
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
[edit] Faster Version
This version uses an array of bits (instead of booleans, that are represented with one byte), and skips even numbers. The output is the same.
import std.stdio, std.math, std.array;
size_t[] sieve(in size_t m) /*pure nothrow*/ {
if (m < 3)
return null;
immutable size_t n = m - 1;
enum size_t bpc = size_t.sizeof * 8;
auto F = new size_t[((n + 2) / 2) / bpc + 1];
F[] = size_t.max;
size_t isSet(in size_t i) nothrow {
immutable size_t offset = i / bpc;
immutable size_t mask = 1 << (i % bpc);
return F[offset] & mask;
}
void resetBit(in size_t i) nothrow {
immutable size_t offset = i / bpc;
immutable size_t mask = 1 << (i % bpc);
if ((F[offset] & mask) != 0)
F[offset] = F[offset] ^ mask;
}
for (size_t i = 3; i <= sqrt(cast(real)n); i += 2)
if (isSet((i - 3) / 2))
for (size_t j = i * i; j <= n; j += 2 * i)
resetBit((j - 3) / 2);
Appender!(size_t[]) result;
result ~= 2;
for (size_t i = 3; i <= n; i += 2)
if (isSet((i - 3) / 2))
result ~= i;
return result.data;
}
void main() {
writeln(sieve(50));
}
[edit] Extensible Version
The output is the same.
import std.algorithm, std.bitmanip;
/// Extensible Sieve of Eratosthenes.
struct IsPrime {
//__gshared static BitArray multiples = [true, true, false];
__gshared static BitArray multiples;
static this() {
multiples.length = 3;
multiples[0] = true;
multiples[1] = true;
}
static bool opCall(in size_t n) /*nothrow*/ {
if (n >= multiples.length) {
// Extend the sieve.
immutable oldLen = multiples.length; // Not nothrow.
immutable newMax = max((oldLen - 1) * 2, n);
multiples.length = newMax + 1;
foreach (immutable p; 2 .. oldLen)
if (!multiples[p])
for (size_t i = p * ((oldLen + p) / p);
i < newMax + 1; i += p)
multiples[i] = true;
foreach (immutable i; oldLen .. newMax + 1)
if (!multiples[i])
for (size_t j = i * 2; j < newMax + 1; j += i)
multiples[j] = true;
}
return !multiples[n];
}
}
version (eratosthenes3_main) {
import std.stdio, std.range;
void main() {
// iota(50).filter!IsPrime().writeln();
iota(50).filter!(i => IsPrime(i))().writeln();
}
}
[edit] Dart
// helper function to pretty print an Iterable
String iterableToString(Iterable seq) {
String str = "[";
Iterator i = seq.iterator();
while(i.hasNext()) {
str += i.next();
if(i.hasNext()) {
str += ",";
}
}
return str + "]";
}
main() {
int limit = 1000;
Set<int> sieve = new Set<int>();
for(int i = 2; i <= limit; i++) {
sieve.add(i);
}
for(int i = 2; i * i <= limit; i++) {
if(sieve.contains(i)) {
for(int j = i * i; j <= limit; j += i) {
sieve.remove(j);
}
}
}
var sortedValues = new List<int>.from(sieve);
sortedValues.sort((int arg1, int arg2) {
return arg1.compareTo(arg2);
});
print(iterableToString(sortedValues));
Expect.equals(168, sieve.length);
}
Output:
[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997]
[edit] DWScript
function Primes(limit : Integer) : array of Integer;
var
n, k : Integer;
sieve := new Boolean[limit+1];
begin
for n := 2 to Round(Sqrt(limit)) do begin
if not sieve[n] then begin
for k := n*n to limit step n do
sieve[k] := True;
end;
end;
for k:=2 to limit do
if not sieve[k] then
Result.Add(k);
end;
var r := Primes(50);
var i : Integer;
for i:=0 to r.High do
PrintLn(r[i]);
[edit] Dylan
With outer to sqrt and inner to p^2 optimizations:
define method primes(n)
let limit = floor(n ^ 0.5) + 1;
let sieve = make(limited(<simple-vector>, of: <boolean>), size: n + 1, fill: #t);
let last-prime = 2;
while (last-prime < limit)
for (x from last-prime ^ 2 to n by last-prime)
sieve[x] := #f;
end for;
block (found-prime)
for (n from last-prime + 1 below limit)
if (sieve[n] = #f)
last-prime := n;
found-prime()
end;
end;
last-prime := limit;
end block;
end while;
for (x from 2 to n)
if (sieve[x]) format-out("Prime: %d\n", x); end;
end;
end;
[edit] E
E's standard library doesn't have a step-by-N numeric range, so we'll define one, implementing the standard iteration protocol.
def rangeFromBelowBy(start, limit, step) {
return def stepper {
to iterate(f) {
var i := start
while (i < limit) {
f(null, i)
i += step
}
}
}
}
The sieve itself:
def eratosthenes(limit :(int > 2), output) {
def composite := [].asSet().diverge()
for i ? (!composite.contains(i)) in 2..!limit {
output(i)
composite.addAll(rangeFromBelowBy(i ** 2, limit, i))
}
}
Example usage:
? eratosthenes(12, println) # stdout: 2 # 3 # 5 # 7 # 11
[edit] eC
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
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] Erlang
Compile the module and make following call to find prime numbers below 100. sieve_module:find_primes_below(100).
%% Task: Implementation of Sieve of Eratosthenes
%% Author: Abhay Jain
-module(sieve_module).
-export([find_primes_below/1]).
find_primes_below(N) ->
NumList = lists:seq(1, N),
determine_primes(NumList, 1, []).
%% Sieve of Eratosthenes algorithm
determine_primes(NumList, Index, Primes) ->
case next_prime(NumList, Index+1, length(NumList)) of
{Prime, PrimeIndex, NewNumList} ->
NewPrimes = lists:append(Primes, [Prime]),
determine_primes(NewNumList, PrimeIndex, NewPrimes);
_ ->
%% All prime numbers have been calculated
Primes
end.
next_prime(NumList, Index, Length) ->
if Index > Length ->
false;
true ->
case lists:nth(Index, NumList) of
0 ->
next_prime(NumList, Index+1, Length);
Prime ->
NewNumList = lists:map(fun(A) ->
if A > Index andalso A rem Index == 0 -> 0;
true -> A
end
end, NumList),
{Prime, Index, NewNumList}
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] Euphoria
constant limit = 1000
sequence flags,primes
flags = repeat(1, limit)
for i = 2 to sqrt(limit) do
if flags[i] then
for k = i*i to limit by i do
flags[k] = 0
end for
end if
end for
primes = {}
for i = 2 to limit do
if flags[i] = 1 then
primes &= i
end if
end for
? primes
Output:
{2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,
97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,
181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,
277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,
383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,
487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,
601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,
709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,
827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,
947,953,967,971,977,983,991,997}
[edit] F#
[edit] Functional
#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
program sieve
implicit none
integer, parameter :: i_max = 100
integer :: i
logical, dimension (i_max) :: is_prime
is_prime = .true.
is_prime (1) = .false.
do i = 2, int (sqrt (real (i_max)))
if (is_prime (i)) is_prime (i * i : i_max : i) = .false.
end do
do i = 1, i_max
if (is_prime (i)) write (*, '(i0, 1x)', advance = 'no') i
end do
write (*, *)
end program sieve
Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
Optimised using a pre-computed wheel based on 2:
program sieve_wheel_2
implicit none
integer, parameter :: i_max = 100
integer :: i
logical, dimension (i_max) :: is_prime
is_prime = .true.
is_prime (1) = .false.
is_prime (4 : i_max : 2) = .false.
do i = 3, int (sqrt (real (i_max))), 2
if (is_prime (i)) is_prime (i * i : i_max : 2 * i) = .false.
end do
do i = 1, i_max
if (is_prime (i)) write (*, '(i0, 1x)', advance = 'no') i
end do
write (*, *)
end program sieve_wheel_2
Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
[edit] GAP
Eratosthenes := function(n)
local sieve, cur, mul;
sieve := [1 .. n]*0;
sieve[1] := 1;
cur := 1;
while cur <= n do
if sieve[cur] = 0 then
mul := cur*cur;
while mul <= n do
sieve[mul] := 1;
mul := mul + cur;
od;
fi;
cur := cur + 1;
od;
return Filtered([1 .. n], x -> sieve[x] = 0);
end;
Eratosthenes(100);
# [ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 ]
[edit] Go
package main
import (
"fmt"
)
func main() {
const limit = 201 // means sieve numbers < 201
// sieve
c := make([]bool, limit) // c for composite. false means prime candidate
c[1] = true // 1 not considered prime
p := 2
for {
// first allowed optimization: outer loop only goes to sqrt(limit)
p2 := p * p
if p2 >= limit {
break
}
// second allowed optimization: inner loop starts at sqr(p)
for i := p2; i < limit; i += p {
c[i] = true // it's a composite
}
// scan to get next prime for outer loop
for {
p++
if !c[p] {
break
}
}
}
// sieve complete. now print a representation.
for n := 1; n < limit; n++ {
if c[n] {
fmt.Print(" .")
} else {
fmt.Printf("%3d", n)
}
if n%20 == 0 {
fmt.Println("")
}
}
}
Output:
. 2 3 . 5 . 7 . . . 11 . 13 . . . 17 . 19 . . . 23 . . . . . 29 . 31 . . . . . 37 . . . 41 . 43 . . . 47 . . . . . 53 . . . . . 59 . 61 . . . . . 67 . . . 71 . 73 . . . . . 79 . . . 83 . . . . . 89 . . . . . . . 97 . . . 101 .103 . . .107 .109 . . .113 . . . . . . . . . . . . .127 . . .131 . . . . .137 .139 . . . . . . . . .149 .151 . . . . .157 . . . . .163 . . .167 . . . . .173 . . . . .179 . 181 . . . . . . . . .191 .193 . . .197 .199 .
A fairly odd sieve tree method:
package main
import "fmt"
type xint uint64
type xgen func()(xint)
func primes() func()(xint) {
pp, psq := make([]xint, 0), xint(25)
var sieve func(xint, xint)xgen
sieve = func(p, n xint) xgen {
m, next := xint(0), xgen(nil)
return func()(r xint) {
if next == nil {
r = n
if r <= psq {
n += p
return
}
next = sieve(pp[0] * 2, psq) // chain in
pp = pp[1:]
psq = pp[0] * pp[0]
m = next()
}
switch {
case n < m: r, n = n, n + p
case n > m: r, m = m, next()
default: r, n, m = n, n + p, next()
}
return
}
}
f := sieve(6, 9)
n, p := f(), xint(0)
return func()(xint) {
switch {
case p < 2: p = 2
case p < 3: p = 3
default:
for p += 2; p == n; {
p += 2
if p > n {
n = f()
}
}
pp = append(pp, p)
}
return p
}
}
func main() {
for i, p := 0, primes(); i < 100000; i++ {
fmt.Println(p())
}
}
[edit] Groovy
This solution uses a BitSet for compactness and speed, but in Groovy, BitSet has full List semantics. It also uses both the "square root of the boundary" shortcut and the "square of the prime" shortcut.
def sievePrimes = { bound ->
def isPrime = new BitSet(bound)
isPrime[0..1] = false
isPrime[2..<bound] = true
(2..(Math.sqrt(bound))).each { pc ->
if (isPrime[pc]) {
((pc**2)..<bound).step(pc) { isPrime[it] = false }
}
}
(0..<bound).findAll { isPrime[it] }
}
Test:
println sievePrimes(100)
Output:
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
[edit] GW-BASIC
10 INPUT "ENTER NUMBER TO SEARCH TO: ";LIMIT
20 DIM FLAGS(LIMIT)
30 FOR N = 2 TO SQR (LIMIT)
40 IF FLAGS(N) < > 0 GOTO 80
50 FOR K = N * N TO LIMIT STEP N
60 FLAGS(K) = 1
70 NEXT K
80 NEXT N
90 REM DISPLAY THE PRIMES
100 FOR N = 2 TO LIMIT
110 IF FLAGS(N) = 0 THEN PRINT N;", ";
120 NEXT N
[edit] Haskell
[edit] Imperative version, with unboxed arrays
Using mutable array of unboxed Bools indexed by Integers:
import Control.Monad (forM_, when)
import Control.Monad.ST
import Data.Array.ST
import Data.Array.Unboxed
primeSieve :: Integer -> UArray Integer Bool
primeSieve top = runSTUArray $ do
a <- newArray (2,top) True -- :: ST s (STUArray s Integer Bool)
let r = ceiling . sqrt $ fromInteger top
forM_ [2..r] $ \i -> do
ai <- readArray a i
when ai $ do
forM_ [i*i,i*i+i..top] $ \j -> do
writeArray a j False
return a
-- Return primes from sieve as list:
primesTo :: Integer -> [Integer]
primesTo top = [p | (p,True) <- assocs $ primeSieve top]
[edit] Unboxed arrays, on odds only
Mutable array of unboxed Bools indexed by Ints representing odds only:
sieveUO :: Int -> UArray Int Bool
sieveUO top = runSTUArray $ do
let m = (top-1) `div` 2
r = floor . sqrt $ fromIntegral top + 1
sieve <- newArray (1,m) True -- :: ST s (STUArray s Int Bool)
forM_ [1..r `div` 2] $ \i -> do
isPrime <- readArray sieve i -- ((2*i+1)^2-1)`div`2 = 2*i*(i+1)
when isPrime $ do
forM_ [2*i*(i+1), 2*i*(i+2)+1..m] $ \j -> do
writeArray sieve j False
return sieve
primesToUO :: Int -> [Int]
primesToUO top | top > 1 = 2 : [2*p + 1 | (p,True) <- assocs $ sieveUO top]
| otherwise = []
Using plain Ints is faster but less general. This represents odds only in the array and is fastest. Empirical time complexity is O(n1.2) in n primes produced, and improving for bigger n s. Memory consumption is low (array seems to be packed) but growing about linearly with n.
[edit] Immutable arrays
This is the most straightforward:
import Data.Array.Unboxed
primesToA m = sieve 3 (array (3,m) [(i,odd i) | i<-[3..m]] :: UArray Int Bool)
where
sieve p a
| p*p > m = 2 : [i | (i,True) <- assocs a]
| a!p = sieve (p+2) $ a//[(i,False) | i <- [p*p, p*p+2*p..m]]
| otherwise = sieve (p+2) a
Its performance sharply depends on compiler optimizations. Compiled with -O2 flag in the presence of the explicit type signature, it is very fast in producing first few million primes. (//) is an array update operator.
[edit] Basic list-based sieve
Straightforward implementation of the sieve of Ertosthenes in its original bounded form:
primesTo m = 2 : eratos [3,5..m] where
eratos (p : xs) | p*p>m = p : xs
| True = p : eratos (xs `minus` [p*p, p*p+2*p..])
minus a@(x:xs) b@(y:ys) = case compare x y of
LT -> x : minus xs b
EQ -> minus xs ys
GT -> minus a ys
minus a b = a
Its time complexity is similar to that of optimal trial division because of limitations of Haskell linked lists, where (minus a b) takes time proportional to length(union a b) and not (length b) as achieved in imperative setting with direct-access memory. Uses ordered list representation of sets.
[edit] List-based tree-merging incremental sieve
Basic version's (((a-b)-c)- ... ) workflow can be rearranged into (a-(b+(c+ ... ))). This is the idea behind Richard Bird's unbounded code presented in M. O'Neill's article. Linear merging structure can be further replaced with indefinitely deepening to the right tree-like structure (a-(b+((c+d)+( ((e+f)+(g+h)) + ... )))).
This finds primes in gaps between composites, and composites as union of the multiples of primes, merging in a tree-like fashion all the multiples streams, running with similar as priority-queue–based version of M. O'Neill's empirical time complexity of just above O(n1.2) (for producing first few million primes at least), and with low space consumption too (not counting the produced sequence space):
primes :: [Integer] -- unbounded merging idea due to Richard Bird
primes = 2 : g (fix g) -- double staged production idea due to Melissa O'Neill
where -- tree merging idea Dave Bayer / simplified formulation Will Ness
g xs = 3 : gaps 5 (unionAll [[p*p, p*p+2*p..] | p <- xs])
fix g = xs where xs = g xs -- global definition to prevent space leak
gaps k s@(x:xs) | k < x = k:gaps (k+2) s -- [k,k+2..]`minus`s, k<=x !
| otherwise = gaps (k+2) xs
unionAll ((x:xs):t) = x : union xs (unionAll (pairs t))
where
pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t
union a@(x:xs) b@(y:ys) = case compare x y of
LT -> x : union xs b
EQ -> x : union xs ys
GT -> y : union a ys
Here's the test entry on Ideone.com, a comparison of more versions, and a similar code with wheel optimization.
[edit] HicEst
REAL :: N=100, sieve(N)
sieve = $ > 1 ! = 0 1 1 1 1 ...
DO i = 1, N^0.5
IF( sieve(i) ) THEN
DO j = i^2, N, i
sieve(j) = 0
ENDDO
ENDIF
ENDDO
DO i = 1, N
IF( sieve(i) ) WRITE() i
ENDDO
[edit] Icon and Unicon
procedure main()
sieve(100)
end
procedure sieve(n)
local p,i,j
p:=list(n, 1)
every i:=2 to sqrt(n) & j:= i+i to n by i & p[i] == 1
do p[j] := 0
every write(i:=2 to n & p[i] == 1 & i)
end
Alternatively using sets
procedure main()
sieve(100)
end
procedure sieve(n)
primes := set()
every insert(primes,1 to n)
every member(primes,i := 2 to n) do
every delete(primes,i + i to n by i)
delete(primes,1)
every write(!sort(primes))
end
[edit] J
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 .
[edit] Java
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 = [];
if (limit >= 2) {
var nums = new Array(limit-1);
for (var i = 2; i <= limit; i++)
nums[i-2] = i;
var last_prime;
var idx = 0;
while ((last_prime = nums[idx]) <= Math.sqrt(limit)) {
if (last_prime != null)
for (var i = idx + last_prime; i < limit - 1; i += last_prime)
nums[i] = null;
idx++;
}
for (var i = 0; i < nums.length; i++)
if (nums[i] != null)
primes.push(nums[i]);
}
return primes;
}
var primes = erasthenes(100);
if (typeof print == "undefined")
print = (typeof WScript != "undefined") ? WScript.Echo : alert;
print(primes);
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] Liberty BASIC
'Notice that arrays are globally visible to functions.
'The sieve() function uses the flags() array.
'This is a Sieve benchmark adapted from BYTE 1985
' May, page 286
size = 7000
dim flags(7001)
start = time$("ms")
print sieve(size); " primes found."
print "End of iteration. Elapsed time in milliseconds: "; time$("ms")-start
end
function sieve(size)
for i = 0 to size
if flags(i) = 0 then
prime = i + i + 3
k = i + prime
while k <= size
flags(k) = 1
k = k + prime
wend
sieve = sieve + 1
end if
next i
end function
[edit] 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] Mathematica
Eratosthenes[n_] := Module[
{numbers = Range[n]},
Do[If[numbers[[i]] != 0,
Do[numbers[[i j]] = 0, {j, 2, n/i}]],
{i, 2, Sqrt[n]}];
Select[numbers, # > 1 &]];
Eratosthenes[100]
[edit] 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
Sample Output:
sieveOfEratosthenes(30)
ans =
2 3 5 7 11 13 17 19 23 29
[edit] Somewhat optimized
function P = erato(x) % Sieve of Eratosthenes: returns all primes between 2 and x
P = [0 2:x] ; % Create vector with all ints between 2 and x where
% position 1 is hard-coded as 0 since 1 is not a prime.
for (n=2:sqrt(x)) % All primes factors lie between 2 and sqrt(x).
if P(n) % If the current value is not 0 (i.e. a prime),
P((2*n):n:x) = 0 ; % then replace all further multiples of it with 0.
end
end % At this point P is a vector with only primes and zeroes.
P = P(P ~= 0) ; % Remove all zeroes from P, leaving only the primes.
return
The optimization lies in fewer steps in the for loop, use of MATLAB's built-in array operations and no modulo calculation.
Limitation: your machine has to be able to allocate enough memory for an array of length x.
[edit] Maxima
sieve(n) := block([u, j, m],
u: makelist(i, i, 1, n),
u[1]: -1,
j: 1,
do (
j: j + 1,
while u[j] < 0 do j: j + 1,
m: j*j,
if m > n then return(sublist(u, nonnegintegerp)),
for k from m step j while k <= n do u[k]: -1
)
)$
[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] MUMPS
ERATO1(HI)
;performs the Sieve of Erotosethenes up to the number passed in.
;This version sets an array containing the primes
SET HI=HI\1
KILL ERATO1 ;Don't make it new - we want it to remain after we quit the function
NEW I,J,P
FOR I=2:1:(HI**.5)\1 FOR J=I*I:I:HI SET P(J)=1
FOR I=2:1:HI S:'$DATA(P(I)) ERATO1(I)=I
KILL I,J,P
QUIT
Example:
USER>SET MAX=100,C=0 DO ERATO1^ROSETTA(MAX) USER>WRITE !,"PRIMES BETWEEN 1 AND ",MAX,! FOR SET I=$ORDER(ERATO1(I)) Q:+I<1 WRITE I,", " PRIMES BETWEEN 1 AND 100 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73,79, 83, 89, 97,
[edit] NetRexx
[edit] Version 1 (slow)
/* NetRexx */
options replace format comments java crossref savelog symbols binary
parse arg loWatermark hiWatermark .
if loWatermark = '' | loWatermark = '.' then loWatermark = 1
if hiWatermark = '' | hiWatermark = '.' then hiWatermark = 200
do
if \loWatermark.datatype('w') | \hiWatermark.datatype('w') then -
signal NumberFormatException('arguments must be whole numbers')
if loWatermark > hiWatermark then -
signal IllegalArgumentException('the start value must be less than the end value')
seive = sieveOfEratosthenes(hiWatermark)
primes = getPrimes(seive, loWatermark, hiWatermark).strip
say 'List of prime numbers from' loWatermark 'to' hiWatermark 'via a "Sieve of Eratosthenes" algorithm:'
say ' 'primes.changestr(' ', ',')
say ' Count of primes:' primes.words
catch ex = Exception
ex.printStackTrace
end
return
method sieveOfEratosthenes(hn = long) public static binary returns Rexx
sv = Rexx(isTrue)
sv[1] = isFalse
ix = long
jx = long
loop ix = 2 while ix * ix <= hn
if sv[ix] then loop jx = ix * ix by ix while jx <= hn
sv[jx] = isFalse
end jx
end ix
return sv
method getPrimes(seive = Rexx, lo = long, hi = long) private constant binary returns Rexx
primes = Rexx('')
loop p_ = lo to hi
if \seive[p_] then iterate p_
primes = primes p_
end p_
return primes
method isTrue public constant binary returns boolean
return 1 == 1
method isFalse public constant binary returns boolean
return \isTrue
- Output
List of prime numbers from 1 to 200 via a "Sieve of Eratosthenes" algorithm: 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199 Count of primes: 46
[edit] Version 2 (significantly, i.e. 10 times faster)
/* NetRexx ************************************************************
* Essential improvements:Use boolean instead of Rexx for sv
* and remove methods isTrue and isFalse
* 24.07.2012 Walter Pachl courtesy Kermit Kiser
**********************************************************************/
options replace format comments java crossref savelog symbols binary
parse arg loWatermark hiWatermark .
if loWatermark = '' | loWatermark = '.' then loWatermark = 1
if hiWatermark = '' | hiWatermark = '.' then hiWatermark = 200000
startdate=Date Date()
do
if \loWatermark.datatype('w') | \hiWatermark.datatype('w') then -
signal NumberFormatException('arguments must be whole numbers')
if loWatermark > hiWatermark then -
signal IllegalArgumentException(-
'the start value must be less than the end value')
sieve = sieveOfEratosthenes(hiWatermark)
primes = getPrimes(sieve, loWatermark, hiWatermark).strip
if hiWatermark = 200 Then do
say 'List of prime numbers from' loWatermark 'to' hiWatermark
say ' 'primes.changestr(' ', ',')
end
catch ex = Exception
ex.printStackTrace
end
enddate=Date Date()
Numeric Digits 20
say (enddate.getTime-startdate.getTime)/1000 'seconds elapsed'
say ' Count of primes:' primes.words
return
method sieveOfEratosthenes(hn = int) -
public static binary returns boolean[]
true = boolean 1
false = boolean 0
sv = boolean[hn+1]
sv[1] = false
ix = int
jx = int
loop ix=2 to hn
sv[ix]=true
end ix
loop ix = 2 while ix * ix <= hn
if sv[ix] then loop jx = ix * ix by ix while jx <= hn
sv[jx] = false
end jx
end ix
return sv
method getPrimes(sieve = boolean[], lo = int, hi = int) -
private constant binary Returns Rexx
p_ = int
primes = Rexx('')
loop p_ = lo to hi
if \sieve[p_] then iterate p_
primes = primes p_
end p_
return primes
[edit] Nial
primes is sublist [ each (2 = sum eachright (0 = mod) [pass,count]), pass ] rest count
Using it
|primes 10 =2 3 5 7
[edit] 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 )
[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
declare
fun {Sieve N}
S = {Array.new 2 N true}
M = {Float.toInt {Sqrt {Int.toFloat N}}}
in
for I in 2..M do
if S.I then
for J in I*I..N;I do
S.J := false
end
end
end
S
end
fun {Primes N}
S = {Sieve N}
in
for I in 2..N collect:C do
if S.I then {C I} end
end
end
in
{Show {Primes 30}}
[edit] PARI/GP
Eratosthenes(lim)={
my(v=Vectorsmall(lim,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", ")))
};
Here is an alternate version which does not use the isprime command:
Sieve(n)=
{
v=vector(n,unused,1);
for(i=2,sqrt(n),
if(v[i],
forstep(j=i^2,n,i,v[j]=0)));
for(i=2,n,if(v[i],print1(i)))
};
[edit] Pascal
Note: Some Pascal implementations put quite low limits on the size of a set (e.g. Turbo Pascal doesn't allow more than 256 members). To compile on such an implementation, reduce the constant PrimeLimit accordingly.
program primes(output)
const
PrimeLimit = 1000;
var
primes: set of 1 .. PrimeLimit;
n, k: integer;
needcomma: boolean;
begin
{ calculate the primes }
primes := [2 .. PrimeLimit];
for n := 1 to trunc(sqrt(PrimeLimit)) do
begin
if n in primes
then
begin
k := n*n;
while k < PrimeLimit do
begin
primes := primes - [k];
k := k + n
end
end
end;
{ output the primes }
needcomma := false;
for n := 1 to PrimeLimit do
if n in primes
then
begin
if needcomma
then
write(', ');
write(n);
needcomma := true
end
end.
[edit] Perl
$max= $ARGV[0];
@primes= ();
@tested= (1);
$j= 1;
while ($j < $max) {
next if $tested[$j++];
push @primes, $j;
for ($k= $j; $k <= $max; $k+=$j) {
$tested[$k-1]= 1;
}
}
print "@primes\n";
A version similar to the above, but much faster:
sub simple_sieve {
my($max) = @_;
return () if $max < 2; return (2) if $max < 3;
my @c;
for(my $t=3; $t*$t<=$max; $t+=2) {
if (!$c[$t]) {
for(my $s=$t*$t; $s<=$max; $s+=$t*2) { $c[$s]++ }
}
}
my @primes = (2);
for(my $t=3; $t<=$max; $t+=2) {
$c[$t] || push @primes, $t;
}
@primes;
}
Using vectors and optimizing away odds. Not much slower than the above, but using 200x less memory:
sub dj_vector {
my($end) = @_;
return () if $end < 2; return (2) if $end < 3;
$end-- if ($end & 1) == 0; # Ensure end is odd
my ($sieve, $n, $limit, $s_end) = ( '', 3, int(sqrt($end)), $end >> 1 );
while ( $n <= $limit ) {
for (my $s = ($n*$n) >> 1; $s <= $s_end; $s += $n) {
vec($sieve, $s, 1) = 1;
}
do { $n += 2 } while vec($sieve, $n >> 1, 1) != 0;
}
my @primes = (2);
do { push @primes, 2*$_+1 if !vec($sieve,$_,1) } for (1..int(($end-1)/2));
@primes;
}
The fastest non-segmented pure Perl sieve, using strings for low memory. About 2x faster than the fast array version with 5.16.0 or later, about 10x faster than the golfed versions below. Uses much less memory than array versions.
sub dj_string {
my($end) = @_;
return () if $end < 2; return (2) if $end < 3;
$end-- if ($end & 1) == 0;
my $s_end = $end >> 1;
my $whole = int( ($end>>1) / 15); # prefill with 3 and 5 marked
my $sieve = '100010010010110' . '011010010010110' x $whole;
substr($sieve, ($end>>1)+1) = '';
my ($n, $limit) = ( 7, int(sqrt($end)) );
while ( $n <= $limit ) {
for (my $s = ($n*$n) >> 1; $s <= $s_end; $s += $n) {
substr($sieve, $s, 1) = '1';
}
do { $n += 2 } while substr($sieve, $n>>1, 1);
}
# If you just want the count, it's very fast:
# my $count = 1 + $sieve =~ tr/0//;
my @primes = (2, 3, 5);
$n = 7-2;
foreach my $s (split("0", substr($sieve, 3), -1)) {
$n += 2 + 2 * length($s);
push @primes, $n if $n <= $end;
}
@primes;
}
Golfing a bit, at the expense of speed:
sub sieve{ my (@s, $i);
grep { not $s[ $i = $_ ] and do
{ $s[ $i += $_ ]++ while $i <= $_[0]; 1 }
} 2..$_[0]
}
print join ", " => sieve 100;
Or with bit strings (much slower than the vector version above):
sub sieve{ my ($s, $i);
grep { not vec $s, $i = $_, 1 and do
{ (vec $s, $i += $_, 1) = 1 while $i <= $_[0]; 1 }
} 2..$_[0]
}
print join ", " => sieve 100;
A short recursive version:
sub erat {Regexp:
my $p = shift;
return $p, $p**2 > $_[$#_] ? @_ : erat(grep $_%$p, @_)
}
print join ', ' => erat 2..100000;
sub sieve {
my ($s, $p) = "." . ("x" x shift);
1 while ++$p
and $s =~ /^(.{$p,}?)x/g
and $p = length($1)
and $s =~ s/(.{$p})./${1}./g
and substr($s, $p, 1) = "x";
$s
}
print sieve(1000);
[edit] Perl 6
sub sieve( Int $limit ) {
my @is-prime = False, False, True xx $limit - 1;
gather for @is-prime.kv -> $number, $is-prime {
if $is-prime {
take $number;
@is-prime[$_] = False if $_ %% $number for $number**2 .. $limit;
}
}
}
(sieve 100).join(",").say;
A recursive version:
multi erat(Int $N) { erat 2 .. $N }
multi erat(@a where @a[0] > sqrt @a[*-1]) { @a }
multi erat(@a) { @a[0], erat(@a.grep: * % @a[0]) }
Â
say erat 100;
Of course, upper limits are for wusses. Here's a version using infinite streams, that just keeps going until you ^C it (works under Niecza):
role Filter[Int $factor] {
method next { repeat until $.value % $factor { callsame } }
}
class Stream {
has Int $.value is rw = 1;
method next { ++$.value }
method filter { self but Filter[$.value] }
}
.next, say .value for Stream.new, *.filter ... *;
[edit] PHP
function iprimes_upto($limit)
{
for ($i = 2; $i < $limit; $i++)
{
$primes[$i] = true;
}
for ($n = 2; $n < $limit; $n++)
{
if ($primes[$n])
{
for ($i = $n*$n; $i < $limit; $i += $n)
{
$primes[$i] = false;
}
}
}
return $primes;
}
[edit] PicoLisp
(de sieve (N)
(let Sieve (range 1 N)
(set Sieve)
(for I (cdr Sieve)
(when I
(for (S (nth Sieve (* I I)) S (nth (cdr S) I))
(set S) ) ) )
(filter bool Sieve) ) )
Output:
: (sieve 100) -> (2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)
[edit] PL/I
eratos: proc options (main) reorder;
dcl i fixed bin (31);
dcl j fixed bin (31);
dcl n fixed bin (31);
dcl sn fixed bin (31);
dcl hbound builtin;
dcl sqrt builtin;
dcl sysin file;
dcl sysprint file;
get list (n);
sn = sqrt(n);
begin;
dcl primes(n) bit (1) aligned init ((*)((1)'1'b));
i = 2;
do while(i <= sn);
do j = i ** 2 by i to hbound(primes, 1);
/* Adding a test would just slow down processing! */
primes(j) = '0'b;
end;
do i = i + 1 to sn until(primes(i));
end;
end;
do i = 2 to hbound(primes, 1);
if primes(i) then
put data(i);
end;
end;
end eratos;
[edit] Pop11
define eratostenes(n);
lvars bits = inits(n), i, j;
for i from 2 to n do
if bits(i) = 0 then
printf('' >< i, '%s\n');
for j from 2*i by i to n do
1 -> bits(j);
endfor;
endif;
endfor;
enddefine;
[edit] PowerShell
[edit] Basic procedure
It outputs immediately so that the number can be used by the pipeline.
function Sieve ( [int] $num )
{
$isprime = @{}
2..$num | Where-Object {
$isprime[$_] -eq $null } | ForEach-Object {
$_
$isprime[$_] = $true
$i=$_*$_
for ( ; $i -le $num; $i += $_ )
{ $isprime[$i] = $false }
}
}
[edit] Processing
Calculate the primes up to 1000000 with Processing, including a visualisation of the process. As an additional visual effect, the layout of the pixel could be changed from the line-by-line layout to a spiral-like layout starting in the middle of the screen.
int maxx,maxy;
int max;
boolean[] sieve;
void plot(int pos, boolean active) {
set(pos%maxx,pos/maxx, active?#000000:#ffffff);
}
void setup() {
size(1000, 1000, P2D);
frameRate(2);
maxx=width;
maxy=height;
max=width*height;
sieve=new boolean[max+1];
sieve[1]=false;
plot(0,false);
plot(1,false);
for(int i=2;i<=max;i++) {
sieve[i]=true;
plot(i,true);
}
}
int i=2;
void draw() {
if(!sieve[i]) {
while(i*i<max && !sieve[i]) {
i++;
}
}
if(sieve[i]) {
print(i+" ");
for(int j=i*i;j<=max;j+=i) {
if(sieve[j]) {
sieve[j]=false;
plot(j,false);
}
}
}
if(i*i<max) {
i++;
} else {
noLoop();
println("finished");
}
}
[edit] Prolog
[edit] using facts as a set of non-prime numbers
Builds a set of non-prime numbers (multiples of already visited numbers) pretty much the same way as imperative languages, by counting up in increments of p for each prime p (or 2p for odd primes). Works with SWI-Prolog.
sieve(N, [2|PS]) :- % PS is list of odd primes up to N
retractall(mult(_)),
sieve_O(3,N,PS).
sieve_O(I,N,PS) :- % sieve odds from I up to N to get PS
I =< N, !, I1 is I+2,
( mult(I) -> sieve_O(I1,N,PS)
; ( I =< N / I ->
ISq is I*I, DI is 2*I, add_mults(DI,ISq,N)
; true
),
PS = [I|T],
sieve_O(I1,N,T)
).
sieve_O(I,N,[]) :- I > N.
add_mults(DI,I,N) :-
I =< N, !,
( mult(I) -> true ; assert(mult(I)) ),
I1 is I+DI,
add_mults(DI,I1,N).
add_mults(_,I,N) :- I > N.
[edit] Using a priority queue
Uses a ariority queue, from the paper "The Genuine Sieve of Eratosthenes" by Melissa O'Neill. Works with YAP (Yet Another Prolog)
?- use_module(library(heaps)).
prime(2).
prime(N) :- prime_heap(N, _).
prime_heap(3, H) :- list_to_heap([9-6], H).
prime_heap(N, H) :-
prime_heap(M, H0), N0 is M + 2,
next_prime(N0, H0, N, H).
next_prime(N0, H0, N, H) :-
\+ min_of_heap(H0, N0, _),
N = N0, Composite is N*N, Skip is N+N,
add_to_heap(H0, Composite, Skip, H).
next_prime(N0, H0, N, H) :-
min_of_heap(H0, N0, _),
adjust_heap(H0, N0, H1), N1 is N0 + 2,
next_prime(N1, H1, N, H).
adjust_heap(H0, N, H) :-
min_of_heap(H0, N, _),
get_from_heap(H0, N, Skip, H1),
Composite is N + Skip, add_to_heap(H1, Composite, Skip, H2),
adjust_heap(H2, N, H).
adjust_heap(H, N, H) :-
\+ min_of_heap(H, N, _).
[edit] PureBasic
[edit] Basic procedure
For n=2 To Sqr(lim)
If Nums(n)=0
m=n*n
While m<=lim
Nums(m)=1
m+n
Wend
EndIf
Next n
[edit] Working example
Dim Nums.i(0)
Define l, n, m, lim
If OpenConsole()
; Ask for the limit to search, get that input and allocate a Array
Print("Enter limit for this search: ")
lim=Val(Input())
ReDim Nums(lim)
; Use a basic Sieve of Eratosthenes
For n=2 To Sqr(lim)
If Nums(n)=#False
m=n*n
While m<=lim
Nums(m)=#True
m+n
Wend
EndIf
Next n
;Present the result to our user
PrintN(#CRLF$+"The Prims up to "+Str(lim)+" are;")
m=0: l=Log10(lim)+1
For n=2 To lim
If Nums(n)=#False
Print(RSet(Str(n),l)+" ")
m+1
If m>72/(l+1)
m=0: PrintN("")
EndIf
EndIf
Next
Print(#CRLF$+#CRLF$+"Press ENTER to exit"): Input()
CloseConsole()
EndIf
Output may look like;
Enter limit for this search: 750 The Prims up to 750 are; 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691 701 709 719 727 733 739 743 Press ENTER to exit
[edit] Python
[edit] Using list lookup
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)
[edit] Using set lookup
The version below uses a set to store the multiples. set objects are much faster (usually O(log n)) than lists (O(n)) for checking if a given object is a member. Using the set.update method avoids explicit iteration in the interpreter, giving a further speed improvement.
def eratosthenes2(n):
multiples = set()
for i in xrange(2, 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
Below code adapted from literateprograms.org using numpy
import numpy
def primes_upto2(limit):
is_prime = numpy.ones(limit + 1, dtype=numpy.bool)
for n in xrange(2, int(limit**0.5 + 1.5)):
if is_prime[n]:
is_prime[n*n::n] = 0
return numpy.nonzero(is_prime)[0][2:]
Performance note: there is no point to add wheels here, due to execution of p[n*n::n] = 0 and nonzero() takes us almost all time.
[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).
import heapq
# generates all prime numbers
def sieve():
# priority queue of the sequences of non-primes
# the priority queue allows us to get the "next" non-prime quickly
nonprimes = []
i = 2
while True:
if nonprimes and i == nonprimes[0][0]: # non-prime
while nonprimes[0][0] == i:
# for each sequence that generates this number,
# have it go to the next number (simply add the prime)
# and re-position it in the priority queue
x = nonprimes[0]
x[0] += x[1]
heapq.heapreplace(nonprimes, x)
else: # prime
# insert a 2-element list into the priority queue:
# [current multiple, prime]
# the first element allows sorting by value of current multiple
# we start with i^2
heapq.heappush(nonprimes, [i*i, i])
yield i
i += 1
Example:
>>> foo = sieve() >>> for i in range(8): ... print(next(foo)) ... 2 3 5 7 11 13 17 19
[edit] Infinite generator with a faster algorithm
The adding of each discovered prime's incremental step info to the mapping should be postponed until the prime's square is seen amongst the candidate numbers, as it is useless before that point. This drastically reduces the space complexity from O(n) to O(sqrt(n/log(n))), in n primes produced, and also lowers the run time complexity quite low (this test entry shows about O(n^1.08) empirical order of growth which is very close to the theoretical value of O(n log(n) log(log(n))), in n primes produced):
def primes():
yield 2 ; yield 3 ; yield 5 ; yield 7 ;
ps = (p for p in primes()) # additional primes supply
p = ps.next() and ps.next() # discard 2, then get 3
q = p*p # 9 - square of next prime to be put
sieve = {} # into sieve dict
n = 9 # the candidate number
while True:
if n not in sieve: # is not a multiple of previously recorded primes
if n < q: # n is prime
yield n
else:
add(sieve, q + 2*p, 2*p) # n==p*p: for prime p, under p*p + 2*p,
p = ps.next() # add 2*p as incremental step
q = p*p
else:
s = sieve.pop(n)
add(sieve, n + s, s)
n += 2 # work on odds only
def add(sieve,next,step):
while next in sieve: # ensure each entry is unique
next += step
sieve[next] = step # next non-marked multiple of a prime
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] Racket
[edit] Imperative versions
#lang racket
;; ugly imperative version
(define (sieve n)
(define non-primes '())
(define primes '())
(for ([i (in-range 2 (add1 n))])
(unless (member i non-primes)
(set! primes (cons i primes))
(for ([j (in-range (* i i) (add1 n) i)])
(set! non-primes (cons j non-primes)))))
(reverse primes))
(sieve 100)
#lang racket
;; a little nicer, but still imperative
(define (sieve n)
(define primes (make-vector (add1 n) #t))
(for* ([i (in-range 2 (add1 n))]
#:when (vector-ref primes i)
[j (in-range (* i i) (add1 n) i)])
(vector-set! primes j #f))
(for/list ([n (in-range 2 (add1 n))]
#:when (vector-ref primes n))
n))
(sieve 100)
[edit] Infinite list using laziness
#lang lazy
;; infinite list using lazy racket (see the language above)
(define nats (cons 1 (map add1 nats)))
(define (sift n l) (filter (λ(x) (not (zero? (modulo x n)))) l))
(define (sieve l) (cons (first l) (sieve (sift (first l) (rest l)))))
(define primes (sieve (rest nats)))
(define (take-upto n l)
(if (<= (car l) n) (cons (car l) (take-upto n (cdr l))) '()))
(!! (take-upto 100 primes))
[edit] Infinite list using threads and channels
#lang racket
;; infinite list using threads and channels (similar to newsqueak)
(define-syntax (bg-loop stx)
(syntax-case stx ()
[(_ expr ...)
(with-syntax ([out! (datum->syntax stx 'out!)])
#'(let ([out (make-channel)])
(define (out! x) (channel-put out x))
(thread (λ() (let loop () expr ... (loop))))
out))]))
(define nats (bg-loop (for ([i (in-naturals 1)]) (out! i))))
(define (filter pred? c)
(bg-loop (define x (channel-get c))
(when (pred? x) (out! x))))
(define (sift n c)
(filter (λ(x) (not (zero? (modulo x n)))) c))
(define (sieve c)
(bg-loop (define x (channel-get c))
(out! x)
(set! c (sift x c))))
(define primes (begin (channel-get nats) (sieve nats)))
(define (take-upto n c)
(let loop ()
(let ([x (channel-get c)]) (if (<= x n) (cons x (loop)) '()))))
(take-upto 100 primes)
[edit] Infinite list using generators
#lang racket
;; yet another variation of the same algorithm, this time using generators
(require racket/generator)
(define nats (generator () (for ([i (in-naturals 1)]) (yield i))))
(define (filter pred g)
(generator () (for ([i (in-producer g #f)] #:when (pred i)) (yield i))))
(define (sift n g) (filter (λ(x) (not (zero? (modulo x n)))) g))
(define (sieve g)
(generator ()
(let loop ([g g]) (let ([x (g)]) (yield x) (loop (sift x g))))))
(define primes (begin (nats) (sieve nats)))
(define (take-upto n p)
(let loop ()
(let ([x (p)]) (if (<= x n) (cons x (loop)) '()))))
(take-upto 100 primes)
[edit] REXX
[edit] no wheel version
The first three REXX versions makes use of a sparse stemmed array: [@.].
As the stemmed array gets heavily populated, the number of entries may slow down the REXX interpreter
substantially, depending upon the efficacy of the hashing technique being used for REXX variables.
/*REXX program generates primes via the sieve of Eratosthenes algorithm.*/
parse arg H .; if H=='' then H=200 /*let the highest # be specified.*/
w = length(H); prime=right('prime',20) /*W is used for formatting output*/
@.=1 /*assume all numbers are prime. */
do j=2 while j*j <= H /*all integers up to √H inclusive*/
if @.j then do m=j*j to H by J /*Prime? Then strike mutliples.*/
@.m=0 /*mark all multiples of J ¬prime.*/
end /*m*/
end /*j*/
#=0 /*count of primes listed so far. */
do n=2 to H /* [↓] if prime, then display it.*/
if @.n then do; # = #+1; say prime right(#,w) " ───► " right(n,w); end
end /*n*/
/*stick a fork in it, we're done.*/
output when using the input default of 200
prime 1 ───► 2
prime 2 ───► 3
prime 3 ───► 5
prime 4 ───► 7
prime 5 ───► 11
prime 6 ───► 13
prime 7 ───► 17
prime 8 ───► 19
prime 9 ───► 23
prime 10 ───► 29
prime 11 ───► 31
prime 12 ───► 37
prime 13 ───► 41
prime 14 ───► 43
prime 15 ───► 47
prime 16 ───► 53
prime 17 ───► 59
prime 18 ───► 61
prime 19 ───► 67
prime 20 ───► 71
prime 21 ───► 73
prime 22 ───► 79
prime 23 ───► 83
prime 24 ───► 89
prime 25 ───► 97
prime 26 ───► 101
prime 27 ───► 103
prime 28 ───► 107
prime 29 ───► 109
prime 30 ───► 113
prime 31 ───► 127
prime 32 ───► 131
prime 33 ───► 137
prime 34 ───► 139
prime 35 ───► 149
prime 36 ───► 151
prime 37 ───► 157
prime 38 ───► 163
prime 39 ───► 167
prime 40 ───► 173
prime 41 ───► 179
prime 42 ───► 181
prime 43 ───► 191
prime 44 ───► 193
prime 45 ───► 197
prime 46 ───► 199
[edit] wheel version with optional prime list suppression
This version skips striking the even numbers (as being not prime).
Also supported is the suppression of listing the primes if the H (high limit) is negative.
Also added is a final message indicating the number of primes found.
/*REXX pgm gens primes via a wheeled sieve of Eratosthenes algorithm.*/
parse arg H .; if H=='' then H=200 /*high# can be specified on C.L. */
oH=H; H=abs(H); w=length(H) /*neg H suppresses prime listing.*/
if 2<=H & oH>0 then say right(1,w+20)'st prime ───► ' right(2,w)
@.=1 /*assume all numbers are prime. */
do j=3 by 2 while j*j <= H /*odd integers up to √H inclusive*/
if @.j then do m=j*j to H by j+j /*Prime? Then strike multiples.*/
@.m = 0 /*mark odd multiples of J ¬prime.*/
end /*m*/
end /*j*/
#= H>1 /*count of primes found so far. */
do n=3 to H by 2 /*count primes & maybe show them.*/
if @.n then do; #=#+1; if oH>0 then say right(#,w+20)th(#)' prime ───► ' right(n,w); end
end /*n*/
say; say right(#,w+20+2) 'prime's(#) "found."
exit /*stick a fork in it, we're done.*/
/*──────────────────────────────────S subroutine────────────────────────*/
s: if arg(1)==1 then return arg(3); return word(arg(2) 's',1) /*plural*/
/*──────────────────────────────────TH subroutine───────────────────────*/
th: procedure;parse arg x;x=abs(x);return word('th st nd rd',1+x//10*(x//100%10\==1)*(x//10<4))
output when using the input default of 200
1st prime ───► 2
2nd prime ───► 3
3rd prime ───► 5
4th prime ───► 7
5th prime ───► 11
6th prime ───► 13
7th prime ───► 17
8th prime ───► 19
9th prime ───► 23
10th prime ───► 29
11th prime ───► 31
12th prime ───► 37
13th prime ───► 41
14th prime ───► 43
15th prime ───► 47
16th prime ───► 53
17th prime ───► 59
18th prime ───► 61
19th prime ───► 67
20th prime ───► 71
21st prime ───► 73
22nd prime ───► 79
23rd prime ───► 83
24th prime ───► 89
25th prime ───► 97
26th prime ───► 101
27th prime ───► 103
28th prime ───► 107
29th prime ───► 109
30th prime ───► 113
31st prime ───► 127
32nd prime ───► 131
33rd prime ───► 137
34th prime ───► 139
35th prime ───► 149
36th prime ───► 151
37th prime ───► 157
38th prime ───► 163
39th prime ───► 167
40th prime ───► 173
41st prime ───► 179
42nd prime ───► 181
43rd prime ───► 191
44th prime ───► 193
45th prime ───► 197
46th prime ───► 199
46 primes found.
output when using the input of: -1000
168 primes found.
output when using the input of: -10000
1229 primes found.
output when using the input of: -100000
9592 primes found.
output when using the input of: -1000000
78498 primes found.
output when using the input of: -10000000
664579 primes found.
output when using the input of: -100000000
9 +++ @.m = 0 /*mark odd multiples of J ¬prime.*/
Error 5 running "C:\sieve_of_Eratosthenes.rex", line 9: System resources exhausted
The above shows one of the weaknesses of this implementation of the Sieve of Eratosthenes: it must keep an array of all (if not most) values which is used to strike out composite numbers.
The System resources exhausted error can be delayed by implementing further optimization to the wheel version.
[edit] wheel version
This version skips striking the even numbers (as being not prime).
/*REXX pgm gens primes via a wheeled sieve of Eratosthenes algorithm.*/
parse arg H .; if H=='' then H=200 /*high# can be specified on C.L. */
w = length(H); prime=right('prime',20) /*W is used for formatting output*/
if 2 <= H then say prime right(1,w) " ───► " right(2,w)
@.=1 /*assume all numbers are prime. */
do j=3 by 2 while j*j <= H /*odd integers up to √H inclusive*/
if @.j then do m=j*j to H by j+j /*Prime? Then strike multiples.*/
@.m = 0 /*mark odd multiples of J ¬prime.*/
end /*m*/
end /*j*/
#=1 /*display a list of odd primes.*/
do n=3 to H by 2 /* [↓] if prime, then display it.*/
if @.n then do; #=#+1; say prime right(#,w) " ───► " right(n,w); end
end /*n*/
/*stick a fork in it, we're done.*/
output is identical to the first (non-wheel) version; program execution is over twice as fast as the non-wheel version.
[edit] Wheel Version restructured
/*REXX program generates primes via sieve of Eratosthenes algorithm.
* 21.07.2012 Walter Pachl derived from above Rexx version
* avoid symbols @ and # (not supported by ooRexx)
* avoid negations (think positive)
**********************************************************************/
highest=200 /*define highest number to use. */
is_prime.=1 /*assume all numbers are prime. */
w=length(highest) /*width of the biggest number, */
/* it's used for aligned output.*/
Do j=3 To highest By 2, /*strike multiples of odd ints. */
While j*j<=highest /* up to sqrt(highest) */
If is_prime.j Then Do
Do jm=j*3 To highest By j+j /*start with next odd mult. of J */
is_prime.jm=0 /*mark odd mult. of J not prime. */
End
End
End
np=0 /*number of primes shown */
Call tell 2
Do n=3 To highest By 2 /*list all the primes found. */
If is_prime.n Then Do
Call tell n
End
End
Exit
tell: Parse Arg prime
np=np+1
Say ' prime number' right(np,w) " --> " right(prime,w)
Return
output is identical to the above versions.
Program execution time is slightly worse than that of the wheel version above
(tested with highest=2000000 and suppressed output: 2.4 seconds vs. 2.3 seconds).
[edit] Ruby
eratosthenes starts with nums = [0, 0, 2, 3, 4, 5, ..., n], then marks (zeroes out) multiples of 2, 3, 5, 7, ... there, then returns all non-zero numbers which are the primes.
def eratosthenes(n)
nums = [0, 0] + (2..n).to_a
(2..Math.sqrt(n).to_i).each do |i|
if nums[i].nonzero?
(i**2..n).step(i) {|m| nums[m] = 0}
end
end
nums.find_all {|m| m.nonzero?}
end
p eratosthenes(100)
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
[edit] With a wheel
eratosthenes2 adds more optimizations, but the code is longer.
- The array
numsonly tracks odd numbers (skips multiples of 2). - The array
numsholds booleans instead of integers, and every multiple of 3 beginsfalse. - The outer loop skips multiples of 2 and 3.
- Both inner loops skip multiples of 2 and 3.
def eratosthenes2(n)
# For odd i, if i is prime, nums[i >> 1] is true.
# Set false for all multiples of 3.
nums = [true, false, true] * ((n + 5) / 6)
nums[0] = false # 1 is not prime.
nums[1] = true # 3 is prime.
# Outer loop and both inner loops are skipping multiples of 2 and 3.
# Outer loop checks i * i > n, same as i > Math.sqrt(n).
i = 5
until(m = i * i
m > n) do
if nums[i >> 1]
i_times_2 = i << 1
i_times_4 = i << 2
while m <= n
nums[m >> 1] = false
m += i_times_2
nums[m >> 1] = false
m += i_times_4 # When i = 5, skip 45, 75, 105, ...
end
end
i += 2
if nums[i >> 1]
m = i * i
i_times_2 = i << 1
i_times_4 = i << 2
while m <= n
nums[m >> 1] = false
m += i_times_4 # When i = 7, skip 63, 105, 147, ...
nums[m >> 1] = false
m += i_times_2
end
end
i += 4
end
primes = [2]
nums.each_index {|i| primes << (i * 2 + 1) if nums[i]}
primes.pop while primes.last > n
return primes
end
p eratosthenes2(100)
This simple benchmark compares eratosthenes with eratosthenes2.
require 'benchmark'
Benchmark.bmbm {|x|
x.report("eratosthenes") { eratosthenes(1_000_000) }
x.report("eratosthenes2") { eratosthenes2(1_000_000) }
}
eratosthenes2 runs about 4 times faster than eratosthenes.
[edit] With the standard library
MRI 1.9.x implements the sieve of Eratosthenes at file prime.rb, class EratosthensesSeive (around line 421). This implementation optimizes for space, by packing the booleans into 16-bit integers. It also hardcodes all primes less than 256.
require 'prime'
p Prime::EratosthenesGenerator.new.take_while {|i| i <= 100}
[edit] Run BASIC
input "Gimme the limit:"; limit
dim flags(limit)
for i = 2 to limit
for k = i*i to limit step i
flags(k) = 1
next k
if flags(i) = 0 then print i;", ";
next i
Gimme the limit:?100 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,
[edit] 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) )
BitSet(2, 3, 5, 7, 11, 13, 17, 19, 23, 29)
Using Streams (the "unfaithful sieve", see the link under the Haskell entry):
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 )
Both println statements give the same results:
List(2, 3, 5, 7, 11, 13, 17, 19, 23, 29)
[edit] Scheme
; Tail-recursive solution :
(define (sieve n)
(define (aux u v)
(let ((p (car v)))
(if (> (* p p) n)
(let rev-append ((u u) (v v))
(if (null? u) v (rev-append (cdr u) (cons (car u) v))))
(aux (cons p u)
(let wheel ((u '()) (v (cdr v)) (a (* p p)))
(cond ((null? v) (reverse u))
((= (car v) a) (wheel u (cdr v) (+ a p)))
((> (car v) a) (wheel u v (+ a p)))
(else (wheel (cons (car v) u) (cdr v) a))))))))
(aux '(2)
(let range ((v '()) (k (if (odd? n) n (- n 1))))
(if (< k 3) v (range (cons k v) (- k 2))))))
; > (sieve 100)
; (2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)
; > (length (sieve 10000000))
; 664579
; Simpler solution, with the penalty that none of 'iota, 'strike or 'sieve is tail-recursive :
(define (iota start stop stride)
(if (> start stop)
(list)
(cons start (iota (+ start stride) stop stride))))
(define (strike lst start stride)
(cond ((null? lst) lst)
((= (car lst) start) (strike (cdr lst) (+ start stride) stride))
((> (car lst) start) (strike lst (+ start stride) stride))
(else (cons (car lst) (strike (cdr lst) start stride)))))
(define (primes limit)
(let ((stop (sqrt limit)))
(define (sieve lst)
(let ((p (car lst)))
(if (> p stop)
lst
(cons p (sieve (strike (cdr lst) (* p p) p))))))
(sieve (iota 2 limit 1))))
(display (primes 100))
(newline)
Output:
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)
Optimised using a pre-computed wheel based on 2 (i.e. odds only):
(define (primes-wheel-2 limit)
(let ((stop (sqrt limit)))
(define (sieve lst)
(let ((p (car lst)))
(if (> p stop)
lst
(cons p (sieve (strike (cdr lst) (* p p) (* 2 p)))))))
(cons 2 (sieve (iota 3 limit 2)))))
(display (primes-wheel-2 100))
(newline)
Output:
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)
Vector-based (faster), works with R5RS:
; initialize v to vector of sequential integers
(define (initialize! v)
(define (iter v n) (if (>= n (vector-length v))
(values)
(begin (vector-set! v n n) (iter v (+ n 1)))))
(iter v 0))
; set every nth element of vector v to 0,
; starting with element m
(define (strike! v m n)
(cond ((>= m (vector-length v)) (values))
(else (begin
(vector-set! v m 0)
(strike! v (+ m n) n)))))
; lowest non-zero index of vector v >= n
(define (nextprime v n)
(if (zero? (vector-ref v n))
(nextprime v (+ n 1))
(vector-ref v n)))
; remove elements satisfying pred? from list lst
(define (remove pred? lst)
(cond
((null? lst) '())
((pred? (car lst))(remove pred? (cdr lst)))
(else (cons (car lst) (remove pred? (cdr lst))))))
; the sieve itself
(define (sieve n)
(define stop (sqrt n))
(define (iter v p)
(cond
((> p stop) v)
(else
(begin
(strike! v (* p p) p)
(iter v (nextprime v (+ p 1)))))))
(let ((v (make-vector (+ n 1))))
(initialize! v)
(vector-set! v 1 0) ; 1 is not a prime
(remove zero? (vector->list (iter v 2)))))
[edit] Seed7
The program below computes the number of primes between 1 and 10000000:
$ include "seed7_05.s7i";
const func set of integer: eratosthenes (in integer: n) is func
result
var set of integer: sieve is EMPTY_SET;
local
var integer: i is 0;
var integer: j is 0;
begin
sieve := {2 .. n};
for i range 2 to sqrt(n) do
if i in sieve then
for j range i ** 2 to n step i do
excl(sieve, j);
end for;
end if;
end for;
end func;
const proc: main is func
begin
writeln(card(eratosthenes(10000000)));
end func;
Original source: [1]
[edit] SNOBOL4
Using strings instead of arrays, and the square/sqrt optimizations.
define('sieve(n)i,j,k,p,str,res') :(sieve_end)
sieve i = lt(i,n - 1) i + 1 :f(sv1)
str = str (i + 1) ' ' :(sieve)
sv1 str break(' ') . j span(' ') = :f(return)
sieve = sieve j ' '
sieve = gt(j ^ 2,n) sieve str :s(return) ;* Opt1
res = ''
str (arb ' ') @p ((j ^ 2) ' ') ;* Opt2
str len(p) . res = ;* Opt2
sv2 str break(' ') . k span(' ') = :f(sv3)
res = ne(remdr(k,j),0) res k ' ' :(sv2)
sv3 str = res :(sv1)
sieve_end
* # Test and display
output = sieve(100)
end
Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
[edit] Tcl
package require Tcl 8.5
proc sieve n {
if {$n < 2} {return {}}
# create a container to hold the sequence of numbers.
# use a dictionary for its speedy access (like an associative array)
# and for its insertion order preservation (like a list)
set nums [dict create]
for {set i 2} {$i <= $n} {incr i} {
# the actual value is never used
dict set nums $i ""
}
set primes [list]
while {[set nextPrime [lindex [dict keys $nums] 0]] <= sqrt($n)} {
dict unset nums $nextPrime
for {set i [expr {$nextPrime ** 2}]} {$i <= $n} {incr i $nextPrime} {
dict unset nums $i
}
lappend primes $nextPrime
}
return [concat $primes [dict keys $nums]]
}
puts [sieve 100] ;# 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
[edit] UNIX Shell
[edit] With array
#!/usr/bin/zsh
function primes() {
typeset -a a
typeset i j
a[1]=""
for (( i = 2; i <= $1; i++ )); do
a[$i]=$i
done
for (( i = 2; i * i <= $1; i++ )); do
if [[ ! -z $a[$i] ]]; then
for (( j = i * i; j <= $1; j += i )); do
a[$j]=""
done
fi
done
print $a
}
primes 1000
function primes {
typeset a i=2 j m=$1
# No for (( ... )) loop in pdksh. Use while loop.
while (( i <= m )); do
a[$i]=$i
(( i++ ))
done
i=2
while (( j = i * i, j <= m )); do
if [[ -n ${a[$i]} ]]; then
while (( j <= m )); do
unset a[$j]
(( j += i ))
done
fi
(( i++ ))
done
# No print command in bash. Use echo command.
echo ${a[*]}
}
primes 1000
Both scripts output a single long line.
2 3 5 7 11 13 17 19 23 ... 971 977 983 991 997
[edit] Using variables as fake array
Bourne Shell and Almquist Shell have no arrays. This script works with bash or dash (standard shell in Ubuntu), but uses no specifics of the shells, so it works with plain Bourne Shell as well.
#! /bin/sh
LIMIT=1000
# As a workaround for missing arrays, we use variables p2, p3, ...,
# p$LIMIT, to represent the primes. Values are true or false.
# eval p$i=true # Set value.
# eval \$p$i # Run command: true or false.
#
# A previous version of this script created a temporary directory and
# used files named 2, 3, ..., $LIMIT to represent the primes. We now use
# variables so that a killed script does not leave extra files. About
# performance, variables are about as slow as files.
i=2
while [ $i -le $LIMIT ]
do
eval p$i=true # was touch $i
i=`expr $i + 1`
done
i=2
while
j=`expr $i '*' $i`
[ $j -le $LIMIT ]
do
if eval \$p$i # was if [ -f $i ]
then
while [ $j -le $LIMIT ]
do
eval p$j=false # was rm -f $j
j=`expr $j + $i`
done
fi
i=`expr $i + 1`
done
# was echo `ls|sort -n`
echo `i=2
while [ $i -le $LIMIT ]; do
eval \\$p$i && echo $i
i=\`expr $i + 1\`
done`
[edit] With piping
This version works by piping 1s and 0s through sed. The string of 1s and 0s represents the array of primes.
# Fill $1 characters with $2.
fill () {
# This pipeline would begin
# head -c $1 /dev/zero | ...
# but some systems have no head -c. Use dd.
dd if=/dev/zero bs=$1 count=1 2>/dev/null | tr '\0' $2
}
filter () {
# Use sed to put an 'x' after each multiple of $1, remove
# first 'x', and mark non-primes with '0'.
sed -e s/$2/\&x/g -e s/x// -e s/.x/0/g | {
if expr $1 '*' $1 '<' $3 > /dev/null; then
filter `expr $1 + 1` .$2 $3
else
cat
fi
}
}
# Generate a sequence of 1s and 0s indicating primality.
oz () {
fill $1 1 | sed s/1/0/ | filter 2 .. $1
}
# Echo prime numbers from 2 to $1.
prime () {
# Escape backslash inside backquotes. sed sees one backslash.
echo `oz $1 | sed 's/./&\\
/g' | grep -n 1 | sed s/:1//`
}
prime 1000
[edit] C Shell
# Sieve of Eratosthenes: Echoes all prime numbers through $limit.
@ limit = 80
if ( ( $limit * $limit ) / $limit != $limit ) then
echo limit is too large, would cause integer overflow.
exit 1
endif
# Use $prime[2], $prime[3], ..., $prime[$limit] as array of booleans.
# Initialize values to 1 => yes it is prime.
set prime=( `repeat $limit echo 1` )
# Find and echo prime numbers.
@ i = 2
while ( $i <= $limit )
if ( $prime[$i] ) then
echo $i
# For each multiple of i, set 0 => no it is not prime.
# Optimization: start at i squared.
@ m = $i * $i
while ( $m <= $limit )
set prime[$m] = 0
@ m += $i
end
endif
@ i += 1
end
[edit] Ursala
with no optimizations
#import nat
sieve = ~<{0,1}&& iota; @NttPX ~&r->lx ^/~&rhPlC remainder@rlX~|@r
test program:
#cast %nL
example = sieve 50
output:
<2,3,5,7,11,13,17,19,23,29,31,37,41,43,47>
[edit] Vala
Without any optimizations:
using Gee;
ArrayList<int> primes(int limit){
var sieve = new ArrayList<bool>();
var prime_list = new ArrayList<int>();
for(int i = 0; i <= limit; i++){
sieve.add(true);
}
sieve[0] = false;
sieve[1] = false;
for (int i = 2; i <= limit/2; i++){
if (sieve[i] != false){
for (int j = 2; i*j <= limit; j++){
sieve[i*j] = false;
}
}
}
for (int i = 0; i < sieve.size; i++){
if (sieve[i] != false){
prime_list.add(i);
}
}
return prime_list;
} // end primes
public static void main(){
var prime_list = primes(50);
foreach(var prime in prime_list)
stdout.printf("%s ", prime.to_string());
stdout.printf("\n");
}
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 .NET
Dim n As Integer, k As Integer, limit As Integer
Console.WriteLine("Enter number to search to: ")
limit = Console.ReadLine
Dim flags(limit) As Integer
For n = 2 To Math.Sqrt(limit)
If flags(n) = 0 Then
For k = n * n To limit Step n
flags(k) = 1
Next k
End If
Next n
' Display the primes
For n = 2 To limit
If flags(n) = 0 Then
Console.WriteLine(n)
End If
Next n
[edit] Vorpal
self.print_primes = method(m){
p = new()
p.fill(0, m, 1, true)
count = 0
i = 2
while(i < m){
if(p[i] == true){
p.fill(i+i, m, i, false)
count = count + 1
}
i = i + 1
}
('primes: ' + count + ' in ' + m).print()
for(i = 2, i < m, i = i + 1){
if(p[i] == true){
('' + i + ', ').put()
}
}
''.print()
}
self.print_primes(100)
Result:
primes: 25 in 100 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,
[edit] XPL0
include c:\cxpl\codes; \intrinsic 'code' declarations
int Size, Prime, I, Kill;
char Flag;
[Size:= IntIn(0);
Flag:= Reserve(Size+1);
for I:= 2 to Size do Flag(I):= true;
for I:= 2 to Size do
if Flag(I) then \found a prime
[Prime:= I;
IntOut(0, Prime); CrLf(0);
Kill:= Prime + Prime; \first multiple to kill
while Kill <= Size do
[Flag(Kill):= false; \zero a non-prime
Kill:= Kill + Prime; \next multiple
];
];
]
Example output:
20 2 3 5 7 11 13 17 19
Categories:
- Programming Tasks
- Prime Numbers
- Clarified and Needing Review
- ACL2
- Bracmat
- GLBasic
- Ada
- ALGOL 68
- AutoHotkey
- AutoIt
- AWK
- BASIC
- Applesoft BASIC
- Locomotive Basic
- ZX Spectrum Basic
- BBC BASIC
- Befunge
- C
- C++
- C sharp
- Clojure
- Lang examples needing attention
- Examples needing attention
- CMake
- Common Lisp
- D
- Dart
- DWScript
- Dylan
- E
- EC
- Eiffel
- Erlang
- Euphoria
- F Sharp
- Forth
- Fortran
- GAP
- Go
- Groovy
- GW-BASIC
- Haskell
- HicEst
- Icon
- Unicon
- J
- Java
- JavaScript
- Liberty BASIC
- Logo
- Lua
- Lucid
- M4
- Mathematica
- MATLAB
- Maxima
- MAXScript
- Modula-3
- MUMPS
- NetRexx
- Nial
- Niue
- Oberon-2
- OCaml
- Oz
- PARI/GP
- Pascal
- Perl
- Perl 6
- PHP
- PicoLisp
- PL/I
- Pop11
- PowerShell
- Processing
- Prolog
- PureBasic
- Python
- Numpy
- R
- Racket
- REXX
- Ruby
- Run BASIC
- Scala
- Scheme
- Seed7
- SNOBOL4
- Tcl
- UNIX Shell
- C Shell
- Ursala
- Vala
- Gee
- Vedit macro language
- Visual Basic .NET
- Vorpal
- XPL0
