rkm0959

쓰기 시작한지는 꽤 되었는데

• 회사 + 연구 + CTF하면서 시간이 잘 안나서 + 게을러져서 (...) 진도가 안빠지고
• 이대로 두었다가는 지금까지 쓴 부분도 공개가 너무 늦어질 것 같고
• 사실 지금까지 쓴 부분만 잘 이해해도 문제 해결에 지장이 없어서

미리 올립니다. 솔직히 머리 속에는 어떻게 글을 써야하는지 아는데 실제로 글을 쓰려면 막막해서 슬프네요 :(

그래도 이거 올리고 사람들이 좀 읽어주면 이 글을 다시 쓸 motivation을 얻을 것 같기도 하고 잘 모르겠습니다

LaTeX 형식은 Evan Chen의 Napkin에서 가져왔습니다. 참고하세요.

11단원의 내용으로 내가 알고 있는 PS/CP 정수론 분야 지식은 (특히, 문제풀이에 직접적인 내용은) 다 정리한 것 같다.

처음에는 그래도 친절하게 설명하려고 노력했는데, 후반부에 가니까 슬슬 다른 일의 우선순위가 높아져서 급하게 작성한 것 같아 아쉽다.

• 초반 부분 : 이미 자료가 있으니까 간략하게 써야지 ㅋㅋㅋ
• 후반 부분 : 내가 다른 일이 있으니까 간략하게 써야지 ㅋㅋㅋ

가 된 것 같아서 조금 그렇긴 한데 ㅋㅋ; 그래도 시간투자를 충분히 하면 이해할 수 있을 정도로는 쓴 것 같다. 나중에 코드, 연습문제도 올릴거고..

일단 이 내용을 다 이해했다면 (또는 관심있는 부분까지 다 이해했다면) 할 수 있는 것은 대강

• 백준 문제풀기. 문제집 제작을 내가 하면 더 제대로 할 수 있을 것이고, 아니어도 solved.ac나 이 블로그 수학 PS일지를 보면서 문제를 고를 수 있다.
• 진짜 고수들이 쓴 자료 읽기. 여러 의미를 담고 있는 말인데, 정리하자면
• 1. 정수론 교과서 읽기. 이유를 설명하지 않고 넘긴게 많으니, 이를 이해하고 싶다면 책을 사서 읽는 것이 제일 빠르다.
• 2. PS/CP 정수론 진짜 고수가 쓴 자료 읽기. 예로, kipa00님의 NTA, HYEA님의 소멤 글, min_25의 블로그, adamant의 블로그 등
• adamant의 블로그는 PS 수학 전반에 대한 다양한 내용이 있다. 특히 FFT 및 다항식 처리 관련 자료가 어마어마하다.
• HYEA님의 글과 min_25의 블로그에는 더 난이도 높은 알고리즘들이 많이 나와있다. 또한, HYEA님의 글 중 min_25의 글 번역본도 있다.
• kipa00님의 NTA는 수학적 깊이가 상당하고, computational number theory의 여러 부분을 다룬다. 수학하는 사람이면 읽는 게 좋을듯.
• 혹시 NTA를 다 읽고 컨텐츠가 부족하다고 느낀다면 (ㅋㅋ;) Victor Shoup의 computational number theory 책을 읽자.
• 3. Project Euler 풀기. 여기에는 고수들이 가득하고, thread 등에서 풀이나 아이디어 공유가 많다. 
• 특히, 10단원, 11단원의 내용을 확장/강화하는 토론들이 많이 이루어졌다. 대표적인 예시로 baihacker의 블로그를 참고하자. 
• 개인적으로 중국/일본에 비해 한국에서 Project Euler를 푸는 사람이 적다는 것에 많은 아쉬움을 느끼고 있다.
• 4. PS/CP scope 벗어나기. 이제 $2^{64}$ 미만의 수만 다루고, 짧은 시간 안에 답을 내야한다는 가정도 버리자.
• 아예 Computational Number Theory를 공부할 수도 있고, (NTA 읽는 것도 방법) 그 대표적인 활용처인 암호학을 공부하는 것도 좋다.
• 개인적으로는 PS/CP 정수론 및 PS/CP에서 겪은 문제풀이 경험이 CTF 암호학으로 넘어가는 것에 정말정말 많은 도움을 주었다.

시간이 없어서 글을 제대로 못 쓰기는 했는데 그래도 질문 받을 시간은 있으니 필요할 때 댓글 달아주세요 :)

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

이 글에서는 다음 내용을 다룬다.

• multiplicative function의 정의, 성질
• Dirichlet Convolution의 정의와 성질
• xudyh's sieve, min_25 sieve를 이용한 multiplicative function의 prefix sum

이 내용을 다루는 튜토리얼 중 가장 접근 가능한 것은 역시 이 Codeforces 블로그. 나도 여기서 배웠다.

multiplicative function

• 기존 정의를 생각하면, multiplicative function이란 $\text{gcd}(a, b)=1$일 때 $f(ab)=f(a)f(b)$인 함수 $f$를 말한다.
• 특히, 우리가 다룰 함수들은 $f(1)=1$을 만족한다. ($f=0$ 같은 특수한 경우를 고려하지 않겠다는 것이다)
• 이제 $n$을 소인수분해하고, $n=p_1^{e_1} \cdots p_k^{e_k}$라고 쓰자. 그러면 $f(n) = f(p_1^{e_1}) \cdots f(p_k^{e_k})$가 성립한다.
• 이는 prime power에 대한 $f$ 값들만 결정하면 모든 자연수에 대한 $f$ 값이 자동 결정이라는 뜻이다.
• 이 방법으로 생각하는 것이 아마 multiplicative function을 이해하는 가장 쉬운 방법일 것이다. 예시를 보자.

• 오일러 파이 함수 : $f(p^k) = p^{k-1}(p-1)$으로 정의된 multiplicative function
• mobius 함수 : $f(p)=-1$, $f(p^k)=  0$ for $k \ge 2$로 정의된 multiplicative function
• 약수의 개수 함수 : $f(p^k) = k+1$로 정의된 multiplicative function
• 약수의 합 함수 : $f(p^k) = 1 + p + \cdots + p^k$로 정의된 multiplicative function
• identity 함수 : $f(p^k) = p^k$로 정의된 multiplicative function
• $[n=1]$ 함수 : $f(p^k) = 0$으로 정의된 multiplicative function
• $f=1$ 함수 : $f(p^k) = 1$으로 정의된 multiplicative function

Dirichlet Convolution

• Dirichlet Convolution -> Dirichlet Series -> Analytic Number Theory의 연결 관계로 중요한 주제인 것으로 안다.
• 하지만 PS/CP에서는 Dirichlet Convolution 그 자체가 중요한 경우는 거의 못 봤고, 아래에서 다룰 xudyh's sieve에서 등장한다.
• 하지만, 실제로 Dirichlet Series와 같은 방식으로 PS/CP 문제풀이를 하는 경우를 Project Euler에서 봤다. 이건 나도 잘 모른다 ㅋㅋ

