miller-rabin (素数判定法)
(src/math/miller-rabin.hpp)

• View this file on GitHub
• Last update: 2023-01-09 22:06:09+09:00
• Include: #include "src/math/miller-rabin.hpp"

概要

正整数nにが素数かどうかを判定する

機能

bool isprime(long long n)

• $n$ が素数かどうかを判定します
• 計算量 $O(\log ^ 2 n)$

計算量

参考

• wikipedia-ja
• wikipedia-en
• C++の実装で参考にしました
• アルゴリズムの理解で参考にしました

Verified with

• test/miller-rabin.test.cpp

Code

#pragma once

#include <vector>

namespace zawa::impl {

long long modpow(__int128_t a, __int128_t b, long long mod) {
	long long res = 1LL;
	a %= mod;
	while (b) {
		if (b & 1) {
			res = ((__int128_t)res * a) % mod;
		}
		a = (a * a) % mod;
		b >>= 1;
	}
	return res;
}

std::vector<long long> cond = { 2, 325, 9375, 28178, 450775, 978504, 1795265022 };

} // namespace zawa::impl

namespace zawa {

bool isprime(long long n) {
	if (n <= 1) {
		return false;
	}
	if (n == 2 or n == 3) {
		return true;
	}
	if (!(n & 1)) {
		return false;
	}
	long long s = 0LL, d = n - 1;
	while (d % 2 == 0) {
		s++;
		d >>= 1;
	}
	for (auto c : impl::cond) {
		if (c % n == 0) {
			return true;
		}
		long long x = impl::modpow(c, d, n);
		if (x != 1) {
			long long t = 0;
			for ( ; t < s ; t++) {
				if (x == n - 1) {
					break;
				}
				x = ((__int128_t)x * x) % n;
			}
			if (t == s) {
				return false;
			}
		}
	}
	return true;
}

} // namespace zawa

#line 2 "src/math/miller-rabin.hpp"

#include <vector>

namespace zawa::impl {

long long modpow(__int128_t a, __int128_t b, long long mod) {
	long long res = 1LL;
	a %= mod;
	while (b) {
		if (b & 1) {
			res = ((__int128_t)res * a) % mod;
		}
		a = (a * a) % mod;
		b >>= 1;
	}
	return res;
}

std::vector<long long> cond = { 2, 325, 9375, 28178, 450775, 978504, 1795265022 };

} // namespace zawa::impl

namespace zawa {

bool isprime(long long n) {
	if (n <= 1) {
		return false;
	}
	if (n == 2 or n == 3) {
		return true;
	}
	if (!(n & 1)) {
		return false;
	}
	long long s = 0LL, d = n - 1;
	while (d % 2 == 0) {
		s++;
		d >>= 1;
	}
	for (auto c : impl::cond) {
		if (c % n == 0) {
			return true;
		}
		long long x = impl::modpow(c, d, n);
		if (x != 1) {
			long long t = 0;
			for ( ; t < s ; t++) {
				if (x == n - 1) {
					break;
				}
				x = ((__int128_t)x * x) % n;
			}
			if (t == s) {
				return false;
			}
		}
	}
	return true;
}

} // namespace zawa

Back to top page