# Sieve Of Eratosthenes In Many Programming Languages

Challenges:

Remember: ConvertSpacesToTabsNotForCode

The code at http://www.quitebasic.com/prj/math/eratosthenes/ marks composite numbers while STEPping through an array, not using MOD.

``` 2>:3g" "-!v\  g30          <
|!`"O":+1_:.:03p>03g+:"O"`|
@               ^  p3\" ":<
2 234567890123456789012345678901234567890123456789012345678901234567890123456789

```

See the SieveOfEratosthenes page for an example in C.

A version where a bit vector is marked and then its indices are listed:

``` (defun sieve (n)
(let ((bit-vector (make-array n :initial-element 1 :element-type 'bit)))
(loop for i from 2 upto (isqrt n) do
(loop for j from i
for index = (* i j)
until (>= index n) do
(setf (sbit bit-vector index) 0)))
(loop for i from 2 below (length bit-vector)
unless (zerop (sbit bit-vector i)) collect i)))

```

The FalseLanguage page has a couple one-liner PrimeNumber finders.

Three fast (no MODs) versions here: http://wiki.forthfreak.net/index.cgi?SieveOfEratosthenes

``` combine op (i1:l1) (i2:l2) = (op i1 i2):(combine op l1 l2)
ns n = n:(ns n)
nmults n = n:(combine (+) (nmults n) (ns n))
diff (i1 : l1) (i2 : l2) | i1 < i2  = i1:(diff     l1 (i2:l2))
| i1 == i2 =    (diff     l1     l2 )
| i1 > i2  =    (diff (i1:l1)    l2 )
filtermults (li:l) = li:(filtermults (diff l (nmults li)))
primes = filtermults [2..]

Main> take 10 primes
[2,3,5,7,11,13,17,19,23,29]
```
This solution uses no multiplication or division, only addition. For each prime, a stream containing all of its multiples is generated, and this stream is subtracted from the integers. The resulting stream contains only the primes.

We must use custom-defined diff because the built-in notElem or the infix \\ "remove" operator won't terminate on an infinite list. -- JonathanTang

The above duplicates a lot of the Prelude functions. Changing variable names we get a one-liner, except for the diff:

``` sieve (p:ns) = p : sieve (diff ns [p, p+p..])
primes = sieve [2..]
```
A true one-liner is DavidTurner's
``` sieve (p:ns) = p : sieve [n | n <- ns, rem n p /= 0]
```
where diff is replaced with the standard ListComprehension syntax and rem divisibility testing. But it is a TrialDivision? algorithm i.e. its theoretical complexity is worse. It is essential for the Sieve of Eratosthenes that it doesn't test for composites but rather directly generates them through iterated addition. This is because the testing of each candidate number must go through all the primes up to its square root, whereas with direct generation each composite is generated only from its prime factors up to its square root.

Both of the above versions start removing the multiples of each produced prime immediately, but actually need to do that only from the prime's square. Thus both of them create many superfluous removal streams. This makes them very slow and of quadratic time complexity, or worse, in n primes produced.

The bounded formulation stops early and thus does not make any unneeded removal steps:

``` primesTo m = sieve [2..] where
sieve (p:ns)
| p*p > m   = p : takeWhile (<= m) ns       -- all ns up to p*p are prime already
| otherwise = p : sieve (diff ns [p*p, p*p+p..])
```
Ideally, each multiple should be removed in O(1) time, or equivalently, O(|b|) for the whole (diff a b) removal step. This is achievable on mutable RandomAccess? arrays. Unfortunately, on SinglyLinkedLists it takes O(|union a b|) time. With BinaryTrees? it would have been O(|b|*log(|a|)). Still, the amount of removal steps is strictly limited to what is necessary. Even though each diff step starts prematurely (right after p instead of at p*p) this code runs much faster than the first one and with EmpiricalOrdersOfGrowth of about ~ n^1.4 in n primes produced.

The unbounded formulation is fixed by postponing the sieving,

``` sieve ns (p:ps) | (h,t) <- span (< p*p) ns = h ++ sieve (diff t [p*p, p*p+p ..]) ps
primes = 2 : sieve [3..] primes
```
In Richard Bird's version, the above (...(((ns-a)-b)-c)-...) workflow is rearranged into (ns-(a+(b+(c+...)))):
``` primes = 2 : diff [3..] (foldr (\p r-> p*p : add [p*p+p, p*p+2*p..] r) [] primes)
| x==y  = x : add    xs    ys
| x>y   = y : add (x:xs)   ys
```
With tree-shaped folding and telescopic multistage recursive production of primes (to minimize memory usage) we get
``` primes = 2 : _Y ((3:) . gaps 5 . _U . map (\p-> [p*p, p*p+2*p..]))
_Y g = g (_Y g)                                  -- non-sharing (recursive) fixpoint combinator
_U ((x:xs):t) = x : (add xs . _U . pairs) t      -- bigU  == sort . concat
pairs (a:b:t) = add a b : pairs t
gaps x (y:ys) | x < y = x : gaps (x+2) (y:ys)    -- gaps 5 == ([5,7..] \\)
| otherwise = gaps (x+2)    ys     -- x==y
```
This runs at about ~ n^1.2 EmpiricalOrdersOfGrowth, in n primes produced, and in practically constant space. WheelFactorization optimization can be further applied. More at http://haskell.org/haskellwiki/Prime_numbers. -- WillNess