• 정의. 두 함수 $f, g$에 대하여, 그 Dirichlet Convolution $(f * g)$를 $(f * g)(n) = \sum_{d|n} f(d)g(n/d)$로 정의한다.
• 사실. $f, g$가 multiplicative 하다면, 그 Dirichlet Convolution $(f * g)$ 역시 multiplicative 하다.
• 증명. 또 $(f * g)(p^k) = f(1)g(p^k) + f(p)g(p^{k-1}) + \cdots + f(p^k)g(1)$로 두고, 위에서 보인 예시처럼 하면 된다.

• $f(n) = n$, $g(n)=1$이면 $(f * g)(n)= \sigma(n) = \sum_{d|n} d$ - 약수의 합
• $f(n)=1$, $g(n)=1$이면 $(f*g)(n) = \tau(n) = \sum_{d|n} 1$ - 약수의 개수
• $f(n) = \mu(n)$, $g(n)=1$이면 $(f*g)(n) = \sum_{d|n} \mu(d) = [n=1]$ - mobius function
• $f(n) = \phi(n)$, $g(n)=1$이면 $(f * g)(n) = \sum_{d|n} \phi(d) = n$
• $f(n) = \mu(n)$, $g(n) = n$이면 $(f * g)(n) = \sum_{d|n} \mu(d) (n/d) = \phi(n)$
• 특히 - mobius inversion을 제대로 표현하는 방법은 역시 $h = (f * 1)$이면 $f = (h * \mu)$라는 것이다.
• $f(n) = n^k \phi(n)$, $g(n) = n^k$면 $(f * g)(n) = n^{2k+1}$ 등.

• multiplicative function $f$가 있다. 우리의 목표는 $\sum_{i=1}^n f(i)$를 $\mathcal{O}(n^{2/3})$에 계산하는 것이다.
• 적당한 multiplicative function $g$가 있어, $g$의 prefix sum과 $(f * g)$의 prefix sum이 계산하기 쉽다고 가정하자.
• 함수 $f$에 대하여, $s_f(n) = \sum_{i=1}^n f(i)$라고 정의한다. 이제 준비 끝. 열심히 식을 조작해보자.

우선 Dirichlet Convolution의 정의를 사용하여 $$s_{f * g}(n) = \sum_{i=1}^n (f * g)(i) = \sum_{i=1}^n \sum_{d|i} f(i/d)g(d)$$ 를 얻는다. 이제 $d$에 대하여 식을 정리하면, 가능한 $i/d$의 값이 $1$부터 $\lfloor n/d \rfloor$이므로 $$s_{f * g}(n) = \sum_{d=1}^n g(d) \sum_{k=1}^{\lfloor n/d \rfloor} f(k) = \sum_{d=1}^n g(d) s_f(\lfloor n/d \rfloor)$$ 이다. $g(1)=1$이니, 이제 이 식을 $s_f(n)$에 대하여 정리하면 결과적으로 $$s_f(n) =  s_{f * g}(n) - \sum_{d=2}^n g(d) s_f(\lfloor n/d \rfloor)$$

• 계속해서 등장하지만, $\lfloor n/d \rfloor$ 형태를 갖는 서로 다른 자연수의 개수는 $\mathcal{O}(\sqrt{n})$개이다.
• 또한, $\lfloor n/d \rfloor$의 특정 값에 대응되는 $d$의 값은 하나의 구간을 이룬다.
• 특히, $g$의 prefix sum이 쉽게 계산될 수 있다고 가정했으므로, 그 구간에 대한 $g$ 값의 합 역시 쉽게 구할 수 있다.
• 여기서는 당연히 $\sum_{i=l}^r g(i) = s_g(r) - s_g(l-1)$이라는 식을 사용한다. 
• 그러니, 우리가 실제로 계산해야 하는 것은 $s_f(\lfloor n/d \rfloor)$ 값이다.

첫 방식은 $\lfloor n/d \rfloor$ 형태의 값들에 대해 $s_f$를 정직하게 DP로 계산하는 것이다.

$s_f(x)$를 계산하기 위해 드는 시간복잡도가 $\mathcal{O}(\sqrt{x})$이니, 우리에게 필요한 시간은 총 $$\sum_{i=1}^{\sqrt{n}} \sqrt{i} + \sum_{i=1}^{\sqrt{n}} \sqrt{n/i} = \mathcal{O}(n^{3/4})$$ 임을 알 수 있다. 이제 지수를 또 $3/4$에서 $2/3$으로 줄여보자. 이를 위해서는 $f$가 multiplicative임을 이용한다.

• 즉, 에라토스테네스의 체나 Linear Sieve를 사용하여, 작은 $i$들에 대한 $f(i)$를 전부 전처리하자.
• 그러면 자연스럽게, 부분합을 계산하여 작은 $i$들에 대한 $s_f(i)$ 역시 전부 전처리할 수 있다.
• 이번에는 $n^{2/3}$ 이하의 $i$에 대한 $s_f(i)$를 전처리하고, $s_f(n)$을 구하기 위한 점화식을 top-down으로 적용하자.
• 이미 전처리된 값은 바로 반환하고, 메모이제이션을 이용하여 중복 계산을 막는다. 그러면 시간복잡도는 $$\mathcal{O}(n^{2/3}) + \sum_{i=1}^{n^{1/3}} \sqrt{n/i} = \mathcal{O}(n^{2/3})$$ 이 된다. 이 방법이 xudyh's sieve이다. 이제 문제는 $s_f$를 구하고 싶을 때, $g$를 어떻게 잡냐는 것이다.

일단 $g$를 잡는 가장 좋은 방법은 위에서 선보인 예시들을 다 참고하는 것이다. 구글 검색하면 예시는 더 나온다.

특히, $1$, $n$, $n^k$ 등 $n$에 대한 다항식으로 나타나는 함수들은 그 prefix sum을 구하기가 상대적으로 쉽다는 점을 이용한다.

어려운 문제들을 풀기 시작하면, $g$를 잡기 위해서 전략적으로 각 prime power들에 대해 $g$의 값을 부여해야 할 수도 있다.

일단 위 문제의 가장 자주 보이는 예시는, $\phi(n)$의 합을 구하기 위하여 $g(n) = 1$, $(f * g) (n) =  n$을 이용하는 것이다.

구현 노트.

• 구현 자체에 크게 어려운 점은 없으나, 역시 $\lfloor n/d \rfloor$가 매우 클 수 있으니 C++의 std::map을 이용해야 한다.
• 하지만 std::map을 사용하면 로그가 하나 더 붙는다. 로그를 붙이지 않고 계산을 할 수는 없을까?
• 애초에 $\lfloor n/d \rfloor \le n^{2/3}$이면 바로 값을 리턴하므로, $\lfloor n/d \rfloor > n^{2/3}$인 경우만 생각하자. 이 경우, $\lfloor n/d \rfloor$ 대신 그냥 $d$를 key 값으로 사용할 수 있다. 
• 그러면 key 값이 $n^{1/3}$ 이하로 매우 작으니, std::map 대신에 그냥 배열을 사용해도 무방하다. 물론, 이 경우 초기화 문제에 주의해야 할 것이다.

min_25 sieve

