rkm0959

관련 코드는 github에서 찾아볼 수 있다.

이 글의 순서는 다음과 같다.

• 원시근의 정의 및 관련 사실들. 핵심 중의 핵심을 제외하고, 증명 대부분은 생략할 것이다.
• 원시근을 찾는 알고리즘과 위수를 계산하는 알고리즘.
• 이산로그 문제와 Baby Step Giant Step. 그리고 그 확장인 Pohlig-Hellman 알고리즘 (optional!)
• 이산제곱근 문제와 그 해결법. 특수한 경우인 modular square root의 경우. 그리고 Adleman-Manders-Miller (매우 optional!) 
• 특히, Adleman-Manders-Miller는 지금까지 나온 경우를 딱 한 번 봤다. (GP of Baltic Sea) 이는 논문 링크 + 구현체 링크로 대체.
• 그리고 마지막으로, FFT와 NTT 사이의 연결고리를 설명한다. 앞서 mobius inversion = 포함-배제를 설명한 것과 같이, FFT = NTT를 설명한다.

원시근이란 무엇인가? 위수란 무엇인가?

• 자연수 $n \ge 2$이 있을 때, $n$과 서로소이며 $n$ 미만인 자연수의 개수는 $\phi(n)$이다.
• 만약 $g, g^2, \cdots , g^{\phi(n)}$을 각각 $n$으로 나눈 나머지가 이 $\phi(n)$개의 자연수 전부가 된다면, $g$를 $n$의 원시근이라 한다.
• 다른 말로 하면, $t$가 $n$과 서로소인 자연수라면 적당한 $1 \le k \le \phi(n)$이 있어 $g^k \equiv t \pmod{n}$이란 것이다.
• 사실 1. 원시근을 갖는 자연수는 $2, 4$와 홀수 소수 $p$와 자연수 $k \ge 1$에 대하여 $p^k$, $2p^k$다.
• 사실 2. $n$이 원시근을 갖는 자연수라면, 원시근은 정확히 $\phi(\phi(n))$개 있다. 다른 말로 하면, "원시근은 꽤 많다".
• 예를 들어, $n=p$가 소수라면, $\phi(\phi(p))=\phi(p-1)$이다. 특히, $\phi(n) > n/(e^{\gamma} \log \log n + 3/\log \log n)$임이 알려져 있다.
• 자연수 $g, n$이 있을 때, (단, $\text{gcd}(g, n)=1$) $g^k \equiv 1 \pmod{n}$을 만족하는 최소 자연수 $k$를 $g$의 위수라고 한다.
• 사실 3. $g$가 $n$의 원시근인 것은, $g$의 위수가 $\phi(n)$인 것과 동치이다. (증명 힌트. $g^i \equiv g^j \pmod{n}$이면 $g^{j-i} \equiv 1 \pmod{n}$) 
• 사실 4. ("학부 대수학의 반의 avatar") $\text{gcd}(g, n)=1$이고 $g$의 위수가 $k$라고 하자. $g^m \equiv 1 \pmod{n}$이 성립할 필요충분조건은 $k|m$이다. 
• 사실 4 증명. $k|m$이면 $g^m \equiv (g^k)^{m/k} \equiv 1^{m/k} \equiv 1 \pmod{n}$이므로 ok. 만약 $g^m \equiv 1 \pmod{n}$인데 $k|m$이 아니라고 하자. 그러면 $m=kt+r$인 $0 < r < k$가 있고, $g^r \equiv g^{m-kt} \equiv 1 \pmod{n}$이다. 그런데 $0 < r < k$이므로 이는 위수의 "최소성"에 모순이다. 즉, 애초에 $k$가 $g^k \equiv 1 \pmod{n}$을 만족하는 "최소의" 자연수였는데, $g^r \equiv 1 \pmod{n}$, $r<k$이므로 모순. 좋은 아이디어.
• 사실 5. $g$의 위수가 $t$라면, $g^d$의 위수는 $t/\text{gcd}(d, t)$가 된다. 사실 4에서 바로 증명할 수 있다.
• 사실 1, 2의 증명은 학부 정수론에서 가장 중요한 부분 중 하나이며, 후에 (사실 1은) 현대대수학에서도 다시 증명하게 된다.

왜 원시근이 강력한가?

• $g$가 원시근이라 하자. $g^k \equiv r \pmod{n}$이라면, $k$를 $r$의 "로그 값"이라고 생각할 수 있겠다.
• 로그가 강력한 이유가 뭐였더라? 거듭제곱을 곱셈으로, 곱셈을 덧셈으로 바꿔주는 것이었다.
• 우리는 거듭제곱은 잘 모르지만, 앞서 2, 3 단원을 통해서 일차합동식 $ax \equiv b \pmod{n}$에 대해서는 안다.
• 그러니 원시근을 사용해서 거듭제곱에 관한 문제를 일차식으로 바꿀 수 있을 것이다. 이 접근을 "이산제곱근"에서 사용한다.

이제부터 원시근을 찾는 방법을 설명한다. $n$이 원시근을 갖는 자연수, 즉 $2, 4, p^k, 2p^k$ 형태라고 가정하자.

• 원시근은 충분히 많으므로, 랜덤하게 $g$를 잡은 뒤 $g$가 원시근인지 판별하는 과정을 원시근을 찾을 때까지 반복하면 된다.
• 원시근이 많다는 것은, 한 번 랜덤한 시도를 했을 때 성공적으로 원시근을 찾을 확률이 나쁘지 않다는 것이다.
• $g$가 원시근인지 판별하려면, 먼저 $\text{gcd}(g, n)=1$인지 확인한다. 그 후, $g$의 위수가 $\phi(n)$인지 확인하면 된다.

그러므로, 원시근을 찾는 방법을 알기 위해서는 위수를 계산하는 알고리즘을 찾으면 된다. 여기서 사실 4가 활용된다.

• $\phi(n)$의 소인수분해가 이미 계산되었다고 가정한다. $\phi(n) = p_1^{e_1} \cdots p_k^{e_k}$라고 가정한다.
• $n=p^k$, $n=2p^k$인 경우 $\phi(n) = p^{k-1}(p-1)$임에 유의하라. 
• 이제 $g$의 위수를 찾는 것이 목표다. $\text{gcd}(g, n)=1$임을 가정한다. 오일러의 정리에 의하여, $g^{\phi(n)} \equiv 1 \pmod{n}$.
• 이는 사실 4에 의하여, 위수의 값이 $\phi(n)$의 약수임을 의미한다. 
• 이제 $g^{\phi(n) / p_1}$의 값을 보자. 만약 이 값이 $1 \pmod{n}$이라면, 사실 4에서 위수가 $\phi(n)/p_1$의 약수임을 얻는다.
• 만약 $g^{\phi(n)/p_1} \equiv 1 \pmod{n}$이었다면, $g^{\phi(n)/p_1^2}$을 보자. (단, $p_1^2 | \phi(n)$) 이 값 역시 $1 \pmod{n}$이라면, 사실 4에서 위수가 $\phi(n) / p_1^2$의 약수.
• $g^{\phi(n)/p_1^2}$은 $1 \pmod{n}$이 아니었다고 하면, 위수가 $\phi(n)/p_1$의 약수였으나 $\phi(n)/p_1^2$의 약수는 아님을 의미한다.
• 이는 다르게 말하면 위수가 $p_1$을 몇 개 가지고 있는지 파악했다는 이야기다. 이제 $p_2$로 넘어가서 같은 작업을 반복, $p_k$까지 하자.
• 즉, $\phi(n)$에서 시작하여, $g^{???} \equiv 1 \pmod{n}$이 성립할 때까지 계속 소인수들을 제거해주면 된다.
• 시간복잡도를 생각해보면, 대강 $(e_1+e_2 + \cdots + e_k) = \mathcal{O}(\log n)$번  $g$의 거듭제곱을 계산하므로, $\mathcal{O}(\log^2 n)$이다.
• 물론, $\phi(n)$의 소인수분해가 필요했음에 주의하라. 필요한 경우, 6단원에서 배운 Pollard-Rho 소인수분해를 활용하자. 

이산로그 문제와 Baby Step Giant Step

• 원시근 $g$와 $\text{gcd}(h, n)=1$을 만족하는 자연수 $h$가 주어졌을 때, $g^x \equiv h \pmod{n}$을 만족하는 $x$를 찾는 문제를 이산로그 문제라 한다.
• 즉, 앞서 언급했던 원시근의 강력함인 "로그"를 실제로 계산하는 문제이다. 이산적인 집합에서 계산하고 있으니 이산로그 문제다.
• 이 문제를 CP/PS에서 보이는 "meet-in-the-middle"을 이용하여 $\mathcal{O}(\sqrt{n} \log n)$ 시간에 해결하는 알고리즘이 있다. 

• $B = \lceil \sqrt{\phi(n)} \rceil$이라고 하자. 그러면 $0$과 $\phi(n)$ 사이의 모든 값은 적당한 $0 \le u, v < B$에 대하여 $uB+v$라고 쓸 수 있다.
• 그러니 우리는 $g^{uB+v} \equiv h \pmod{n}$을 만족하는 $0 \le u, v < B$를 찾으면 된다. 
• 그런데 이 식을 변형하면 $(g^B)^{u} \equiv h \cdot g^{-v} \pmod{n}$을 얻는다. 이제 meet-in-the-middle 각이 보인다.
• $h \cdot g^{-0}$부터 $h \cdot g^{-(B-1)}$를 C++의 std::map과 같은 자료구조에 넣어놓자. 이때, key 값을 $h \cdot g^{-v}$로 하고, value 값을 $v$로 하도록 한다. 
• 이제 $u$에 대하여 순회하면서, $(g^B)^{u} \pmod{n}$이 해당 자료구조의 key 값에 있는지 확인하자. 있다면, value 값에서 $v$를 얻는다.
• 이 방식으로 $u, v$를 찾을 수 있으며, 시간복잡도는 $\mathcal{O}(\sqrt{n} \log n)$이다. 로그는 C++의 std::map 때문에 붙는다.
• 이 알고리즘은 이렇게 $\pmod{n}$에 대한 이산로그를 구하는 것에서만 쓰일 수 있는 것이 아니라, finite abelian group 구조를 갖는 모든 문제에서 적용될 수 있다. 애초에 이산로그 문제 자체가 finite abelian group 구조에서 정의되는 문제이며, 지금 우리가 다루고 있는 문제는 그 중 일부인 경우이다. 
• 말은 어렵게 했으나, 실제로 문제를 풀면 다른 생각이 들 수도 있다. BOJ의 xorshift32 문제를 참고하라.
• Note. $g$가 원시근이 아니어도, $g$의 위수를 $k$라고 하면 같은 방법으로 $\mathcal{O}(\sqrt{k} \log k)$ 알고리즘을 생각할 수 있다. 위에서 $\phi(n)$을 $k$로 바꾸면 된다.

Pohlig-Hellman 알고리즘 (optional)

• 위 Baby-Step-Giant-Step 알고리즘을 강화한 것이 Pohlig-Hellman 알고리즘이다. $\phi(n) = p_1^{e_1} \cdots p_k^{e_k}$라고 하자. 
• $g^x \equiv h \pmod{n}$인 $x$를 찾는다는 것은, 사실 $x \pmod{\phi(n)}$을 찾는다는 것이다.
• 그러니, 중국인의 나머지 정리를 생각하면 $x \pmod{p_i^{e_i}}$들을 계산한 다음 합치면 충분하다.

여기서 $\mathcal{O}(\sqrt{p_i} \log p_i)$ 시간에 $x \pmod{p_i}$를 찾는 방법을 제시한다. 

$g^x \equiv h \pmod{n}$이면, $(g^{\phi(n)/p_i})^x \equiv h^{\phi(n)/p_i} \pmod{n}$ 역시 성립한다. 여기서 $g^{\phi(n)/p_i}$의 위수가 $p_i$임에 주목하자.

이 새로운 이산로그 문제를 풀면 $x \pmod{p_i}$를 얻고, 이를 위해서 소모되는 시간은 $\mathcal{O}(\sqrt{p_i} \log p_i)$가 된다.

이제 이를 확장하여, $x \pmod{p_i^2}$을 찾는 방법을 고안하자. 물론, $p_i^2 | \phi(n)$이라고 가정한다.

마찬가지로, $(g^{\phi(n)/p_i^2})^x \equiv h^{\phi(n)/p_i^2} \pmod{n}$을 생각한다. $g^{\phi(n)/p_i^2}$의 위수는 $p_i^2$이다.

그러니 이 이산로그 문제를 풀면 $x \pmod{p_i^2}$을 얻는데, 중요한 점은 우리가 이미 $x \pmod{p_i}$를 안다는 것이다.

이미 계산한 값인 $x \pmod{p_i}$를 $r$이라고 하자. 이제 $x \pmod{p_i^2} = p_ik + r$이라고 쓸 수 있고, 목표는 $k$가 된다. 그러면 $$ (g^{\phi(n)/p_i})^k \equiv h^{\phi(n)/p_i^2} \cdot (g^{\phi(n)/p_i^2})^{-r} \pmod{n}$$을 얻는다. 이제 $k$를 찾기 위해 이산로그 문제를 풀면 $\mathcal{O}(\sqrt{p_i} \log p_i)$ 시간에 $x \pmod{p_i^2}$을 얻는다. 

이를 $p_i^{e_i}$까지 반복하고, 이를 각 prime power에 대해 반복하면 끝. 시간복잡도는 $\mathcal{O}(\sum e_i \sqrt{p_i} \log p_i )$다.

관련 문제로는 위에 언급한 xorshift32를 강화한 xorshift64가 있다. 

그리고 필자도 굉장히 최근에 배운 트릭을 하나 소개한다. 매우 optional인 건 맞지만, 흥미롭고 Pohlig-Hellman 보다 쉽다.

• $g^x \equiv h \pmod{p^2}$를 푼다고 가정하자. 단 $g$는 $p^2$에 대한 원시근이다. 
• 먼저 $g^{a_1} \equiv h \pmod{p}$를 푼다. 그러면 여기서 얻는 결과는 $x \equiv a_1 \pmod{p-1}$이 될 것이다.
• 이제 $g^{a_2} \equiv h \pmod{p^2}$을 푼다. 이제 $a_2 \equiv a_1 \pmod{p-1}$을 알고 있으므로, $a_2 \pmod{p}$를 알면 된다.
• 이를 위해서, $g^{(p-1)a_2} \equiv h^{p-1} \pmod{p^2}$을 푼다. 그런데, 놀라운 점은 $g^{p-1}$, $h^{p-1}$ 모두 $1 \pmod{p}$라는 것이다.
• 그래서 사실 $(1+pu)^{a_2} \equiv (1+pv) \pmod{p^2}$ 형태의 식을 풀게 되며, 특히 $g$가 원시근이면 $u$는 $p$의 배수가 아니다.
• 또한, 이항정리에서 $(1+pu)^{a_2} \equiv 1+pua_2 \pmod{p^2}$이다. 그러니 결국 $ua_2 \equiv v \pmod{p}$를 풀면 ok.
• 그러니 Baby-Step-Giant-Step을 두 번 돌릴 필요가 없으며, 한 번만 돌려도 된다. 이는 확장이 가능한 아이디어다. 이 링크 참고.

이산제곱근 문제

