Sieve of Atkin C++ Implementation does not tag some Primes

Question

I implemented the Sieve of Atkin in C++ to return a Vector of type bool, but it does not tag some Primes.

// Example program
#include <iostream>
#include <vector>

std::vector<bool> listPrimes(int limit){
  std::vector<bool> primes(limit);

  primes[2] = primes[3] = true;

  for(int i=1; i*i < limit; ++i){
    for(int j=1; j*j < limit; ++j){

      int n = (4*i*i) + (j*j);
      if (n <= limit && (n % 12 == 0 || n % 12 == 5 ))
        primes[n] = !primes[n];

      n = (3*i*i) + (j*j);
      if (n <= limit && n % 12 == 7 )
        primes[n] = !primes[n];

      n = (3*i*i) - (j*j);
      if ( i > j && n <= limit && n % 12 == 11 )
        primes[n] = !primes[n];
    }
  }

  for(int i=5; i*i < limit; ++i ){
    if(primes[i])
      for(int j=i*i; j < limit; j+=i*i)
          primes[i] = false;
  }

  return primes;
}

int main()
{
  std::vector<bool> primes = listPrimes(100);

  for(int i=0; i < 100; ++i)
    if(primes[i])
        std::cout << i << ", ";

    return 0;
}

Here is the output of the Given Code. 2, 3, 11, 17, 19, 23, 29, 31, 41, 43, 47, 53, 59, 67, 71, 72, 79, 83, 89,

What am I doing wrong ?

Answer 1

You have 3 mistakes (typos) in your code:

1. Replace n <= limit with n < limit everywhere.

2. Replace n % 12 == 0 with n % 12 == 1.

3. Replace primes[i] = false; with primes[j] = false;

After fixing all 3 mistakes above your algorithm works totally correctly! Fixed code is at the end of my post, copy-paste listPrimes() function only from it.

In order to check correctness of your algorithm, I wrote a special (quite advanced) C++ program located at the end of my post.

I compare your results with alternative primes computation algorithm Sieve of Eratosthenes that I implemented in a separate function.

More than that I check not only one limit value but ALL limit values till 2^14, and partially (with some step) limit values till 2^27. So I test really a lot of examples!

And all of them (test examples) are passing after those 3 fixes above. So I can assume then your algorithm is correct.

As I understand your code implements alternative way of computing Atkin. The one Atkin Wiki provides totally different way of computation. Instead of n % 12 it uses set of reminders modulus 60 (n % 60).

But surprisingly your version works too! Probably it is a simplified (more slow) version of Atkin.

My code outputs to console progress of checking limit till some power of 2 and show if all limits are checked or with some step. See example of console output located after the code below.

Try it online!

#include <iostream>
#include <iomanip>
#include <vector>
#include <stdexcept>
#include <string>
#include <sstream>
#include <cstdint>
#include <cmath>

#define ASSERT_MSG(cond, msg) { if (!(cond)) throw std::runtime_error("Assertion (" #cond ") failed at line " + std::to_string(__LINE__) + "! Msg: '" + std::string(msg) + "'"); }
#define ASSERT(cond) ASSERT_MSG(cond, "")

using std::size_t;

std::vector<bool> listPrimes(int limit){
  std::vector<bool> primes(limit);

  primes[2] = primes[3] = true;

  for(int i=1; i*i < limit; ++i){
    for(int j=1; j*j < limit; ++j){

      int n = (4*i*i) + (j*j);
      if (n < limit && (n % 12 == 1 || n % 12 == 5 ))
        primes[n] = !primes[n];

      n = (3*i*i) + (j*j);
      if (n < limit && n % 12 == 7 )
        primes[n] = !primes[n];

      n = (3*i*i) - (j*j);
      if ( i > j && n < limit && n % 12 == 11 )
        primes[n] = !primes[n];
    }
  }
  
  for(int i=5; i*i < limit; ++i ){
    if(primes[i])
      for(int j=i*i; j < limit; j+=i*i)
          primes[j] = false;
  }

    //std::cout << "Limit " << primes.size() << ": "; for (size_t i = 0; i < primes.size(); ++i) if (primes[i]) std::cout << i << ", "; std::cout << std::endl;
  return primes;
}

std::vector<bool> SieveOfEratosthenes(size_t end) {
    std::vector<bool> primes(end, true);
    primes[0] = primes[1] = false;
    for (size_t i = 0;; ++i) {
        if (i * i >= end)
            break;
        if (!primes[i])
            continue;
        for (size_t j = i * i; j < end; j += i)
            primes[j] = false;
    }
    //std::cout << "Limit " << primes.size() << ": "; for (size_t i = 0; i < primes.size(); ++i) if (primes[i]) std::cout << i << ", "; std::cout << std::endl;
    return primes;
}

int main() {
    try {
        for (size_t end = 4; end <= (1ULL << 27); ++end) {
            bool const is_all = end <= (1 << 14);
            size_t const step_bits = end <= (1 << 18) ? 8 : end <= (1 << 22) ? 16 : end <= (1 << 24) ? 21 : 25;
            if ((end & (end - 1)) == 0)
                std::cout << "Checked " << (is_all ? "all" : "step 2^" + std::to_string(step_bits))
                    << " till 2^" << std::llround(std::log2(end)) << std::endl << std::flush;
            if (!is_all && end % (1 << step_bits) != 0)
                continue;
            auto const primes_ref = SieveOfEratosthenes(end);
            auto const primes_op = listPrimes(end);
            ASSERT_MSG(primes_ref == primes_op, "Failed for limit " + std::to_string(end) +
                (primes_ref.size() != primes_op.size() ? ", size difference " + std::to_string(primes_ref.size()) +
                    " (ref) vs " + std::to_string(primes_op.size()) + " (OP)" : "") + ", diff: " +
                [&]{
                    std::ostringstream ss;
                    for (size_t j = 0; j < std::min(primes_op.size(), primes_ref.size()); ++j)
                        if (primes_op[j] != primes_ref[j]) {
                            if (primes_op[j])
                                ss << "extra " << j << ", ";
                            else
                                ss << "absent " << j << ", ";
                        }
                    return std::string(ss.str());
                }());
        }
        return 0;
    } catch (std::exception const & ex) {
        std::cout << "Exception: " << ex.what() << std::endl;
        return -1;
    }
}

Console output:

Checked all till 2^2
Checked all till 2^3
Checked all till 2^4
Checked all till 2^5
Checked all till 2^6
Checked all till 2^7
Checked all till 2^8
Checked all till 2^9
Checked all till 2^10
Checked all till 2^11
Checked all till 2^12
Checked all till 2^13
Checked all till 2^14
Checked step 2^8 till 2^15
Checked step 2^8 till 2^16
Checked step 2^8 till 2^17
Checked step 2^8 till 2^18
Checked step 2^16 till 2^19
Checked step 2^16 till 2^20
Checked step 2^16 till 2^21
Checked step 2^16 till 2^22
Checked step 2^21 till 2^23
Checked step 2^21 till 2^24
Checked step 2^25 till 2^25
Checked step 2^25 till 2^26
Checked step 2^25 till 2^27