# Compiling Eratosthenes Sieve

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]] Commented Mar 6, 2021 at 18:57
• See the answer by @ChipHurst for several examples of compiled prime sieves. Commented Mar 6, 2021 at 20:25

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}

• 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! Commented Mar 9, 2021 at 13:15

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.