• 우선 필요한 정의부터 나열하고 시작하도록 하자. 모든 설명은 이 글을 따른다.
• 이제부터 모든 나눗셈은 C++ 방식의 정수 나눗셈이다. 자연수만 다루니 부호는 신경쓰지 말자. (10장에서도 이럴 걸 ㅋㅋ)
• $p_k$는 $k$번째 소수. 1-index를 사용한다. 즉, $p_1, p_2, p_3, \cdots = 2, 3, 5, \cdots$.
• $lpf(n)$은 $n$이 갖는 최소의 소인수. 단, $n=1$인 경우 $lpf(1)=1$로 정의한다.
• $F_{prime}(n)$은 $\sum_{2 \le p \le n} f(p)$이며, 이 합은 소수 $p$에 대한 것이다.
• $F_k(n)$은 $\sum_{p_k \le lpf(i), 2 \le i \le n} f(i)$이다. 특히, $\sum_{i=1}^n f(i) = F_1(n) + f(1) = F_1(n)+1$이다.
• 즉, $F_k(n)$은 최소 소인수가 $p_k$ 이상인 자연수에 대한 $f$ 값의 합이다. 이제 본격적으로 시작.

식 조작을 하자. 각 자연수의 최소 소인수가 무엇인가, 그리고 그 소인수를 몇 개 가지고 있는가로 분류하면, $$F_k(n) = \sum_{p_k \le lpf(i), 2 \le i \le n} f(i) = \sum_{k \le i, p_i \le n} \sum_{c \ge 1, p_i^c \le n} f(p_i^c) \left( 1 + F_{i+1}(n/p_i^c) \right)$$ 과 같다. 즉, $p_i^c$와 $1$의 곱, 또는 최소 소인수가 $p_{i+1}$이면서 $n/p_i^c$ 이하인 자연수의 곱으로 분류한다. 

한편, $n < p_k$면 $F_k(n) = 0$임을 사용하고, 적당히 식을 정리하면 $$F_k(n) = \left( \sum_{k \le i, p_i^2 \le n} \sum_{c \ge 1, p_i^{c+1} \le n} \left( f(p_i^{c+1}) + f(p_i^c) F_{i+1}(n/p_i^c) \right) \right) + \sum_{k \le i, p_i \le n} f(p_i)$$ 를 얻는다. 여기서는 기존 식에서 $p_i^2 > n$이면 나오는 항이 $f(p_i)$ 하나임을 이용했다. $F_{prime}$ 표현을 이용하면, $$F_k(n) = \left( \sum_{k \le i, p_i^2 \le n} \sum_{c \ge 1, p_i^{c+1} \le n} \left( f(p_i^{c+1}) + f(p_i^c) F_{i+1}(n/p_i^c) \right) \right) + F_{prime}(n) - F_{prime}(p_{k-1})$$ 을 얻는다. 이제 이 점화식을 써먹는 방법을 고민해보자.

• 위 점화식을 top-down으로 계산한다고 가정하자.
• 또 $\lfloor n/i \rfloor$ 형태의 값들만 등장한다. 이제 $\mathcal{O}(\sqrt{n})$개의 값들만 등장함은 익숙하다.
• 재귀 $k$ 값들을 보면, $p_i^2 \le n$을 아는 상태에서 $F_{i+1}(*)$ 값들을 호출함을 볼 수 있다. 
• 그러니 $F_k(n)$을 구해야 한다는 재귀호출이 나왔다면, $p_{k-1}$이 큰 값은 아니라는 것이다. 
• 왜 이게 중요하냐면, $F_{prime}(p_{k-1})$의 값이 전처리의 영역에 있다는 것을 설명하기 때문이다. 
• 즉, $n^{1/2}$ 정도의 값까지 $F_{prime}$의 값을 전처리하면, $F_{prime}(p_{k-1})$ 값은 $\mathcal{O}(1)$에 정리된다.
• 이제 $\lfloor n/i \rfloor$ 형태의 값들에 대해 $F_{prime}$ 값을 전부 구해야 한다는 문제가 있다.
• 만약 $f(p)$가 $p$에 대한 다항식이라면, 앞에서 배운 Lucy-Hedgehog 알고리즘으로 이를 전부 전처리할 수 있다.

시간복잡도 분석은 자료에 따르면 $\mathcal{O}(n^{1-\epsilon})$이다. 하지만 시간복잡도에 비해 매우 좋은 성능을 보인다.

이 알고리즘의 강점은 xudyh's sieve처럼 $g$를 고민해서 잡을 필요가 없다는 점, 그리고 $f$에 대해 가정하는 점이 적다는 점이다.

• 이 방법을 제작한 것으로 알려진 min_25가 추후 블로그에 $\mathcal{O}(n^{2/3})$ 시간복잡도를 갖는 강화판을 올렸다.
• 위 점화식을 살짝 변형 + 메모이제이션을 추가하면 $\mathcal{O}(n^{3/4})$ 알고리즘이 얻어진다고 한다. baihacker의 블로그 참고. 

구현 노트.

• 중요한 사실인데, $p_i < \sqrt{n} < p_{i+1}$이라고 해서 $n/p_i < p_{i+1}$인 것은 아니다.
• 그러니 전처리를 하는 소수/체의 범위를 조금 조심해서 고를 필요가 있다. 마음 편하게 $\sqrt{n} + 500$ 정도로 두자.

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

이 글에서는 다음 알고리즘을 다룬다.

• Lucy-Hedgehog의 Algorithm 소개 및 최적화.
• Meissel-Lehmer Algorithm 소개 및 최적화.

Lucy-Hedgehog Algorithm

• 글을 읽기 전에, 경건한 마음을 가지고 Project Euler에 가입하자. 그 다음, (쉬운) 10번 문제를 풀자.
• 그러면 여러 사람이 자신의 접근과 풀이, 코드를 공유하는 thread에 접근할 수 있다.
• 이제 https://projecteuler.net/thread=10;page=5#111677로 가면, 이 알고리즘의 기원을 볼 수 있다.

• Note. 동적계획법만 알고 있어도 이 알고리즘을 충분히 이해할 수 있다.

$n$ 이하의 소수의 개수를 $\mathcal{O}(n^{3/4})$에 계산하는 방법을 소개한다. 먼저, $n^{1/2}$ 이하의 소수들을 에라토스테네스의 체를 이용하여 전처리한다. 

이제 $dp[v][m]$을 에라토스테네스 체에서 $m$까지 보았을 때, $1$ 이상 $v$ 이하 자연수 중 지워지지 않은 것의 개수라고 하자.

• 즉, $dp[v][m]$은 소수거나, 모든 소인수가 $m$ 초과인 $1$ 이상 $v$ 이하 자연수의 개수이다.
• 특히, $(m+1)^2 > v$라면 모든 소인수가 $m$ 초과인 $v$ 이하 자연수는 무조건 소수이므로, $dp[v][m]$은 $v$ 이하 소수의 개수다.

이제 $dp$ 식을 생각해보자.

