C++20とBoost C++ライブラリで素数を求めるコード

my_primes.cpp

// コンパイルオプション：g++ -std=c++20 -Wall -O3 -march=native -mtune=native my_primes.cpp
#include <chrono>    // for std::chrono
#include <cmath>     // for std::log
#include <cstdint>   // for std::int32_t
#include <format>    // for std::format
#include <iostream>  // for std::iostream
#include <vector>    // for std::vector
#include <boost/math/quadrature/gauss.hpp> // boost::math::quadrature::gauss

namespace {
static auto constexpr MAX = 1000000000;
static auto constexpr NGAUSS = 100;
double li(double x);
void seek_prime();
}  // namespace

int main()
{
    std::cout << std::format("max = {}\n", MAX);
    seek_prime();

    return 0;
}

namespace {
double li(double x)
{
    auto const f = [](auto x) { return 1.0 / std::log(x); };
    
    return boost::math::quadrature::gauss<double, NGAUSS>::integrate(f, 2.0, x);
}

void seek_prime()
{
    using namespace std::chrono;

    auto const start = system_clock::now();

    auto const max = MAX + 1;
    auto const pmax = static_cast<std::int32_t>(li(static_cast<double>(max))) + 4000;

    std::vector<bool>         sieve(max / 2, true);
    std::vector<std::int32_t> prime(pmax);

    for (auto i = 3; i * i < max; i += 2) {
        if (sieve[i / 2]) {
            for (auto j = i * i / 2; j < max / 2; j += i) {
                sieve[j] = false;
            }
        }
    }

    prime[0] = 2;
    auto pnum = 1;
    for (auto i = 1; i < max / 2; i++) {
        if (sieve[i]) {
            prime[pnum] = i * 2 + 1;
            pnum++;
        }
    }

    auto const end = system_clock::now();
    auto const elapsed = std::chrono::duration<double>(end - start).count();

    std::cout << std::format("{:.3f} (sec)\n素数が{}個見つかりました\n", elapsed, pnum);
}
}  // namespace