素数計数

math/prime_count.hpp · GitHub で見る

共通ページ

#include "lib/math/prime_count.hpp"

auto a = prime_count(10);             // 4: 2, 3, 5, 7
auto b = prime_count(100000000000ULL); // 4118054813

操作	内容・返り値	計算量
prime_count(n)	n 以下の素数の個数を std::uint64_t で返す	\(\displaystyle O\!\left(\frac{n^{3/4}}{\log n}\right)\)

n は非負の std::uint64_t です。上限自身を含み、prime_count(0) と prime_count(1) は 0、prime_count(2) は 1 です。素数篩の prime_sieve(n) は n 未満を列挙するため、上限の扱いが異なります。

全ての素数を列挙せず、奇数だけの部分篩と整数商の再利用で計算します。小さい素因数を除いた候補だけを残し、最後に素数の積による補正を行います。繰り返し現れる除算は逆数乗算と整数の補正に置き換えます。整数平方根も積を使って補正するため、浮動小数点の丸めによる数え違いはありません。

前計算表や共有する可変状態は持たず、呼び出しごとに必要な領域を確保して終了時に解放します。64 bit 全域を整数演算で扱いますが、計算時間とメモリは入力に応じて増えます。特に最大値付近では数十 GB の作業領域を必要とするため、実用上の上限は環境の資源に依存します。\(n=10^{11}\) での主要配列は約 4 MiB です。

時間は \(\displaystyle O\!\left(\frac{n^{3/4}}{\log n}\right)\)、追加空間は \(\displaystyle O(\sqrt n)\) です。n < 2 では定数時間で終了します。

検証コード

• test/prime_count.test.cpp

コード

#ifndef PRIME_COUNT_HPP
#define PRIME_COUNT_HPP

#include <cmath>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <vector>

namespace prime_count_detail {
template<bool FastDivide> std::uint64_t compute(std::uint64_t n) {
	if (n < 2) return 0;
	std::uint64_t root = std::sqrt(static_cast<long double>(n));
	while (__uint128_t(root) * root > n) --root;
	while (__uint128_t(root + 1) * (root + 1) <= n) ++root;
	std::size_t size = (root + 1) / 2;
	std::vector<std::uint32_t> small(size), rough(size);
	std::vector<std::uint64_t> large(size), quotient(size);
	std::vector<std::uint8_t> removed(root + 1);
	for (std::size_t i = 0; i < size; ++i) {
		small[i] = i;
		rough[i] = 2 * i + 1;
		if constexpr (FastDivide) quotient[i] = double(n) / rough[i];
		else quotient[i] = n / rough[i];
		large[i] = (quotient[i] - 1) / 2;
	}
	std::uint32_t count = 0;
	for (std::uint64_t p = 3; p * p <= root; p += 2) {
		if (removed[p]) continue;
		removed[p] = 1;
		for (std::uint64_t i = p * p; i <= root; i += 2 * p) removed[i] = 1;
		const std::uint64_t reciprocal = UINT64_MAX / p + FastDivide;
		auto divide = [=](std::uint64_t x) {
			std::uint64_t q = (__uint128_t(x) * reciprocal) >> 64;
			if constexpr (FastDivide) return q;
			else return q + (x - q * p >= p);
		};
		std::size_t next = 0, k = 0;
		for (; k < size && std::uint64_t(rough[k]) * p <= root; ++k) {
			const std::uint32_t r = rough[k];
			if (removed[r]) continue;
			large[next] = large[k] - large[small[r * p / 2] - count] + count;
			quotient[next] = quotient[k];
			rough[next++] = r;
		}
		for (; k < size; ++k) {
			const std::uint32_t r = rough[k];
			if (removed[r]) continue;
			large[next] = large[k] - small[(divide(quotient[k]) - 1) / 2] + count;
			quotient[next] = quotient[k];
			rough[next++] = r;
		}
		size = next;
		std::size_t i = (root - 1) / 2;
		for (std::uint64_t j = (root / p - 1) | 1; j >= p; j -= 2) {
			const std::uint32_t delta = small[j / 2] - count;
			const std::size_t end = j * p / 2;
			for (; i >= end; --i) small[i] -= delta;
		}
		++count;
	}
	__int128_t result = 1 + __int128_t(large[0]) + __int128_t(size) * (size - 1) / 2;
	result += (__int128_t(count) - 1) * (size - 1);
	for (std::size_t k = 1; k < size; ++k) result -= large[k];
	for (std::size_t k = 1; k < size; ++k) {
		const std::uint64_t p = rough[k], reciprocal = UINT64_MAX / p + FastDivide;
		auto divide = [=](std::uint64_t x) {
			std::uint64_t q = (__uint128_t(x) * reciprocal) >> 64;
			if constexpr (FastDivide) return q;
			else return q + (x - q * p >= p);
		};
		const std::size_t end = small[(divide(quotient[k]) - 1) / 2] - count;
		if (end <= k) break;
		for (std::size_t j = k + 1; j <= end; ++j)
			result += small[(divide(quotient[j]) - 1) / 2];
		result -= __int128_t(end - k) * (count + k - 1);
	}
	return static_cast<std::uint64_t>(result);
}
}

inline std::uint64_t prime_count(std::uint64_t n) {
	if constexpr (std::numeric_limits<double>::is_iec559 && std::numeric_limits<double>::digits == 53)
		if (n <= (std::uint64_t(1) << 42)) return prime_count_detail::compute<true>(n);
	return prime_count_detail::compute<false>(n);
}

#endif

Back to top page