Sieve of Eratosthenes: Difference between revisions
→{{header|C}}: add C++ version |
m →{{header|C++}}: Boost link, syntax highlighting |
||
Line 246: | Line 246: | ||
Possible optimizations include sieving only odd numbers (or more complex wheels), packing the sieve into bits to improve locality (and allow larger sieves), etc. |
Possible optimizations include sieving only odd numbers (or more complex wheels), packing the sieve into bits to improve locality (and allow larger sieves), etc. |
||
=={{header|C++}}== |
=={{header|C++}}== |
||
{{works with|Boost}} |
|||
< |
<cpp>// yield all prime numbers less then limit. |
||
template<class UnaryFunction> |
template<class UnaryFunction> |
||
void primesupto(int limit, UnaryFunction yield) |
void primesupto(int limit, UnaryFunction yield) |
||
Line 258: | Line 259: | ||
for (unsigned k = n*n, ulim = static_cast<unsigned>(limit); k < ulim; k += n) |
for (unsigned k = n*n, ulim = static_cast<unsigned>(limit); k < ulim; k += n) |
||
//NOTE: |
//NOTE: "unsigned" is used to avoid an overflow in `k+=n` for `limit` near INT_MAX |
||
is_prime[k] = false; |
is_prime[k] = false; |
||
} |
} |
||
Line 265: | Line 266: | ||
if (is_prime[n]) |
if (is_prime[n]) |
||
yield(n); |
yield(n); |
||
}</ |
}</cpp> |
||
Full program: |
Full program: |
||
< |
<cpp>/** |
||
$ g++ -I/path/to/boost sieve.cpp -o sieve && sieve 10000000 |
$ g++ -I/path/to/boost sieve.cpp -o sieve && sieve 10000000 |
||
*/ |
*/ |
||
Line 310: | Line 311: | ||
cout << "limit sum pi(n)\n" |
cout << "limit sum pi(n)\n" |
||
<< limit << " " << sum << " " << count << endl; |
<< limit << " " << sum << " " << count << endl; |
||
}</ |
}</cpp> |
||
=={{header|Common Lisp}}== |
=={{header|Common Lisp}}== |
Revision as of 18:29, 6 October 2008
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
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. 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
[[Category:{{{1}}} examples needing attention]]Property "Example requires attention" (as page type) with input value "{{{1}}}" contains invalid characters or is incomplete and therefore can cause unexpected results during a query or annotation process.
This solution uses an array of bits representing boolean values to identify the prime numbers. The limit on the outer loop is set to the square root of the number received from the command line. This limit minimizes the number of unneeded iterations in the algorithm.
The lowest index of the array of boolean is set to 2 so that all calculations begin at 2.
with Ada.Text_Io; use Ada.Text_Io; with Ada.Command_Line; use Ada.Command_Line; with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions; procedure Sieve is package Nat_Io is new Ada.Text_Io.Integer_Io(Natural); use Nat_Io; type Bool_Array is array(Positive range <>) of Boolean; pragma pack (Bool_Array); procedure Get_Sieve (A : out Bool_Array) is Test_Num : Positive; Limit : Positive := Positive(sqrt(float(A'Last))); begin A := (Others => True); -- set all values to true for Num in A'First..Limit loop if A(Num) then Test_Num := Num * Num; while Test_Num <= A'Last loop A(Test_Num) := False; Test_Num := Test_Num + Num; end loop; end if; end loop; end Get_Sieve; pragma Inline(Get_Sieve); N : Positive := 2; begin if Argument_Count > 0 then N := Positive'Value(Argument(1)); end if; declare Vals : Bool_Array(2..N); begin Get_Sieve(Vals); for I in Vals'range loop if Vals(I) then Put_Line(Positive'Image(I)); end if; end loop; end; end Sieve;
This version uses tasking to try to improve the performance on multiprocessor systems which are becoming common. It is still an implementation of the Sieve of Eratosthenes: the list of numbers is in the main task, and the list of primes is in the list of Siever tasks.
with Ada.Command_Line; with Ada.Text_IO; procedure Sieve is task type Siever is entry Check (Value : in Positive); entry Print; end Siever; subtype Siever_Task is Siever; type Sieve_Ptr is access Siever; task body Siever is Number : Positive; Candidate : Positive; Next : Sieve_Ptr; begin -- Siever accept Check (Value : in Positive) do Number := Value; end Check; All_Values : loop select accept Check (Value : in Positive) do Candidate := Value; end Check; if Candidate rem Number /= 0 then if Next = null then Next := new Siever_Task; end if; Next.Check (Value => Candidate); end if; or accept Print do Ada.Text_IO.Put (Item => Integer'Image (Number) ); if Next /= null then Next.Print; end if; end Print; exit All_Values; end select; end loop All_Values; end Siever; Head : Sieve_Ptr := new Siever; Max : Positive := 20; begin -- Sieve if Ada.Command_Line.Argument_Count > 0 then Max := Integer'Value (Ada.Command_Line.Argument (1) ); end if; All_Values : for I in 2 .. Max loop Head.Check (Value => I); end loop All_Values; Head.Print; end Sieve;
This version is the same as the previous version, but with a wheel of 2.
with Ada.Command_Line; with Ada.Text_IO; procedure Sieve_With_Wheel is task type Siever is entry Check (Value : in Positive); entry Print; end Siever; subtype Siever_Task is Siever; type Sieve_Ptr is access Siever; task body Siever is Number : Positive; Candidate : Positive; Next : Sieve_Ptr; begin -- Siever accept Check (Value : in Positive) do Number := Value; end Check; All_Values : loop select accept Check (Value : in Positive) do Candidate := Value; end Check; if Candidate rem Number /= 0 then if Next = null then Next := new Siever_Task; end if; Next.Check (Value => Candidate); end if; or accept Print do Ada.Text_IO.Put (Item => Integer'Image (Number) ); if Next /= null then Next.Print; end if; end Print; exit All_Values; end select; end loop All_Values; end Siever; Head : Sieve_Ptr := new Siever; Max : Positive := 20; Value : Positive := 3; begin -- Sieve_With_Wheel if Ada.Command_Line.Argument_Count > 0 then Max := Integer'Value (Ada.Command_Line.Argument (1) ); end if; All_Values : loop exit All_Values when Value > Max; Head.Check (Value => Value); Value := Value + 2; end loop All_Values; Ada.Text_IO.Put (Item => " 2"); Head.Print; end Sieve_With_Wheel;
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
Befunge
2>:3g" "-!v\ g30 < |!`"O":+1_:.:03p>03g+:"O"`| @ ^ p3\" ":< 2 234567890123456789012345678901234567890123456789012345678901234567890123456789
C
Plain sieve, without any optimizations:
#include <stdlib.h> #include <math.h> char* eratosthenes (int n) { char* sieve; int i, j, m; /* calloc initializes to zero */ sieve = calloc (n+1, sizeof(char)); m = (int) sqrt ((double) n); sieve[0] = 1; sieve[1] = 1; for (i = 2; i <= m; i++) { if (sieve[i] == 0) { for (j = i*i; j <= n; j += i) { sieve[j] = 1; } } } return sieve; }
Possible optimizations include sieving only odd numbers (or more complex wheels), packing the sieve into bits to improve locality (and allow larger sieves), etc.
C++
<cpp>// yield all prime numbers less then limit. template<class UnaryFunction> void primesupto(int limit, UnaryFunction yield) {
std::vector<bool> is_prime(limit, true); const int sqrt_limit = static_cast<int>(std::sqrt(limit)); for (int n = 2; n <= sqrt_limit; ++n) if (is_prime[n]) {
yield(n);
for (unsigned k = n*n, ulim = static_cast<unsigned>(limit); k < ulim; k += n)
//NOTE: "unsigned" is used to avoid an overflow in `k+=n` for `limit` near INT_MAX
is_prime[k] = false;
}
for (int n = sqrt_limit + 1; n < limit; ++n) if (is_prime[n])
yield(n); }</cpp>
Full program:
<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;
}</cpp>
Common Lisp
(defun primer () ((lambda (&rest q) (lambda () ((lambda (u n) (funcall u u n)) (lambda (u n) (if (some (lambda (x y) (and (/= 1 y) (= (mod x y) 0))) ((lambda (x &aux (y (list x))) (setf (cdr y) y)) n) q) (funcall u u (1+ n)) (setf q (cons n q) n n))) (1+ (car q))))) 1)) ;; Example use: (defparameter *png* (primer)) (loop repeat 50 collect (funcall *png*))
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
Forth
: sieve ( n -- ) here over erase 2 begin 2dup dup * > while here over + c@ 0= if 2dup dup * do 1 here i + c! dup +loop then 1+ repeat drop ." Primes: " 2 do here i + c@ 0= if i . then loop ; 100 sieve
Fortran
Fortran 90 and later
PROGRAM SIEVE IMPLICIT NONE INTEGER :: i, j, number LOGICAL, ALLOCATABLE :: pflags(:) WRITE(*,*) "Enter number to search to" READ (*,*) number ALLOCATE (pflags(2:number)) pflags = .TRUE. DO i = 2, SQRT(REAL(number)) IF (pflags(i)) THEN DO j = i+i, number, i pflags(j) = .FALSE. END DO END IF END DO DO i = 2, number IF (pflags(i)) WRITE(*,*) i END DO DEALLOCATE (pflags) END PROGRAM SIEVE
Since 2 is trivially the only even prime number this version just sieves odd numbers
ALLOCATE (pflags((number-1)/2)) pflags = .TRUE. DO i = 1, SQRT(REAL(number-1)/2) IF (pflags(i)) THEN DO j = 3*i+1, (number-1)/2, 2*i+1 pflags(j) = .FALSE. END DO END IF END DO DO i = 1, (number-1)/2 IF (pflags(i)) WRITE(*,*) 2*i + 1 END DO
Haskell
[[Category:{{{1}}} examples needing attention]]Property "Example requires attention" (as page type) with input value "{{{1}}}" contains invalid characters or is incomplete and therefore can cause unexpected results during a query or annotation process.
Imperative version, with unboxed arrays:
import Control.Monad import Control.Monad.ST import Data.Array.ST import Data.Array.Unboxed primeSieve n = runST sieve where newSieve :: ST s (STUArray s Integer Bool) newSieve = newArray (2,n) True sieve :: ST s (UArray Integer Bool) sieve = do s <- newSieve 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 unsafeFreeze s
Return primes from sieve as list:
primesTo n = [p | p <- [2..n], s ! p] where s = primeSieve n
One often sees a purely functional version similar to
primes :: [Integer] primes = sieve [2..] where sieve (p:xs) = p : sieve [x | x <- xs, x `mod` p > 0]
As Melissa O'Neill demonstrates in The Genuine Sieve of Eratosthenes, this version is not better than trial division. She also presents an actual purely functional version that uses priority queues, and is much easier to combine with wheels (pre-calculated sieves for small primes) than the imperative version. Even combined with a simple wheel, the speed of this version is similar to the imperative version.
J
[[Category:{{{1}}} examples needing attention]]Property "Example requires attention" (as page type) with input value "{{{1}}}" contains invalid characters or is incomplete and therefore can cause unexpected results during a query or annotation process.
sieve=: 3 : 0 m=. <.%:y z=. $0 b=. y{.1 while. m>j=. 1+b i. 0 do. b=. b+.y$(-j){.1 z=. z,j end. z,1+I.-.b )
"Wheels" are readily implemented as follows:
sieve1=: 3 : 0 m=. <.%:y z=. y (>:#]) 2 3 5 7 b=. 1,}.y$+./(*/z)$&>(-z){.&.>1 while. m>j=. 1+b i. 0 do. b=. b+.y$(-j){.1 z=. z,j end. z,1+I.-.b )
The use of 2 3 5 7 as wheels provides a 20% time improvement for n=1000 and 2% for n=1e6 .
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); }
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
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
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
Nial
primes is sublist [ each (2 = sum eachright (0 = mod) [pass,count]), pass ] rest count
Using it
|primes 10 =2 3 5 7
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";
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;
Python
Using list lookup
<python> def eratosthenes(n):
multiples = [] print 2 for i in xrange(3, n+1): if not i in multiples: print i for j in xrange(i*i, n+1, i): multiples.append(j) eratosthenes(100)</python>
Using set lookup
The version below uses a set
to store the multiples. set
objects are much faster (usually O(log n)) than list
s (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.
<python> def eratosthenes2(n):
multiples = set() print 2 for i in xrange(3, n+1): if not i in multiples: print i multiples.update(xrange(i*i, n+1, i)) eratosthenes2(100)</python>
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
.
<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]</python>
Using generator
<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
</python> Example: <python> >>> list(iprimes_upto(15)) [2, 3, 5, 7, 11, 13]</python>
Using numpy
Belowed code adapted from literateprograms.org using numpy <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:]</python>
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: <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]</python>
Examples: <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])</python>
The above-mentioned examples demonstrate that the given wheel based optimization does not show significant performance gain.
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