Compiling Eratosthenes Sieve

Question

My best implementation of Eratosthenes sieve is here.

EratosthenesSieve[n_Integer?Positive] := Module[{lst = Range[2, n]},
  Do[
    If[Positive@lst[[i]],Part[lst, Range[i^2 + 2 i, n - 1, lst[[i]]]] = 0],
    {i, Floor[Sqrt[n + 0.0]]}];
  Pick[lst, Unitize@lst, 1]]

Now we can make a list of primes < 120.

EratosthenesSieve[120]
(*{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}*)

I was able to make a compiled version of the above, and you will see it doesn't have MainEvaluate in the CompiledFunction. My compiled version is here.

Needs["CompiledFunctionTools`"];
EratosthenesSieveC=Compile[{{n,_Integer}},Module[{lst=Range[2,n]},
  Do[If[Positive@lst[[i]],
    Do[Part[lst,j]=0,{j,Range[i^2+2i,n-1,lst[[i]]]}]],{i,Floor[Sqrt[n+0.0]]}
  ];
  Complement[lst,{0}]
]];
CompilePrint[EratosthenesSieveC]

However, the compiled version is slower than my version that isn't compiled. Is there a way to make EratosthenesSieve faster?

Here's a functional version - it's slow though maybe somebody can figure out how to speed it up: EratosthenesSieve[n_Integer] := Reap[ With[{sv = {lst, x} |-> Select[lst, ! Divisible[#, x] &]}, NestWhile[(Sow[#[[1]]]; sv[#, #[[1]]]) &, Range[2, n], # != {} &]]][[2, 1]] — flinty, Mar 06 '21 at 18:57
See the answer by @ChipHurst for several examples of compiled prime sieves. — KennyColnago, Mar 06 '21 at 20:25

score 7 · Accepted Answer · answered Mar 07 '21 at 13:58

7

Clear[primes];
primes[n_] :=
  Module[{p = Range[1, n, 2]},
   p[[1]] = 2;
   Do[If[p[[(k + 1)/2]] != 0, p[[(k^2 + 1)/2 ;; ;; k]] = 0], {k, 3, n^.5, 2}];
   SparseArray[p]["NonzeroValues"]
   ];

primes[10^8] // Length // AbsoluteTiming

{0.756629, 5761455}

answered Mar 07 '21 at 13:58

chyanog

15,542
3
40
78

I tried many variations of this and others, and I can't find anything faster. The fact that it is so concise, makes it all the more elegant. Wonderful code! – Ted Ersek Mar 09 '21 at 13:15

Henrik Schumacher · Answer 2 · 2021-03-06T18:13:21.673

I do not know how to make the compiled version faster than the your high level version, but I know how to improve the latter. By using Span (;; ;;) instead of a Range, you save the time for biulding the range first. (However, one has to fiddle a bit because Span does not like it if its first argument is greater than its second.

EratosthenesSieve[n_Integer?Positive] := Module[{lst, stop},
  lst = Range[2, n];
  stop = Floor[Sqrt[n]];
  If[stop^2 + 2 stop > n - 1, stop = stop - 1];
  Do[
   If[lst[[i]] > 0, lst[[i^2 + 2 i ;; n - 1 ;; lst[[i]]]] = 0]
   , {i, 1, stop}];
  Pick[lst, Unitize@lst, 1]
  ]

This is a bit more than twice as fast on my machine.

I think the issue with Compile is that Mathematics creates a bound check for each write operations into an array. This is to prevent segfaults that might. force the kernel to quit. But these checks come at a cost. These checks might also prevent auto vectorization by the compiler.

I have just tested the algorithm in plain C++ with parallelization by OpenMP (Quad Core), and I have not been able to make more than 20% faster than the above version of EratosthenesSieve. That's understandable: the bulk of work lies in lst[[i^2 + 2 i ;; n - 1 ;; lst[[i]]]] = 0 and Mathematica uses already a very efficient library for this in the background.

Compiling Eratosthenes Sieve

2 Answers2