Miller-Rabin 소수 판별법

추상대수학 시간에 배운 소수판정법인데, 아주 유용하여 기록해놓는다.

Miller-Rabin 소수 판별법은 확률에 의존하는 소수판정법인데, 라그랑주 정리를 이용한다.

라그랑주 정리는 p가 소수이면, $f(x) \equiv 0 ( \text{mod}\:p)$는 최대 $deg(f(x))$개의 incongruent solution을 가진다는 정리이다.

라그랑주 정리에 따르면, $x^2 \equiv 1 ( \text{mod}\: p)$는 최대 두 개의 해를 가진다.

그런데 우리는 두 가지의 자명해를 알고있다. 바로 $\pm1$이다. (단, p > 2라 하자)

따라서, 만약 자명해가 아닌 $x^2 \equiv 1 ( \text{mod}\: p)$를 만족하는 해를 발견한다면, p는 소수가 아니라는 것을 알 수 있다.

이 성질과, 페르마의 소정리를 교묘하게 이용한 것이 바로 Miller-Rabin 소수 판별법이다.

먼저, $p - 1 = 2^s \times d$ 꼴로 놓자.

p가 소수라면 페르마의 소정리에 의해 $a^{p - 1}  = a^{2^s \times d} \equiv 1 ( \text{mod}\: p)$이 성립하므로,

$a^{2^{s - 1} \times d} ( \text{mod}\: p)$를 확인한다. 이 값이 $\pm 1$이 아니면, 자명해가 아닌 $x^2 \equiv 1 ( \text{mod}\: p)$의 해를 발견했으므로 p는 소수가 아니다.

값이 1이면 $a^{2^{s - 2} \times d} ( \text{mod}\: p)$에 대해 반복하면 되고, 값이 -1이면 더이상 확인할 수 있는 방법이 없으므로 그만둔다.

대충 이런 원리로 체크를 한다.

실제로 알고리즘을 구현할 때에는 편의상 $a^d ( \text{mod}\: p)$부터 시작하여 $a^{2^s \times d} ( \text{mod}\: p)$까지 가면서 확인하는 것이 편리하다.

$a^d ( \text{mod}\: p) \equiv 1$일 경우나 중간에 -1이 튀어나왔을 경우, 그 뒤는 항상 1이 되므로 검사를 중단해도 된다.

합성수로 판정되는 경우는 ($\pm 1$이 아닌 값) -> 1로 가는 경우이다.

앞서 언급한 바 있지만, 이 판정법은 확률에 의존하므로 정확한 결과를 주지는 못한다.

때문에 보통 x를 바꾸어가며 k번 테스트를 거친다. k번 판정을 통과할 경우, 정확도는 $4^{-k}$ 정도라고 한다.

하지만, 비교적 작은 수들에 대해서는 몇몇 a에 대해서만 판정을 해보면 된다는 것이 알려져 있다.

우리가 보통 사용하는 unsigned int는 a = 2, 7, 61에 대해서만 판정을 진행해보면 되고,

unsigned long long의 경우에는 37까지의 소수에 대해서만 판정을 진행해보면 된다. (다만, unsigned long long은 제곱하고 mod를 구하는 것을 구현할 때 어려움이 따른다)

아래는 구현해놓은 코드이다.

bool isPrime(unsigned int n) {
	if (n < 2) return 0;
	if (n % 2 == 0) {
		if (n == 2) return 1;
		else return 0;
	}
	int s = __builtin_ctz(n - 1);
	for (unsigned long long a: {2, 7, 61}) {
		if (a >= n) continue;
		int d = (n - 1) >> s;
		unsigned long long now = 1;
		while (d != 0) {
			if (d & 1) now = (now * a) % n;
			a = (a * a) % n;
			d >>= 1;
		}
		if (now == 1) goto success;
		for (int i = 0; i < s; ++i) {
			if (now == n - 1) goto success;
			now = (now * now) % n;
		}
		return 0;
		success:;
	}
	return 1;
}

조금 깔끔해보이려고 goto문을 썼는데, 의미 전달에 문제가 없기를 바란다.