Miller-Rabin 소수 판정법

p가 소수일 때, p의 배수가 아닌 모든 a에 대해 ${ a }^{ p-1 }\equiv 1 \mod p$이다.

홀수 d에 대해 $p-1={ 2 }^{ s }\cdot d$이라고 하자. 이제 이렇게 쓸 수 있다.

$$\Large{a^{2^{s} \cdot d} \equiv 1 \mod p}$$

s가 1 이상이라면 2를 하나 떼어낼 수 있다.

$$\Large{
a^{2^{s} \cdot d} \equiv 1 \mod p \\
a^{2^{(s-1)} \cdot 2 \cdot d} \equiv 1 \mod p \\
a^{(2^{(s-1)} \cdot d) \cdot 2} \equiv 1 \mod p \\
\left(a^{2^{(s-1)} \cdot d}\right)^{2} \equiv 1 \mod p
}$$

마지막 식은 수학적으로 이렇게 쓸 수 있다.

$$\Large{
a^{2^{(s-1)} \cdot d} \equiv \pm 1 \mod p
}$$

만일 우변이 +1이 나왔다면 왼쪽 지수의 (2^{s-1} \cdot d)를 다시 한 번 쪼갤 수 있다.

이렇게 계속 쪼개나가다보면, 반드시 다음 중 하나가 성립한다는 결론에 이르를 수 있다.

$$\begin{cases}\Large{
a^{2^{s-1} \cdot d} \equiv -1 \mod p \\
a^{2^{s-2} \cdot d} \equiv -1 \mod p \\
\vdots \\
a^{2^{1} \cdot d} \equiv -1 \mod p \\
a^{2^{0} \cdot d} \equiv -1 \mod p \\
a^{2^{0} \cdot d} \equiv +1 \mod p \\
}\end{cases}
$$

이렇게도 쓸 수 있다.

$$\large{\begin{cases}
a^{(p-1)} & \equiv -1 \mod p \\
a^{(p-1)/2} & \equiv -1 \mod p \\
a^{(p-1)/2/2} & \equiv -1 \mod p \\
\vdots \\
a^{(p-1)/2/2\cdots /2} & \equiv -1 \mod p \\
a^{(p-1)/2/2\cdots /2} & \equiv +1 \mod p
\end{cases}}
$$

(p-1)이 더 이상 2로 나뉘지 않을 때까지 나눈다.

Miller-Rabin primality test는 이 명제의 대우 명제를 이용한다.
어떤 n에 대해 저 모든 식을 만족하지 않는 a를 찾을 수 있다면 n은 소수가 아니라는 것.
즉 몇 가지 a에 대해서 위의 조건을 확인해서, 하나라도 거짓이 나온다면 "소수가 아님"을 출력,
모든 a에 대해서 통과한다면 "소수임"을 출력하는 것이 Miller-Rabin 소수 판정법이다.

"소수가 아님"을 출력할 때는 확실한 수학적인 증거가 있으므로 괜찮다.
그런데 넣어본 a에 대해 통과할 때 사실은 "아직 어느쪽이든 될 수 있다"가 맞는 말이지만, 이때 "소수임"을 반환하기로 했으므로 이 부분이 오류의 여지가 있다.

그럼 a를 어떤 걸 해야 하느냐?
Miller는 $O(\left(\log n)^2\right))$ 이하의 모든 a를 확인하면 됨을 알았다. (단, 이는 일반화된 리만 가설이 참일 때를 가정한다고 한다. 얘는 보통 참으로 믿는단다.) 나중엔 저 big-O notation의 계수가 2 이하임이 알려졌다. 즉 $2\left(\log n\right)^{2}$ 이하의 모든 a를 확인하면 된단다.

조금 더 실용적으로, unsigned long long int 안에 들어가는(2^64 이하의) 수들은 a를 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 (첫 12개의 소수)인 경우만 시도해보면 된다고 한다.

typedef unsigned long long T;
T ipow(T b, T e, T mod){
    T ret = 1; b %= mod;
    while(e){
        if(e&1) (ret *= b) %= mod;
        (b *= b) %= mod; e >>= 1;
    }
    return ret;
}

// returns 1 when n is not prime
bool miller_counter(T n,T a)
{
    if(a % n == 0) return 0;
    for(T e=n-1;;){
        T tmp = ipow(a, e, n);
        if(tmp == n-1) return 0;
        if(e&1){ return tmp != 1; }
        e >>= 1;
    }
}

int miller_base[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37,0};

bool isprime(T x) {
    for(int i=0; i<12; ++i) if(miller_counter(x,miller_base[i])) return 0;
    return 1;
}