7
$\begingroup$

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?

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

2 Answers 2

7
$\begingroup$
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}

$\endgroup$
1
  • $\begingroup$ 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! $\endgroup$
    – Ted Ersek
    Commented Mar 9, 2021 at 13:15
6
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.