• 이번에는 $x, b$가 주어졌을 때, (단, $\text{gcd}(b, n)=1$) $a^x \equiv b \pmod{n}$을 만족하는 $a$를 구하는 것이 목표다.
• 중국인의 나머지 정리의 아이디어를 가져오자. $n = p_1^{e_1} \cdots p_k^{e_k}$라 하면, $a \pmod{p_i^{e_i}}$를 각각 구한 뒤 합치는 전략을 세울 수 있다.
• 그러니 $n$이 prime power인 경우에서만 해결하면 충분하다. 
• 여기서는 $n$이 홀수 소수의 거듭제곱인 경우만 다루겠다. 이제 원시근 $g$의 존재를 가정해도 된다.

$a = g^u$, $b = g^v$인 $u, v$를 이산로그 알고리즘으로 찾자. 목표는 $g^{ux} \equiv g^v \pmod{n}$인 $x$를 구하는 것이다.

이는 $g^{ux-v} \equiv 1 \pmod{n}$과 같으니, 사실 4에 의하여 $ux \equiv v \pmod{\phi(n)}$을 푸는 것과 같다.

이는 일차합동식이므로, 2단원에서 배운 내용을 통해 해결할 수 있다. $u$를 알면, 이를 가지고 $a \pmod{n}$ 역시 구할 수 있다.

특히, $\text{gcd}(x, \phi(n))=1$이라면, $tx \equiv 1 \pmod{\phi(n)}$인 $t$를 찾을 수 있다.

그러면 $b^t \equiv a^{tx} \equiv a \pmod{n}$이 성립함을 알 수 있고, (오일러 정리) $a \pmod{n}$을 쉽게 찾을 수 있다.

위 예시에서 $2$의 거듭제곱이 여러모로 까다롭다는 것을 알 수 있다. 이렇게 $2$의 거듭제곱을 가지고 귀찮게하는 문제는 적다.

• 하지만, 이 사실이 도움이 될 수 있다 : 모든 $x \pmod{2^n}$은 (단, $x$ 홀수) $5^k$ 또는 $-5^k$로 나타낼 수 있다.
• 즉, $x$가 $1 \pmod{4}$인 경우, 전부 $5$의 거듭제곱으로 표현이 되며, $3 \pmod{4}$인 경우, $-x$가 $5$의 거듭제곱으로 표현이 된다.
• 그러니 $5$를 $2^n$의 일종의 원시근과 비슷한 무언가로 생각할 수 있겠다. 이를 통해 이산로그 등 문제를 접근할 수 있다.

modular square root

• 이 파트에서는 특수한 경우인, $x^2 \equiv a \pmod{n}$을 다룬다. 편의상 $\text{gcd}(a, n) = 1$을 가정한다.
• 역시 중국인의 나머지 정리 아이디어를 사용하여, $x^2 \equiv a \pmod{p^k}$만 해결하면 된다.

먼저 $p$가 홀수 소수인 경우부터 보자. 특히, $k=1$인 경우, 즉 $x^2 \equiv a \pmod{p}$인 경우를 보자.

• 이차잉여. 위 식이 해를 가질 조건은 $a^{(p-1)/2} \equiv 1 \pmod{p}$이다. 
• 증명. 스케치만 하자면, 해를 가질 조건이 $a$가 $g^{even}$인 것과 동치이며 이건 다시 $a^{(p-1)/2} \equiv 1 \pmod{p}$와 동치이다.

해를 갖는 경우, $\mathcal{O}(\log^2 p)$ 시간에 해를 찾는 알고리즘이 있다. Tonelli-Shanks 알고리즘이며, 설명은 이 블로그에 있다.

특히, $x$가 해면 $-x$ 역시 해임에 주목하자. 해는 정확히 2개 있을 것이다. (이유가 궁금하다면, 이 정리 참고) 이제 문제는 $k \ge 2$인 경우이다. 

이제부터 각 $x^2 \equiv a \pmod{p^k}$의 해를 $x^2 \equiv a \pmod{p^{k+1}}$의 해로 "lift"하는 방법을 생각해보자.

이미 $r^2 \equiv a \pmod{p^k}$인 $r$의 값을 찾았다고 생각하자. 해는 정확히 $r, -r$ 두 개가 있을 것이다. 이제 $x^2 \equiv a \pmod{p^{k+1}}$을 풀자.

$x^2 \equiv a \pmod{p^k}$ 역시 성립하니, $x \equiv r \pmod{p^k}$ 또는 $x \equiv -r \pmod{p^k}$가 성립해야 한다.

• 이제 $x \equiv r \pmod{p^k}$인 경우를 가정하자. 그러면 $x \pmod{p^{k+1}}$은 $r + t p^k$ 형태로 쓸 수 있다. 단, $0 \le t < p$. 
• 이제 $x^2 \pmod{p^{k+1}}$을 생각하면, $r^2 + 2rt p^k + t^2 p^{2k} \equiv r^2 + 2rt p^k \pmod{p^{k+1}}$과 같다.
• 그러니, 문제는 결국 $r^2 + 2rt p^k \equiv a \pmod{p^{k+1}}$이고, 이를 변형하면 $2rtp^k \equiv (a-r^2) \pmod{p^{k+1}}$이다. 
• 여기서 $a-r^2 \pmod{p^{k+1}}$이 $p^k$의 배수임에 주목하자. 
• 그러니 원하는 것은 $(2rt - (a-r^2)/p^k) \cdot p^k$가 $p^{k+1}$의 배수가 되는 것, 즉 $2rt \equiv (a-r^2)/p^k \pmod{p}$이다. 
• 2와 $r$이 전부 $p$와 서로소이므로, 이 합동식을 풀어 $t \pmod{p}$를 얻는다. 이제 $x = r + tp^k$가 해이다.
• 시작을 $x \equiv -r \pmod{p^k}$로 얻었다면, $-r-tp^k$를 얻었을 것이다. 그러니 여전히 해는 $x, -x$ 정확히 두 개다.

이러한 과정을 잘 나타내는 정리가 Hensel's Lemma다. 이제 $p=2$인 경우를 생각해보자.

$k=1, 2$인 경우, 즉 $\pmod{2}$, $\pmod{4}$에서 식을 푸는 것은 brute force로 가능할 것이다. 이 경우는 생략.

이제부터 $k \ge 3$을 가정한다. 핵심적인 사실은, $x^2 \equiv 1 \pmod{8}$이 모든 홀수 $x$에 대해서 성립한다는 것이다.

• $x^2 \equiv a \pmod{2^k}$가 해를 가질 필요충분조건은 $a \equiv 1 \pmod{8}$.
• 이때, 해는 정확히 $4$개 존재한다. ($\pmod{2^k}$에서) $x$가 해라면 $x, -x, x+2^{k-1}, -x+2^{k-1}$이 해가 된다.

역시 "lift"할 각을 재보자. $x^2 \equiv a \pmod{2^3}$의 해는 (애초에 $a \equiv 1 \pmod{8}$이니) $x = 1, 3, 5, 7 \pmod{8}$이다.

해가 4개가 있고, 하나만 알면 나머지 세 개의 해를 쉽게 얻을 수 있다. 그러니 해를 하나만 관리해도 충분할 것이다. $x = 1$로 시작하자.

$x^2 \equiv a \pmod{2^k}$의 해를 하나 찾았다고 하고, 이를 $r$이라 하자. 그러면 $(r + t2^{k-1})^2 \equiv r^2 + rt 2^k \pmod{2^{k+1}}$이다. 

만약 $r^2 \equiv a \pmod{2^{k+1}}$이라면, $r$이 그대로 $x^2 \equiv a \pmod{2^{k+1}}$의 해가 된다.

그렇지 않다면, $r^2 \equiv a + 2^k \pmod{2^{k+1}}$일 것이다. 그러면 $x = r + 2^{k-1}$이 해가 된다. 이제 끝.

Adleman-Manders-Miller (매우 optional)

• 이산제곱근을 더욱 효율적으로 계산하는 알고리즘으로, Tonelli-Shanks의 확장이다.
• 이 글, 이 논문, 그리고 KAIST의 ICPC팀(이었던) 더불어민규당의 팀노트를 (discrete kth root) 참고하자.
• PS/CP 환경에서는 쓰이는 것은 저 오픈컵을 제외하고 본 적이 없다. 개인적으로도 암호학 대회 환경을 제외하면 쓴 적이 없다.

FFT와 NTT

• FFT와 FFT의 원리를 이미 알고 있다고 가정한다. $2^n = N$이 FFT의 "크기"라고 가정하자.
• FFT에서 사용되는 사실은 $\omega = e^{2\pi i / N}$이라 할 때, $\omega^{N}=1$이고 $1, \omega, \cdots , \omega^{N-1}$이 모두 다르다는 것이다.
• 또한, inverse FFT 과정에서는 계산 후에 $N$을 나누어야 하므로, 복소수에서 $N$을 나누는 연산이 잘 정의됨을 이용한다.
• NTT도 다르지 않다. $p-1 = 2^ab$이고, $a$가 충분히 큰 소수 $p$와 $p$에 대한 원시근 $g$를 준비하자. $2^a = N$이라 하자.
• $ \omega = g^{b}$라 하면, $\omega^{N} = 1$이고, $1, \omega, \cdots , \omega^{N-1}$이 모두 다르다. (모든 연산은 $p$로 나눈 나머지로 본다)
• 마찬가지로, $N$이 $p$와 서로소이므로, $N$으로 나누는 연산이 잘 정의된다. (modular inverse가 존재하므로)
• 그러니, FFT와 NTT는 정확히 같은 알고리즘이다. 실제 코드도 별로 다를 것이 없다.

독자들은 정수론은 잘 모르더라도 포함-배제의 원리는 어느 정도 숙지하고 있을 것이다. 

포함-배제의 원리는 애초에 교육과정의 경우의 수 세는 문제에서도 매우 많이 등장하며, "mainstream CP/PS"에서도 꽤 등장하기 때문이다.

이 파트에서 이야기하고자 하는 것은 mobius function/inversion이 포함-배제와 다를 게 없다는 것이다. 

(사실, 이 사실은 제곱 ㄴㄴ 문제를 해결한 사람이라면 이미 느꼈을 것이라고 생각한다. 여기서 더 자세하게 보자.)

이제부터 포함-배제의 원리 이야기로 넘어간다. 

• 포함-배제의 원리의 가장 익숙한 형태는, 유한집합 $A_1, A_2, \cdots ,A_n$에 대하여 $$\displaystyle |\cup_{i} A_i| = \sum_{i} |A_i| - \sum_{i< j} |A_i \cap A_j| + \sum_{i < j < k} |A_i \cap A_j \cap A_k| - \cdots  + (-1)^{n-1} | A_1 \cap A_2 \cap \cdots \cap A_n|$$이 성립한다는 명제일 것이다. 특히 $n=2, 3$인 경우는 고등학교 수학. 
• 위 식을 이해하는 방법은 다양할 것이다. 하지만 개인적인 생각으로 가장 쉬운 방법은 "원소별 기여"를 보는 것이다.
• $x \in \cup_i A_i$를 아무거나 하나 잡자. 이 원소는 $A_1, \cdots , A_n$ 중 정확히 $k$개의 집합에 속한다고 하자. 물론, $1 \le k \le n$.
• $x$의 존재로 위 식의 좌변은 $1$만큼 커진다. 그러니, $x$가 우변에 "기여하는 값"이 $1$임을 증명하면 충분하다.
• $\sum_i |A_i|$에 $x$가 기여하는 횟수는 $k$번, $\sum_{i < j} |A_i \cap A_j|$에 기여하는 횟수는 $x$가 $A_i, A_j$ 각각에 속해야 하므로 $\displaystyle \binom{k}{2}$번.
• 이제 느낌을 알 것 같다. $x$가 우변에 기여하는 총 값은 $\displaystyle \sum_{i=1}^n (-1)^{i-1} \binom{k}{i}$이다. 단, $n < m$면 $\displaystyle \binom{n}{m} = 0$이라 정의한다.
• 그런데 이항정리에서 $\displaystyle \sum_{i=0}^n (-1)^i \binom{k}{i} = \sum_{i=0}^k (-1)^i \binom{k}{i} = (1-1)^k =  0$이다. 이를 통해서, $x$가 기여하는 값이 $1$임을 얻는다.

이제 본격적으로 작업을 시작하자. 내용은 이 글을 따라간다.

• 포함-배제의 원리의 다른 형태를 소개한다. 전체집합 $U$가 유한집합이라 하자. $U$의 각 부분집합 $S$에 대하여 값 $f(S) \in \mathbb{R}$을 주자. 
• 이제 $U$의 각 부분집합 $A$에 대하여 $g(A) = \sum_{S \subset A} f(S)$를 정의하자. 우리의 목표는 $g$의 값을 가지고 $f$의 값을 복원하는 것이다.
• 이때, $f(A) = \sum_{S \subset A} (-1)^{|A|-|S|} g(S)$가 성립한다. 이를 위와 똑같은 방법으로 증명해보자.
• 우변은 $\sum_{S \subset A} (-1)^{|A|-|S|} \sum_{T \subset S} f(T) = \sum_{T \subset S \subset A} (-1)^{|A|-|S|} f(T)$다. 이제 $f(T)$가 "기여하는 값"을 생각해보자. 
• $T = A$인 경우, $\sum_{A \subset S \subset A} (-1)^{|A|-|S|} = (-1)^0 = 1$이다. 그러니 우변에 $f(A)$라는 값이 추가된다.
• $T \subsetneq A$인 경우, $\sum_{T \subset S \subset A} (-1)^{|A|-|S|}$를 계산해보자. $|A|=m+n$, $|T|=m$이라고 하면, $T \subset S \subset A$를 만족하는 크기 $m+k$의 집합 $S$는 총 $\displaystyle \binom{n}{k}$개 존재할 것이다. 그러므로, $\displaystyle \sum_{T \subset S \subset A} (-1)^{|A|-|S|} = \sum_{k=0}^n \binom{n}{k} \cdot (-1)^{n-k} = (1-1)^n = 0$이다. 여기서 $n \neq 0$이 사용되었다. 이는 $f(T)$가 사실 우변에 없음을 (0이 곱해져서 더해짐을) 의미한다. 그러니 사실 우변을 정리하면 남는 것은 $f(A)$ 뿐이며, 증명 끝.

위 "새로운 형태"는 포함-배제 느낌이 나지만, 뭔가 갑갑함을 느끼는 독자들이 있을 것이다.

• 우리에게 익숙한 형태로 돌아가서, $A_1, A_2, \cdots , A_n$이란 유한집합이 있다고 가정하자.
• 이제 다른 형태로 간다. $U = \{1, 2, \cdots , n\}$이라 하자. 각 $U$의 부분집합 $S$에 대해, $f(S)$를 다음과 같이 정의한다.
• $f(U) = 0$이다. 만약 $S$가 $U$가 아니라면, $f(S)$는 "속하는 집합의 index 모음이 정확히 $U - S$인 원소의 개수"다.
• 이게 무슨 뜻이냐? $n=3$이라 하자. $f(\{1\})$은, 정확히 $A_2, A_3$에만 속하는 원소의 개수, 즉 $|A_1^C \cap A_2 \cap A_3|$이다.
• 이제 $g(A)$를 "포함-배제의 새로운 형태"에서 다룬 방식으로 정의하고, $0 = f(U) = \sum_{S \subset U} (-1)^{|U|-|S|} g(S)$라는 식을 써보자.
• 이 식이 우리에게 익숙한 포함-배제의 형태와 동치임을 알 수 있을 것이다. 좋은 Exercise이며, CP/PS에서도 충분히 활용될 수 있는 연습이다.
• 힌트를 주자면, $g(A)$는 "속하는 집합의 index 모음이 $U-A$를 포함하는 원소의 개수"다. 앞서 다룬 $n=3$ 예시에서 $g(\{1\}) = |A_2 \cap A_3|$다.
• 특히, $g(U) = |\cup_i A_i|$가 됨을 확인할 수 있다. 이제 식을 전부 대입하고 약간의 고민을 하면 원하는 결과를 얻을 것이다.

