检查埃拉托色尼移动筛的唯一性

Question

我的代码如下：

int main(){
    vector<int> primes;
    vector<int> primesSum;
    int tally = 1;
    bool primeFound = false;
    while (true) {
        int i = 0;
        while (!primeFound) {
            primeFound = true;
            tally++;
            for (i = 0; i < (int)primes.size(); i++) {
                if (tally == primesSum[i]) {
                    primeFound = false;
                    primesSum[i] += primes[i];
                }
            }
        }
        primeFound = false;
        primesSum.push_back(tally*2);
        primes.push_back(tally);
    }
}

它似乎可以很好地生成素数，但我想知道我是否可以通过实施所有素数在 2 和 3 之后都是 6n +/- 1 的规则来提高它的效率。但这似乎会牺牲我最初空向量的成就。

我见过使用 6n +/- i 规则的质数验证器，但实际上不是在筛选程序中，可能是因为它破坏了 2 过滤器。

编辑：顺便说一句，我宁愿不使用模数，因为这是该程序的另一项成就，除非模数比我目前拥有的更有效（时间方面而不是空间方面）。

Answer 1

这里是我编写的一些代码，用于使用仅考虑 +-1 mod 6 整数的筛子有效生成素数。它被编写为使用协程库的生成器，但逻辑是等效的。它使用

vector

的

uint64_t

作为位掩码。

没有与计算素数相关的模数，但它们用于访问位图。

更新

我从生成器中提取了我的原始代码，并创建了一个简单的函数，将素数放入一个向量中。一种方法使用 +-1 mod 2 搜索，另一种方法使用 +-1 mod 6 搜索。我还将 OP 的代码提取到一个函数中。

我在 M1 Max (arm64) 上运行了一些简单的时序测量。

功能	质数<1M	质数<10M	质数<100M
筛分指数	12904 毫秒	仍在运行...
sieve_2n	2 毫秒	30 毫秒	222 毫秒
sieve_6n.	1 毫秒	15 毫秒	132 毫秒

我将创建一个 GitHub 存储库来构建代码。

代码

auto sieve_mod6_prime_seq(int max = int{1} << 20) {
    std::vector<int> primes;
    primes.push_back(2);
    primes.push_back(3);

    auto max_index = max / 3;
    auto bits_per = sizeof(uint64_t) * CHAR_BIT;
    auto nwords = (bits_per + max_index - 1) / bits_per;
    std::vector<uint64_t> words(nwords);

    words[0] |= 1;
    size_t wdx = 0;
    while (wdx < nwords) {
        auto b = std::countr_one(words[wdx]);
        auto p = 3 * (64 * wdx + b) + 1 + (b bitand 1);
        if (b < 64 and p < max) {
            primes.push_back(p);

            for (auto j = p; j < max; j += 6 * p) {
                auto idx = j / 3;
                auto jdx = idx / 64;
                auto jmask = uint64_t{1} << (idx % 64);
                words[jdx] |= jmask;
            }

            for (auto j = 5 * p; j < max; j += 6 * p) {
                auto idx = j / 3;
                auto jdx = idx / 64;
                auto jmask = uint64_t{1} << (idx % 64);
                words[jdx] |= jmask;
            }
        }
        else {
            ++wdx;
        }
    }
    return primes;
}

auto sieve_mod2_prime_seq(int max = int{1} << 20) {
    std::vector<int> primes;

    auto bits_per = 2 * sizeof(uint64_t) * CHAR_BIT;
    auto nwords = (bits_per + max - 1) / bits_per;
    std::vector<uint64_t> words(nwords);
    if (max > 2)
        primes.push_back(2);
    words[0] |= 1;
    size_t idx = 0;
    while (idx < nwords) {
        auto bdx = 2 * std::countr_one(words[idx]) + 1;
        auto p = 128 * idx + bdx;
        if (bdx < 128 and p < max) {
            primes.push_back(p);
            for (auto j = p; j < max; j += 2 * p) {
                auto jdx = j / 128;
                auto jmask = uint64_t{1} << ((j % 128) / 2);
                words[jdx] |= jmask;
            }
        }
        else {
            ++idx;
        }
    }
    return primes;
}

auto sieve_index(int max = int{1} << 20) {
    std::vector<int> primes;
    std::vector<int> primesSum;
    int tally = 1;
    bool primeFound = false;
    while (tally < max) {
        int i = 0;
        while (!primeFound) {
            primeFound = true;
            tally++;
            for (i = 0; i < (int)primes.size(); i++) {
                if (tally == primesSum[i]) {
                    primeFound = false;
                    primesSum[i] += primes[i];
                }
            }
        }
        primeFound = false;
        primesSum.push_back(tally*2);
        primes.push_back(tally);
    }
    return primes;
}

template<class Work>
void measure(std::ostream& os, std::string_view desc, Work&& work) {
    chron::StopWatch timer;
    timer.mark();
    work();
    auto millis = timer.elapsed_duration<std::chrono::milliseconds>().count();
    os << fmt::format("{:>30s}: {} ms", desc, millis) << endl;
}

int tool_main(int argc, const char *argv[]) {
    ArgParse opts
        (
         argValue<'n'>("number", 1000, "Number of primes"),
         argFlag<'v'>("verbose", "Verbose diagnostics")
         );
    opts.parse(argc, argv);
    auto n = opts.get<'n'>();
    // auto verbose = opts.get<'v'>();

    measure(cout, "sieve_index", [&]() { sieve_index(n); });
    measure(cout, "sieve_2n", [&]() { sieve_mod2_prime_seq(n); });
    measure(cout, "sieve_6n", [&]() { sieve_mod6_prime_seq(n); });
    return 0;
}

结束更新

您可以通过将每个

vector

替换为

co_yield

来轻松地调整它以构建一个

vec.push_back

素数。

coro::Generator<uint64_t> sieve_mod6_prime_seq(uint64_t max = uint64_t{1} << 20) {
    co_yield 2;
    co_yield 3;

    auto max_index = max / 3;
    auto bits_per = sizeof(uint64_t) * CHAR_BIT;
    auto nwords = (bits_per + max_index - 1) / bits_per;
    std::vector<uint64_t> words(nwords);

    words[0] |= 1;
    size_t wdx = 0;
    while (wdx < nwords) {
        auto b = std::countr_one(words[wdx]);
        auto p = 3 * (64 * wdx + b) + 1 + (b bitand 1);
        if (b < 64 and p < max) {
            co_yield p;

            for (auto j = p; j < max; j += 6 * p) {
                auto idx = j / 3;
                auto jdx = idx / 64;
                auto jmask = uint64_t{1} << (idx % 64);
                words[jdx] |= jmask;
            }

            for (auto j = 5 * p; j < max; j += 6 * p) {
                auto idx = j / 3;
                auto jdx = idx / 64;
                auto jmask = uint64_t{1} << (idx % 64);
                words[jdx] |= jmask;
            }
        }
        else {
            ++wdx;
        }
    }
    co_return;
}