## Sieve of Eratosthenes

### February 19, 2009

Over two millenia ago, Eratosthenes, who calculated the circumference of the earth, the distance to the Sun and the tilt of the Earth’s axis, developed a system of latitude and longitude, and invented the leap day, created a systematic method to enumerate the prime numbers that is still in use today. Eratosthenes was born in Cyrene (present-day Libya), lived from 276 B.C. to 194 B.C., and spent most of his life in Alexandria, Egypt, where he was the second Chief Librarian of the Great Library, succeeding Apollonius of Rhodes; he was a good friend of Archimedes.

The Sieve of Eratosthenes starts by making a list of all the numbers up to a desired maximum; we’ll illustrate the method by calculating the prime numbers through thirty:

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

Now take the first number on the list, 2, and cross off every second number:

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

(Although it may not be obvious, the number 4 is crossed off the list; in some fonts, the cross-bar of the 4 coincides with the strike-through bar.) Next, take the next number on the list that isn’t crossed off, 3, and cross off every third number; some of them have already been crossed off:

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

Repeat that last step for the next un-crossed number on the list, 5:

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

And so on, each time crossing off all multiples of the next un-crossed number on the list. The list of prime numbers are all those that haven’t been crossed off:

2 3 5 7 11 13 17 19 23 29

This method is called a sieve because it sweeps through a range of numbers, with each prime number, as it is discovered, blocking all its multiples from falling through as prime numbers. The sieve admits several optimizations. First, only odd numbers are considered, since the initial sifting crosses off all the even numbers except 2, which is handled separately. Second, crossing off starts at the square of the number being sifted, since all smaller primes have already been crossed off by previous steps of the sieve; for instance, sifting by 3 starts at 9, since 6 was already crossed off when sifting by 2. Third, sifting stops at the square root of the maximum number in the sieve, since any non-primes larger than the square root must have already been crossed off at previous levels of the sieve; thus, in the above example there is no need to sieve on the prime number 7, or any larger prime number, since the square of 7 is greater than 30, which is the largest number in the list.

Write a function that takes a single argument n and returns a list of prime numbers less than or equal to n using the optimized sieving algorithm described above. Apply the function to the argument 15485863 and count the number of primes returned.

Pages: 1 2

### 60 Responses to “Sieve of Eratosthenes”

1. Dr Frankenstein said
2. […] this is the second Programming Praxis exercise I’ve done, and this one was slightly more difficult to get right, or at least test. […]

3. Stuart said

I’m confused by one aspect of the above solution. I don’t understand why startj is defined as (+ (* 2 i i) (* 6 i) 3).

I understand why p is (+ i i 3) – that’s the sequence of odd numbers, starting at 3. But what relationship does that have with startj?

Can anyone enlighten me?

4. programmingpraxis said

That’s the implementation of the second optimization.

Suppose you have already sieved 2, 3, and 5 and are now beginning to sieve 7. Start by adding 7 to the list of primes. Then 14, but that has already been sifted out by 2; likewise, 21 has been sifted out by 3, 28 has been sifted out by 2, 35 has been sifted out by 5, and 42 has been sifted out by 2 (and also 3, but that doesn’t matter). The first multiple of 7 that gets sifted out is 49, which is 7 times 7. In general, when sifting by n, sifting starts at n-squared because all the previous multiples of n have already been sieved.

The rest of the expression has to do with the cross-reference between numbers and sieve indexes. There’s a 2 in the expression because we eliminated all the even numbers before we ever started. There’s a 3 in the expression because Scheme vectors are zero-based, and the numbers 0, 1 and 2 aren’t part of the sieve. I think the 6 is actually a combination of the 2 and the 3, but it’s been a while since I looked at the code, so I’ll leave it to you to figure out.

Good question!

5. programmingpraxis said

The number 2 is handled specially, so the vector v looks this:

```index  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 number 3  5  7  9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49```

When you want to begin sieving 7, the first number you sift out is 7 × 7 = 49. The number 7 is at index i = 2 in the vector, and the number 49 is at index startj = (+ (* 2 i i) (* 6 i) 3) = 23 in the vector.

6. Stuart said

So, it seems we can alternatively define startj as:

`(quotient (- (expt p 2) 3) 2)`

That is,

