Sieve Of Eratosthenes In Many Programming Languages

Challenges:

Remember: ConvertSpacesToTabsNotForCode

BasicLanguage:

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

BefungeLanguage: 

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

CeeLanguage:

See the SieveOfEratosthenes page for an example in C.

CommonLisp:

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.

ForthLanguage:

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

GameOfLife:

http://www.radicaleye.com/lifepage/patterns/primes.html

HaskellLanguage: 

 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)
 add (x:xs) (y:ys)  | x<y   = x : add    xs (y:ys)
                    | 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

JavaScript: 

 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)

PerlLanguage:

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) 

 $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;'.

PythonLanguage:

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()

RubyLanguage: 

 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)
   (let ((h1 (head s1)) (h2 (head 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)
   (let ((h1 (head s1)) (h2 (head 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) 
(let ((p (head 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-cons (head (head sts))
              (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))

SmalltalkLanguage:

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

CategoryInManyProgrammingLanguages

EditText of this page (last edited February 16, 2014) or FindPage with title or text search