이제 집합 $X$에 대하여 $\mu(X) = (-1)^{|X|}$라 정의하자.

• 우리의 결론은 $g(A) = \sum_{S \subset A} f(S)$면, $f(A) = \sum_{S \subset A} (-1)^{|A|-|S|} g(S) = \sum_{S \subset A} \mu(A-S) g(S)$이 성립한다는 것이다.
• $\mu$를 정의한 상태에서, "새로운 형태"에 대한 증명을 다시 보자. 증명의 핵심은 $$\sum_{T \subset S \subset A} (-1)^{|A|-|S|} = 0 $$이 $T \subsetneq A$일 때 성립함을 보인 것이다. $(-1)^{|A|-|S|} = \mu(A-S)$이고, $T \subset S \subset A$는 $A-S \subset A-T$와 동치이므로, 이 결과는 사실 $$ \sum_{X \subset (A-T)} \mu(X) = 0$$과 같다. 특히, $T \subsetneq A$이므로 $A-T$는 nonempty다. 결과적으로, 증명의 핵심은 모든 nonempty가 아닌 집합 $U$에 대하여 $$\sum_{X \subset U} \mu(X) = 0$$인 것이다. 한편, $U$가 공집합일 때는 좌변의 값이 $1$임 역시 자명하다.

집합이 아니라, 중복 원소가 허용되는 multiset이라면 어떻게 될까? 이때는 $\mu$를 다음과 같이 정의하자.

• $X$가 중복 원소를 갖는 multiset이라면, $\mu(X) = 0$.
• $X$가 중복 원소를 가지지 않는다면, 기존의 경우와 같이 $\mu(X) = (-1)^{|X|}$.

이제 여전히 $g(A) = \sum_{S \subset A} f(S)$이면, $f(A) = \sum_{S \subset A} \mu(A-S) g(S)$임을 보이자.

• 증명의 핵심만 보이면 충분하다. 즉, $\sum_{X \subset U} \mu(X) = 0$인 것만 보이면 충분하다. (물론, $U$가 공집합이면 $1$)
• $U$의 부분 multiset 중 중복 원소가 존재하는 것은 좌변에 기여하지 않는다. 그러니 좌변에 기여하는 것은 중복 원소가 없는 경우 뿐.
• 그러니 단순히 $U$에서 중복 원소를 제거한 집합을 보고, 기존 결과를 그대로 적용하면 증명이 끝난다. 

이제 정수론에서의 mobius function을 이야기 할 수 있다. 

• 자연수 $n$을 소인수분해하여, $p_1^{e_1} p_2^{e_2} \cdots p_k^{e_k}$라 쓰자. 이제 $n$을 multiset $\{p_1, \cdots , p_1, p_2, \cdots , p_2, \cdots p_k \}$로 대응시킨다. 단, $p_i$는 $e_i$번 등장.
• 이제 "집합의 포함 관계"는 "자연수의 나누어떨어짐 관계"로 대응이 되며, mobius function 역시 그대로 대응된다. 
• $\sum_{d|n} \mu(d) = [n=1]$ 역시 $\sum_{S \subset A} \mu(S) = [A = \emptyset] $과 그대로 대응된다.
• 그러니 포함-배제와 mobius inversion, mobius function에 관한 문제풀이는 다를 것이 없다. 
• 즉, 독자가 포함-배제에 관해서 가지고 있는 문제풀이 아이디어, 전략이나 팁들 대부분이 이번 7단원 내용에도 적용될 수 있다.

잠깐 정수론과 무관한 이야기.

• $g(A)$의 값이 주어졌을 때 $f(A)$의 값을 전부 구하는 문제를 다시 생각해보자.
• 순서대로 $f(A) = g(A) - \sum_{S \subsetneq A} f(S)$를 적용했다면, subset을 전부 iterate 해야하므로 $\mathcal{O}(3^n)$의 시간이 걸린다.
• 하지만 $f(A) = \sum_{S \subseteq A} (-1)^{|A|-|S|} g(S)$임을 이용하면, $(-1)^{|A|} f(A) = \sum_{S \subseteq A} (-1)^{|S|} g(S)$임을 얻는다.
• 이제 SOS DP를 이용하여 $f$ 전부를 $\mathcal{O}(2^n n)$ 시간에 계산할 수 있다. 시간 절약하기 좋은 테크닉.  

이제부터는 문제풀이에는 별로 중요하지 않은 이야기.

• "나누어떨어짐 관계" $a|b$와 "포함 관계" $S \subset A$ 사이에는 무언가 비슷한 점이 있는 것 같다.
• $a|a$이며, $a|b$, $b|a$면 $a=b$ (자연수 기준), $a|b$, $b|c$면 $a|c$가 성립한다.
• $S \subset S$이며, $S \subset T$, $T \subset S$면 $S = T$, $S \subset T$, $T \subset U$면 $S \subset U$가 성립한다.
• 이러한 조건을 만족시키는 relation을 partial order라고 한다. mobius function은 이러한 partial order를 갖는 poset에서 전부 정의할 수 있다.

관련 코드는 github에서 찾아볼 수 있다.

이 단원의 내용은 특정 값을 계산하는 알고리즘과 그 template code가 명확하게 있지 않다.

즉, 이 단원에 관한 문제 풀이는 template을 복사/붙여넣기하는 방식으로 풀기 어렵다. 

대신, 문제풀이를 위해서 알아야 하는 이론의 양은 상대적으로 적다. 적은 이론 + 많은 식 정리로 풀리는 문제가 대부분이다.

하지만 "이 단원의 내용이 포함-배제의 원리와 같다"는 점을 제대로 이해하면, 문제풀이가 더 쉬워질 수도 있다. 

그러나 이 점을 제대로 이해하지 않아도 문제는 풀 수 있으니, 이 부분은 부록으로 옮긴다. 그래도 부록을 한 번 읽는 것을 추천.

mobius function $\mu(n)$은 각 자연수에 대하여 다음과 같이 정의된다.

• $n = 1$인 경우, $\mu(n) = 1$이다.
• 만약 소수 $p$가 존재하여 $p^2|n$이라면, $\mu(n) = 0$.
• $n = p_1p_2 \cdots p_k$라면 (단, $p_i$는 서로 다른 소수) $\mu(n) = (-1)^k$

이때, mobius function은 다음과 같은 성질을 갖는다.

• $\mu$는 multiplicative. 즉, $m, n$이 서로소인 자연수라면 $\mu(mn) = \mu(m)\mu(n)$. 이는 정의에서 바로 증명할 수 있다.
• 각 자연수 $n$에 대하여 $\sum_{d|n} \mu(d)$의 값은 $n=1$인 경우 $1$이고, $n \ge 2$인 경우에는 $0$이다. 증명은 뒤로 미룬다.
• mobius inversion. $f(n) = \sum_{d|n} g(d)$라면, $g(n) = \sum_{d|n} \mu(n/d) f(d)$가 성립한다.
• 예시. 오일러 정리 단원에서 $n = \sum_{d|n} \phi(d)$라고 했다. 그러니 $\phi(n) = \sum_{d|n}  d \mu(n/d)$다.
• $\sum_{d|n} d \mu(n/d)$라는 식은, $n$의 소인수분해를 생각한 후 계산하면 $n(1-1/p_1)(1-1/p_2) \cdots (1-1/p_k)$와 같다.
• $\mu(1), \mu(2), \cdots , \mu(n)$을 전부 계산하는 것은 에라토스테네스의 체로 할 수 있다. github 코드 참조.

사실 이게 CP/PS에서 쓰는 성질 전부다. 이걸 잘 조합해서 문제를 풀어내는 것이 어려운 것이다. 

• mobius function이란 말을 들으면 mobius inversion을 생각하는 경우를 봤다. 하지만 실제로 mobius inversion 식 자체는 많이 쓰이지 않는다.

• 사실 mobius function을 이용한 문제풀이의 핵심은 $\sum_{d|n} \mu(d)$가 $n=1$일 때만 $1$, 아니면 $0$이라는 점이다. 

• 계속 강조하지만, mobius inversion은 포함-배제와 같다. 부록 참고.

이제부터 시범을 보인다. 식 조작의 편의를 위하여, Iverson bracket을 도입한다.

• 명제 $P$에 대하여, $[P]$의 값은 $P$가 참이면 $1$, 아니면 $0$이다. 

• 이제 $\sum_{d|n} \mu(d) = [n=1]$이라고 쓸 수 있다. 이 식이 문제풀이의 핵심이다.

BOJ 16409 : Coprime Integers

문제에서 요구하는 값은 $$\sum_{x=a}^b \sum_{y=c}^d [\text{gcd}(x, y)=1]$$이다. $\sum_{d|n} \mu(d) = [n=1]$을 사용하여 Iverson bracket을 없애면, $$ \sum_{x=a}^b \sum_{y=c}^d \sum_{g|x, g|y} \mu(g) $$를 얻는다. 이제 식을 바라보는 관점을 바꾸자. $x, y$에 대해서 iterate 한 다음 $g|x, g|y$를 만족하는 $g$를 보는 것이 아니라, $g$를 고정하고 $g$의 배수가 되는 $x, y$에 대해서 iterate 하자. 이렇게 식을 변형하는 방식을 보통 "double counting"이라고 부른다. 그러면 우리의 식이 $$ \sum_{g=1}^{\text{min}(b, d)} \mu(g) \cdot ([a, b] \text{ 에 존재하는 } g \text{의 배수 개수}) \cdot ([c, d] \text{ 에 존재하는 } g \text{의 배수 개수}) $$가 된다. 물론, $[a, b]$에 존재하는 $g$의 배수 개수는 $\lfloor b/g \rfloor - \lfloor (a-1)/g\rfloor$이고, $[c, d]$에 대해서도 마찬가지다. 즉, 답은 $$ \sum_{g=1}^{\min(b,d)} \mu(g) \cdot (\lfloor b/g \rfloor - \lfloor (a-1)/g \rfloor) \cdot (\lfloor d/g \rfloor - \lfloor (c-1)/g \rfloor)$$가 된다. 이제 $\mu$를 전처리하면 답을 쉽게 계산할 수 있다.

여기서도 볼 수 있지만, 약수에 대해 iterate 하는 것보다 배수에 대해 iterate 하는 것이 더 쉽다는 점을 많이 이용한다. (아래에 예시를 더 넣어두겠다)

시범을 한 번 더 보이겠다. 이번에는 $$\sum_{x=1}^n \sum_{y=1}^n \text{gcd}(x, y)$$를 구해보자. 이를 위해서 $\text{gcd}(x, y) = d$인 $(x, y)$의 개수를 각 $d$에 대해서 구해보자. 이 값을 구하려면, $$ \sum_{x=1}^n \sum_{y=1}^n [\text{gcd}(x, y)=d]$$를 구하면 된다. 그러니, 기존 식의 값은 $$ \sum_{d=1}^n d \cdot \sum_{x=1}^n \sum_{y=1}^n [\text{gcd}(x, y)=d]$$와 같다. $\text{gcd}(x, y)=d$려면, $x, y$ 역시 $d$의 배수여야 하니 $x=da$, $y=db$라 쓸 수 있고, 이때 $1 \le a, b \le \lfloor n/d \rfloor$와 $\text{gcd}(a, b)=1$이다. 즉, 식은 $$ \sum_{d=1}^n d \cdot \sum_{a=1}^{\lfloor n/d \rfloor} \sum_{b=1}^{\lfloor n/d \rfloor} [\text{gcd}(a, b)=1]$$과 같다. 이제 바로 앞의 문제의 결과를 적용하면 $$ \sum_{d=1}^n \sum_{g=1}^{\lfloor n/d \rfloor} d \mu(g) \cdot \lfloor \frac{n}{gd} \rfloor^2$$을 얻는다. 여기서 자연수 $n, a, b$에 대하여 $\lfloor n/(ab) \rfloor = \lfloor \lfloor n/a \rfloor / b \rfloor$임을 이용하였다.

• Finish 1. $\mu$를 전처리하자. 위 합을 for 문으로 계산하면, 시간복잡도는 $\sum_{d=1}^n \lfloor n/d \rfloor = \mathcal{O}(n \log n)$ (by Harmonic Sequence)

• Finish 2. $gd$라는 값이 심상치 않다. $gd = T$라고 치환을 한 다음, $T$를 기준으로 계산을 해보자. 그러면 $$ \sum_{T=1}^n \lfloor n/T \rfloor^2 \cdot \sum_{d|T} (d \cdot \mu(T/d))$$라는 식을 얻는다. 그런데 $\phi(T) = \sum_{d|T} (d \cdot \mu(T/d))$임을 앞에서 언급했다. 그러니 식은 $ \sum_{T=1}^n \lfloor n/T \rfloor^2 \phi(T)$와 같다. 

위 문제의 mobius function을 사용하지 않는 풀이를 소개한다. $cnt[d]$를 $\text{gcd}(x, y) = d$인 $x, y$의 개수라고 하자. 사실 우리는 $\text{gcd}(x, y) = d$인 $x, y$의 개수는 세기 어려워해도, $\text{gcd}(x, y)$가 $d$의 배수인 것은 쉽게 셀 수 있다. $x, y$가 각각 $d$의 배수이면 충분하므로, 정확히 $\lfloor n/d \rfloor^2$개 존재할 것이다. 이는 사실 $$ cnt[d] + cnt[2d] + \cdots + cnt[d \cdot \lfloor n/d \rfloor]  = \lfloor n/d \rfloor^2$$이 성립한다는 것이다. 이를 통해서, $cnt$ 배열을 모두 구할 수 있다. 다음과 같은 코드를 생각하면 이해가 될 것이다.

더욱 놀라운 점은, Harmonic Sequence의 성질 때문에 이 코드 역시 시간복잡도가 $\mathcal{O}(n \log n)$이란 것이다.

• 이 풀이를 보면서, "배수에 대해 iterate 하는 것이 효율적이다"라는 말의 의미를 느낄 수 있다.

• 또한, 이 풀이가 굉장히 포함-배제의 원리 느낌이 난다는 것을 느낄 수 있을 것이다. 

물론, 실제로 문제를 해결하려면 $\sum d \cdot cnt[d]$를 계산하면 끝이다. 

이제부터 "배수에 대해 iterate 하는 것"과 double counting을 사용한 간단하지만 중요한 예제를 몇 개 제시한다.

• $\tau(n)$을 $n$의 약수의 개수라고 정의한다. $\sum_{i=1}^n \tau(i)$를 어떻게 계산할까?  

• $\sigma(n)$을 $n$의 약수의 합이라고 정의한다. $\sum_{i=1}^n \sigma(i)$를 어떻게 계산할까?