– p-squared, to give us the number to start sifting at
– shifted by 3 because we ignore the first 3 numbers (0, 1 and 2)
– divided by 2 to get the index, because we ignore every second number.

Is that correct?

7. Jebb said

The straightforward implementation in C, storing all the odd numbers in an array and then zeroing all non-primes. Not being much of a programmer, I’ll admit I’m amazed my laptop can compute 1,000,000 primes in under half a second… And I suppose the code could even be made to run twice as fast on my c2duo by having two threads sifting through the array.

```#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define index(A) (A - 3) / 2    //index of A in the array of odd numbers
#define vector(A) 2 * A + 3     //inverse operation of index()

long *fill_num(long n);         //Generates array of odd numbers > 2
int seek_primes(long *numbers, long n);     //Zeroes all non-primes

int main()
{
long ulimit;
long *nombres;
int premiers;
char c;
printf("Upper limit for the calculation of primes?\n");
scanf("%ld", &ulimit);
ulimit = index(ulimit) + 1;     //+ 1: the array is indexed from 0
nombres = fill_num(ulimit);
premiers = seek_primes(nombres, ulimit) + 1;    //2 is also prime
printf("Found %d primes\n", premiers);
free(nombres);
return 0;
}

int seek_primes(long *numbers, long n)
{
long i, j, primes = 0, max;
//Third optimisation: sift until sqrt of upper limit
printf("sift until %ld\n", max = (long)sqrt(vector(n)));
for (i = 0; i < n; i++)
if (numbers[i] != 0) {
primes++;
if (i < max) {
j = numbers[i] * numbers[i];
while (j < vector(n)) {
numbers[index(j)] = 0;
j = j + 2 * numbers[i];
}
}
}
return primes;
}

long *fill_num(long n)
//n is the size of the array, vector(n - 1) the highest odd number
{
long i;
long *tmp;
tmp = (long *)malloc(n * sizeof(long));
for (i = 0; i < n; i++)
tmp[i] = vector(i);
return tmp;
}
```
8. Jebb said

…turns out the multi-threaded version does not work, at least the way I’ve implemented it. I’ve effectively split the seek_primes in two threads, one sifting for (i = 0; i < n; i+=2), and the other one for (i = 1; i < n; i+=2).

It appears my first thread finds non-zero elements in the numbers[] array in seek_primes, but there is no guarantee these same elements are not going to be zeroed by the second thread later.

In fact up to 100 I find 3 primes too many, and up to 1000 about 38 (and of course the actual number found might even be different from one run to the next).

9. Robert Berman said

My python solution is at http://pastebin.com/Qup0TuFu.
I generated all primes from the integer list of 0-100000. Since the last prime in that list is 99991, using the square of the primes method, it is relatively easy and fast deriving the prime sequence of most numbers. I think I began having problems around the number of primes for the integer 9990000000.

10. neveu said

Here’s my Clojure version:
(require ‘[clojure.contrib.math :as math])

;; faster for dropping a few leading nils than (into [] (drop-while nil? sieve))
(defn drop-while-nil? [sieve] (if (first sieve) sieve (recur (subvec sieve 1))))

;; returns sieve with all multiples of (first sieve) set to nil,
;; leading nils removed
(defn remove-every-nth [sieve]
(let [count-sieve (count sieve), n (first sieve)]
(loop [s sieve, i 1, ni n]
(if (> ni count-sieve)
(drop-while-nil? (subvec s 1))
(recur (assoc s ni nil), (inc i), (* n (inc i)))))))

(defn eratosthenes [r]
(loop [s (into [] (range 2 r)), primes [1] ]
(if (> (peek primes) (math/sqrt r))
(concat primes (remove nil? s))
(recur (remove-every-nth s), (conj primes (first s))))))