• 만약 $m$이 소수가 아니라면, 애초에 에라토스테네스 체에서 하는 일이 없다. 이때는 $dp[v][m]=dp[v][m-1]$.
• 만약 $m^2 > v$라면, $m$이 소수여도 이번에 새로 지워지는 값은 없을 것이다. 여전히 $dp[v][m]=dp[v][m-1]$.
• $m$이 소수라고 가정하자. 이번에 제거되는 값들은 최소 소인수가 $m$인 값들이다.
• 즉, $m$과 "최소 소인수가 $m-1$ 초과인 자연수"을 곱한 값들이 지워진다. 
• 이 값의 개수를 구하려면, $\lfloor v/m \rfloor$ 이하 자연수 중 최소 소인수가 $m-1$ 초과인 자연수의 개수를 구하면 된다.
• 그런데 $dp[\lfloor v/m \rfloor][m-1]$에 소수이거나, 최소 소인수가 $m-1$ 초과인 값의 개수가 저장되어 있다.
• 이 값에 $m-1$ 이하 소수의 개수를 빼면 원하는 값이 나올 것이다. 그 값은 $dp[m-1][m-1]$에 저장되어 있다.
• 그러므로, $dp[v][m] = dp[v][m-1] - (dp[\lfloor v/m \rfloor][m-1] - dp[m-1][m-1])$이라는 결과가 나온다.

이제 필요한 $dp$ 상태를 생각해보자. 우리는 $dp[n][\lfloor \sqrt{n} \rfloor]$가 필요하다.

• 먼저 두 번째 인자 $m$을 보자. 우리의 관찰에 의해, $m^2 > v$인 경우는 추가적으로 고려할 필요가 없다. 
• 특히, 우리가 다룰 모든 두 번째 인자 $m$ 값은 $\sqrt{n}$ 이하이다.
• 이제 첫 인자인 $n$을 생각하자. 필요한 첫 번째 인자의 값은 $\lfloor n/i \rfloor$ 형태의 값을 가짐을 알 수 있다.
• 즉, 필요한 첫 번째 인자의 값은 $1, 2, \cdots, \lfloor \sqrt{n} \rfloor$과 $\lfloor n/1 \rfloor, \lfloor n/2 \rfloor, \cdots , \lfloor n/\lfloor \sqrt{n} \rfloor \rfloor$다.
• 중요한 점은, 다루어야 할 첫 번째 인자의 값의 종류가 $\mathcal{O}(\sqrt{n})$개라는 점이다. 이제 준비 완료.

이제 본격적으로 알고리즘을 구체화한다.

• DP에 익숙한 사람들은 바로 눈치를 챘겠지만, 위 DP 식은 정말 토글링을 하기 좋은 식이다. 
• 그러니 $dp$를 1차원 배열로 정의하고, 순서대로 $m$ 값을 늘려나가보자. $dp$ 배열은 $dp[v][m]$을 갖고 있을 것이다.
• 그러기 위해서는 사실 $dp[v][1]$을 정의해야 하는데, 이는 시작하자마자 체에서 1을 제거했다고 치고 $dp[v][1]=v-1$로 둔다.
• 이제 큰 $v$값부터 순서대로 $dp[v][m]=dp[v][m-1]-(dp[\lfloor v/m \rfloor][m-1] - dp[m-1][m-1])$을 적용한다.
• 그러다가, $v \le m^2$인 시점이 오면 업데이트를 중단하고 다음 $m$으로 넘어가도록 한다. 
• 물론, 우리가 실제로 업데이트를 해야 하는 $v$의 값들은 앞서 언급한 $\mathcal{O}(\sqrt{n})$ 개의 값들이다.
• 마지막으로 얻은 $dp[n]$의 값이 우리가 원하는 $n$ 이하의 소수 개수다. 이게 알고리즘 끝이다. 예상보다 간단했을 것.
• 특히, 실제로 관심이 있던 각 값 $v$에 대하여, $dp[v]$는 $v$ 이하의 소수 개수다. 즉, 예를 들어, $n/2$, $n/3$ 이하의 소수 개수도 동시에 구해진다.

시간복잡도 분석을 하자.

• 이 알고리즘은 "실제로 계산해야 하는" $dp[v][m]$의 개수만큼 계산을 한다.
• 특히, 우리는 $m^2 \le v$인 경우에만 추가적인 계산을 하니, 우리가 다루는 각 $v$값에 대해 $\sqrt{v}$를 더한만큼 계산을 한다. 이는 $$\sum_{i=1}^{\sqrt{n}} \sqrt{i} + \sum_{i=1}^{\sqrt{n}} \sqrt{n/i} = \mathcal{O}(\sqrt{n} \cdot \sqrt{n}^{1/2}) + \mathcal{O}(\sqrt{n} \cdot \sqrt{n}^{1/2}) = \mathcal{O}(n^{3/4})$$ 다. 여기서 $1+1/\sqrt{2} + 1/\sqrt{3} + \cdots + 1/\sqrt{n} = \mathcal{O}(\sqrt{n})$을 이용하였다.

구현 노트.

• $dp$ 배열을 관리하기 위하여 가장 쉽게 생각할 수 있는 방법은 역시 C++의 std::map일 것이다. 하지만 로그가 붙는다.
• 로그를 붙이지 않고 구현을 하려면, $dp$를 말 그대로 배열에 저장해야 한다. 이를 위해서는, $v = 1, 2, \cdots, \lfloor \sqrt{n} \rfloor$에 대응되는 길이 $\lfloor \sqrt{n} \rfloor$의 배열과, $v = \lfloor n/1 \rfloor, \lfloor n/2 \rfloor, \cdots ,  \lfloor n/\lfloor \sqrt{n} \rfloor \rfloor$에 대응되는 길이 $\lfloor \sqrt{n} \rfloor$의 배열을 따로 두면 된다. 자세한 디테일은 github 참고.

추가 코멘트.

• 필자는 이 알고리즘을 많이 써먹었다. 이 알고리즘의 장점은, 확장이 쉽게 가능하다는 것이다.
• 예를 들어, 소수의 개수가 아니라 소수의 합, 소수의 제곱의 합 등을 원하는 경우도 $dp$ 식만 살짝 바꿔주면 해결된다.
• 심지어, $n$ 이하 $4k+1$ 꼴 소수의 개수/합, $4k+3$ 꼴 소수의 개수/합 역시 약간의 아이디어를 추가하면 할 수 있다.
• 확장이 쉽게 가능한 이유는, 역시 기본 원리인 dp 자체가 간단한 편이기 때문이다. 이해하기도 쉬운 편인 알고리즘이라 생각. 
• 사실 시간복잡도 분석을 빡빡하게 하지 않았다. 실제로는 $m$이 소수인 경우만 중요하다는 사실도 역시 써먹을 수 있다.
• 빡빡한 시간복잡도 분석을 위해서는, $n$ 이하 소수의 개수가 약 $n/\ln n$개 존재한다는 Prime Number Theorem을 참고하라.

이 알고리즘을 $\mathcal{O}(n^{2/3})$ 수준으로 최적화시킬 수 있다. 여기서 "로그가 시간복잡도에 얼마나 붙는지"는 따지지 않겠다.

• 지수를 3/4에서 2/3으로 내리는 것은 11단원에서도 등장한다. 보통은, "에라토스테네스의 체를 더 써먹자"의 방식으로 최적화가 된다.
• 사실 우리는 기존 알고리즘에서 에라토스테네스의 체로 $1$부터 $n^{1/2}$까지의 소수 판별만 했다. 체를 더 써먹자.

이제부터, 기존 점화식 $dp[v][m] = dp[v][m-1] - (dp[\lfloor v/m \rfloor][m-1] - dp[m-1][m-1])$을 top-down으로 돌려보자.

