筛子阿特金C++实现没有标记某些Primes

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

本文关键字:Primes 实现 C++      更新时间:2023-10-16

我在C++中实现了 Sieve of Atkin 以返回布尔类型的向量,但它没有标记一些素数。

// 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;
}

这是给定代码的输出。 2, 3, 11, 17, 19, 23, 29, 31, 41, 43, 47, 53, 59, 67, 71, 72, 79, 83, 89,

我做错了什么?

代码中有 3 个错误(拼写错误):

  1. 随处用n < limit替换n <= limit

  2. n % 12 == 0替换为n % 12 == 1

  3. primes[i] = false;替换为primes[j] = false;

修复上述所有 3 个错误后,您的算法完全正确运行!固定代码位于我帖子的末尾,复制粘贴listPrimes()仅从中起作用。

为了检查您的算法的正确性,我编写了一个特殊的(相当高级的)C++程序,位于我的帖子末尾。

我将您的结果与我在单独的函数中实现的替代素数计算算法埃拉托色尼筛进行比较。

不仅如此,我不仅检查一个limit值,而且检查所有极限值直到 2^14,以及部分(带有一些步骤)极限值直到 2^27。所以我测试了很多例子!

所有这些(测试示例)都在上述 3 个修复之后通过。所以我可以假设你的算法是正确的。

据我了解,您的代码实现了计算 Atkin 的替代方法。一个Atkin维基提供了完全不同的计算方式。它不是n % 12而是使用一组提醒模数 60 (n % 60)。

但令人惊讶的是,您的版本也可以工作!可能是阿特金的简化(更慢)版本。

我的代码输出到控制台检查limit直到 2 的幂的进度,并显示是否检查了所有限制或具有某个步骤。请参阅位于以下代码后面的控制台输出示例。

在线试用!

#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;
}
}

控制台输出:

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