$\sum_{i=1}^n \tau(i)$는 결국 $g|i$, $1 \le i \le n$을 만족하는 순서쌍 $(g, i)$의 개수이다. 순서쌍의 개수를 $i$ 기준으로 세면, 각 $i$에 대해 조건을 만족하는 $g$가 $\tau(i)$개 있기 때문이다. 순서쌍의 개수를 $g$ 기준으로 세면, 조건을 만족하는 $i$가 $\lfloor n/g \rfloor$개 있음을 알 수 있다. 그러니 순서쌍의 개수는 $\sum_{g=1}^n \lfloor n/g \rfloor$이다. 그러므로, $$ \sum_{i=1}^n \tau(i) = \sum_{i=1}^n \lfloor n/i \rfloor $$이다. 이는 결국 "각 정수에 대해 약수의 개수를 센다"에서 "각 정수에 대해 배수의 개수를 센다"로 문제를 바꿔 생각해서 얻은 결과이다. 위 식의 우변은 자명하게 $\mathcal{O}(n)$에 계산할 수 있다. 또한, $\lfloor n/i \rfloor$로 가능한 값의 개수가 $\mathcal{O}(\sqrt{n})$개이며, 그 값을 갖게 하는 $i$의 범위 역시 계산할 수 있음은 이미 1단원에서 언급하였다. 이를 사용하면, 우변의 값을 $\mathcal{O}(\sqrt{n})$에 계산할 수 있다. 이러한 방식은 후에 multiplicative function을 다룰 때 더 나온다.

$\sum_{i=1}^n \sigma(i)$ 역시, 같은 방법으로 생각하면 $\sum_{g=1}^n g \lfloor n/g \rfloor$임을 얻는다. 이는 $\mathcal{O}(n)$에 자명하게 계산할 수 있다. 특정 구간에 속한 정수의 합, 즉 $\sum_{i=a}^b i$를 쉽게 계산할 수 있으므로, $\lfloor n/g \rfloor$로 가능한 값이 $\mathcal{O}(\sqrt{n})$개임을 이용하여 우변의 값을 $\mathcal{O}(\sqrt{n})$에 계산할 수 있다. github 코드 참고.  

관련 코드는 github에서 찾아볼 수 있다.

이 글에서 다룰 내용의 많은 부분에서 완벽한 수학적 논의는 생략할 것이다.

Miller-Rabin 소수 판별법

이 부분에 대해서는 이미 좋은 자료가 있습니다. https://casterian.net/archives/396

우리는 주어진 $n$이 소수임을 확인하고자 한다. 먼저 $n-1 = d \cdot 2^s$인 정수 $s, d$를 찾는다.

만약 $n$이 소수라면, 각 자연수 $1 \le a < n$에 대하여 다음이 성립한다.

• [P] $a^d \equiv 1 \pmod{n}$이거나, $r=0, 1, \cdots , s-1$ 중 적어도 하나에 대해 $a^{2^rd} \equiv -1 \pmod{n}$

증명을 해보자. $a^{n-1}-1 = a^{d \cdot 2^s}-1 = (a^d-1)(a^d+1)(a^{2d}+1)(a^{2^2d}+1) \cdots (a^{2^{s-1}d}+1)$이다.

그런데 페르마 소정리에 의하여 좌변은 $0 \pmod{n}$이므로, 우변 역시 $n$의 배수임을 알 수 있다.

$n$이 소수이므로, 여러 자연수의 곱이 $n$의 배수라는 것은 그 중 하나 이상이 $n$의 배수임을 의미한다. 

그러므로, $a^d \equiv 1 \pmod{n}$이거나 적당한 $0 \le r <s$가 있어 $a^{2^rd} + 1 \equiv 0 \pmod{n}$이 성립한다.

$n$이 주어졌을 때, $s, d$의 값을 계산하는 것은 쉽게 할 수 있다. 

또한, $1 \le a < n$을 만족하는 $a$가 주어졌을 때, [P]가 성립하는지 여부도 역시 빠르게 확인할 수 있다. 

이제 다음과 같은 알고리즘을 하나 고안할 수 있다.

• 테스트 횟수 $k$를 하나 정한다. $k$를 어떻게 잡는가에 대한 이야기는 나중에 한다.
• $1 \le a < n$을 만족하는 $a$를 랜덤하게 하나 잡자.
• [P]가 성립하는지 확인하자. 만약 성립하지 않는다면, $n$이 소수가 아니라는 뜻이고 알고리즘 종료.
• [P]가 성립한다면, 다시 처음으로 돌아가서 알고리즘을 반복한다. $k$번에 걸쳐 항상 [P]가 성립한다면 소수라고 판단하고 종료.

만약 $n$이 소수라면, $a$를 어떻게 잡더라도 [P]가 성립할 것이니, 무조건 소수라고 판단할 것이다.

$n$이 합성수라면 어떨까? $n$이 소수라면 [P]가 성립한다고 했지, $n$이 합성수면 무조건 [P]가 거짓이라고 하지는 않았다.

위 알고리즘이 정당한 소수 판별 알고리즘이 되려면, $n$이 합성수일 경우 높을 확률로 $n$이 소수가 아니라는 결론을 내려야 한다.

즉, $k$번의 테스트를 거쳤을 때 [P]가 거짓이 되는 $a$의 값이 매우 높은 확률로 한 번 이상 뽑혀야 한다. 이는 결국

• $n$이 합성수인 경우, [P]가 거짓이 되는 $a$가 얼마나 많은가? 

라는 문제와 동치이다. 이제 중요한 결론을 증명하지 않고 제시한다. 증명은 구글 ㄱㄱ

• $n$이 합성수인 경우, [P]가 거짓이 되는 $a$는 적어도 $3n/4$개 존재한다.
• 즉, $n$이 합성수인 경우, 한 번의 테스트를 통과할 확률은 $1/4$ 이하이다. 

이는 $k$번의 독립적인 테스트를 거쳤을 때, 합성수가 모든 테스트를 통과할 확률이 $4^{-k}$ 이하라는 뜻이다. 

그러니 $k=20$ 정도로 잡아주면 매우 높은 확률로 정답을 반환하는 좋은 소수 판별 알고리즘이 될 수 있다.

하지만 랜덤이 싫을 수 있다. 랜덤을 제거한 Miller-Rabin 알고리즘의 여러 variant가 존재한다.

• Generalized Riemann Hypothesis를 가정하면, $2 \le a < 2(\ln n)^2$를 전부 시도하면 정확하게 소수 판별을 할 수 있다.
• $n < 2^{32}$라면, $a = 2, 7, 61$인 경우만 따져보면 된다. 이때, $n = 2, 7, 61$인 경우를 따로 처리해야 함에 주의하라.
• $n < 2^{64}$라면, $a$로 $2$ 이상 $40$ 이하의 소수만 따져보면 된다. 이때, $n$ 역시 $2$ 이상 $40$ 이하의 소수인 경우에 주의하라.    

랜덤을 제거한 variant가 아마 실제 성능도 좋을 것이다. $n < 2^{64}$의 variant에서도 테스트를 12번만 하기 때문이다.

구현 노트

• 구현이 까다로운 알고리즘은 아니지만, C++의 int 범위를 넘어가는 modulus에서 여러 연산을 해야한다. 
• 특히 까다로운 것은 곱셈이다. $a, b, n$이 모두 64bit 정수일 때, $ab \pmod{n}$을 계산해야 한다.
• 이를 위해서 쓸 수 있는 첫 번째 방법은 이집트 곱셈, 즉 이진법을 사용한 곱셈이다. 이를 위해서 로그 시간이 필요하다.
• 두 번째 방법은 __int128을 동원하는 것이다. $ab$의 값을 __int128에서 계산한 뒤, $n$에 대한 나머지를 구한 후 64bit 정수로 반환한다. 
• $n$이 작으면, $\mathcal{O}(\sqrt{n})$ 알고리즘이 Miller-Rabin 알고리즘보다 빠를 수 있다. 최적화가 필요한 경우 참고하자.
• 더 적은 테스트 횟수로 $2^{64}$ 미만 자연수의 소수 여부를 판별할 수 있다. 이때, $a$ 값이 소수가 아닌 경우 "$a$의 약수인 소수"들을 특별히 고려해야 한다. 
• 위 두 가지 포인트를 추가로 고려하여 작성한 Miller-Rabin 구현체를 원한다면, 맨 위에 소개한 casterian님의 블로그를 참고하라.

Pollard-Rho 소인수분해

이 부분 역시 좋은 자료가 있습니다. https://aruz.tistory.com/140

핵심 아이디어 1. Floyd's Cycle Detection Algorithm

잠깐 전형적인 CP/PS 토픽인 그래프로 넘어오자. functional graph는 각 정점이 outdegree 1을 갖는 directed graph다.

정점 $v$가 있으면 그 정점에서 나가는 간선이 정확히 하나이니, 그 간선을 따라 도달하는 정점을 $f(v)$라 정의할 수 있다. 

즉, $f(v)$는 $v$의 "다음 정점"이다. 이제 $v \rightarrow f(v) \rightarrow f(f(v)) \rightarrow \cdots$라는 path를 생각하자.

결국 이 경로는 하나의 일직선을 그리다가, 이미 방문한 정점을 다시 방문하면서 cycle에 들어가게 될 것이다. 이 cycle을 어떻게 찾을까?

물론 여러 방법이 있겠으나, Floyd의 "토끼와 거북이" 알고리즘이 매우 효율적이다. 

• 잠시 $f^n(v) = f(f^{n-1}(v))$로 정의하자. 즉, $f^n(v)$는 $v$의 "$n$번째 뒤"의 정점이다.
• 토끼 한 마리와 거북이 한 마리를 $v$에 두자. 각 동물들의 위치를 변수 $r, t$로 명명하고 관리하겠다.
• 각 단계마다, 토끼는 두 칸, 거북이는 한 칸 이동시킨다. 즉, $r \leftarrow f(f(r))$, $t \leftarrow f(t)$를 반복한다.
• 만약 $r = t$인 시점이 오면, 사이클을 찾은 것이다. 물론, 맨 처음에 $r=t=v$인 경우는 제외하고 생각한다.

• 왜 사이클을 찾은 것인가? $k$번의 단계를 거친 후에 $r = t$임을 확인했다고 가정하자.
• 이는 $f^k(v) = f^{2k}(v)$임을 의미한다. 즉, $k$번째 뒤 정점과 $2k$번째 뒤 정점이 같다는 것을 의미한다.
• 이제 $k$번째 정점인 $f^k(v)$부터 시작하여, 다음 정점으로 이동하는 것을 반복한다. 
• 다시 시작한 $f^k(v)$에 도착할 때까지 이를 반복하면, 지금까지의 경로가 완벽하게 찾고자 하는 사이클이 된다.

• 왜 사이클이 찾아지는가? $i < j$이고 $f^i(v) = f^j(v)$인 $i, j$가 존재할 때, (즉 사이클이 있을 때) 자연수 $k$가 있어 $f^k(v) = f^{2k}(v)$임을 보이면 된다.
• "한 칸씩 격차가 벌어지는 걸 결국 사이클 길이만큼 격차가 벌어져서 다시 만났다"고 생각하면 직관적이며, 증명도 똑같다.
• 속도 차이가 나는 두 사람이 400m 트랙에서 달렸을 때, 실제로 달린 거리가 400m 차이가 나서 만날 수 있는 것과 같은 원리.
• 자세한 증명을 원한다면, $t(j-i) \ge i$인 자연수 $t$를 하나 잡고, $k=t(j-i)$라고 정의해보면 증명됨을 확인할 수 있다.

• 왜 효율적인가? 일단 시간복잡도가 (사이클 도달에 필요한 거리) + (사이클의 길이) 정도로 매우 효율적이다.
• 공간복잡도도 효율적이다. 각 정점 $v$에 대해서 $f(v)$를 계산하는데 드는 메모리를 제외하면, $\mathcal{O}(1)$만큼의 메모리가 추가적으로 필요하다.
• 특히, 여러 경우에서 $f(v)$는 진짜 간선으로 표현되는 "다음 정점"이 아니라, 말 그대로 $v$에 대한 함수가 된다. 
• 이는 아래의 경우도 마찬가지. 애초에 이름이 functional graph인 이유가 있다. 이 경우, $f(v)$를 계산하는데 드는 메모리도 $\mathcal{O}(1)$.

핵심 아이디어 2. Pseudorandom Sequence. Birthday Paradox

Birthday Paradox는 $\{1, 2, \cdots , n\}$에서 $\sqrt{n}$급 개수의 원소를 독립적으로 뽑으면, 같은 원소를 두 번 이상 뽑을 확률이 상당히 높다는 결과이다.

이는 $n=365$일 때, 대충 20~30명의 학급에서 생일이 겹치는 경우가 많다는 내용으로 알려진 결과일 것이다.

이제 우리의 목표는 주어진 합성수 $n$을 소인수분해 하는 것이다. $p$를 $n$의 가장 작은 소인수라고 하자. 우리는 $n$은 알고 있으나 $p$는 모른다.

이제 다음과 같은 방식으로 랜덤한 수열을 잡는다. $x, c$ 각각을 $1$ 이상 $n$ 미만 랜덤한 자연수라고 두자.

$x_0 = x$, $x_i = (x_{i-1}^2 + c) \pmod{n}$으로 수열을 잡자. 여기서 짚고 넘어가야 할 부분은

• 수열 $\langle x_k \rangle$은 pseudorandom. 즉, 완전 랜덤한 것은 아니지만 랜덤에 가깝다.
• $x_i$의 값은 $x_{i-1}$의 값에 의하여 완전히 결정된다. 그러니, $x_{i-1} \rightarrow x_i$를 간선으로 생각하는 functional graph를 생각할 수 있다.
• 즉, 이 수열의 값을 정점으로 갖는 functional graph를 생각할 수 있고, 마찬가지로 cycle 개념도 생각할 수 있다.
• $x_i \equiv (x_{i-1}^2 + c) \pmod{p}$ 역시 성립한다. 그러니, $y_i = x_i \pmod{p}$로 정의하면, $y_i \equiv (y_{i-1}^2 + c) \pmod{p}$이다.
• 현재 계산할 수 있는 것은 $x_i$이며, $y_i$는 실제로는 $p$를 모르기 때문에 계산할 수 없는 값이다.

이제 본격적으로 알고리즘을 설계해보자.

• $\text{gcd}(x_{2k} - x_k, n)$이라는 값을 생각해보자. 이 값을 계산하려면 $x_{2k}, x_k$의 값이 필요하며, 이는 "토끼와 거북이"의 방식으로 관리할 수 있다.
• $y_k = y_{2k}$라면, (즉 수열 $y$에서 cycle이 발견되었다면) $\text{gcd}(x_{2k} - x_k, n)$은 $p$의 배수가 된다. 
• Birthday Paradox에 의해, 이렇게 수열 $y$에서 cycle이 발견되려면 평균적으로 $\mathcal{O}(\sqrt{p})$ 정도 index가 되어야 한다. 
• 비슷하게, 수열 $x$에 대해서 cycle이 발견되려면 평균적으로 $\mathcal{O}(\sqrt{n})$ 정도 index가 되어야 한다.
• 즉, 평균적으로 $\mathcal{O}(\sqrt{p})$ 정도의 계산량 이후로 $\text{gcd}(x_{2k}-x_k, n)$이 $1$이 아닌 시점이 오며, 이때 $\text{gcd}(x_{2k}-x_k , n)$은 높은 확률로 $n$이 아니다.
• 이는 $\text{gcd}(x_{2k}-x_k, n)$이 $n$의 "nontrivial divisor"가 된다는 뜻이다. 이를 반복하면 소인수분해를 얻는다.

