title: Porting optimized Sieve of Eratosthenes from Python to C++
body: <p>Some time ago I used the (blazing fast) primesieve in python that I found here: <a href="https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188">Fastest way to list all primes below N</a></p>
  
  <p>To be precise, this implementation:</p>
  
  <pre><code>def primes2(n):
      """ Input n&gt;=6, Returns a list of primes, 2 &lt;= p &lt; n """
      n, correction = n-n%6+6, 2-(n%6&gt;1)
      sieve = [True] * (n/3)
      for i in xrange(1,int(n**0.5)/3+1):
        if sieve[i]:
          k=3*i+1|1
          sieve[      k*k/3      ::2*k] = [False] * ((n/6-k*k/6-1)/k+1)
          sieve[k*(k-2*(i&amp;1)+4)/3::2*k] = [False] * ((n/6-k*(k-2*(i&amp;1)+4)/6-1)/k+1)
      return [2,3] + [3*i+1|1 for i in xrange(1,n/3-correction) if sieve[i]]
  </code></pre>
  
  <p>Now I can slightly grasp the idea of the optimizing by automaticly skipping multiples of 2, 3 and so on, but when it comes to porting this algorithm to C++ I get stuck (I have a good understanding of python and a reasonable/bad understanding of C++, but good enough for rock 'n roll).</p>
  
  <p>What I currently have rolled myself is this (<code>isqrt()</code> is just a simple integer square root function):</p>
  
  <pre><code>template &lt;class T&gt;
  void primesbelow(T N, std::vector&lt;T&gt; &amp;primes) {
      T sievemax = (N-3 + (1-(N % 2))) / 2;
      T i;
      T sievemaxroot = isqrt(sievemax) + 1;
  
      boost::dynamic_bitset&lt;&gt; sieve(sievemax);
      sieve.set();
  
      primes.push_back(2);
  
      for (i = 0; i &lt;= sievemaxroot; i++) {
          if (sieve[i]) {
              primes.push_back(2*i+3);
              for (T j = 3*i+3; j &lt;= sievemax; j += 2*i+3) sieve[j] = 0; // filter multiples
          }
      }
  
      for (; i &lt;= sievemax; i++) {
          if (sieve[i]) primes.push_back(2*i+3);
      }
  }
  </code></pre>
  
  <p>This implementation is decent and automatically skips multiples of 2, but if I could port the Python implementation I think it could be much faster (50%-30% or so).</p>
  
  <p>To compare the results (in the hope this question will be successfully answered), the current execution time with <code>N=100000000</code>, <code>g++ -O3</code> on a Q6600 Ubuntu 10.10 is 1230ms.</p>
  
  <p>Now I would love some help with either understanding what the above Python implementation does or that you would port it for me (not as helpful though).</p>
  
  <p><strong>EDIT</strong></p>
  
  <p>Some extra information about what I find difficult.</p>
  
  <p>I have trouble with the techniques used like the correction variable and in general how it comes together. A link to a site explaining different Eratosthenes optimizations (apart from the simple sites that say "well you just skip multiples of 2, 3 and 5" and then get slam you with a 1000 line C file) would be awesome.</p>
  
  <p>I don't think I would have issues with a 100% direct and literal port, but since after all this is for learning that would be utterly useless.</p>
  
  <p><strong>EDIT</strong></p>
  
  <p>After looking at the code in the original numpy version, it actually is pretty easy to implement and with some thinking not too hard to understand. This is the C++ version I came up with. I'm posting it here in full version to help further readers in case they need a pretty efficient primesieve that is not two million lines of code. This primesieve does all primes under 100000000 in about 415 ms on the same machine as above. That's a 3x speedup, better then I expected!</p>
  
  <pre><code>#include &lt;vector&gt;
  #include &lt;boost/dynamic_bitset.hpp&gt;
  
  // http://vault.embedded.com/98/9802fe2.htm - integer square root
  unsigned short isqrt(unsigned long a) {
      unsigned long rem = 0;
      unsigned long root = 0;
  
      for (short i = 0; i &lt; 16; i++) {
          root &lt;&lt;= 1;
          rem = ((rem &lt;&lt; 2) + (a &gt;&gt; 30));
          a &lt;&lt;= 2;
          root++;
  
          if (root &lt;= rem) {
              rem -= root;
              root++;
          } else root--;
  
      }
  
      return static_cast&lt;unsigned short&gt; (root &gt;&gt; 1);
  }
  
  // https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
  // https://stackoverflow.com/questions/5293238/porting-optimized-sieve-of-eratosthenes-from-python-to-c/5293492
  template &lt;class T&gt;
  void primesbelow(T N, std::vector&lt;T&gt; &amp;primes) {
      T i, j, k, l, sievemax, sievemaxroot;
  
      sievemax = N/3;
      if ((N % 6) == 2) sievemax++;
  
      sievemaxroot = isqrt(N)/3;
  
      boost::dynamic_bitset&lt;&gt; sieve(sievemax);
      sieve.set();
  
      primes.push_back(2);
      primes.push_back(3);
  
      for (i = 1; i &lt;= sievemaxroot; i++) {
          if (sieve[i]) {
              k = (3*i + 1) | 1;
              l = (4*k-2*k*(i&amp;1)) / 3;
  
              for (j = k*k/3; j &lt; sievemax; j += 2*k) {
                  sieve[j] = 0;
                  sieve[j+l] = 0;
              }
  
              primes.push_back(k);
          }
      }
  
      for (i = sievemaxroot + 1; i &lt; sievemax; i++) {
          if (sieve[i]) primes.push_back((3*i+1)|1);
      }
  }
  </code></pre>