素数判定（高速版）

math/is_prime.hpp · GitHub で見る

全体サイト（256 KiB向け） · 高速実装

軽量実装（64 KiB）はこちら

#include "lib/math/is_prime.hpp"

bool a = is_prime(97); // true
bool b = is_prime(1); // false

操作	内容・返り値	計算量
is_prime(n)	std::uint64_t 型の n が素数なら true、それ以外は false を返す	\(O(\log n)\)

非負整数を 1 個ずつ判定するときに使います。n は std::uint64_t で、64 bit 全域を決定的に判定します。0 と 1 は false、2 は true です。ある上限までの素数を全て列挙する場合は 素数篩 を使えます。

Bradley Berg の基数表を使います。ヘッダーをコピー・再配布するときは、表に付属する著作権表示と利用許諾を残してください。

時間は \(O(\log n)\)、追加空間は \(O(1)\) です。

依存ライブラリ

• 素数判定の基数表
• 素数判定の共通演算
• Montgomery 64
• Montgomery 32

検証コード

• test/discrete_log.test.cpp
• test/factorize.test.cpp
• test/is_prime.test.cpp
• test/mod_pow_log.test.cpp
• test/multiplicative_convolution.test.cpp

コード

#ifndef IS_PRIME_HPP
#define IS_PRIME_HPP

#include "../components/is_prime_core.hpp"
#include "../components/is_prime_bases.hpp"

namespace is_prime_detail {
inline bool fast_prime(const Montgomery32 &m) {
	using V = Montgomery32::U;
	V d = m.mod - 1;
	int s = std::countr_zero(d);
	d >>= s;
	const std::array<V, 1> bases{bases32[(m.mod * 0xad625b89u) >> 24]};
	if (m.mod < (V(1) << 30)) return witnesses(m, d, s, bases, [&](V a, V b) {
		std::uint64_t x = std::uint64_t(a) * b;
		return V(x >> 32) + m.mod - V((std::uint64_t(V(x) * m.inverse) * m.mod) >> 32);
	});
	return witnesses(m, d, s, bases, [&](V a, V b) { return m.mul(a, b); });
}

inline bool fast_prime(const Montgomery64 &m) {
	if (m.mod < (U(1) << 32)) return fast_prime(Montgomery32(m.mod));
	U d = m.mod - 1;
	int s = std::countr_zero(d);
	d >>= s;
	U base = bases64[(std::uint32_t(m.mod) * 0xad625b89u) >> 18];
	constexpr std::array<U, 8> second{15, 135, 13, 60, 15, 117, 65, 29};
	auto check = [&](auto multiply) {
		if (m.mod < (U(1) << 49)) return witnesses(m, d, s, std::array<U, 2>{2, base}, multiply);
		return witnesses(m, d, s, std::array<U, 3>{2, base, second[base >> 13]}, multiply);
	};
	if (m.mod < (U(1) << 62)) {
		U inverse = -m.inverse;
		return check([&](U a, U b) {
			__uint128_t x = __uint128_t(a) * b;
			return U((x + __uint128_t(U(x) * inverse) * m.mod) >> 64);
		});
	}
	return check([&](U a, U b) { return m.mul(a, b); });
}
}

inline bool is_prime(std::uint64_t n) {
	return is_prime_detail::check(n, [](const auto &m) { return is_prime_detail::fast_prime(m); });
}

#endif

Back to top page