$n$의 non-trivial divisor를 얻기까지 평균적으로 $\mathcal{O}(\sqrt{p}) = \mathcal{O}(n^{1/4})$번의 GCD 연산이 필요하다.

구현 노트.

• 이 알고리즘에서는 random을 사용한다. C++의 rand()를 사용하기 보다는, mt19937 등을 사용하는 것을 추천한다.
• 앞에서도 언급했지만 $n$이 합성수라는 조건이 붙는다. 미리 Miller-Rabin 소수 판별을 통해 $n$이 소수인지 판별을 해야한다.
• $\text{gcd}(x_{2k}-x_k, n)=n$일 확률이 낮다고 했지, 불가능하다고 하지는 않았다. 이 경우, 시작하는 값 $x, c$를 다시 잡고 알고리즘을 다시 돌리자.
• $\text{gcd}(x_{2k}-x_k, n)=d$가 $1, n$이 아니라면, 문제를 $d$와 $n/d$를 각각 소인수분해 하는 것으로 바꿀 수 있다. 이제 재귀적으로 문제를 풀자.
• 앞에서 언급한 64bit 정수의 곱셈에 대한 내용이 여기서도 적용된다. 이집트 곱셈을 쓰거나, __int128을 사용하자.

관련 코드는 github에서 찾아볼 수 있다.

이 글에서는 팩토리얼과 이항계수를 $n$으로 나눈 나머지를 계산하는 방법에 대하여 다룬다.

• 중국인의 나머지 정리를 기반으로 한 아이디어는 여기서도 등장한다.
• $n = p_1^{e_1} p_2^{e_2} \cdots p_k^{e_k}$에 대한 나머지를 구하려면, 각 $p_i^{e_i}$에 대한 나머지를 구하면 된다.
• 그 후, 단순히 중국인의 나머지 정리를 사용하여 합쳐주면 된다. 그러니 문제를 prime power에 대해서만 해결하면 된다.

이 글에서 다룰 내용은 상당한 난이도를 갖는 강화판들이 있다. 이에 대한 것은 결과만 언급하고 링크를 걸겠다.

먼저 팩토리얼, 이항계수를 $p$로 나눈 나머지를 계산하는 방법을 알아보자. $p$는 소수.

팩토리얼

주제 1. 주어진 $n, p$에 대하여 $n! \pmod{p}$ 계산하기. 

단순히 $n! \pmod{p}$를 계산한다고 가정하자. $n \ge p$인 경우, 값은 $0$이다. 

그러니, $n < p$인 경우만 고려해도 충분하다. $n$이 작으면, $\mathcal{O}(n)$ 알고리즘이 자명하다.

$n, p$가 모두 $10^{9}$-scale로 큰 경우는, $\mathcal{O}(\sqrt{p}\log p)$ 시간에 해결할 수 있으나 상당히 어려운 문제다. 이 두 문제와 이 글 참고. 

주제 2. $1!, 2!, \cdots , n!$ 각각을 $p$로 나눈 나머지와 각각의 $\pmod{p}$에 대한 modular inverse 계산하기.

modular inverse를 계산하는 이유는 역시 이항계수를 계산하기 위해서이다. 이 주제는 수많은 카운팅 문제들에서 사용된다.

$1!, 2!, \cdots n!$ 각각을 $p$로 나눈 나머지는 점화식 $k! \equiv (k \cdot (k-1)!) \pmod{p}$를 이용하여 $\mathcal{O}(n)$ 시간에 계산할 수 있다.

팩토리얼들의 modular inverse를 구하는 것은 정말 여러 방법이 있다.

• 방법 1. 우선 각각의 modular inverse를 그냥 구하는 방법이 있다. 확장 유클리드 알고리즘을 쓰면 된다. 이 경우 $\mathcal{O}(n \log p)$의 시간 소요.
• 방법 2. $1, 2, \cdots, n$ 각각의 modular inverse를 $\mathcal{O}(n)$에 구하는 방법이 있다. 이를 사용하면 팩토리얼의 modular inverse도 똑같이 $\mathcal{O}(n)$에 ok.
• 알고리즘의 설명을 위해 프로그래밍 같은 notation을 잠시 도입한다. $i$의 modular inverse를 $inv[i]$라 하고, $a$를 $b$로 나눈 나머지를 $a \% b$라 하자.
• Fact. $2 \le i<p$이면, $inv[i] \equiv inv[p \% i] \cdot (p - \lfloor p/i \rfloor) \pmod{p}$가 성립한다. ($p$가 소수임이 필요하다)
• 증명. $p = qi + r$이라 하자. $p \% i = r$이며, $\lfloor p/i \rfloor = q$이다.
• 이제 $i \cdot inv[p \% i] \cdot (p - \lfloor p/i \rfloor) \equiv i \cdot inv[r] \cdot (p-q) \equiv inv[r] \cdot (-qi) \equiv inv[r] \cdot r \equiv 1 \pmod{p}$이므로 증명 끝.
• 여기서 $p \% i < i$이며, $inv[1]=1$임을 생각하자. 위 Fact는 결국 $inv$ 배열을 채우기 위한 DP 식이 된다. 
• $k!$의 modular inverse는 $(k-1)!$의 modular inverse와 $k$의 modular inverse를 곱해서 얻을 수 있다.
• 방법 3. 매우 쉽고 똑똑한 아이디어다. $n!$의 modular inverse를 먼저 확장 유클리드 알고리즘으로 구하자.
• 그 후, $k!$의 modular inverse가 $(k+1)!$의 modular inverse와 $(k+1)$의 곱임을 이용한다. 역순으로 계산하는 것이다.
• 이 방법으로 $1, 2, \cdots n$ 각각의 modular inverse도 $\mathcal{O}(n)$ 시간에 구할 수 있다. 
• $k$의 modular inverse를 구하려면, $k!$의 modular inverse와 $(k-1)!$을 곱하면 된다.

방법 2, 3은 둘 다 $\mathcal{O}(n)$ 알고리즘이다. 하나만 알아도 충분하지만, 극한의 시간 최적화가 필요한 경우에는 둘 다 알아야 한다.

• 만약 $1, 2, \cdots , n$ 각각의 modular inverse를 얻는 것이 주요 목적이라면, 방법 2를 사용하는 것이 맞다.
• 만약 $1!, 2!, \cdots, n!$ 각각의 modular inverse를 얻는 것이 주요 목적이라면, 방법 3을 사용하는 것이 맞다.
• 위에서 언급한 modular inverse들을 전부 구하는 것이 목적이라면, 어느 것을 사용해도 문제가 없다.
• 이유를 파악하기 위해서는 각 방법들이 언급한 값들을 계산하기 위해서 몇 번의 나머지 연산을 하는지 생각해보자.

개인적인 추천은 그냥 방법 2, 3을 다 알고 있는 것이다. 방법 2의 원리까지는 이해할 필요가 없으니, 코드만 알아두면 된다.

주제 3. $n! = p^e m$라고 쓸 경우 (단, $m$은 $p$의 배수가 아니다) $e$와 $m \pmod{p}$를 계산하기.

이는 일종의 "$n!$의 0이 아닌 마지막 자리 구하기"와 같은 부류의 문제라고 볼 수 있다. 이 문제를 $p$진법에서 푸는 것이기 때문.

이 문제가 필요한 이유는 이항계수의 계산이다. $n!$이 $p$의 배수라고 해서 $\binom{n}{k} =n!/(k!(n-k)!)$이 $p$의 배수일 이유는 없다.

$n < p$인 경우는 $e=0$이며, $m \pmod{p}$는 단순히 $n! \pmod{p}$이다. 그러니 앞선 논의를 그대로 사용할 수 있다.

사실 이 말은 주제 3의 문제가 주제 1의 문제를 포함하고, 더 어려운 문제라는 이야기다. 그러니 $p$가 크면 이 문제도 풀기 힘들 것이라고 예측할 수 있다.

여기서 다룰 경우는 $p$가 $\mathcal{O}(p)$ 알고리즘이 돌아갈 정도로 작고, $n$은 매우 클 수 있는 경우이다.

여기서도 문제를 축소하는 방법을 사용할 것이다. $n=pk+r$, $0 \le r < p$라고 하자. 아래 그림과 함께 보자.