(map #(count (time (eratosthenes %))) [10 100 1000 10000 100000 ])

11. Jim Wise said

Simple and simple tail-recursive solutions, but missing the “start eliminating at (* p p)” optimization:

```#lang racket

;; .. : num num -> list-of-nums
;; given first and last, return a list of the range first .. last, inclusive
(define (.. first last)
(if (= first last)
(list first)
(cons first (.. (+ first 1) last))))

(define (erat l n)
(let ([p (car l)])
(if (> (sqr p) n)
l
(cons p (erat (filter (lambda (n) (not (= (modulo n p) 0))) (cdr l)) n)))))

(define (tail-erat l n accum)
(let ([p (car l)])
(if (> (sqr p) n)
(append (reverse accum) l)
(tail-erat (filter (lambda (n) (not (= (modulo n p) 0))) (cdr l)) n (cons p accum)))))

; main program

(let* ([args (vector->list (current-command-line-arguments))]
[n (string->number (car args))]
[l (.. 2 n)])
;  (display (erat l n))
(display (tail-erat l n '()))
(newline))
```
12. programmingpraxis said

Did you calculate the number of primes less than 15485863? I think it will take a while.

13. Jim Wise said

Yes, it takes about fifty times longer than the vector-based solution above, mostly in GC.

Using Petite Chez Scheme for profiling:

(time (display (length (…))))
6 collections
3021 ms elapsed cpu time, including 248 ms collecting
3022 ms elapsed real time, including 247 ms collecting
110325408 bytes allocated, including 40188320 bytes reclaimed
(time (length (tail-erat (…) …)))
2767 collections
167705 ms elapsed cpu time, including 103201 ms collecting
167817 ms elapsed real time, including 103268 ms collecting
23602922896 bytes allocated, including 24373515936 bytes reclaimed

14. Jim Wise said

Moving to an in-place filtering of the list (see below) brought this time down by a factor of ten, though it’s still slower than the vector version:

(time (length (tail-erat (…) …)))
105 collections
11498 ms elapsed cpu time, including 5577 ms collecting
11503 ms elapsed real time, including 5590 ms collecting
1135547536 bytes allocated, including 3928896 bytes reclaimed

This was achieved by replacing filter with the following in-place filter! routine (and staying in Chez Scheme, where set-cdr! is available).

A more general filter!, implemented as syntax! so it could consistently alter the list in place but not depend on returning a list would be interesting, but isn’t needed for this — this implementation returns the result of modifying the list in place, but the original value passed in is only modified from the first value not matching pred? on…

```(define (filter! pred? l)
(cond [(null? l) '()]
[(pred? (car l)) (filter! pred? (cdr l))]
[else
(let loop ([follower l]
(cond
(set-cdr! follower (cdr leader))
(loop follower (cdr leader))]
```
15. programmingpraxis said

You’re missing the point. The Sieve is fast because it is based on addition. You are using division (the modulo function). Your algorithm is not the Sieve of Eratosthenes.

16. Jim Wise said

Oh, I see… let me revisit…

17. Jim Wise said

Naive attempt at moving to vector and using addition over vector indices gives a slight speedup, looking further:

(time (length (primes 15485863)))
4 collections
10296 ms elapsed cpu time, including 366 ms collecting
10300 ms elapsed real time, including 366 ms collecting
155888768 bytes allocated, including 124098560 bytes reclaimed

```(define (primes n)
(define (zeroing-pass! vec start step)
(cond
[(> start n) vec]
[(not (vector-ref vec (sub1 step))) vec]
[else
(vector-set! vec (sub1 start) #f)
(zeroing-pass! vec (+ start step) step)]))

(define (erat! vec)
(let* ([stop (div n 2)])
(let loop ([p 2])
(if (= p stop)
vec
(begin
(zeroing-pass! vec (+ p p) p)
(loop (+ p 1)))))))

(let ([vec (make-vector n #t)])
(erat! vec)
(reverse
(let loop ([i 1]
[l '()])
(if (= i n)
l
(if (vector-ref vec i)
(cons (add1 i) l)
l)))))))
```
18. Jim Wise said

I got curious as to how big a difference the additon-vs-module distinction is on modern hardware (historically, it has been very large, as you point out).

I used the following code to run 10^8 integer additions, modulo operations, and modulo operations with a power-of-two modulus (another traditional distinction in performance).

I may well be missing something big, but the output shows about a 25% difference in speed between additions and modulo operations. This obviously adds up, but looking at the above code, it seems to be small compared to other differences (particularly, the various operations which control how many of whichever operation are carried out, and which reduce the infrastructural cost of the algorithm).

(In the code below, I use three vectors in an attempt to minimize the effect of cache read-ahead on subsequent vector-maps during earlier testing with a smaller n.)

```(let* ([n (expt 10 6)]
[times 100]
[time-reps (lambda (f v)
(time (do ([i 0 (+ i 1)])
([= i times])
(vector-map f v))))]
[v1 (make-vector n 1203006)]
[v2 (make-vector n 1203006)]
[v3 (make-vector n 1203006)])

(time-reps (lambda (n) (+ n 57365)) v1)
(newline)

(display "Timing modulo (power of 2)\n")
(time-reps (lambda (n) (mod n 65536)) v2)
(newline)

(display "Timing modulo (non-power of 2)\n")
(time-reps (lambda (n) (mod n 57365)) v3)
(newline))

```
```Timing addition
(time (do ((...)) ...))
299 collections
6471 ms elapsed cpu time, including 2642 ms collecting
6474 ms elapsed real time, including 2641 ms collecting
2833256224 bytes allocated, including 2762557728 bytes reclaimed

Timing modulo (power of 2)
(time (do ((...)) ...))
299 collections
8536 ms elapsed cpu time, including 2586 ms collecting
8545 ms elapsed real time, including 2597 ms collecting
2833185520 bytes allocated, including 2838706016 bytes reclaimed

Timing modulo (non-power of 2)
(time (do ((...)) ...))
299 collections
8544 ms elapsed cpu time, including 2601 ms collecting
8548 ms elapsed real time, including 2596 ms collecting
2833185520 bytes allocated, including 2821785344 bytes reclaimed

```
19. Grégory LEOCADIE said

my OCaml version. Because 15000000 is bigger than the bound for basic array, i had to use Bigarray.

```open Bigarray
open Sys
open Format

let timeit f =
let tb = time () in
let res = Lazy.force_val f in
printf "time: %f s\n" ((time()) -. tb);
res

let ba_fold_left f init arr =
let adim = Array1.dim arr in
let rec ba_fold_left' i res =
match i >= 0 with
| true -> ba_fold_left' (i - 1) (f res i (Array1.get arr i))
| _ -> res
in
ba_fold_left' (adim - 1) init
;;

let sieve n =
let adim = (n - 1)/ 2 in
let arr = Array1.create int8_unsigned c_layout adim in
Array1.fill arr 0;
let rec erase i mo =
if i < adim
then
begin
if (2 * i + 3) mod mo = 0
then Array1.set arr i 1;
erase (i + mo) mo
end
in
let rec sieve' i =
let e = 2 * i + 3 in
match e*e with
| p when p > n -> ()
| p when (Array1.get arr i) = 0
-> erase (p / 2 - 1) e; sieve' (i + 1)
| _ -> sieve' (i + 1)
in
sieve' 0;
2 :: (ba_fold_left (fun x k v ->
if v = 0 then (k * 2 + 3) :: x
else x) [] arr)
;;

let _ =
Printf.printf " %d "
(List.length (timeit (Lazy.lazy_from_fun (fun () -> (sieve 15485863)))))
```
20. […] ancient Sieve of Eratosthenes that computes the list of prime numbers is inefficient in the sense that some composite numbers are […]

21. David said

This was my first Factor program I wrote a few months ago, to learn the language. Updated it with the optimization to start sieving at the square of the prime being sieved. (Didn’t know about that trick before :)

```USING: kernel math sequences math.ranges bit-sets locals formatting io ;
FROM: sets => adjoin delete in? members ;
IN: sieve

CONSTANT: items/line 10

: .row ( seq -- )
[ "%5d " printf ] each nl ;

: .vec ( seq -- )
[ dup length items/line > ]
[ items/line cut-slice swap .row ]
while
.row ;

: init-sieve  ( max-prime -- bitset )
dup <bit-set>
2 rot [a,b) [ over adjoin ] each ;

: filter-composite ( bitset start end -- )
[ dup sq ] dip  rot  <range>
[ over delete ] each drop ;

:: primes ( maxprime -- vec )
[let
maxprime 1 + init-sieve :> sieve

2
[ dup sq maxprime > ]
[ dup sieve in?
[ dup maxprime sieve -rot filter-composite ]
when
1 +
]
until drop
sieve members
] ;

: .primes ( maxprime -- )
primes .vec ;
```

Session:

```( scratchpad ) 15485863 primes length .
1000000
( scratchpad ) 127 .primes
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
```
22. […] As mentioned in my last article, I started doing some challenges from the Programming Praxis website. And here comes my PHP solution to the second challenge. […]

23. Graham said

I’m working my way through old exercises; here’s my submission.

24. Graham said

Python was a bit slow for me, and I’m trying to learn Common Lisp, so here’s
another version. I was only able to figure out how to get 2/3 of the optimizations
working, but it’s still absurdly fast compared to my last submission.

25. ftt said
```#!/usr/bin/env python

import itertools
import sys

# note: the 'modulo' version with the same cut-offs uses much less memory in Python
# but is also way slower
def sieve(n):
# optimization 1: only odd numbers
primes = range(3, n + 1, 2)
upper_bound = n ** 0.5
for base in xrange(len(primes)):
if not primes[base]:
continue
# optimization 3: stop at the square root of n
if primes[base] >= upper_bound:
break
# optimization 2: start with the square
for i in xrange(base + (base + 1) * primes[base], len(primes), primes[base]):
primes[i] = None
primes.insert(0, 2)
return filter(None, primes)

def main():
print 'Enter an upper bound for the sieve.'
result = sieve(n)
print 'Number of primes: {0}'.format(len(result))

if __name__ == '__main__':
main()
```
26. Rens said

I just started learning ruby. Here is my first try.

def prime_numbers(n)
return [] if n < 2 # no prime numbers before 2
max = (n**0.5).to_i; # define until when we need to loop
primes = []
# optimization 1: checking only odd numbers
primes = 3.step( n, 2 ).to_a
primes.each do |p|
next unless p
# optimization 3: checking until the root of the maximum number
break unless p <= max
# optimization 2: checking starts at the square
((p**2)/2-1).step( primes.length-1, p) do |i|
primes[i] = nil
end
end
primes.insert(0, 2).delete(nil) # add the prime number 2 and remove false numbers
primes
end

puts prime_numbers(15485863).length
[\sourcecode]

27. kawas said

Another clojure version,
because the first one die on me with OutOfMemoryError to compute (eratosthenes 15485863)

```; we work only with odd numbers till N
(defn odds [n] (vec (range 3 (inc n) 2)))

; it's easy to calculate index of a value in our vector of odds
(defn odds-idx [v] (/ (- v 3) 2))

; just sieve the vector from start=(index of v²)
(defn sieve [max-idx nums v]
(let [start (odds-idx (* v v))]
(reduce #(assoc %1 %2 nil) nums (range start max-idx v))))

; get primes equal or less than N
(defn primes [n]
(let [max-i (Math/floor (odds-idx n))]
(loop [nums (odds n) i 0]
(if-let [v (nums i)]
(if (> (* v v) n)
(cons 2 (remove nil? nums))
(recur (sieve max-i nums v) (inc i)))
(recur nums (inc i))))))

(count (primes 15485863))
1000000
```
28. An imperative OCaml Batteries `BitSet`-based version, without multiplications, only additions, subtractions and shifts:

```let sieve_of_eratosthenes (proc : int -> unit) n =
let limit = (n - 3) asr 1 in
let sieve = BitSet.create_full (limit + 1) in
let i = ref 0
and p = ref 3
and q = ref 9 in
while !q < n do
if BitSet.is_set sieve !i then begin
let j = ref ((!q - 3) asr 1) in
while !j <= limit do
BitSet.unset sieve !j;
j := !j + !p
done
end;
incr i;
p := !p + 2;
q := !q + (!i + 1) lsl 3
done;
proc 2;
foreach (BitSet.enum sieve) (fun i -> if i <= limit then proc (i lsl 1 + 3))
```
29. CyberSpace17 said

My C++ solution where a “crossing out” means “made 0”
http://ideone.com/fdmux

30. moink said

Here’s a Matlab solution:

```function primes=sieve(maxnum)
%Sieve of Erasthenes for finding primes less than or equal to the input
%argument
for i=3:2:floor(sqrt(maxnum))
end
end
```

It’s a tiny bit faster than the nearly-identical Matlab implementation, called primes.

31. A rather simple version of a Python solution.
``` def primes_sieve(limit): limit += 1 not_prime = [False] * limit primes = []```

``` for i in xrange(2, limit): if not_prime[i]: continue for j in xrange(i*2, limit, i): not_prime[j] = True primes.append(i) return primes ```

```if __name__ == "__main__": print len(primes_sieve(15485863)) ```

32. Matthias Altmann said

I tried another solution in C++ using 2 threads, one running upwards on the numbers, the other running downwards. They stop if the reach another.
The solution is correct, but the performance decreases severly.

Here is my code:

#include
#include
#include
#include
#include
using namespace std;

const unsigned long MAX_NUMBER = 15485863;
vector numbers(MAX_NUMBER+1);
unsigned long l1,l2;

// thread specific arguments
{
unsigned long n; // data to sieve
bool bw; // backward or forward search
};

void init_numbers(vector * numbers);

void *prime_sieve(void *arg);

void init_numbers(vector * numbers)
{
unsigned long i;
(*numbers)[2]=true;
for (i=3;ibw))
{
while (l1 <= l2)
{
if (numbers[l1] != false)
{
c = l1*l1;
pos_move = (l1 << 1);
while (c n)
{
numbers[c]=false;
c += pos_move;
}
}
l1++;
}
}
else
{
while(l2 >= l1)
{
if (numbers[l2] != false)
{
c = l2*l2;
pos_move = (l2 << 1);
while (c n)
{
numbers[c]=false;
c += pos_move;
}
}
l2–;
}
}
}