목표는 실제로 계산을 전부 하는 것이 아니라, 어떤 $dp$ 값들이 필요한지를 탐색하는 것이다. 시작은 물론 $dp[n][\lfloor \sqrt{n} \rfloor]$을 호출하는 것이다.

• $dp[m-1][m-1]$은 $m$ 미만의 소수 개수와 같다. 이는 에라토스테네스의 체로 처리할 수 있으니, $dp$ 값으로 간주하지 않는다.
• 만약 $m=1$이라면, $dp[v][m]=v-1$임을 이미 알고 있으니, 이를 필요한 $dp$ 값으로 간주하지 않는다.
• 만약 $v \le n^{2/3}$이라면, "$dp[v][m]$을 계산해야 한다"라는 정보를 저장하고 리턴한다.
• 재귀 과정에서 이미 방문한 위치 $v, m$에 다시 도달했다면, 이미 $dp[v][m]$을 계산하기 위해 필요한 $dp$ 위치를 알고 있는 상태이므로 바로 리턴.
• 즉, $(v,m)$의 자손이 $(v,m-1)$과 $(\lfloor v/m \rfloor ,m-1)$인 이진트리를 생각하고, $m=1$이거나 $v \le n^{2/3}$인 경우를 리프 노드로 둔다.

시간복잡도를 분석하자.

• 사실 1. 위 재귀 루틴의 시간복잡도는 $\mathcal{O}(n^{2/3})$이다. 
• 증명. $dp$ 상태 중, $v \ge n^{2/3}$인 것의 개수를 구하면 충분하다. 실제로 관심 있는 $v$만 고려하면, 이는 $$\sum_{i=1}^{n^{1/3}} \sqrt{n/i} = \mathcal{O}(\sqrt{n} \cdot \sqrt{n^{1/3}}) = \mathcal{O}(n^{2/3})$$
• 사실 2. 결과적으로 우리가 "계산해야 하는 $dp$ 값"들의 위치는 $\mathcal{O}(n^{2/3})$개이다.
• 증명. 물론, 이는 사실 1에서 바로 따라온다. 애초에 시간복잡도가 저러니 저장한 위치의 개수도 마찬가지.

이제 우리는 $\mathcal{O}(n^{2/3})$개의 $dp$ 값만 구해주면 된다. 이 값들을 다 구하면, 다시 재귀를 돌려 답을 얻을 수 있다.

즉, 이전에 "$dp[v][m]$을 계산해야 한다"고 하고 리턴한 경우, 이번에는 실제로 계산된 $dp[v][m]$의 값을 리턴하면 된다!

이제부터 우리는 알고리즘의 영역으로 온다. 처음으로 이 가이드에서 C++ STL에 없는 자료구조를 사용한다.

• 오프라인으로 $dp$ 값들을 계산한다. $m$의 값이 증가하는 순서로 필요한 값들을 계산할 것이다.
• $dp[v][m]$을 계산하는 것은, 결국 정의상 $m$까지 체를 돌렸을 때 $1$부터 $v$까지 살아남은 자연수의 개수를 구하는 것이다.
• 현재 자연수가 남아있는지 여부를 0/1로 표현하면, 이는 결국 $[1, v]$에서 구간합을 구하는 것과 마찬가지다.
• 에라토스테네스의 체에서 값을 지우는 것은 결국 해당 자연수에 대응되는 값을 1에서 0으로 바꾸는 것과 같다.
• 그러니 이 문제는 사실 Point Update + Segment Query. 세그먼트 트리로 해결할 수 있다. (속도를 위해 펜윅을 쓰는 것을 추천한다.)

이제 전체 문제가 $\mathcal{O}(n^{2/3})$ 시간에 해결된다. 필요한 에라토스테네스의 체 크기가 $n^{2/3}$이었음에 주목하자.

• 이 최적화 방법은 앞서 언급한 Lucy-Hedgehog 알고리즘의 여러 확장에도 바로 적용된다.
• 예를 들어, 소수의 합을 구하고 싶다면 자연수 $k$가 남아있는지 여부를 0/$k$로 나타내면 된다.

Meissel-Lehmer Algorithm

• 사실 Lucy-Hedgehog와 큰 차이는 없다. 개인적으로 많이 사용하지 않아, 얼마나 확장이 쉽게 되는지는 잘 모르겠다.
• 하지만 알고리즘의 원리가 근본적으로 비슷하니, 아마 크게 다르지는 않을 것이라고 생각한다. 독자에게 맡긴다.

이번에는 notation을 살짝 바꾼다. 그래도 Lucy-Hedgehog 알고리즘과 비슷한 점이 많으니, 비교해가며 읽자.

• $p_k$를 $k$번째 소수로 정의한다. 단, $1$-index를 이용한다. 즉, $p_1, p_2, p_3 = 2, 3, 5$.
• $\pi(x)$를 $x$ 이하의 소수의 개수로 정의한다. 이는 매우 standard 한 notation이다. 예를 들어, $\pi(x) \approx x / \ln x$.
• $\phi(x, a)$를 모든 소인수가 $p_{a+1}$ 이상인 $x$ 이하 자연수 개수라 하자. 중요한 점은, $p_1, p_2, \cdots ,p_a$ 역시 제거했다는 것이다.
• 또한, $\phi(x, a)$를 생각할 때 $1$도 포함한다. 그러니, 정확하게 말하자면 $p_1, p_2, \cdots , p_a$의 배수를 제거한 결과를 보는 것이다.
• 이어서, $P_k(x, a)$를 모든 소인수가 $p_{a+1}$ 이상인 $x$ 이하 자연수 중, 소인수 개수가 정확히 $k$개인 것의 개수라 하자.
• 이때, "소인수 개수"란, $n=p_1^{e_1} \cdots p_k^{e_k}$라 할 때 $e_1 + e_2 + \cdots + e_k$로 정의된다. 
• 특히, $p_{a+1}^k > x$라면 $P_k(x, a) = 0$이 성립한다. 이제 필요한 용어의 정의는 끝났다.

$\pi(x)$를 계산하려면, 필요한 값은

• 우선 억울하게 제거당한 소수들 $p_1, p_2, \cdots , p_a$가 있다. 총 $a$개.
• 이제 모든 소인수가 $p_{a+1}$ 이상인 $x$ 이하 자연수가 필요하다. 총 $\phi(x, a)$개.
• 그런데 $1$은 소수가 아니니, $1$을 빼주어야 할 것이다.
• 또한, 소인수 개수가 $2$ 이상인 것에는 관심이 없으니, $P_2(x, a) + P_3(x,a) + \cdots$를 빼야한다. 즉, $$\pi(x) = a + \phi(x, a) - 1 - \sum_{k=2}^\infty P_k(x, a)$$
• 물론, $p_{a+1}^k > x$면 $P_k(x, a)=0$이므로, 실제로는 무한합을 계산할 필요가 없다. 

특히, 우리는 $a = \pi(x^{1/3})$을 선택할 것이다. 그러면 $p_{a+1}^3 > x$이므로, $P_3$부터 쓸모가 없다.

그러므로, $\pi(x) = a + \phi(x, a) - 1 - P_2(x, a)$를 계산하면 된다. 이제 각 항을 계산할 전략을 설계해보자.

• $a$의 계산. 단순히 $x^{1/3}$ 이하의 소수의 개수이므로, 에라토스테네스의 체로 정리된다.

