10. 소수의 개수 빠르게 세기

관련 코드는 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})$를 선택하는 것도 방법이다. 방법은 정말 많다..

실제로 어떻게 구현되었는지 보려면, 이 링크를 타고 여러 코드를 탐색해보자. 상당히 빠르다 :) 

$\mathcal{O}(n^{2/3})$ 알고리즘들이 메모리를 많이 잡아먹는데 (물론, 적게 먹는 variant가 있는 것으로 안다) 주의해야 한다.

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