unsigned long get_primes(vector numbers, unsigned long n)
{
unsigned long primes=0;
for(unsigned int i=0;i<=n;i++)
{
primes+=numbers[i];
}
return primes;
}

int main(int argc, char *argv[])
{
int rc;
unsigned long n = MAX_NUMBER;
init_numbers(&numbers);

l1 = 2;
l2 = (int)(sqrt(n))+1;

void* retval;

cout << get_primes(numbers, n) << endl;
}

33. Ikenna said

Using python:

from math import *

def sieve2(n):
index = 0
numberList = range(3,n+1,2)
while numberList[index] = sqrt(n):
break
return [2] + numberList

34. Ikenna said

^
|
|

Didn’t post correctly

35. j0sejuan said
```package main

import "fmt"

func PrimesLessEqual(n int) [] int {
c, p, w, l := 0, 2, make([] bool, n + 1), make([] int, n)
for ; p <= n; p++ {
for p <= n && w[p] { p++ }
if p <= n {
l[c] = p; c++
for k := p; k <= n; k += p { w[k] = true }
}
}
return l[0:c]
}

func main() {
fmt.Println(len(PrimesLessEqual(15485863)))
}
```
36. dchild said

compact ruby implementation:

```list = (3..ARGV[0].to_i).step(2).to_a

(0..(Math.sqrt(list.last)/2)-1).each do |i|
(((list[i]**2)/2) - 1..list.size).step(list[i]) { |j| list[j] = 0 } if (list[i] > 0)
end

puts ([2] + list.select {|x| x > 0}).size
```
37. A Scala solution:

```import scala.math._
import scala.collection.mutable._

object Main extends Application {
val limit = 15485863 //20
val limitSquareRoot = pow(limit,0.5) toInt

def sieve(factor: Int, xs: Array[Int]): Array[Int] = {
for (i <- (factor * factor) until (limit + 1) by factor) {
xs.update((i-2), -1)
}
xs
}

def sieveOfEtratosthenes: Array[Int] = {
def numbersArray(start: Int, end: Int): Array[Int] = start until end toArray
val xs = numbersArray(2, limit + 1)
for (i  x != -1}
}
println(sieveOfEtratosthenes.length)// foreach println)
}

Main
```
38. Previous post did not copy properly:

```import scala.math._
import scala.collection.mutable._

object Main extends Application {
val limit = 15485863 //20
val limitSquareRoot = pow(limit,0.5) toInt

def sieve(factor: Int, xs: Array[Int]): Array[Int] = {
for (i <- (factor * factor) until (limit + 1) by factor) {
xs.update((i-2), -1)
}
xs
}

def sieveOfEtratosthenes: Array[Int] = {
def numbersArray(start: Int, end: Int): Array[Int] = start until end toArray
val xs = numbersArray(2, limit + 1)
for (i  x != -1}
}
println(sieveOfEtratosthenes.length)// foreach println)
}

Main
```
39. Ramapriya said