• $\phi(x, a)$의 계산. 이 부분은 Lucy-Hedgehog 알고리즘의 dp와 크게 다르지 않다.
• $\phi(x, a)$를 계산해보자. 일단 $\phi(x, a-1)$에서 시작하면, $a-1$에서 $a$로 넘어갈 때 새로 제거되는 수의 개수를 세면 된다.
• 이 값을 계산하려면, "최소 소인수가 정확히 $p_a$인 $x$ 이하의 자연수"의 개수를 세면 된다.
• 이는 결국 "모든 소인수가 $p_{a}$ 이상인 $\lfloor x/p_a \rfloor$ 이하의 자연수"의 개수와 같다. 이는 $\phi(\lfloor x/p_a \rfloor, a-1)$.
• 그러므로, dp 식 $\phi(x, a) = \phi(x, a-1) - \phi(\lfloor x/p_a \rfloor, a-1)$이 성립함을 알 수 있다.

• 우선 Lucy-Hedgehog에서 본 DP와 여러 성질이 겹친다. $\lfloor x/i \rfloor$ 형태의 값만 필요하다는 점 역시 똑같다.
• 토글링이 가능한 DP 형태라는 점도 똑같다. 그러니 이미 살펴본 접근 방식과 거의 동일한 방법이 여기서도 먹힌다.
• 하지만 이 방법의 중요한 점은, $p_a^2 > x$여도 $\phi(x, a)$와 $\phi(x, a-1)$이 다르다는 것이다. 
• 에라토스테네스의 체와 다르게 여기서는 $p_a$도 제거하므로, $\phi(x, a) = \phi(x, a-1) - 1$이 성립하겠다.
• 물론, $x < p_a$면 $\phi(x, a) = \phi(x, a-1)$이 성립함은 당연하다.
• 그래서 기존 $\mathcal{O}(n^{3/4})$ 알고리즘을 적용하기가 조금 귀찮아진다. $1$을 진짜 정직하게 빼주면 당연히 시간복잡도가 유지되지 않는다.
• 그러니 $p_a^2 > x$면 $\phi(x, a)$를 업데이트 하지 않고, 대신 "실제로 $1$이 몇 번 더 빠졌어야 하는가"를 필요할 때 계산해주자.

• 재밌는 점은 $\mathcal{O}(n^{2/3})$ 최적화 알고리즘은 그대로 적용된다는 것이다.
• 언급한 사실 1, 사실 2 역시 그대로 유지되며, 오프라인 쿼리 + 세그트리/펜윅트리로 $\phi$를 계산할 수 있음도 동일하다.

• $P_2(x, a)$의 계산. 이 값은 다행히도, 정직하게 계산할 수 있다.
• $P_2(x, a)$는 결국 $a<b \le c$를 만족하면서 $p_bp_c \le x$인 $(b, c)$의 순서쌍의 개수를 구하는 것이다.
• $b$를 고정하고, $c$의 개수를 세자. 그러면 $p_c \le x/p_b$와 $b \le c$가 동시에 성립해야 하니, $\pi(\lfloor x/p_b \rfloor) - (b-1)$이다.
• 가능한 $b$의 범위는 물론 $a+1$부터 $\pi(x^{1/2})$까지다. 즉, 이 결과를 합치면 우리는 $$P_2(x, a) = \sum_{b=a+1}^{\pi(x^{1/2})} \left( \pi(\lfloor x/p_b \rfloor) - (b-1) \right)$$ 를 얻는다. 이제 이 값을 빠르게 계산하면 된다. 여기서 핵심은 주어진 범위의 $b$에 대하여, $x/p_b < x^{2/3}$이 성립한다는 것이다. 
• 그러니 실제로 필요한 것은 $x^{2/3}$까지의 에라토스테네스의 체 결과다. 그러니, 모든 계산을 $\mathcal{O}(n^{2/3})$ 시간에 할 수 있다.

이를 종합하면 우리가 원하는 알고리즘이 나온다. 최적화로 $\mathcal{O}(n^{2/3})$까지 가는 방법도 설명하였다. 그런데,

• 소수 카운팅 알고리즘의 속도를 측정할 수 있는 가장 쉬운 방법은, Codeforces의 Four Divisors 문제를 해결하는 것이다.
• 당장 Editorial에서도 이 글에서 다룬 $\mathcal{O}(n^{2/3})$ 테크닉을 언급한다. 자세한 설명이 있으니 읽어보는 것도 좋다.
• 그런데, 정작 submission을 실행 속도로 정렬한 다음 몇 개의 제출을 읽어보면 오프라인이고 세그트리고 다 없다.
• 즉, 시간복잡도로 표현하는 속도와 실제 실행 속도와의 괴리감이 약간 있는 부분이라고 볼 수 있겠다. 이제 실제로 빠른 알고리즘을 간략하게 설명한다.
• (물론, 효율적인 코드를 매우 잘 작성하는 분들이라면 $\mathcal{O}(n^{2/3})$ 알고리즘으로 저 분들을 이길 수도 있을 것이다.)

• 먼저, 에라토스테네스의 체를 많이 돌린다. $n=10^{11}$ 기준으로, 500만 정도 크기로 돌린 것 같다.
• 앞서 우리는 $\pi(x^{1/3})$에서 식을 자르는 방식으로 $P_3(x, a)$ 항을 날렸다. 이번에는 $\pi(x^{1/4})$에서 식을 잘라서 $P_3(x, a)$ 항을 살린다.
• 즉, $a = \pi(x^{1/4})$를 잡고, $\pi(x) = a + \phi(x, a) - 1 - P_2(x, a) - P_3(x, a)$를 계산한다.
• 이제 추가적으로 $P_3$에 대한 전개식을 계산해야 한다. $p_bp_cp_d$ 형태를 보면 되는데, 진짜 정직하게 $b, c$에 대한 이중 for 문을 돌린다.

• 이제 $\phi(x, a)$를 계산해야 한다. 이 값을 계산하기 위해서, 상당한 최적화를 한다.
• $\phi(x, a)$의 context에서, 특정 값이 살아있는지 여부는 그 값을 $p_1p_2 \cdots p_a$로 나눈 나머지에 의해 결정된다.
• 그러니, $a$가 작은 경우 ($a \le 7$을 이용한다) 아예 어떤 나머지가 살아남는지를 통째로 전처리를 한다. 크기 $p_1 \cdots p_a$짜리 전처리다.
• 이러면 $a \le 7$인 입력이 들어온 경우, 아예 $\mathcal{O}(1)$에 $\phi(x, a)$를 반환할 수 있다. 이것이 첫 최적화.

• 만약 $x < p_a^2$이라면, $\phi(x, a)$는 사실 $p_{a+1}$ 이상 $x$ 이하 소수의 개수이다. 
• 그러니, 만약 $x$ 이하의 소수 개수가 체로 전처리된 상황이라면, 바로 답을 반환할 수 있다. 이것이 두 번째 최적화.

• 만약 $x < p_a^3$이라면, $\phi(x, a) = \pi(x) - a + 1 + P_2(x, a)$다. 이때, $\pi(x)$의 값이 이미 체로 전처리된 상황이라 하자.
• 이 경우에는, $P_2(x, a)$를 기존 방식과 동일한 방법으로 구해 $\phi(x, a)$를 얻는다. 이것이 세 번째 최적화.

