3

Just for fun, I'm trying to implement the pseudocode from this StackOverflow answer for the highly factorized Sieve of Eratosthenes in C++. I can't figure out why my code returns both prime and non-prime numbers. Am I implementing these for loops incorrectly? Should I be using while loops instead? I suspect that I'm not incrementing the for loops properly. Any help would be greatly appreciated. I've spent several hours trying to hunt down the flaw.

GordonBGood's pseudocode is inserted as comments, and I've used all the same variable names.

#include <iostream>
#include <vector>
#include <cmath>

const int limit = 1000000000;
const std::vector<int> r {23,29,31,37,41,43,47,53, 59,61,67,71,73,79,83, //positions + 19
                      89,97,101,103,107,109,113,121,127, 131,137,139,
                      143,149,151,157,163,167,169,173,179,181,187,191,193, 
                      197,199,209,211,221,223,227,229};

int main()
{
// an array of length 11 times 13 times 17 times 19 = 46189 wheels initialized
// so that it doesn't contain multiples of the large wheel primes
// for n where n ← 210 × w + x where w ∈ {0,...46189}, x in r: // already
//    if (n mod cp) not equal to 0 where cp ∈ {11,13,17,19}: // no 2,3,5,7
//        mstr(n) ← true else mstr(n) ← false                // factors
  std::vector<bool> mstr(limit);
  int n;
  for (int w=0; w <= 46189; ++w) {
    for (auto x = begin(r); x != end(r); ++x) {
      n = 210*w + *x;
      if (n % 11 != 0 && n % 13 != 0 && n % 17 != 0 && n % 19 != 0)
        mstr[n]=true;
      else
        mstr[n]=false; 
    }
  }

// Initialize the sieve as an array of the smaller wheels with
// enough wheels to include the representation for limit
// for n where n ← 210 × w + x, w ∈ {0,...(limit - 19) ÷ 210}, x in r:
//    sieve(n) ← mstr(n mod (210 × 46189))    // init pre-culled primes.
  std::vector<bool> sieve(limit+1000);
  for (int w=0; w <= (limit-19)/210; ++w) {
    for (auto x = begin(r); x != end(r); ++x) {
      n = 210*w + *x;
      sieve[n] = mstr[(n % (210*46189))];
    }
  }

// Eliminate composites by sieving, only for those occurrences on the
// wheel using wheel factorization version of the Sieve of Eratosthenes
// for n² ≤ limit when n ← 210 × k + x where k ∈ {0..}, x in r
//    if sieve(n):
//        // n is prime, cull its multiples
//        s ← n² - n × (x - 23)     // zero'th modulo cull start position
//        while c0 ≤ limit when c0 ← s + n × m where m in r:
//            c0d ← (c0 - 23) ÷ 210, cm ← (c0 - 23) mod 210 + 23 //means cm in r
//            while c ≤ limit for c ← 210 × (c0d + n × j) + cm
//                    where j ∈ {0,...}:
//                sieve(c) ← false       // cull composites of this prime  
  int s, c, c0, c0d, cm, j;
  for ( auto x = begin(r); x != end(r); ++x) {
    for ( int k=0; (n=210*k + (*x)) <= sqrt(limit); ++k){
      if (sieve[n]) {
        s = n*n - n*((*x)-23);
        for ( auto m = begin(r); (c0=s+n*(*m)) <= limit &&  m != end(r); ++m) {
      c0d = (c0-23)/210;
      cm  = (c0-23)%210 + 23;
      for ( int j=0; (c=210*(c0d+n*j)+cm) <= limit; ++j) {
        sieve[c] = false;
          }
        } 
      } 
    } 
  } 

// output 2, 3, 5, 7, 11, 13, 17, 19,
// for n ≤ limit when n ← 210 × k + x where k ∈ {0..}, x in r:
//     if sieve(n): output n  
  std::cout << "2\n3\n5\n7\n11\n13\n17\n19\n";
  for ( auto x = begin(r); x != end(r); ++x) {
    for ( int k = 0; (n=210*k + (*x)) <= limit; ++k) {
      if (sieve[n]);
    std::cout << n << std::endl;
    }
  }
  std::cout << std::endl;

  return 0;
}
Community
  • 1
  • 1
tn3rt
  • 240
  • 2
  • 13
  • Knowing how to use a debugger is a very important skill when programming, especially in languages like C++. Using a debugger, you can step through the code, line by line, while being able to see values of variables (including all values in a non-dynamic array) as well as test out sub-expressions and even full function calls. It's a lot of work to go through code with large number of loops, especially if they are nested, or recursive code, but it will help you understand what happens in your code and what the problem might be. – Some programmer dude Jun 14 '15 at 14:13
  • 3
    Always compile with a high warning level. If you do so, it should warn you about the extra semicolon you have at the end of `if (sieve[n]);`, which makes the `if` useless. – ooga Jun 14 '15 at 14:36
  • Wow. I can't believe I didn't catch that semicolon on the if statement. Thank you for spotting that. For the record, the -Wall flag didn't catch it but the -Wextra flag did. There's still something inherently wrong with the code though. It's leaving out prime numbers and then printing some numbers that aren't prime. I'll have to take another look at the code with gdb. – tn3rt Jun 14 '15 at 15:09
  • [here's a short Haskell wheel-based sieve code](http://pastebin.com/tahtbqWW), maybe it could serve as a pseudocode for you. – Will Ness Jun 17 '15 at 17:39
  • Thank you. I will look over it and study it carefully. – tn3rt Jun 17 '15 at 17:43
  • HTH. here's [a newer version](http://pastebin.com/ZUxxuLV9) of it. – Will Ness Jun 18 '15 at 03:48

0 Answers0