#include
#include
#include
#include
#include
unsigned LIMIT 1000000000;
int main()
{
int all[LIMIT],i,j=2,k=0,lim;
float start,stop;
printf(“Enter the limit(excluding 0)”);
scanf(“%d”,&lim);
for(i=0;i<lim;i++)
{
all[i]=i;
}
all[0]=0;
all[1]=0;
stop=(unsigned)sqrt(lim);
while(j<stop)
{
start=0;
if(all[j]!=0)
{
start=all[j]*all[j];
if(all[start]%all[j]==0)
{
all[start]=0;
start++;
}
}
j++;
}
for(i=2;i<lim-2;i++)
{
if(all[i]!=0)
{
printf(" %d",all[i]);
}
}
getch();
return 0;
}

40. Ramapriya said

Corrected solution :) very simple works like a charm

#include “stdafx.h”
#include
#include
#include
#include
#include
int main()
{
int all[100],i,j=2,k=0,lim,start;
float stop;
printf(“Enter the limit(excluding 0)”);
scanf(“%d”,&lim);
lim=lim+1;
for(i=0;i<=lim;i++)
{
all[i]=i;
}
all[0]=0;
all[1]=0;
stop=sqrt((float)lim);
while(j<stop)
{
start=0;
if(all[j]!=0)
{

start=all[j]*all[j];
while(start<=lim)
{
all[start]=0;
start=start+j;
}
}
j++;
}
for(i=2;i<=lim-2;i++)
{
if(all[i]!=0)
{
printf(" %d",all[i]);
}
}
getch();
return 0;
}