• 만약 이 모든 경우에 해당하지 않는다면, 그제서야 점화식 $\phi(x, a) = \phi(x, a-1) - \phi(\lfloor x/p_a \rfloor, a-1)$을 호출한다.
• 사실 기존 Meissel-Lehmer 알고리즘에서 $a = \pi (x^{1/2})$를 선택하는 것도 방법이다. 방법은 정말 많다..

필자는 가능하면 Lucy-Hedgehog를 쓰고 (특히 Project Euler와 같은 시간제한이 없는 경우에는 더더욱) 만약 시간제한이 빡빡하다면 Meissel-Lehmer의 최적화된 구현체, 즉 위 Four Divisors 문제의 링크에 있는 알고리즘을 사용하여 문제를 푼다. 모두 흥미로운 알고리즘들이다. 이 정도로 이번 내용은 끝.

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

이 글에서는 다음 4가지 문제를 다룬다. 모두 "유클리드 알고리즘에서 등장한 문제 축소 방식"을 이용한다.

• Lattice Point Counting. $ax+by \le c$, $x, y \ge 1$을 만족하는 격자점 $(x, y)$의 개수를 세는 법.
• Linear Inequality with Modulo. $L \le Ax \pmod{M} \le R$을 만족하는 최소의 음이 아닌 정수 $x$를 구하는 법.
• Modulo Maximization. $1 \le x \le C$이며 $Ax \pmod{M}$이 최댓값이 되도록 하는 $x$를 구하는 법.
• Minimal Fraction. $A/B < x/y < C/D$이며, $y$가 최소인 분수 $x/y$를 찾는 법. 단 $A, C, D \ge 0$, $B \ge 1$.
• 그리고 사실 이 내용에서 중요한 사례라고 볼 수 있는, continued fraction, 연분수를 매우 간단하게 소개한다. 연분수에 대한 내용은 정수론 책을 찾아보거나, 링크에 걸린 글을 읽는 것을 추천. 필자가 연분수에 대한 내용을 직접 소개하기에는 지식이나 엄밀한 이해/직관이 부족한 것 같다. 

Lattice Point Counting

• $ax + by \le c$, $x, y \ge 1$을 만족하는 격자점 $(x, y)$의 개수를 세자. 이를 $(a, b, c)$에 대한 문제라고 부르겠다.
• 일반성을 잃지 않고, $a \ge b$를 가정한다. $a = qb+r$인 $q \ge 1$과 $0 \le r < b$가 존재한다.
• $ax + by \le c$는 $(qb+r)x + by \le c$, 즉 $b(qx+y) + rx \le c$와 같다.
• 그러므로, $bu+rv \le c$이며 $u \ge qv+1$, $v \ge 1$를 만족하는 격자점 $(u, v)$의 개수를 세면 된다. ($x=v$, $y=u-qv$에 대응된다)
• $u$를 고정했을 때 $v \le (c-bu)/r$과 $v \le (u-1)/q$라는 두 부등식이 있다. 둘 중 어느 것이 더 "강한 부등식"인지를 기준으로 경우를 나눈다.
• 계산을 열심히 하면, 그 기준이 $u = \lfloor cq/a \rfloor$임을 알 수 있다. 이제 경우를 나눠보자.

• $u > \lfloor cq/a \rfloor$인 경우, $v \le (c-bu)/r$만 만족하면 된다. 그러니 $bu+rv \le c$만 생각하면 된다.
• $b(u-\lfloor cq/a \rfloor) + rv \le c - b \lfloor cq/a \rfloor$이어야 하며, 원하는 조건은 $u-\lfloor cq/a \rfloor , v \ge 1$이다.
• 이 경우는, 문제를 $(b, r, c - b \lfloor cq/a \rfloor)$에 대한 문제로 축소시켰다.

• $u \le \lfloor cq/a \rfloor$인 경우, $v \le (u-1)/q$만 만족하면 된다. 이 식은 $qv + 1 \le u \le \lfloor cq/a \rfloor$과 같다.
• $v=1$이라면, 가능한 $u$의 개수는 $\lfloor cq/a \rfloor - q$개다. $v$가 $1$ 증가할 때마다, $u$로 가능한 개수가 $q$개 감소한다.
• 결국 이 경우에는, 등차수열의 합으로 가능한 경우의 수를 표현할 수 있다. $\mathcal{O}(1)$ 선으로 정리할 수 있다.

그러니 문제를 유클리드 호제법과 동일한 시간복잡도로 해결할 수 있다. 

$\sum_{k=1}^n \lfloor (a + bk) / c \rfloor$ 형태의 값 역시 같은 알고리즘으로 계산할 수 있다.

• 먼저 $a = qc + r$이면, $\lfloor (a+bk) / c \rfloor = q + \lfloor (r+bk)/c \rfloor$이다. 
• 그러니, $qn$을 값에 더한다음 $a$를 $r$로 바꿀 수 있다. 이제부터 $0 \le a < c$인 경우만 봐도 무방.
• 위 식은 $1 \le x \le n$, $1 \le y$이며 $a+bx \ge cy$인 격자점 $(x, y)$의 개수를 센다.
• 이 식을 변형하여, $cy+b(n+1-x) \le (n+1)b+a$ 형태를 얻는다. 
• $t = n+1-x$라 하면, 목표는 $1 \le t \le n$과 $cy+bt \le (n+1)b + a$다.
• 애초에 $t \ge n+1$일 수 없는게, 그러면 $cy+bt \ge c + b(n+1) > a + b(n+1)$이 성립한다.
• 그러니, 단순히 $cy + bt \le (n+1)b + a$, $y, t \ge 1$인 $(y, t)$의 개수를 세면 끝이다. 즉, $(c, b, (n+1)b+a)$ 문제를 풀자.

Linear Inequality with Modulo

• $L \le Ax \pmod{M} \le R$을 만족하는 최소의 음이 아닌 정수 $x$를 구하는 법을 찾자. 이를 $(A, M, L, R)$ 문제라 하자.
• 먼저 해가 존재하는지 찾자. $Ax \pmod{M}$은 $\text{gcd}(A, M)$의 배수 전체를 표현할 수 있다. (2단원 내용)
• 그러니 $L$과 $R$ 사이에 $\text{gcd}(A, M)$의 배수가 존재하는지만 확인하면 된다. 
• 만약 $L > R$이라면, 구간이 "한 바퀴를 돌았다"고 생각하여, $0$을 문제의 답으로 본다.
• 즉, $[L, R]$이라는 구간을 수평선의 구간으로 보지 말고, 둘레가 $M$인 원의 한 구간으로 보는 것이다.

• 이제 본격적으로 문제를 해결하자. 먼저 $0$이 조건을 만족하는지 확인하자. $L > R$이나 $L=0, R=0$인지 보면 된다. 
• 만약 $2A > M$이라면, $L \le Ax \pmod{M} \le R$과 $M-R \le (M-A)x \pmod{M} \le M-L$이 동치임을 이용하자.
• 이를 사용하면, 문제를 $(M-A, M, M-R, M-L)$ 문제로 옮길 수 있다. 이제 $2(M-A) < M$이다. 이제부터 $2A \le M$을 가정한다.
• 만약 $L \le Ax \le R$인 (단, $L, R$ 모두 $[0, M)$에 속한다. $Ax$는 $M$으로 나눈 나머지로 보지 않는다.) $x$가 있다면, 그 중 가장 작은 것을 반환하자.

