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 ?
You have 3 mistakes (typos) in your code:
Replace n <= limit with n < limit everywhere.
Replace n % 12 == 0 with n % 12 == 1.
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.
#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