41. Ramapriya said

now if any of you guys take a look here how can you make this work for a billionth using the constant did not give me result,are linked list only solution???

42. Evgeni said

C# solution:

public class Primer
{
public int[] GetPrimes(int max)
{
if (max < 0) max = max * -1;
if (max 0)
{
for (int i = start * start; i < primeIndexes.Length; i += start)
{
primeIndexes[i] = true;
}
start++;
}

List primes = new List(max / 10);
for (int i = 1; i primeIndexes.Length)
return 0;

for (int i = minIndex; i < primeIndexes.Length; i++)
{
if (primeIndexes[i] == false)
return i;
}
return 0;
}
}

Just over a second on Core2 Duo.

43. Evgeni said

Hm, it didn’t format properly, here’s a gist:

44. Andri Juanda said

Hi, just learn python and still beginner level, so I don’t know

def sieve(number):
array = [False,False] + [True] * (number – 1)
result = []

prime_found = 0
index = 2
while index <= number:
if array[index]:
prime_found = index
result.append(index)
index *= index
while index <= number:
array[index] = False
index += prime_found
index = prime_found + 1
else:
index += 1
return result

But I remember in java, an array of boolean is really cheap, only 1 bit per element. and if I switch so False is Boolean while True is not or not examined yet, this will probably safe some space and array creation time.

45. Andri Juanda said

Sorry, Just disappointed with the way my code shown here, so I just try this to see how to format my code appropriately:

def just_a_test(‘my apology’):
print ‘Again, really sorry’

46. Colin Williams said

Java solution:

```import java.util.*;

class Sieve {
public static void main(String[] args) {
for (Integer prime : sieveFor(Integer.parseInt(args[0]))) {
System.out.println(prime.toString());
}
}

public static Integer[] sieveFor(int maximum) {
return (new Sieve(maximum)).primes();
}

private int maximum;
public Sieve(int maximum) {
this.maximum = maximum;
}

public Integer[] primes() {
TreeSet<Integer> primes = allNumbers();

Integer prime = primes.first();
while (prime <= largestSievableNumber()) {
prime = primes.higher(prime);

Integer potential_composite = primes.higher(firstUnsievedComposite(prime) - 1);
while (potential_composite != null) {
if (potential_composite % prime == 0) {
primes.remove(potential_composite);
}
potential_composite = primes.higher(potential_composite);
}
}

return primes.toArray(new Integer[0]);
}

private TreeSet<Integer> allNumbers() {
TreeSet<Integer> numbers = new TreeSet<Integer>();

for (int i = 3;i <= this.maximum;i += 2) {
}

return numbers;
}

private Integer largestSievableNumber() {
return new Integer((int) Math.pow(this.maximum, 0.5));
}

private ArrayList<Integer> compositesFor(int prime, int maximum) {
ArrayList<Integer> composites = new ArrayList<Integer>();
for (int i = firstUnsievedComposite(prime);i <= this.maximum;i += prime) {
}

return composites;
}

private Integer firstUnsievedComposite(int prime) {
return new Integer((int) Math.pow(prime, 2));
}
}```

