11. more on multiplicative functions

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

등등, prime power만 보고 모든 값을 얻을 수 있는 함수들을 나열할 수 있을 것이다.

이제 multiplicative function의 여러 성질을 어렵지 않게 증명할 수 있다. 예시를 하나 보인다.

사실. $f(n)$이 multiplicative라면, $g(n) = \sum_{d|n} f(d)$ 역시 multiplicative.

증명. $g(p^k) = 1 + f(p) + \cdots + f(p^k)$인 multiplicative function이 됨을 보일 수 있다. 예를 들어, $$ g(p_1^{e_1} \cdots p_k^{e_k}) = \sum_{d_i \le e_i} f(p_1^{d_1} \cdots p_k^{d_k}) = \sum_{d_i \le e_i} f(p_1^{d_1}) \cdots f(p_k^{d_k}) = \sum_{d_1 \le e_1} f(p_1^{d_1}) \cdot \sum_{d_2 \le e_2} f(p_2^{d_2}) \cdots \sum_{d_k \le e_k} f(p_k^{d_k})$$이고 이는 정의상 $g(p_1^{e_1}) \cdots g(p_k^{e_k})$와 동일하다. 증명 끝. 이제 본격적인 주제로 들어가자.

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의 prefix sum을 계산하는 것으로 넘어가자.

xudyh's sieve

• 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$ 정도로 두자.