``` function sieve(max) {
var D = [], primes = [];
for (var q=2; q<max; q++) {
if (D[q]) {
for (var i=0; i<D[q].length; i++) {
var p = D[q][i];
if (D[p+q]) D[p+q].push(p);
else D[p+q]=[p];
}
delete D[q];
} else {
primes.push(q);
if (q*q<max) D[q*q]=[q];
}
}
return primes;
}
sieve(100)

```

Haskell's a fun language, but this seems like a poor choice for comparing it to other languages. Take one of the worst: perl; this takes just under 1 second on my athlon-xp to generate the first 78498 primes; how does your Haskell program compare? I suspect that its performance is "pathetic".

Although I've timed neither the haskell or your code, I'd disagree. The haskell code (to me, at least) indicates very clearly what it is doing while being quite concise. This, combined with the fact that I find it unlikely that anybody is too concerned with quickly generating the first few prime numbers makes it seem to me that in this context this haskell code is a fine example, especially since my understanding of this page is that we're looking for clear and idiomatic implementations instead of optimized. (Of course, the performance of the haskell code could be easily doubled without a major sacrifice in clarity, but the fact that it was not indicates to me that its purpose is clarity, which should be the first intent in most code)

• Well, naively, I tried timing the code above using ghc -O foo.hs, and after more than thirty minutes on a 2.8GHz, dual-core box using GHC 6.8.5, it's still waiting to show the result of main = do putStr \$ show \$ sieve [2..] !! 78498. Even doubling the performance isn't sufficient to come anywhere close to the Perl developer's claimed performance. You're going to need a monadic implementation of this algorithm, even naively, to even come close to the Perl implementation in run-time. I also question the code's legibility as well. It reads too much like an APL one-liner to me. That being said, I do still find Haskell's one-line brute-force method to be much easier to read and understand than even the multi-lined Perl copy below.
• Yes, that first Haskell code version above is very slow, has quadratic complexity in n primes produced, and should be taken as specification only. Tested ghc-compiled at IdeOne (http://ideone.com/P0E81) it can't even produce 20,000 primes in 15 seconds; yet the 2nd, bounded and guarded version primesTo has complexity below O(n^1.5) and takes just 0.37 seconds on Ideone.com to calculate all the 78498 primes below 1,000,000, as I write this in April 2011. -- WillNess (small copyedits, Jan 2012).
• Strongly agreed. This example is a great fit for a functional approach.

``` \$n=1000000;for(\$t=3;\$t*\$t<\$n;\$t+=2){if(!\$c[\$t]){for(\$s=\$t*\$t;\$s<\$n;\$s+=\$t*2){\$c[\$s]++}}}
print "2\n";for(\$t=3;\$t<\$n;\$t+=2){\$c[\$t]||print "\$t\n"}
```
or

``` \$n=1000000;
for(\$t=3; \$t*\$t<\$n; \$t+=2) {
if (!\$c[\$t]) {
for(\$s=\$t*\$t; \$s<\$n; \$s+=\$t*2) { \$c[\$s]++ }
}
}
print "2\n";
for(\$t=3; \$t<\$n; \$t+=2) {
\$c[\$t] || print "\$t\n";
}
```
I'm confused. Where does \$c come from? Where is @c defined?

What you're seeing is AutoVivification and the lack of 'use strict;'.

Based on a recipe in the PythonCookbook:

``` from __future__ import generators

def eratosthenes(maxnum):
D = {}  # map composite integers to primes witnessing their compositeness
q = 2   # first integer to test for primality
while q <= maxnum:
if q not in D:
yield q       # not marked composite, must be prime
D[q*q] = [q]  # first multiple of q not already marked
else:
for p in D[q]:  # move each witness to its next multiple
D.setdefault(p+q,[]).append(p)
del D[q]        # no longer need D[q], free memory
q += 1

for p in eratosthenes(19): print p,
# outputs 2 3 5 7 11 13 17 19

```
Animated version, using the GuidoVanRobot toolkit:

``` import powerMode

def main():
for num in range(2, 101):
x = num % 10
y = num / 10
if x == 0:
x += 10
y -= 1
world.setBeepers(x, y+1, num)
world.positionRobot(1,1, 'E')
setDelayAmount(0.05)
beginDisplay(11, 11)
clearMultiples(2)
clearMultiples(3)
clearMultiples(5)
clearMultiples(7)
celebrate()

def clearMultiples(factor):
num = 1
while num <= 100:
if num > factor and num % factor == 0:
world.beepers[world.robot] = 0
pause(0.4)
move()
if num % 10 == 0:
turnaround()
zipTenSpaces()
turnright()
move()
turnright()
num += 1
turnright()
zipTenSpaces()
turnleft()

def turnaround():
turnleft()
turnleft()

def zipTenSpaces():
for i in range(10):
setDelayAmount(0.01)
move()
setDelayAmount(0.05)

def turnright():
turnaround()
turnleft()

def spin():
turnaround()
turnaround()

def celebrate():
while 1:
turnleft()

main()

```

``` def sieve(upper)
nums = (2..upper).to_a
iend = nums.size - 1
sqrt = Math.sqrt(upper).to_i
(0..(nums.index(sqrt))).each do |i|
n = nums[i]    # n is either prime or nil
(i + n).step(iend, n) { |j| nums[j] = nil } unless n.nil?
end
nums.compact     # eliminate nil entries
end

sieve(19)          # returns [2, 3, 5, 7, 11, 13, 17, 19]

```

SchemeLanguage: (Using SiCp style streams)

``` ;;;; Stream Implementation
(define (head s) (car s))          ;; _odd_ non-memoized streams,
(define (tail s) ((cdr s)))        ;;  per SRFI-41
(define-syntax s-cons
(syntax-rules () ((s-cons h t) (cons h (lambda () t)))))

;;;; Stream Utility Functions
(define (from-By x s)
(s-cons x (from-By (+ x s) s)))
(define (take n s)
(cond
((> n 1) (cons (head s) (take (- n 1) (tail s))))
((= n 1) (list (head s)))      ;; don't force it too soon
(else ())))     ;; so (take 4 (s-map / (from-By 4 -1))) works
(define (drop n s)
(cond
((> n 0) (drop (- n 1) (tail s)))
(else s)))
(define (s-map f s)
(s-cons (f (head s)) (s-map f (tail s))))
(define (s-diff s1 s2)
(cond
((< h1 h2) (s-cons h1 (s-diff  (tail s1)       s2 )))
((< h2 h1)            (s-diff        s1  (tail s2)))
(else                 (s-diff  (tail s1) (tail s2))))))
(define (s-union s1 s2)
(cond
((< h1 h2) (s-cons h1 (s-union (tail s1)       s2 )))
((< h2 h1) (s-cons h2 (s-union       s1  (tail s2))))
(else      (s-cons h1 (s-union (tail s1) (tail s2)))))))

;;;; odd multiples of an odd prime
(define (mults p) (from-By (* p p) (* 2 p)))

;;;; The Sieve itself, bounded, ~ O(n^1.4) in n primes produced
;;;;   (unbounded version runs at ~ O(n^2.2), and growing worse)
;;;;   **only valid up to m**, includes composites above it
(define (primes-To m)
(define (sieve s)
(cond ((> (* p p) m) s)
(else (s-cons p
(sieve (s-diff (tail s) (mults p))))))))
(s-cons 2 (sieve (from-By 3 2))))

;;;; all the primes' multiples, tree-merged, removed;
;;;;    ~O(n^1.17..1.15) time in producing 100K .. 1M primes
;;;;    ~O(1) space (O(pi(sqrt(m))) probably)
(define (primes-TM)
(define (no-mults-From from)
(s-diff (from-By from 2) (s-tree-join (s-map mults odd-primes))))
(define odd-primes
(s-cons 3 (no-mults-From 5)))
(s-cons 2 (no-mults-From 3)))

;;;; join an ordered stream of streams (here, of primes' multiples)
;;;; into one ordered stream, via an infinite right-deepening tree
(define (s-tree-join sts)                               ;; sts -> s
(define (join-With of-Tail sts)                       ;; sts -> s
(s-union (tail (head sts)) (of-Tail (tail sts)))))
(define (pairs sts)                                   ;; sts -> sts
(s-cons (join-With head sts) (pairs (tail (tail sts)))))
(join-With (lambda (t) (s-tree-join (pairs t))) sts))

;;;; Print 10 last primes from the first thousand primes
(begin
(display (take 10 (drop 990 (primes-To 7919)))) (newline)
(display (take 10 (drop 990 (primes-TM)))) (newline))

```

Didn't realize it till afterwards, but I ended up with something pretty close to the Lisp solution above, but using nil instead of 0 and 1. -- RamonLeon

``` Integer>>sieveOfEratosthenes
| primes |
primes := Array new: self.
2 to: self do: [:each |
(primes at: each)
ifNil: [each + each to: self
by: each do: [:notPrime | primes at: notPrime put: notPrime]]].
^(2 to: self) select: [:each | (primes at: each) isNil]

19 sieveOfEratosthenes -> (2 3 5 7 11 13 17 19)

```

MATLAB:

function [ primes ] = sieveOfEratosthenes( N )

```  sieve=1:N;
p=2;
while p^2<=N
sieve(p^2:p:N)=0;
p=find(sieve>p,1);
end
primes = sieve(sieve>1);
```
end EditText of this page (last edited February 16, 2014) or FindPage with title or text search