• 그렇지 않다는 것은, 일단 $At < L \le R < A(t+1)$인 음이 아닌 정수 $t$가 있음을 의미하고, $R-L+1 < A$임을 얻는다. 이제 진짜 시작. 
• $L \le Ax \pmod{M} \le R$인 최소 $x$를 찾는 건, $L + My \le Ax \le R + My$인 정수 $x$가 있도록 하는 최소 음이 아닌 정수 $y$를 찾는 것이다.
• 이는 다시 $L \le Ax- My \le R$과 같은데, 이걸 $\pmod{A}$로 보면 $L \pmod{A} \le (-My) \pmod{A} \le R \pmod{A}$이다.
• 그러니 $y$의 값을 구하기 위해서는 $(- M \pmod{A}, A, L \pmod{A}, R \pmod{A})$ 문제를 풀어야한다.
• $y$의 값을 구하면, $x$의 값을 구하기 위해서는 단순히 $L + My \le Ax$인 최소의 $x$를 찾기만 하면 된다.

우리는 modulus에 해당하는 값을 $M$에서 $A$로 줄였으며, $2A \le M$이니 이는 최소 반 줄인 것이다. 그러니 로그 시간.

• 최소의 음이 아닌 정수 $x$가 아니라, $T$ 이상의 정수 $x$로 조건을 바꿔도 해결할 수 있다.
• 이 경우, $x = y + T$라고 쓰고 조건을 $(L-AT) \pmod{M} \le Ay \pmod{M} \le (R-AT) \pmod{M}$이라 바꾸자.
• 이제 위 조건을 만족하는 최소의 음이 아닌 정수 $y$를 구한 뒤, $x = y+T$를 대입하면 끝이다.
• 이를 반복하면, $L \le Ax \pmod{M} \le R$을 만족하는 자연수 $x$를 순서대로 나열할 수도 있다.

Modulo Maximization/Minimization

• $1 \le x \le C$를 만족시키면서 $Ax \pmod{M}$이 최댓값이 되도록 하는 $x$를 구하자. 이를 MAX-$(A, M, C)$ 문제라 하자.
• $1 \le x \le C$를 만족시키면서 $Ax \pmod{M}$이 최솟값이 되도록 하는 $x$를 구하자. 이를 MIN-$(A, M, C)$ 문제라 하자.
• MAX 문제나 MIN 문제나 사실 같은 문제다. $A$를 $-A$로 바꿔서 생각하면 두 문제가 같은 문제가 되기 때문이다.
• 이 문제는 "Linear Inequality with Modulo" 문제에 이분탐색을 추가해서 풀 수 있다.
• adamant의 블로그 글이 시각화도 좋고, 설명도 좋고, 추가할 점도 특별히 없는 것 같다. 링크를 거는 것으로 대체한다.
• 그런 줄 알았는데, 실제로 코드를 비교해보니 약간의 문제가 있는 것으로 보인다. 일단은 링크의 아이디어만 알아두자.

Minimal Fraction

• $A/B < x/y < C/D$를 만족시키면서, $y$가 최소인 분수 $x/y$를 찾는 경우. 이를 $(A, B, C, D)$ 문제라 하자.
• 위에서도 언급했지만, $A, C, D \ge 0$, $B \ge 1$을 가정한다. 우리의 목표 $x, y$ 역시 음이 아닌 정수여야 한다.
• 사실 $y$를 최소화하는 것이나 $x$를 최소화하는 것이나 큰 차이가 없다. $y$가 결정되었다면, $Ay/B < x < Cy/B$인 최소의 $x$를 잡으면 된다. 
• $y$가 커지면, 자연스럽게 가능한 $x$의 범위도 "오른쪽으로" 이동하게 될 것이다.
• $D=0$인 경우, $C/D$를 $\infty$로 간주하도록 하겠다. 이 경우까지 처리해야 문제가 풀린다.
• 이 문제는 아래에 언급할 Stern-Brocot Tree를 이용한 풀이가 있다. 링크가 걸린 myungwoo님의 글을 참고하라.

• $A/B < C/D$가 성립하지 않는 경우, 답은 당연히 없고, 이 결과를 반환한다.
• $A=0$인 경우. 이때는 $0 < x/y < C/D$인 것만 생각하면 된다. $x=1$로 고정하고, $1/y < C/D$인 최소의 $y$를 잡으면 ok.
• $A/B$와 $C/D$ 사이에 (양 끝 제외) 정수가 존재하는 경우. $y=1$을 정하고 $x$를 존재하는 정수 중 가장 작은 것으로 잡으면 ok.
• $C/D$가 정수인 경우. 이때, $C/D=n$이라 하자. 이 경우, $A/B < n-1/y$인 최소 $y$를 잡고, $Ay/B < x < ny$인 최소 $x$를 잡으면 ok.
• 이제 $n \le A/B < C/D < n+1$인 자연수 $n$이 존재하는 경우만 남았다. 이 단계가 핵심.
• 목표는 $A/B < x/y < C/D$인 최소 $x, y$를 찾는 것이다. $n$을 빼준 다음, 전부 역수를 취해주자.
• 이제 목표는 $D/(C-nD) < y/(x-ny) < B/(A-nB)$인 최소의 $x,  y$를 찾는 것이다.
• 애초에 분모를 최소화하는 것이나 분자를 최소화하는 것이나 같은 문제이므로, $(D, C-nD, B, A-nB)$ 문제를 풀면 된다.
• 여기서 $n$이 $A, C$를 $B, D$로 나눈 몫임을 생각하면, 이 과정이 유클리드 호제법과 같음을 알 수 있다. 끝. 

관련 문제로 Google Code Jam 2019 Round 2 C번 문제가 있다. 아래 myungwoo님의 글도 그 당시 나왔다.

Continued Fraction

• $[a_0 ; a_1, \cdots , a_k] = a_0 + (1 / (a_1 + (1 / (a_2 + \cdots ) ) ) )$라고 정의하자. 이를 연분수라고 부른다.
• 모든 유리수는 유한한 연분수로 나타낼 수 있다. $a=qb+r$이라 하면, $a/b = q + r/b = q + 1/(b/r)$이다.
• 이는 $(a, b)$에 대한 문제를 $(b, r)$에 대한 문제로 줄인 것과 같다. 유클리드 호제법 그 자체...
• 연분수는 정말 다양한 성질들이 있다. 이는 학부 정수론 책을 참고하는 것을 추천.
• 또한, 위에서도 언급한 adamant의 Codeforces 블로그를 참고하는 것을 추천한다. 좋은 자료가 많다.
• 연분수와 밀접한 연관이 있는 개념으로 Stern-Brocot Tree가 있다. 이에 대한 글은 myungwoo님의 이 글을 참고하라. 
• 위 링크에서도 나와있지만, Stern-Brocot Tree로 유리수에서 효율적인 이분탐색을 하는 등 매우 흥미로운 알고리즘들이 많다.

관련 코드는 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 코드 참고.