jUnit Tests:

```import static org.junit.Assert.*;
import org.junit.Test;

public class SieveTest {
@Test
public void even() {
assertArrayEquals(new Integer[]{2}, sieveFor(2));
}

@Test
public void firstOdd() {
assertArrayEquals(new Integer[]{2, 3}, sieveFor(3));
assertArrayEquals(new Integer[]{2, 3}, sieveFor(4));
}

@Test
public void compositeOdd() {
assertArrayEquals(new Integer[]{2, 3, 5, 7}, sieveFor(9));
}

@Test
public void largeList() {
assertArrayEquals(new Integer[]{2, 3, 5, 7, 11, 13, 17, 19, 23, 29}, sieveFor(30));
}

private Integer[] sieveFor(int maximum) {
return Sieve.sieveFor(maximum);
}
}```

Output on 15485863 took a while, although I didn’t time it. Does anybody have some performance suggestions, my java’s a bit rusty.

47. johnw said

My solution in C using all optimizations. It prints primes as it finds them, and then once it hits the square root mark prints any remaining untouched variables in the array. Chews through N=15485863 in about 3.3 seconds.

48. blutorange said

In lua, not efficient, but short:

function sieve(n,s)
for i=2,math.ceil(math.sqrt(n)) do
for j=2*i,n,i do s[j]=false end
end
return s
end
s=sieve(15485863,{false,false})
for i=1,15485863 do if s[i]==nil then print(i) end end

Well, it “only” takes only 5.5s for 15485863 (and 262MB RAM)…

49. blutorange said

function sieve(n,s)
for i=2,math.ceil(math.sqrt(n)) do
if (i==2) or (i%2==1) then for j=i*i,n,(i%2+1)*i do s[j]=false end end
end
return s
end
s=sieve(15485863,{false,false})
for i=1,15485863 do if s[i]==nil then print(i) end end

(a few improvements, now it only takes 2.3s [and 1.5s for printing…])

50. Chris said

Just stumbled across this page. I have written a multi-threaded nth-prime algorithm in C# (it finds the millionth prime in 7 milliseconds using 16 threads on a 4 core – 8 hyper-threaded core – i7 3770s running at 3.1 GHz while using only 405 KB for the sieve) and was looking for information on how to convert it to racket. I’m new to racket and still trying to wrap my head around the functional approach to the sieve.

The various examples in functional languages should give me what I was looking for. Thanks to everyone for sharing.

If anyone is interested in the multi-threaded algorithm, let me know (I will be notified of follow up comments). Alternatively, as it appears that this page is all about coding challenges, how about a multi-threaded sieve challenge? I also want to convert my algo to c and compile it with clang to see what kind of boost clang’s optimers might provide vs. compiling with gcc (and vs. my visual studio compiled c#).

Note that my algorithm is an nth-prime algorithm – the inverse of what was requested on this page. In other words, given the input 1000000, my algo outputs 15485863.

There are three reasons I want to convert it to racket:

1. To use racket’s native big integers (the C# version begins to return clipped answers somewhere in the neighborhood of the 240 millionth prime).
2. To see what kind of performance penalty racket’s native big integers incurr (a simple tail recursive factorial in wicked fast in racket – faster than any big integer factorial I can find in c#).
3. It’s fun to learn new languages.

51. This is my PHP version that calculates half a million before the terminal crashes, don’t have a local copy of PHP to test with at the minute so using a remote low powered box, had to chunk the number sequences as a result of the low performing box.

52. Mike said

Python / Numpy solution

```import itertools as it
import numpy as np

def sieve(n):
if n > 2:
yield 2

if n > 3:
limit = int(n ** 0.5)
size = (n - 1) // 2
flags = np.ones([size], dtype=bool)
primes = it.compress(it.count(3, 2), flags)

for p in primes:
yield p
if p > limit:
yield from primes
break
flags[(p*p - 3) // 2::p] = False

len(list(sieve(15485863)))
```
53. […] start at p^2 instead of p and the outer loop can stop at the square root of n. I’ll leave the optimized version for you to work […]