Millor-Rabin(ミラー・ラビン素数判定法) (math/millor_rabin.cpp)

• View this file on GitHub
• View document part on GitHub
• Last update: 2023-10-07 01:33:17+09:00

なにこれ

ミラー・ラビン素数判定法を用いた素数判定．

関数

• isprime(n)：$n$ が素数であれば true を，そうでなければ false を返す．$n \leq $ LLONG_MAX

計算量

• $O(\log^3 n)$

補足

定数倍 $(k=7)$ がついている．

以下，ループさせたときの実行時間（Local / AtCoder $[\textrm{ms}]$）

• $n \leq 2 \times 10^9$ で乱数
  • loop = $10^5$ : $0.024$ / $0.015$
  • loop = $10^6$ : $0.237$ / $0.150$
• $n \leq 9 \times 10^{18}$ で乱数
  • loop = $10^5$ : $0.064$ / $0.029$
  • loop = $10^6$ : $0.639$ / $0.292$

テスト用コード

参考

• Miller–Rabin primality test - Wikipedia
• Deterministic variants of the Miller-Rabin primality test

Verified with

• test/aoj_ALDS1_1_C_millorRabin.cpp
• test/yosupo_primality_test.cpp

Code

#pragma once

#ifndef call_include
#define call_include
#include <bits/stdc++.h>
using namespace std;
#endif

bool isprime(long long n) {
	if(n == 2) return true;
	if(n < 2 || (n & 1) == 0) return false;

	if(n < 200000) {
		for(long long i = 3; i * i <= n; i += 2)
			if(n % i == 0) return false;
		return true;
	}

	long long d = n - 1;
	int s = 0;
	while(!(d & 1)) {
		s++;
		d >>= 1;
	}

	vector<int> bases = (n < 4759123141ll) ? vector({2, 7, 61}) :
											 vector({2, 325, 9375, 28178, 450775, 9780504, 1795265022});

	for(int k : bases) {
		if(k >= n) break;

		__int128_t r = 1, q = k;
		long long t = d;
		while(t > 0) {
			if(t & 1) (r *= q) %= n;
			t >>= 1;
			(q *= q) %= n;
		}
		if(r == 1 || r == n - 1) continue;

		for(int i = 1; i < s; i++) {
			(r *= r) %= n;
			if(r == n - 1) break;
		}
		if(r != n - 1) return false;
	}

	return true;
}

#line 2 "math/millor_rabin.cpp"

#ifndef call_include
#define call_include
#include <bits/stdc++.h>
using namespace std;
#endif

bool isprime(long long n) {
	if(n == 2) return true;
	if(n < 2 || (n & 1) == 0) return false;

	if(n < 200000) {
		for(long long i = 3; i * i <= n; i += 2)
			if(n % i == 0) return false;
		return true;
	}

	long long d = n - 1;
	int s = 0;
	while(!(d & 1)) {
		s++;
		d >>= 1;
	}

	vector<int> bases = (n < 4759123141ll) ? vector({2, 7, 61}) :
											 vector({2, 325, 9375, 28178, 450775, 9780504, 1795265022});

	for(int k : bases) {
		if(k >= n) break;

		__int128_t r = 1, q = k;
		long long t = d;
		while(t > 0) {
			if(t & 1) (r *= q) %= n;
			t >>= 1;
			(q *= q) %= n;
		}
		if(r == 1 || r == n - 1) continue;

		for(int i = 1; i < s; i++) {
			(r *= r) %= n;
			if(r == n - 1) break;
		}
		if(r != n - 1) return false;
	}

	return true;
}

Back to top page