• $1$부터 $n$까지의 수 중, $p$의 배수가 아닌 것에 먼저 집중하자.
• $[1, p-1], [p+1, 2p-1], \cdots , [p(k-1)+1, pk-1], [pk+1, pk+r]$ 형태의 구간으로 나누어 볼 수 있을 것이다.
• 첫 $k$개의 구간을 보자. 각 구간에 있는 정수들을 곱한 값은 $(p-1)! \pmod{p}$이다. (사실 이는 $-1 \pmod{p}$다. 윌슨의 정리)
• $[pk+1, pk+r]$의 구간에 있는 정수들을 곱한 값은 $r! \pmod{p}$이다. 
• 이제 $p$의 배수들인 것에 집중하자. $p, 2p, \cdots , kp$가 있다. 여기서 $p$를 $k$개 뽑아낼 수 있다.
• 그러면 남는 것은 $1, 2, \cdots , k$고, 이는 정확히 $k!$과 같다. 여기서 우리는 문제를 $k!$에 대한 문제로 축소시켰다.
• $k!$을 $p^{e'} m'$으로 쓸 때, $e'$와 $m' \pmod{p}$의 값을 알아냈다고 가정하자. 이제 $n! = p^e m$이라 쓸 때 $e$와 $m \pmod{p}$를 구할 수 있을까?
• 가능하다. 위 논의를 정리하면, $e = e' + k$, $m \equiv m' \cdot (p-1)!^k \cdot r! \pmod{p}$를 반환하면 된다.
• 앞에서도 언급했지만 $(p-1)! \equiv -1 \pmod{p}$이므로, $m$을 계산할 때는 $k$의 홀짝성만 봐도 충분하다.
• 당연히 구현은 재귀적으로 하는 것이 편하다. 또한, 알고리즘을 돌리기 전 $1!, 2!, \cdots , (p-1)! \pmod{p}$를 전처리하는 것이 좋다. 
• 문제가 얼마나 축소되는가? : $n$에서 $k$로 축소시켰고, $k \le n/p \le n/2$이다. 그러니 반 이상 감소하고, 로그 시간에 ok.
• 여러분도 눈치챘겠지만, 여기서 $e = \lfloor n/p \rfloor + \lfloor n/p^2 \rfloor + \cdots $임을 얻는다.

이항계수

이제 이항계수 $\binom{n}{k}$를 $\pmod{p}$에 대하여 계산하는 방법을 보자. $\binom{n}{k} = \binom{n}{n-k}$이니, $k \le n/2$만 보자. 여러 경우가 있다.

Case 1. $k<p$이며 $k$가 작은 경우. 

이 경우, $\displaystyle \binom{n}{k} = \frac{n \cdot (n-1) \cdots (n-k+1)}{1 \cdot 2 \cdots k}$임을 이용하여 직접 계산을 할 수 있다.

분자, 분모를 모두 $\pmod{p}$로 구한 후, 분모의 modular inverse를 구해 나눗셈을 곱셈으로 대체한다. 

물론, $k<p$이므로 $\text{gcd}(k!,p)=1$이며, 그러니 확장 유클리드 알고리즘으로 modular inverse를 구할 수 있다.

사실, $\displaystyle \binom{n}{k} = \frac{n-k+1}{k} \binom{n}{k-1}$이 성립함을 이용하여 더 많은 것을 할 수 있다.

위 식을 차례대로 적용하여, 단순히 $\binom{n}{k} \pmod{p}$를 구하는 것이 아니라 $\binom{n}{0}, \binom{n}{1}, \cdots , \binom{n}{k}$ 각각을 $p$로 나눈 나머지를 구할 수 있다.

Case 2. 작은 $n, k$에 대하여 전부 전처리를 하고 싶은 경우.

이 경우, 파스칼의 삼각형을 이용하여 DP와 같은 방식으로 계산을 할 수 있다. 정수론이 필요없는 영역이다.

Case 3. $p$가 작은 경우, 즉 $\mathcal{O}(p)$ 알고리즘이 동작할 수 있는 경우.

먼저 $1!, 2!, \cdots , (p-1)!$ 각각을 $p$로 나눈 나머지와, 그들의 modular inverse를 전처리하자.

방법 1. Lucas 정리 사용

많이 사용되는 방법 중 하나다. 다음 결과를 증명하지 않고 제시한다. 증명은 위 링크에 걸려있다. 

정리. $m = m_km_{k-1} \cdots m_0$, $n=n_k n_{k-1} \cdots n_0$가 $m, n$의 $p$진법 전개라면, $\displaystyle \binom{m}{n} \equiv \binom{m_k}{n_k} \cdot \binom{m_{k-1}}{n_{k-1}} \cdots \binom{m_0}{n_0} \pmod{p}$다.

증명의 아이디어가 (Frobenius Endomorphism) 문제 풀이에 사용된 걸 본 적이 없는 것은 아니나, 모두 매우 어려운 문제들이다.

이 결과는 문제를 단순히 $m, n$이 $p$ 미만일 경우에 $\binom{m}{n} \pmod{p}$를 계산하는 문제로 만들어준다.

$m<n$인 경우 결과는 $0$이며, 그렇지 않은 경우 $m!/(n!(m-n)!)$을 계산하여 답을 얻을 수 있다.

이제는 어느 정도 당연한 이야기지만, 나눗셈은 modular inverse를 곱하는 것으로 계산한다.

전처리를 했기 때문에, 각 $\binom{m_i}{n_i} \pmod{p}$를 구하는 것에는 $\mathcal{O}(1)$의 시간이 필요하다.

그러니 시간복잡도는 $\mathcal{O}(k) = \mathcal{O}(\log_p m)$이 되겠다. 매우 빠르고 간단한 방법이다. 

방법 2. 앞서 다룬 주제 3 사용

주제 3을 이용하여도 문제를 해결할 수 있다. $\displaystyle \binom{n}{k} = \frac{n!}{k!(n-k)!}$이라 하자.

이제 $n! = p^{e_1} m_1$, $k! = p^{e_2} m_2$, $(n-k)!=p^{e_3} m_3$라고 쓰자. $m_1, m_2, m_3$는 $p$의 배수가 아니다.

이때, 주제 3에서 다룬 알고리즘으로 $e_1, e_2, e_3$와 $m_1, m_2, m_3 \pmod{p}$를 로그 시간에 계산할 수 있다.

그러면 이제 $\displaystyle \binom{n}{k} =  p^e m$이라 쓸 때, $e$와 $m \pmod{p}$를 구할 수 있을까?

있다. $e = e_1 - e_2 - e_3$, $m \equiv m_1 (m_2m_3)^{-1} \pmod{p}$가 성립하기 때문이다. $a^{-1}$은 $a$의 잉여역수.

만약 $e>0$이라면 주어진 이항계수는 $p$의 배수가 된다. 그렇지 않다면, 이항계수를 $p$로 나눈 나머지는 $m \pmod{p}$와 같다. 끝.

이제 $\pmod{p}$의 논의를 끝냈다. 이제 본격적으로 prime power에 대하여 논의해보자. 

팩토리얼

주제 1. 주어진 $n, p$에 대하여 $n! \pmod{p^e}$ 계산하기.

$n \ge pe$이면 $n!$에서 $p$의 배수가 $e$개 이상 들어가므로, $n! \pmod{p^e}$는 $0$이 된다. 

$n$이 작으면 단순한 $\mathcal{O}(n)$ 알고리즘이 작동한다는 점은 기존 주제 1의 논의에서와 동일하다. 

주제 2. $1!, 2!, \cdots , n!$ 각각을 $p^e$로 나눈 나머지와 각각의 $\pmod{p^e}$에 대한 modular inverse 계산하기.

문제가 잘 정의되려면 역시 $n < p$를 가정해야 한다. 이 경우, 방법 1과 방법 3이 잘 적용된다. 

주제 3. $n! = p^t m$라고 쓸 경우 (단, $m$은 $p$의 배수가 아니다) $t$와 $m \pmod{p^e}$를 계산하기.

$t$를 계산하는 것은 기존 주제 3과 애초에 목표가 동일하다. $m \pmod{p^e}$를 계산하는 것이 목표. 이 역시 기존 주제 3과 비슷한 방법으로 할 수 있다.

• 쉽게 가자. $n=pk+r$, $0 \le r < p$라고 쓰자. 
• $1$부터 $n$까지의 수 중, 이번에는 $p$의 배수인 것에 먼저 집중해보자.
• 여기서는 $p$가 $k$개 추가되고, $k!$이 튀어나온다. 이는 문제의 "축소"된 형태이다. 여기까진 기존 주제 3과 동일.
• 이제 $p$의 배수가 아닌 것을 봐야 한다. 잠시 $n=p^ek' + r'$이라고 하자. $0 \le r' < p^e$. 
• $\mathcal{O}(p^e)$의 시간을 투자하여, DP를 이용해 다음 배열 $val[i]$를 ($0 \le i < p^e$) 채워놓자. 미리 전처리하는 것이다.
• $val[i]$ : $1$ 이상 $i$ 이하 자연수 중 $p$의 배수가 아닌 것을 곱한 값을 $p^e$로 나눈 나머지.
• $[1, p^e-1]$ 중 $p$의 배수가 아닌 것부터 $[p^e(k'-1)+1, p^ek'-1]$ 중 $p$의 배수가 아닌 것이 먼저 있을 것이다.
• 그 후 $[p^ek'+1, p^ek'+r']$ 중 $p$의 배수가 아닌 것이 여기에 있다. 이를 모두 종합하면?
• $val[p^e-1]^{k'} \cdot val[r'] \pmod{p^e}$이라는 값이 튀어나온다. 결국 기존 주제 3과 거의 다를 바가 없다. 
• 또한, $val[p^e-1]$이 $\pmod{p^e}$로 $1$ 또는 $-1$이라는 결과 역시 알려져 있다. 그러니 $k'$의 홀짝성만 봐도 ok.

이 알고리즘은 $\mathcal{O}(p^e)$ 급의 계산이 필요하다. 하지만 실제로는 $\mathcal{O}(pe)$ 급으로 계산이 가능하다. 매우 어렵다.

이 강화된 알고리즘을 공부하고 싶다면, HYEA님의 이 글과 원본이라고 볼 수 있는 min_25의 이 글을 보자.

이항계수

여전히 $k \le n/2$만 보는 아이디어는 유효하다.

Case 1. $k<p$이며 $k$가 작은 경우. 

이 경우, $1, 2, \cdots , k$ 각각의 $\pmod{p^e}$에 대한 modular inverse가 존재하므로 아이디어가 그대로 적용된다.

Case 2. 작은 $n, k$에 대하여 전부 전처리를 하고 싶은 경우.

이 경우, 파스칼의 삼각형을 이용하여 DP와 같은 방식으로 계산을 할 수 있다. 정수론이 필요없는 영역이다.

Case 3. $p$가 작은 경우, 즉 $\mathcal{O}(p^e)$ 또는 $\mathcal{O}(pe)$ 알고리즘이 동작할 수 있는 경우.

이 문제를 Lucas 정리로 풀려고 하면 (즉, Lucas 정리의 확장판을 찾으려 한다면) 상당히 고생할 것이다. 적어도 내 기억은 즐겁지 않았다.

하지만 주제 3을 이용한 풀이를 생각한다면, 기존 Case 3과 완전히 동일한 풀이가 적용된다. 주제 3을 여기에 수록한 이유다.

이제 $\pmod{p^e}$에 대한 이야기를 끝냈으니, 이론적으로 모든 내용을 다 했다. 하지만 몇 가지 더 이야기거리가 남았다.

• 중국인의 나머지 정리는 항상 강조하지만 수학에서나 정수론 알고리즘에서나 매우 중요한 결과고 쓸모도 많다.
• 하지만 이걸 제대로 써먹으려면 (즉, prime power로 문제를 쪼갠 후 다시 합치는 방식을 적용) 소인수분해를 해야 한다.
• 기본적으로 소인수분해는 매우 어려운 문제다. CP/PS 환경에서는 보통 커봐야 $2^{64}$인 수를 다루니 Pollard-Rho 알고리즘을 배우면 조금 낫다.
• 하지만 Pollard-Rho가 엄청 빠른 알고리즘인 것도 아니다. 즉, 일반적인 $\pmod{M}$에 대해서 사용할 수 있는 알고리즘을 더 다룰 필요가 있다.
• 당연한 것이지만 아래에서 다룰 "일반적인 알고리즘"은 위에서 다룬 $\pmod{p}$, $\pmod{p^e}$의 경우에서도 사용할 수 있다. 
• 아래에서 다룰 알고리즘들은 modulus가 갖는 소인수 $p$나 prime power $p^e$의 크기 등에 영향을 받지 않는다.

팩토리얼

주제 1. $\mathcal{O}(n)$ 알고리즘은 당연히 아직도 유효하다.

주제 2. 팩토리얼만 구하는 $\mathcal{O}(n)$ 알고리즘은 당연히 아직도 유효하다.

modular inverse를 전부 구하는 문제가 잘 정의되려면 $1, 2, \cdots , n$이 전부 $M$과 서로소여야 한다.

이 경우, 방법 1, 3이 잘 적용된다. (소수, 소수의 거듭제곱 여부를 사용한 적이 없다)

주제 3. 일반적인 modulus $M$에 대해서 논의하기 애매한 주제다. 이 때문에 이항계수를 계산하기가 어렵다.

이항계수

Case 1. $k$가 작고, $1, 2, \cdots , k$ 전부가 $M$과 서로소인 경우.

이 경우에는 기존 알고리즘을 그대로 적용할 수 있다. $\pmod{M}$에 대한 modular inverse를 취할 수 있기 때문이다.

Case 2. 파스칼의 삼각형을 그대로 적용할 수 있다. 이미 말했지만, 이 부분은 정수론이 아니라 DP에 가깝다.

Case 3. 하기 힘들다. 애초에 각 소수가 작다는 가정이 전제되는 만큼, 소인수분해가 필요하다.

Case 4. 특수한 경우를 하나 더 소개한다. 바로 $n, k$가 둘 다 적당히 작은 경우다.

이 경우, 앞서 소개한 에라토스테네스의 체를 이용한 로그 시간 소인수분해를 하여 문제를 해결할 수 있다.

즉, $1, 2, \cdots , n$ 전부를 소인수분해 할 수 있으니, $\displaystyle \binom{n}{k} = \frac{n \cdot (n-1) \cdots \cdot ( n-k+1)}{1  \cdot 2 \cdots \cdot k}$ 역시 그 분모/분자를 전부 소인수분해 할 수 있다.

이를 통해서 아예 $\displaystyle \binom{n}{k}$의 소인수분해를 얻을 수 있다. 이제 나눗셈에 대한 걱정이 아예 없으니, 단순 계산이 가능해진다. 끝.

다음 정리들은 정수론의 기초적이고 중요한 결과들이다. 증명은 찾아보면 쏟아져 나온다. 혹시 필요하면 댓글 :)

• 페르마의 소정리. 소수 $p$와 자연수 $a$가 $\text{gcd}(a, p)=1$을 (즉, $a$는 $p$의 배수가 아니다) 만족하면 $a^{p-1} \equiv 1 \pmod{p}$. 
• 오일러 정리. 앞서 오일러 파이 함수를 정의했다. 사실 $\phi(n)$은 $n$ 이하이며 $n$과 서로소인 자연수의 개수이다. 
• 오일러의 정리는, 자연수 $n, a$가 $\text{gcd}(a, n)=1$을 만족할 경우, $a^{\phi(n)} \equiv 1 \pmod{n}$임을 말한다. 페르마의 소정리는 특수 케이스.

잠깐 오일러 파이 함수에 대한 기본 사실들을 빠르게 짚고 넘어가자. 증명은 검색하면 충분히 나올듯. 필요하면 댓글 :)

• 앞서 언급했듯이, $n$의 서로 다른 소인수가 $p_1, p_2, \cdots , p_k$라면 $\phi(n) = n(1-1/p_1)(1-1/p_2) \cdots (1-1/p_k)$이다. 
• $\phi$는 multiplicative function이다. 이는 $n, m$이 서로소인 자연수일 때, $\phi(nm) = \phi(n)\phi(m)$임을 의미한다.
• $m|n$이면 $\phi(m)| \phi(n)$ 역시 성립한다. 이 사실은 이 단원에서 사용될 것이다.
• 이 외에도, $\sum_{d|n} \phi(d) = n$ 등 재밌는 성질이 많다. $\sum_{d|n}$이란, $n$의 약수 $d$에 대하여 더한다는 뜻.

그런데 생각보다 이 정리들을 CP/PS 환경에서 사용할 곳을 떠올리기는 쉽지 않다. 

•  보통 많이 쓰이는 방식은 $p$가 소수일 때, $\pmod{p}$에서 $a$의 ($a, p$는 서로소) modular inverse를 구하는 것이다.
•  즉, $a^{p-1} \equiv 1 \pmod{p}$에서 $a^{p-2} \cdot a \equiv 1 \pmod{p}$를 얻고, $a^{p-2}$가 modular inverse임을 얻는 것이다.
• 하지만 우리가 앞에서 살펴보았듯이, modular inverse를 구하는 것은 소수가 아닌 경우에도 쓸 수 있는 일반적이고 빠른 방법이 있다.
• 페르마의 소정리는 그렇다 쳐도, 오일러 정리는? modular inverse를 오일러 정리로 구하려면 $\phi(n)$ 값이 필요하다.
• $\phi(n)$을 구하려면 기본적으로 $n$의 소인수분해가 필요하다. $n$이 $10^9$-scale로 큰 경우에는 에라토스테네스의 체도 쓸 수 없으니, 진짜 소인수분해를 해야한다. 이때 시간복잡도는 현재 지식상으로는 $\mathcal{O}(\sqrt{n})$. 로그 시간을 보장하는 확장 유클리드와 비교된다. 그러면 어디서 쓸까?
• 물론, 오일러 파이 함수 자체는 CP/PS 환경에서 많이 쓰인다. 카운팅 문제에서 Burnside's Lemma 등을 활용할 때도 사용되고, 후에 다룰 내용에도 등장한다. 애초에 쓸모가 엄청 많은 함수다. 여기서 사용할 곳이 적다고 말하는/주장하는 것은 오일러 정리의 활용이다.
• (물론, 이렇게 말하지만 오일러 정리도 PS에서 아래에서 설명하는 예시 외의 문제에서 쓰이지 않는 것은 당연히 아니다)

보통 CP/PS에서 (그리고 사실 예전 KMO 1차에서) 이러한 정리들을 사용하는 환경은, "Power Tower"를 계산하기 위한 것이다. 대충 이런 문제다. 

모두의 편의를 위하여 $[a_1, a_2, \cdots , a_k] = a_1^{a_2^{ \cdots ^{a_k}}}$라고 정의하자. 우리의 목표는 $[a_1, a_2, \cdots , a_k] \pmod{n}$을 효율적으로 계산하는 것이다.

또한, 편의상 Power Tower $[a_1, a_2, \cdots, a_k]$의 길이를 등장하는 자연수의 개수 $k$라고 정의하겠다.

원리에 특별히 관심이 없고 알고리즘 및 시간복잡도만 필요한 분들은, 마지막 문단으로 바로 넘어가도 된다. 그래도 한 번 다 읽어보는 것을 추천. 

상대적으로 간단한 경우부터 보자. $\text{gcd}(a_1, n) = 1$을 가정한 상태에서 생각을 해보자. 오일러 정리에 의하여, 우리는 $a_1^{\phi(n)} \equiv 1 \pmod{n}$을 알고 있다. 

그러므로, $a_1$을 밑으로 할 때 지수에 $\phi(n)$에 대한 나머지를 취해도 결과값의 $n$에 대한 나머지가 변하지 않는다.

즉, $[a_1, \cdots , a_k] \equiv a_1^{[a_2, \cdots , a_k]} \equiv a_1^{[a_2, \cdots , a_k] \pmod{\phi(n)}} \pmod{n}$가 성립할 것이다. 이 식의 의미하는 바는 결국

• $a_1, \cdots , a_k, n$에 대한 문제를 $a_2, \cdots, a_k, \phi(n)$으로 간단하게 "축소", "환원"할 수 있다.

라는 것이다. 왜 "축소"일까? 이를 바라보는 것에는 두 가지 관점이 가능하다. 

• Power Tower의 "길이"가 $k$에서 $k-1$로 축소되었다고 볼 수 있다.
• $\phi(n) < n$은 ($n \neq 1$이면) 앞서 언급한 $\phi$의 성질/정의에서 자명하다. 그러니 modulus가 $n$에서 $\phi(n)$으로 축소되었다고 볼 수 있다.

어느 관점으로 보는 게 좋을까? 첫 번째 관점은 단순히 $\mathcal{O}(k)$의 횟수로 끝남을 의미할 뿐이다. 두 번째를 보자.

• 사실 1. $\phi(n)$은 $n \ge 3$에 대하여 모두 짝수이다. $\phi(n)$은 $1$ 이상 $n$ 이하 자연수 중 $n$과 서로소인 것의 개수라고 했었다. 그런데 $\text{gcd}(a, n)=1$이면 $\text{gcd}(n-a, n) = 1$도 성립한다. (갑자기 이해가 안된다면, 유클리드 알고리즘 부분을 다시 한 번 읽고 오는 것을 추천한다) 그러니 $a$와 $n-a$를 (합이 $n$인) 한 쌍으로 묶을 생각을 할 수 있다. 이렇게 묶이지 않는 경우는 $n$이 짝수일 때 $n/2$ 뿐이다. 그런데 $n \ge 3$이므로 $n$이 짝수이면 $\text{gcd}(n/2, n) = n/2 > 1$이 되고, 이 경우는 고려할 필요가 없음을 알 수 있다. 그러니 쌍을 전부 맺어줄 수 있고, 개수는 짝수.
• 사실 2. $n$이 짝수이면, $\phi(n) \le n/2$이다. $2$ 이상 $n$ 이하 모든 짝수들이 $n$과 공통 소인수 $2$를 가지므로, $n$과 서로소가 될 수 없다. 그러니 서로소인 것의 개수는 많아봐야 $n-n/2 = n/2$개이며, 증명 끝. 이 두 사실을 이용하여, 우리가 문제를 "축소"하는 방식이 매우 강력함을 보여보자.

우리의 두 번째 관점은 modulus $n$의 감소에 집중하고 있다. $n$이 작으면 문제가 확실히 간단해진다.

$n=2$라면, $a_1$의 홀짝성만 판단하여 답을 바로 결정할 수 있다. $n=1$이면 그냥 답이 $0$이다. 그러니, 목표는 결국

• $n \rightarrow \phi(n) \rightarrow \phi(\phi(n)) \rightarrow \cdots $를 거칠 때, 언제 $1$이 될 것인가?

가 된다. 이러한 과정이 많아봐야 $\mathcal{O}(\log n)$번 필요함을 보이자.

• $n \rightarrow \phi(n)$ 과정에서는 별 논의를 할 수 없지만, 적어도 $\phi(n) < n$임은 안다.
• 그 이후부터는 재밌어진다. 사실 1에 의하여, 현재 $\phi(n)$은 짝수거나 $\phi(1)=\phi(2)=1$이다.
• $\phi(n)=1$인 경우에는 그냥 끝이다. $\phi(n)$이 짝수라 가정하자. 사실 2에 의하여, $\phi(\phi(n)) \le \phi(n)/2$이다.
• 즉, 한 번 과정을 거치면 modulus가 무조건 반 이상 감소한다. 이제 많아봐야 로그 횟수의 과정이 필요함은 당연하다.

물론, 위 논의는 문제를 축소하는 과정에서 계속 밑과 modulus가 서로소라는 가정이 성립할 때 가능한 것이다. 

즉, 처음에는 $\text{gcd}(a_1, n)=1$이어야 하며, 그 뒤에는 $\text{gcd}(a_2, \phi(n))=1$, $\text{gcd}(a_3, \phi(\phi(n)))=1$ 등이 필요하다.

이 가정을 제거하고 일반적인 문제를 해결하려면, 많은 노력이 필요하다. 우선 빠르게 전처리를 하나 하자.

• 이제부터 $a_i \ge 2$가 각 $i$에 대하여 성립함을 가정한다. $1$은 몇 번 곱해도 $1$이므로, $[a_1, a_2, \cdots, a_k , 1, b_1, b_2, \cdots, b_t] = [a_1, a_2, \cdots, a_k]$가 항상 성립한다. 즉, $1$이 등장하면 그 $1$과 뒤에 있는 tower를 다 날려도 된다. 이렇게 쓸모없는 부분을 날리는 과정을 거치면, $a_i \ge 2$인 상태가 된다. 

이제 다시 문제를 축소할 준비를 하자. $[a_1, \cdots, a_k] \equiv a_1^{[a_2, \cdots,  a_k]} \pmod{n}$는 정의상 여전히 성립.

문제를 해결할 각을 천천히 재보자. 개인적인 생각으로, 핵심 아이디어는 다음과 같다.

• 아이디어 1. 중국인의 나머지 정리와 소인수분해 : 사실 앞 단원에서 해도 되는 이야기지만, 여기서 한다. 
• 소인수분해는 자연수 $n$을 적당한 서로 다른 소수 $p_1, \cdots , p_k$에 대하여 $p_1^{e_1} p_2^{e_2} \cdots  p_k^{e_k}$ 형태로 유일하게 쓸 수 있다는 결과이다.
• 중국인의 나머지 정리는 $x \pmod{p_1^{e_1}}, x \pmod{p_2^{e_2}}, \cdots , x \pmod{p_k^{e_k}}$를 알면 이를 합쳐서 $x \pmod{n}$을 얻을 수 있다는 결과이다.
• 물론 반대 방향도 마찬가지. 즉, $n$에 대한 나머지를 아는 것과 $p_1^{e_1}, p_2^{e_2}, \cdots , p_k^{e_k}$에 대한 나머지를 아는 것은 완벽하게 동치이다.
• 그러니 $n$을 소인수분해하고, $n$이 갖는 각 prime power $p^e$에 대하여 $a_1^{[a_2, \cdots, a_k]} \pmod{p^e}$를 구할 생각을 할 수 있다.
• 이 생각의 흐름은 순수하게 수학적으로 보나, 지금/나중에 다룰 알고리즘의 중간 과정으로 보나 매우 중요하다. 

• 아이디어 2. 소수에 대한 케이스 분류 : 이제 목표를 $a_1^{[a_2, \cdots , a_k]} \pmod{p^e}$를 계산하는 것으로 정했다.
• Case 1. $a_1$이 $p$의 배수가 아닌 경우. 이 경우에는 $\text{gcd}(a_1, p^e)=1$이므로 기존에 했던 것과 똑같이 할 수 있다. 
• 즉, $[a_1, a_2, \cdots, a_k] \equiv a_1^{[a_2, \cdots , a_k] \pmod{\phi(p^e)}} \pmod{p^e}$가 성립한다.
• Case 2. $a_1$이 $p$의 배수인 경우. 이 경우가 우리가 새롭게 다루어야 하는 경우라고 볼 수 있겠다.
• 이 경우, $[a_2, \cdots , a_k]$의 값이 중요하다. 당장 이 지수가 $e$ 이상이라면, $a_1^{[a_2, \cdots, a_k]}$는 $p^e$의 배수가 된다.
• 그러니, 우리가 고안해야 하는 것은 $[a_2, \cdots, a_k]$가 $e$ 이상인지 판별하고, 아니라면 실제 값이 얼마인지 구하는 알고리즘이다.
• 길이 $4$의 Power Tower만 해도 $[2, 2, 2, 2] = 2^{16}$이므로, 이는 (CP/PS 환경에서) 커봐야 $60$ 언저리일 $e$를 훌쩍 넘는다.
• 즉, 길이 $3$ 이하의 Power Tower에 대해서만 고려해도 된다. 이제부터는 정수론이 아닌 단순 알고리즘의 영역이니, 직접 고민해보는 것 추천.
• 어쨌든 이 과정에서 $a_1^{[a_2, \cdots , a_k]} \pmod{p^e}$를 구할 수 있다. 각 경우에서 해답을 얻었으니, 중국인의 나머지 정리로 합치면 된다.

위 알고리즘은 어쨌든 문제를 환원하는 알고리즘이고, 실제로도 잘 작동할 것이다. 하지만 좀 복잡하고 더러우니, 깔끔하게 바꿔보자.

• 보조정리. $m|n$이면, $a \equiv (a \pmod{n}) \pmod{m}$이 성립한다. 
• 증명. $a \pmod{n}$은 어쨌든 적당한 정수 $k$에 대하여 $a+kn$으로 쓸 수 있다. $(a+kn)-a=kn$은 $m$의 배수이므로 증명 끝.

다시 아이디어 2의 케이스 분류 중 첫 번째 케이스로 돌아가보자. 

• Case 1. $a_1$이 $p$의 배수가 아닌 경우.
• $p^e|n$이므로, $\phi(p^e)|\phi(n)$이다. 보조정리에서 $([a_2, \cdots , a_k] \pmod{ \phi(n)}) \equiv [a_2, \cdots , a_k] \pmod{\phi(p^e)}$다.
• 그러므로, 오일러 정리에서 $[a_1, \cdots, a_k] \equiv a_1^{[a_2, \cdots , a_k] \pmod{\phi(n)}} \pmod{p^e}$가 성립한다.
• 이는 문제를 $n \rightarrow \phi(p_1^{e_1}), \phi(p_2^{e_2}), \cdots , \phi(p_k^{e_k})$로 복잡하게 문제를 환원하지 않고, 기존처럼 $n \rightarrow \phi(n)$ 형태로 문제를 환원할 수 있음을 의미한다.
• 즉, $[a_2, \cdots , a_k] \pmod{\phi(n)}$만 구하면 별다른 추가적인 처리 없이 원하는 계산을 할 수 있다는 의미다. 

이제부터 $[a_2, \cdots , a_k]$의 크기에 대한 케이스를 나눈다. 

• $[a_2, \cdots , a_k]$가 $100$ 이하인 경우. 이 경우에는 애초에 Power Tower의 길이 $k-1$이 $3$ 이하임을 앞에서 설명하였다. 그러니 $[a_2, \cdots,  a_k]$를 직접 계산한 뒤, $a_1^{[a_2, \cdots , a_k]} \pmod{n}$을 계산하면 된다. 즉, 이 경우는 중국인의 나머지 정리, 소인수분해, 앞선 케이스 분류 등의 아이디어가 필요없다.
•  $[a_2, \cdots , a_k]$가 $100$ 초과인 경우. 이 경우에는 $p|a_1$인 소수 $p$들에 대하여, $p^e$가 $n$의 prime power인 경우 ($e \le 100$일테니) $[a_1, a_2, \cdots , a_k] \equiv a_1^{[a_2, \cdots , a_k]} \equiv 0 \pmod{p^e}$가 될 것이다. 물론, $a_1$이 $p$의 배수가 아닌 경우는 위에서 다루었다.

이제 다음을 보일 준비가 완료되었다.

• Claim. $[a_2, \cdots , a_k] > 100$인 경우, $[a_1, \cdots, a_k] \equiv a_1^{[a_2, \cdots , a_k] \pmod{\phi(n)} + 100 \cdot \phi(n)} \pmod{n}$이다.
• 중국인의 나머지 정리를 생각하면, 결국 $n$의 소인수분해에 등장하는 각 prime power $p^e$에 대해 양변이 같은 나머지를 가짐을 보이면 된다.
• Case 1. $a_1$이 $p$의 배수가 아닌 경우. $\phi(p^e)|\phi(n)$이 성립하므로, 오일러의 정리에서 (즉, $a_1^{\phi(p^e)} \equiv 1 \pmod{p^e}$)
• $[a_1, \cdots , a_k] \equiv a_1^{[a_2, \cdots , a_k] \pmod{\phi(n)}} \equiv a_1^{[a_2, \cdots , a_k] \pmod{\phi(n)} + 100 \cdot \phi(n)} \pmod{p^e}$다. 증명 끝.
• Case 2. $a_1$이 $p$의 배수인 경우. $e \le 100$이고, $[a_2, \cdots , a_k] \ge 100$이고, $100 \cdot \phi(n) \ge 100$이다. 
• 그러니, $[a_1, \cdots , a_k] \equiv a_1^{[a_2, \cdots, a_k]} \equiv  0 \equiv a_1^{[a_2, \cdots , a_k] \pmod{\phi(n)} + 100 \cdot \phi(n)} \pmod{p^e}$다. 증명 끝.

위의 $100$이라는 값은 CP/PS 환경을 가정한 것이다. 실제로는 modulus가 $n$일 때, $\log_2 n$ 초과인 값을 아무거나 잡으면 된다. (너무 크게 잡으면 고생한다)

CP/PS 환경에서 웬만하면 $2^{100}$을 넘는 modulus는 볼 수 없을 것이므로, 그냥 $100$이라고 잡았다. 이렇게 하는 것이 짜기도 편하다.

이제 알고리즘을 정리해보자. 

• 입력은 $a_1, \cdots , a_k$와 $n$. $[a_1, \cdots , a_k] \pmod{n}$을 푸는 것이 목표다. 단, $a_i \ge 2$를 가정한다.
• 만약 $k \le 4$라면, $[a_2, \cdots , a_k]$가 $100$ 이하인지 판별하는 과정을 거친다.
• 만약 $100$ 이하라면, $[a_2, \cdots , a_k]$를 직접 계산하고 $a_1^{[a_2, \cdots , a_k]} \pmod{n}$을 반환한다.
• 그렇지 않다면, $\phi(n)$을 계산해준 뒤, $[a_2, \cdots , a_k] \pmod{\phi(n)}$을 계산한다. (즉, $a_2, \cdots , a_k$와 $\phi(n)$에 대한 문제를 푼다.)
• 그 후 $a_1^{[a_2, \cdots , a_k] \pmod{\phi(n)} + 100 \cdot \phi(n)} \pmod{n}$을 반환한다. 끝.

$a_i$가 $1$이 되어도 빠르게 잘 작동하는 알고리즘을 구축하는 것은 여러분의 몫이다. (정수론적 내용/테크닉이 추가되지 않는다) 

이 알고리즘은 맨 처음 여러 가정을 추가한 경우의 알고리즘처럼, modulus를 $n$에서 $\phi(n)$으로 축소한다.

이 과정이 로그 번 반복된 후에 modulus가 $1$이 된다는 사실은 그 때와 달라지지 않는다. 이제 마지막으로 복잡도 분석을 하자. 

각 단계에서 가장 핵심적인 비용은 역시 $\phi(n)$의 계산이다. 에라토스테네스의 체 등을 사용하지 않았다고 가정하면, $\phi(n)$의 계산은 $n$의 소인수분해가 필요하며 이는 $\mathcal{O}(\sqrt{n})$에 할 수 있다. 최대 $\mathcal{O}(\log n)$ 번의 단계가 필요하므로, 전체 시간복잡도는 $\mathcal{O}(\sqrt{n} \log n)$이다. 끝.

관련 코드는 github에서 찾아볼 수 있다.

이 주제도 이미 좋은 자료가 있습니다. https://casterian.net/archives/609

사실 3단원까지만 잘 이해해도 CP/PS 하는데는 문제가 없습니다. 레드까지는 문제 없을듯 ㅋㅋ 

소개가 늦었지만, 앞서 다룬 $ax \equiv b \pmod{c}$와 같이 mod 연산을 다루는 방정식을 합동식이라 한다.

중국인의 나머지 정리는 자연수 $a_1, n_1, a_2, n_2, \cdots , a_k , n_k$가 주어졌을 때, 연립합동식 $$ x \equiv a_i \pmod{n_i} \quad i = 1, 2, \cdots k $$를 푸는 알고리즘이다. 우리의 첫 목표는 $k=2$인 경우를 푸는 것이다. 나중에 가면 알겠지만, $k=2$만 풀면 끝이다.

$x \equiv a_1 \pmod{n_1}$, $x \equiv a_2 \pmod{n_2}$가 성립한다고 가정하자. 정의상, $x = n_1 y + a_1$인 정수 $y$가 있다.

그러니 우리는 $n_1y + a_1 \equiv a_2 \pmod{n_2}$, 또는 $n_1 y \equiv a_2 - a_1 \pmod{n_2}$를 풀면 된다.

그런데 이러한 형태의 합동식은 이미 푸는 방법을 알고 있다. 바로 앞에서 배웠기 때문.

• $g = \text{gcd}(n_1, n_2)$라 하자. 만약 $a_2-a_1$이 $g$의 배수가 아니라면 해는 존재하지 않는다.
• 만약 $a_1 \equiv a_2 \pmod{g}$라면, 앞선 결과에서 우리는 $y \equiv t \pmod{n_2/g}$ 형태의 결과를 얻는다.
• 이는 다른 말로, 적당한 자연수 $m$이 있어서 $y = (n_2/g) m  + t$이라는 것이다.
• $x = n_1y + a_1$에 그대로 대입하면, $x = (n_1n_2/g) m + (n_1t+ a_1)$이다. 즉, $x \equiv (n_1t+ a_1) \pmod{n_1n_2/g}$.
• 여기서 $n_1n_2 / g$가 $\text{lcm}(n_1, n_2)$임을 확인하자. $\text{gcd}(a, b) \times \text{lcm}(a, b) = a \times b$라는 결과 때문.

결국 $k=2$일 때의 결론은, 연립합동식이 해가 없거나 아니면 $x \equiv \text{something} \pmod{\text{lcm}(n_1, n_2)}$라는 결론을 얻는다는 것이다.

이는 연립합동식의 해가 존재하는 경우, 두 합동식을 하나의 합동식으로 합칠 수 있다는 의미가 된다.

• Example. $x \equiv 17 \pmod{42}$, $x \equiv 23 \pmod{60}$을 풀어보자.
• $x \equiv 17 \pmod{42}$에서, $x = 42m + 17$이라고 쓸 수 있다. 
• 이제 $42m + 17 \equiv 23 \pmod{60}$을 풀자. 이는 $42m \equiv 6 \pmod{60}$을 푸는 것과 같다.
• 앞에서 배운 방식을 생각하면, $42m+60n = 6$을 풀면 된다. $\text{gcd}(42, 60) = 6$이므로, $6$으로 식을 나누자.
• $7m+10n = 1$을 풀자. 이는 확장 유클리드 알고리즘으로 가능하며, $m=3$, $n=-2$를 얻는다.
• 여기서 우리는 $m \equiv 3 \pmod{10}$을 얻는다. $m=10t+3$이라 쓰고, $x=42m+17$에 대입하자.
• 그러면 $x = 42(10t+3)+17 = 420t + 143$을 얻고, 이는 $x \equiv 143 \pmod{420}$을 의미한다. 

이제 $k \ge 3$인 경우도 간단하게 해결할 수 있다. 두 합동식을 하나의 합동식으로 합치는 것을 반복하기만 하면 된다.

만약 중간 과정에서 "해가 없다"라는 결론이 나온다면, 전체 연립합동식 역시 해가 없다는 것을 바로 파악할 수 있다.

중국인의 나머지 정리, Chinese Remainder Theorem을 줄여서 CRT로 많이 부른다. 

중국인의 나머지 정리의 심화/응용 내용으로 Garner's Algorithm이 있다. CRT로 큰 수를 다루는 방법이다. 관심이 있다면 링크 참조.

관련 코드는 github에서 찾아볼 수 있다.

이 주제는 이미 좋은 자료가 있습니다. https://casterian.net/archives/601

이 글은 유클리드 알고리즘은 이미 어느 정도 숙지한 독자를 위해 작성하였다. 유클리드 알고리즘 관련 자료 역시 많다.

글의 핵심 목적은 유클리드 알고리즘의 핵심 아이디어 2개를 짚고, 확장 유클리드 알고리즘도 사실 별반 다를 거 없다는 것을 주장하는 것이다.

앞쪽 내용인만큼 제대로 이해해야 뒤 내용을 이해할 수 있으니, 시간을 많이 투자해보자. 

유클리드 알고리즘은 자연수 $a, b$가 주어졌을 때, $\text{gcd}(a, b)$를 구하는 방법이다. 결국 핵심 아이디어는

• 아이디어 1. 몫, 나머지를 구하여 $a = bq + r$로 쓰자. 이때 $\text{gcd}(a, b) = \text{gcd}(b, r)$이다. 
• Why? $g|a$, $g|b$가 성립한다고 가정하면 $g|(a-bq)=r$이므로, $g|b$, $g|r$.
• 반대로 $g|b$, $g|r$이 성립한다고 가정하면 $g|(bq+r)=a$이므로, $g|a$, $g|b$. 
• 즉, $g$가 $a, b$의 공약수인 것과 $g$가 $b, r$의 공약수인 것은 완전히 동치이다. 그러니 최대공약수 역시 동일.

• 아이디어 2. 위 사실은 문제를 $(a, b)$에서 $(b, r)$로 "축소"하는 방법을 제시해준다. 즉, $\text{gcd}(a, b)$를 구하는 문제를 $\text{gcd}(b, r)$을 구하는 문제로 축소.
• 얼마나 줄어들까? $a \ge b$를 가정하자. 그러면 $q \ge 1$이며 $2r \le (q+1)r = rq + r \le bq + r = a$이므로 $2r \le a$.
• 그러므로, $br \le (1/2) ab$임을 얻는다. 즉, 문제를 정의하는 두 인자의 곱이 적어도 반 이상은 감소한다.

두 아이디어를 합치면, 로그 시간을 거쳐서 최대공약수를 얻는 방법이 완성된다.

$(a, b)$를 $(b, r)$로 변환하는 과정을 거치면, 로그 시간 이내에 $(g, 0)$ 형태의 순서쌍을 얻는 것이 보장되기 때문이다.

이 시점에서 $\text{gcd}(a, b) = g$를 얻고, 알고리즘을 종료시키면 된다. 여기서부터는 구현의 문제다.

• Example. $\text{gcd}(576, 204)$를 구한다고 가정하자. 다음 과정을 거친다.
• $576 = 204 \times 2 + 168$이므로, $(576, 204)$에서 $(204, 168)$로 문제를 축소
• $204 = 168 \times 1 + 36$이므로, $(204, 168)$에서 $(168, 36)$으로 문제를 축소
• $168 = 36 \times 4 + 24$이므로, $(168, 36)$에서 $(36, 24)$로 문제를 축소
• $36 = 24 \times 1 + 12$이므로, $(36, 24)$에서 $(24, 12)$로 문제를 축소
• $24 = 12 \times 2 + 0$이므로, $(24, 12)$에서 $(12, 0)$으로 문제를 축소. 답은 이제 $12$. 

아이디어 1, 2는 CP/PS의 많은 정수론 문제에서 사용되는 아이디어다. 이를 유클리드 알고리즘의 활용에서 다룬다.

확장 유클리드 알고리즘은 자연수 $a, n$이 주어졌고 $\text{gcd}(a, n) = 1$일 때, $ax \equiv 1 \pmod{n}$인 $x$를 찾는 알고리즘이다.

(엄밀하게 말하자면, 자연수 $a, b$에 대하여 $ax + by = \text{gcd}(a, b)$인 $x, y$를 찾는 알고리즘이다. 사실상 똑같은 말)

이 문제를 바라보는 가장 정석적인 방법은 변수를 하나 추가하는 것이다.

• $ax \equiv 1 \pmod{n}$이라는 것은, 정수 $y$가 있어 $ax + ny = 1$이라는 것과 동치이다.
• 그러니 사실 우리는 $ax + ny = 1$을 만족하는 정수 순서쌍 $(x, y)$를 찾으면 된다. 
• 여기서 우리가 $\text{gcd}(a, n) = 1$이라는 조건을 추가한 이유도 알 수 있다. 위 식의 좌변이 무조건 $\text{gcd}(a, n)$의 배수이기 때문.

목표는 이제 $\text{gcd}(a, b) = 1$을 가정하고 $ax + by = 1$을 푸는 것이다. 문제를 바라보는 새로운 시각을 얻었으니, 위 아이디어 1, 2를 적용해보자.

• 아이디어 1. 몫, 나머지를 구하여 $a = bq+r$로 쓰자. $a, b$에 대한 문제를 $b, r$에 대한 문제로 바꿀 수 있는가?
• 가능하다. $ax + by = 1$은 $(bq+r)x + by = 1$과 동치이며, 정리하면 $b(qx+y) + rx = 1$로 쓸 수 있다.
• 그러니, $bx'+ry'=1$인 정수 순서쌍 $(x', y')$을 구하면, $qx + y = x'$, $x = y'$을 풀어 $ax + by = 1$인 $(x, y)$를 얻는다.
• 위 식을 다시 정리하면, 결국 $x = y'$, $y = x' - qy'$이므로, 이 순서쌍 $(x, y)$는 정수 순서쌍이다. 
• 이는 $(b, r)$에 대한 문제를 풀면, 빠르게 $(a, b)$에 대한 문제를 풀 수 있음을 의미한다.

• 아이디어 2. 역시 아이디어 1은 문제 축소의 의미를 담는다. 기존 유클리드 알고리즘에서 한 논의가 그대로 성립한다.
• 결과적으로, 위 과정을 로그번 거치면 결국 $(1, 0)$에 대한 문제를 해결하는 것으로 문제가 환원된다. 이는 $x=1, y=0$이라는 해가 있다.

두 아이디어를 합치면, 로그 시간을 거쳐서 조건을 만족하는 $x$를 찾을 수 있다. 이를 modular inverse라고 부르며, 이는 유일하게 존재한다.

알고리즘을 설명하는 과정에서 $ax + by = 1$을 만족하는 정수 순서쌍 $(x, y)$가 존재함까지 증명했음을 참고하자. 

나아가면, $ax + by = \text{gcd}(a, b)$인 $x, y$가 존재함을 얻는다. 이를 Bezout's Identity라고 한다.

• Example. $3x \equiv 1 \pmod{7}$을 만족하는 $x$를 구해보자.
• $3x + 7y = 1$인 정수 순서쌍 $(x, y)$를 찾으면 충분하다. 
• 이는 $3(x+2y) + y = 1$과 같다. 그러니 이제 $3a + b = 1$인 정수 순서쌍 $(a, b)$를 찾자.
• 이는 $1(3a+b) + 0a = 1$과 같다. 이제 $3a + b = 1$, $a = 0$을 풀어 $a = 0, b = 1$을 얻는다.
• 다시 위로 돌아가서, $3(x+2y) + y = 1$을 $3a + b = 1$로 변환한 것을 생각해보자.
• 우리의 목표는 이제 $x + 2y = 0$, $y=1$. 이를 풀면 $x= -2$와 $y = 1$을 얻는다.
• 실제로 $3(-2) \equiv 1 \pmod{7}$임을 쉽게 확인할 수 있다.

• $ax+by=1$이면, 각 정수 $t$에 대하여 $a(x+bt) + b(y-at) = 1$이 성립한다. 즉 $(x+bt, y-at)$ 역시 해.
• $t$를 잘 잡아서, $x+bt$가 $0$ 이상 $b$ 미만이 되도록 할 수 있다. overflow 방지를 위해 필요한 작업이다.
• 즉, 해를 재귀적으로 계산하는 과정에서, 해가 지나치게 크거나 작아지지 않도록 하는 과정이 필수적이다. 
• 위 Example을 보면 느낌이 오겠지만, 이 알고리즘은 재귀적으로 작성하는 것이 정신건강에 좋다.

더 나아가면, 아무런 가정 없이도 $ax \equiv b \pmod{n}$과 같은 식을 해결할 수 있다.

• 마찬가지 아이디어로, $ax + ny = b$라는 식을 해결하는 것으로 문제를 바꾸자.
• $g = \text{gcd}(a, n)$을 유클리드 알고리즘으로 구하자. 위 식의 좌변 $ax+ny$는 무조건 $g$의 배수임을 알 수 있다.
• 그러니, $b$가 $g$의 배수가 아니라면 해가 없음을 바로 알 수 있다. 
• 이제 $g|b$를 가정하고, 식을 $(a/g)x + (n/g)y = (b/g)$로 써보자. 이제 $\text{gcd}(a/g, n/g) = 1$이다.
• 확장 유클리드 알고리즘으로, $(a/g)x' + (n/g)y' = 1$인 정수 순서쌍 $(x', y')$을 찾을 수 있다.
• 이제 $x = (b/g)x'$, $y = (b/g) y'$을 선택하면 $ax + ny = b$임을 알 수 있다.
• 또한, $(x, y)$가 해라면 $(x + n/g \cdot t, y - a/g \cdot t)$ 역시 각 정수 $t$에 대하여 해임을 알 수 있다.
• 즉, $x \equiv \text{something} \pmod{n/g}$ 형태의 해를 얻는다. 자세한 디테일은 github 참고.

• Example. $12x \equiv 18 \pmod{54}$를 해결해보자.
• $12x + 54y = 18$을 해결하면 충분하다. $\text{gcd}(12, 54) = 6$임을 유클리드 알고리즘으로 계산한다.
• $18$은 $6$의 배수이니, 해가 존재한다. 식을 $6$으로 나누어주자. $2x + 9y = 3$을 얻을 수 있다.
• 이제 $\text{gcd}(2, 9) = 1$이니, $2x' + 9y' = 1$인 $x', y'$를 확장 유클리드 알고리즘으로 찾을 수 있다.
• $x'=5$, $y'=-1$을 찾았다고 하자. 이제 이 값을 $3$배 해주면 $x = 15$, $y = -3$을 얻는다.
• 이제 $x \equiv 15 \pmod{9}$, 또는 $x \equiv 6 \pmod{9}$가 답이 된다. 끝!  

사실 위 논의에서 일부 엄밀한 수학적 논의를 생략했다. 내용을 추가할 겸 간략하게 설명한다.

• 우리는 이미 확장 유클리드 알고리즘으로 $\text{gcd}(x, y) = 1$이면 $ax+by=1$인 정수 $a, b$가 존재함을 안다. (엄밀하게 증명한 것은 아니나, 나름 합리적으로 유도했다. "수학자들의 증명"을 원하면 당연히 정수론 책 참고.) 이를 통해 $\text{gcd}(a, b)=1$, $a|c$, $b|c$가 성립하면 $ab|c$가 성립함도 증명할 수 있다. 또한, $a|bc$이고 $\text{gcd}(a,b)=1$이면 $a|c$가 성립함도 증명할 수 있다. 특히, $g=\text{gcd}(a, b)$라면 $\text{gcd}(a/g, b/g)=1$이니 $a|bc$에서 $(a/g)|(b/g)c$를 얻고, 즉 $(a/g)|c$임 역시 얻을 수 있다. 즉, $a/\text{gcd}(a,b)$가 $c$의 약수가 되어야 한다. 이 "증명할 수 있다"라는 부분을 자세히 알고 싶다면 정수론 책 참고 or 댓글로 질문. 이 정도 내용은 알아야 정수론 문제 해결에 지장이 없을 것이라고 생각하여 추가한다.

• modular inverse는 유일한가? $x, y$가 모두 $\pmod{n}$에 대한 $a$의 modular inverse라고 하자. 그러면 가정 상 $\text{gcd}(a, n)=1$이고 $n|(ax-1)$, $n|(ay-1)$이니 두 식을 빼면 $n|a(x-y)$를 얻는다. 한편, $a|bc$이고 $\text{gcd}(a, b)=1$이면 $a|c$임이 알려져 있다. 그러므로 $n|(x-y)$를 얻고 $x \equiv y \pmod{n}$이 성립한다. 즉, modular inverse는 ($n$으로 나눈 나머지를 보았을 때) 유일함을 얻는다. 증명 끝.
• $ax \equiv b \pmod{n}$의 해는 정말 저게 다인가? $ax+ny=b$를 생각할 때, $(x, y)$가 해라면 $(x + n/g \cdot t, y - a/g \cdot t)$가 해임은 직접 대입을 통해서 자명하게 얻을 수 있다. 이것들만이 해라는 것 역시 보여야 한다. 만약 $ax+ny = ax'+ny'=b$라면, $a(x-x')=n(y'-y)$가 된다. 여기서 $x-x'$이 $n/g$의 배수, $y-y'$이 $a/g$의 배수임을 증명할 수 있다. 즉, 해가 되려면 $(x + n/g \cdot t, y - a/g \cdot t)$ 형태를 가져야 한다. 증명 끝.