Mertens Trick (xudyh's sieve)

모두 Happy New Year!

 

 

 

▣ 정수론 시리즈 ─ 순서대로 읽으시는 걸 추천합니다.

1. 뫼비우스 함수와 디리클레 합성곱

2. Double Counting과 Harmonic Lemma

3. 곱셈적 함수 테이블 만들기

4. Mertens Trick (xudyh's sieve) (현재 글)

 

Intro

  이번에도 문제 하나를 소개하면서 시작한다.

Problem. 오일러 토션트 함수의 합
양의 정수 $N$이 주어질 때, 다음 식의 값을 1초 안에 구하는 프로그램을 작성하시오. (단, $N \le 10^9$이다.) $$S_{\phi}(N) = \sum_{n=1}^N \phi(n)$$

 

  이전 글에서 소개했듯이, 오일러 토션트 함수의 테이블을 $\mathcal O(N \log N)$에 채울 수 있으므로 $\phi(n)$의 부분합 $S_{\phi}(N)$ 또한 $\mathcal O(N \log N)$에 구할 수 있다(Linear Sieve를 사용하면 $\mathcal O(N)$까지도 가능하다). 하지만 이 시간복잡도로는 $10^9$이라는 범위를 감당할 수 없다. 그리고 생각해 보면 부분합 테이블을 채우는 시간복잡도와 하나의 부분합만을 구하는 시간복잡도가 같은 건 좀 의아하지 않은가? 어딘가 시간복잡도를 개선할 수 있는 방법이 있을 것 같은데….

 

$\mathcal O(\sqrt N)$ Solution?

  뭐가 걱정인가? 우리는 이미 전전 게시글에서 $S_{f \ast g}(N)$의 값을 $\mathcal O(\sqrt N)$에 구하는 방법을 배우지 않았던가?

$$\begin{aligned} S_{f \ast g}(N) &= \sum_{d=1}^N f(d) S_g {\left( \left\lfloor \frac{N}{d} \right\rfloor \right)} \\ &= \sum_{i=1}^t S_g(a_i) (S_f(e_i) - S_f(s_i - 1)) &\cdots\ (\star) \end{aligned}​$$

이고($t$, $a_i$, $s_i$, $e_i$는 해당 게시글의 정의와 동일하며, 이후에도 마찬가지다), $\phi = {\rm Id} \ast \mu$임을 알고 있으므로

$$\begin{aligned} S_{\phi}(N) &= \sum_{d=1}^N d \cdot S_{\mu} {\left( \left\lfloor \frac{N}{d} \right\rfloor \right)} \\ &= \sum_{i=1}^t S_{\mu}(a_i) (e_i - s_i + 1) \end{aligned}​$$

가 된다. 이제 $S_{\mu}(n)$만 구할 수 있으면 문제를 풀 수 있다!

 

  ……어떻게?

 

메르텐스 함수(Mertens Function)

  사실 뫼비우스 함수의 부분합 $S_{\mu}(N)$은 특별히 메르텐스 함수라는 이름이 붙어 있을 정도로 수학계에서 유명한(?) 함수다.

Definition. 메르텐스 함수(Mertens Function)
양의 정수 $N$에 대해 메르텐스 함수 $M(N)$을 다음과 같이 정의한다.
$$M(N) := \sum_{n=1}^N \mu(n)$$

 

그러니 이제부터 뫼비우스 함수의 부분합을 $S_{\mu}(n)$ 대신 메르텐스 함수 $M(n)$으로 표기하기로 한다. 이 함수가 유명한 이유는 그 악명 높은 난제인 리만 가설과도 관련이 있기 때문이라고 하는데, 일단은 본 글의 주제와는 크게 관련이 없으므로 관심이 있는 분들은 알아서(...) 찾아보기 바란다. 여튼 $S_{\phi}(N)$을 구하려면 $M(n)$을 빨리 구할 수 있으면 되지만, 이것도 만만치 않아 보인다. 일단 메르텐스 함수에 대해서는 $(\star)$ 식을 적용할 수 없다! 뫼비우스 함수를 두 함수의 디리클레 합성곱으로 나타낼 방법을 찾지 못했기 때문이다. 게다가 설사 찾는다고 하더라도, 그 두 함수의 부분합을 빨리 계산할 수 있느냐도 문제다. 다시 말해, 위의 $\mathcal O(\sqrt N)$짜리 방법으로 $\phi(n)$의 부분합을 구하는 것은 힘들다는 소리다.

 

  그런데 뫼비우스 함수를 두 함수의 합성곱으로 나타내진 못해도, 뫼비우스 함수가 들어간 디리클레 합성곱 식은 있다. $$\mu \ast 1 = \varepsilon$$가 그 예시인데, 잘 보면 $1(n)$의 부분합 $S_1(N)$과 $\varepsilon(n)$의 부분합 $S_{\varepsilon}(N)$은 각각 $S_1(N) = N$, $S_{\varepsilon}(N) = 1$로 아주 쉽게 구할 수 있다! 혹시 「Double Counting과 Harmonic Lemma」 게시글의 후반부에서 던졌던 질문을 기억하는가? 그때

 

합성하는 두 함수 중 한 함수의 부분합과 합성곱의 부분합을 빠르게 구할 수 있을 때,

나머지 하나의 함수의 부분합도 좀 더 효율적으로 구할 수 있는가?

 

라는 질문을 던졌었다. 그러고 보니 오일러 토션트 함수도 ${\rm Id} = \phi \ast 1$라는 관계식이 성립하고, 거기에 $S_{{\rm Id}}(N) = \frac{N(N+1)}{2}$와 $S_1(N) = N$로 나머지 두 함수의 부분합을 쉽게 구할 수 있지 않은가? 따라서 위의 질문을 해결한다면, $S_{\phi}(N)$과 $M(N)$의 값도 빠르게 구할 수 있게 되는 것이다!

 

$\mathcal O(N^{3/4})$ Solution: xudyh's sieve

  본격적으로 위 질문에 대한 답을 찾아 보자. 우선 문제 상황을 정리해 보면,

 

• 우리는 어떤 수론적 함수 $f(n)$의 부분합 $S_f(N)$을 계산하려 한다.
• 일단 $S_f(N)$을 선형 시간(혹은 선형로그 시간) 정도에 구할 수는 있지만(앞의 게시글에서 보았듯이 보통 $f(n)$이 곱셈적 함수일 경우 이를 만족하게 된다), 그것보다 더 빠르게 구하고 싶다.
• 그런데 또 다른 수론적 함수 $g(n)$이 있어서, $S_g(N)$과 $S_{f \ast g}(N)$의 값을 빠르게(조금 더 명확히 이야기하자면 $\mathcal O(1)$에 가깝게) 구할 수 있다!
• 이때 $S_f(N)$을 어떻게 빠르게 구할 수 있는가?

정도가 될 것이다. 일단 $S_{f \ast g}$ 이야기가 나오니까, 위 $(\star)$ 식에서부터 출발하는 게 좋을 것 같다. 편의상 $f$와 $g$의 자리를 바꿔 보면(디리클레 합성곱은 교환법칙이 성립하기에 상관없다), $(\star)$ 식은

$$\begin{aligned} S_{g \ast f}(N) &= \sum_{d=1}^N g(d) S_f {\left(\left\lfloor \frac{N}{d} \right\rfloor \right)} \\ &= g(1)S_f(N) + g(2)S_f {\left(\left\lfloor \frac{N}{2} \right\rfloor \right)} + \cdots + g(N)S_f(1) \end{aligned}$$

가 되는데 여기서 우리가 관심 있는 것은 $S_f(N)$이므로, 그에 관해 정리해 보면 다음과 같다.

$$g(1)S_f(N) = S_{g \ast f}(N) - \sum_{d=2}^N g(d) S_f {\left(\left\lfloor \frac{N}{d} \right\rfloor \right)}$$

$$\color{red}{\therefore\, S_f(N) = \frac{1}{g(1)} \left( S_{g \ast f}(N) - \sum_{d=2}^N g(d) S_f {\left(\left\lfloor \frac{N}{d} \right\rfloor \right)} \right) \quad \cdots\ (1)}$$

 

  식 (1)을 잘 보니, $S_f(N)$을 구하는 데 더 작은 $n$들의 $S_f(n)$ 값이 필요한 재귀적인 식이다. 그렇다면 메모이제이션을 통해 위 식을 빠르게 계산할 수 있지 않을까? 여기서 중요해지는 건 실질적으로 $S_f(n)$을 계산하는 횟수가 얼마나 되느냐다. 일단 Harmonic Lemma에 의해, 하나의 $n$에 대해 $S_f(n)$을 구할 때는 더 작은 $S_f$를 $\mathcal O(\sqrt n)$번 계산하게 된다. 따라서

$$\begin{aligned} S_f(n) &= \frac{1}{g(1)} \left( S_{g \ast f}(n) - \sum_{d=2}^n g(d) S_f {\left(\left\lfloor \frac{n}{d} \right\rfloor \right)} \right) \\ &= \frac{1}{g(1)} \left( S_{g \ast f}(n) - \sum_{i=2}^t S_f(a_i) (S_g(e_i) - S_g(s_i - 1)) \right) \quad \cdots\ (2) \end{aligned}$$

로 정리할 수 있다. 그런데 식 (2)에서 $S_f(a_i)$를 구하는 과정을 또 생각해 보면 $S_f(\lfloor \frac{a_i}{d} \rfloor)$들을 계산해야 하고, 이들을 계산할 때 또 더 작은 $S_f$들을 계산해야 하고… 와 같이 이어질 텐데, 그렇다면 모든 재귀적인 과정을 통틀어 $S_f$를 너무 많이 계산하게 되는 게 아닐까? 다행히도, 이와 관련해서 또 하나의 보조정리가 준비되어 있다.

Lemma. Integer Division Lemma
임의의 양의 정수 $p$, $q$, $r$에 대해 다음이 성립한다.
$$\left\lfloor \frac{\left\lfloor \frac{p}{q} \right\rfloor}{r} \right\rfloor = \left\lfloor \frac{p}{qr}\right\rfloor$$

 

Proof.  정수 $k_1$, $r_1$, $k_2$, $r_2$에 대해 $p = qk_1 +r_1\,(0 \le r_1 \le q-1)$, $k_1 = rk_2 + r_2\,(0 \le r_2 \le r-1)$라고 하면, (좌변) $= \left\lfloor \frac{k_1}{r} \right\rfloor = k_2$이고 (우변) $= \left\lfloor \frac{qk_1 + r_1}{qr} \right\rfloor = \left\lfloor \frac{qrk_2 + qr_2 + r_1}{qr} \right\rfloor$이다. 이때 $0 \le qr_2 + r_1 \le qr_2 + q-1 \le qr-1$이므로 $k_2 = \left\lfloor \frac{qrk_2}{qr} \right\rfloor \le$ (우변) $\le \left\lfloor \frac{qrk_2 + qr - 1}{qr} \right\rfloor = k_2$가 되어 (우변) $= k_2 =$ (좌변)이 된다.     $\square$

 

  이 보조정리는 정수 나눗셈 연산을 계속해서 적용하는 건, 하나의 정수 나눗셈 연산으로 표현 가능하다는 것과 같다. 한마디로, 식 (1)을 재귀적으로 계산하는 과정에서 나오는 모든 $S_f$들은 전부 $S_f(\lfloor \frac{N}{k} \rfloor)$ 꼴로 나타낼 수 있다는 거다! 따라서 Harmonic Lemma에 의해, $S_f(N)$을 계산할 때 필요한 모든 서로 다른 $S_f(n)$은 $S_f(a_1)$, $S_f(a_2)$, $\cdots$, $S_f(a_t)$들뿐이고 이는 많아야 $2\lfloor \sqrt N \rfloor$개다!

 

  자, 우리에게 필요한 열쇠는 모두 주어졌다. 이제 식 (1)에 메모이제이션을 적용했을 때의 시간복잡도를 분석해 보자. $a_i$가 작은 $S_f(a_i)$부터(즉 $S_f(a_t)$부터 역순으로) 구하고 메모이제이션한다고 생각하면, 위의 논의에 의해 현재 $a_i$에서 $S_f(a_i)$를 계산할 때 필요한 $S_f(n)$들은 이미 메모이제이션이 되어 있을 것이다. 즉 단순히 $2\left\lfloor \sqrt{a_i} \right\rfloor$개의 값만 더하면 되어, $S_f(a_i)$를 구하는 데 $\mathcal O(\sqrt{a_i})$번의 연산이 필요하다! 따라서 $S_f(N)$을 계산할 때의 총 연산량 $T(N)$은

$$T(N) = \sum_{i=1}^t \mathcal O( \sqrt{a_i}) = \mathcal O {\left( \sum_{i=1}^t \sqrt{a_i} \right)}$$

이고, Harmonic Lemma 증명 과정에서 알 수 있다시피

$$\begin{aligned} \sum_{i=1}^t \sqrt{a_i} &\le \sum_{i=1}^{\lfloor \sqrt N \rfloor} \sqrt{\left\lfloor \frac{N}{i} \right\rfloor} + \sum_{i=1}^{\lfloor \sqrt N \rfloor} \sqrt i \\ &\cong \int_0^{\sqrt N} \sqrt{\frac{N}{x}}\,dx + \int_0^{\sqrt N} \sqrt x\,dx = \color{magenta}{\frac{8}{3}N^{\frac{3}{4}}} \end{aligned}$$

가 되어, 결국 식 (1)에 메모이제이션을 버무렸을 때의 시간복잡도는 $T(N) = \mathcal O(N^{3/4})$가 된다! (경이로운 시간복잡도이지 않은가?) 이로써 우리는 우리의 희망 사항이었던 $S_f(N)$을 $\mathcal O(N)$보다 빠르게 구하기에 성공했다! 이렇게 $g$와 $f \ast g$의 부분합을 이용해서 $f$의 부분합을 선형 시간보다 빠르게 구하는 알고리즘을 xudyh's sieve(혹은 메르텐스 함수의 이름을 따서 Mertens Trick이라고도 한다)이라 부른다.

잠깐! xudyh's sieve에서 'xudyh'는 대체 무엇을 뜻하는 것인가?
'xudyh'라는 명칭이 하도 괴상해서 대체 무슨 뜻인지, 발음은 어떻게 하는 것인지(...)에 대해 질문하는 것을 꽤 봤다. 예전에 내가 타 블로그에서도 밝혔듯이, 'xudyh'는 codeforces의 랭커 Retired_MiFaFaOvO의 예전 코포 핸들이자 Project Euler 포럼 내의 닉네임이다. (이 사실은 이 댓글에서 확인할 수 있다.) 덧붙여, 나는 '수디'와 '주디'의 중간 정도로 발음하는 편이다.

 

$\mathcal O(N^{2/3})$ Solution: Precomputation

  그런데 $\mathcal O(N^{3/4})$라는 시간복잡도를 도출한 건 좋았지만, 조금 애매하다. 앞에 $\frac{8}{3}$이라는 상수까지 붙어 있어서, $N = 10^9$이면 $\frac{8}{3} N^{3/4} \fallingdotseq 1.5 \times 10^7$이 되어 재귀를 돌기에는 상당히 부담스러운 규모이고, 시간제한도 아슬아슬할 수도 있다. 시간복잡도를 조금 더 줄여볼 수는 없을까? 생각해 보니 우리가 아직 이용하지 않은 조건이 있었다! 그건 바로 $S_f(N)$을 선형 시간(혹은 선형로그 시간)에 구할 수 있다는 조건이다. 만약 작은 값에 대해서는 미리 함수 테이블을 선형 시간에 만들어서 전처리해 놓는다면 어떨까?

 

  구체적으로, 어떤 $K \ge \sqrt N$에 대해 $K$ 이하인 모든 $n$의 $S_f(n)$ 테이블을 미리 $\mathcal O(K)$에 전처리하고 $S_f(N)$을 구한다고 했을 때의 시간복잡도를 분석해 보면(편의상 Big-O Notation은 생략하였다),

$$\begin{aligned} T(N) &\cong K + \sum_{i=1}^{\frac{N}{K}} \sqrt{\left\lfloor \frac{N}{i} \right\rfloor} \quad(\because K \ge \sqrt N) \\ &\cong K + \int_0^{\frac{N}{K}} \sqrt{\frac{N}{x}}\,dx = K + \frac{2N}{\sqrt K} \end{aligned}$$

가 되고, 이 식은 산술-기하평균 부등식에 의해

$$K = \frac{2N}{\sqrt K} \quad \Rightarrow \quad K = (2N)^{2/3} \,(\,\ge \sqrt N)$$

일 때 최솟값 $T(N) = 2(2N)^{2/3}$을 가진다. 즉 시간복잡도는 $\color{royalblue}{T(N) = \mathcal O(N^{2/3})}$이 되어, 차수가 $\frac{3}{4}$에서 $\frac{2}{3}$으로 줄었다! 이 정도의 시간복잡도면 $N = 10^9$일 때도 1초 내에 충분히 돌아간다. 이를 오일러 토션트 함수에 적용하여 코드를 작성하면 글의 서두에서 소개했던 문제가 해결된다! 참고로 xudyh's sieve라는 이름을 위에서 한발 먼저 소개했지만, 보통 xudyh's sieve라 하면 전처리를 추가한 버전인 이 $\mathcal O(N^{2/3})$짜리 알고리즘을 가리키는 경우가 많으니 알아두자.

 

구현 시 참고할 점 및 코드

• $S_f(n)$을 메모이제이션해야 하는데, $n = \left\lfloor \frac{N}{d} \right\rfloor$이 $N$까지 가능하기에 단순히 일반 배열에 저장하는 건 안 되고, C++의 std::map(혹은 hashmap) 등의 자료구조를 사용해야 한다. (사실 배열만으로도 해결할 수 있으나, 자료구조를 이용하는 편이 코드가 간단하다.)
• 그런데 std::map을 사용하면 시간복잡도에 로그가 추가로 붙고 상수도 상당히 크므로 웬만하면 해시맵을 이용하도록 하자(C++에 std::unordered_map으로 구현되어 있다). 어차피 인덱스로 가능한 $n$의 값은 $2\lfloor \sqrt N \rfloor$개로 정해져 있으므로, 해시맵 코드를 터뜨리기로 유명한 코드포스 등에서도 자유롭게 이용할 수 있다!
• 위에서 최적의 전처리 범위 $K$는 $(2N)^{2/3}$이라고 했지만 애초에 시간복잡도를 대략적으로만 분석한 데다가, 구현을 어떻게 하느냐에 따라 실제 최적의 $K$값이 달라질 수 있다. 그러므로 문제를 실제로 풀 때에는, $K$를 $N^{2/3}$ 근방에서 조금씩 변경해 가며 여러 번 제출해서 최적의 $K$를 찾는 걸 추천한다.
• xudyh's sieve가 굉장히 놀라운 알고리즘이긴 하나, 함수 $g$의 조건이 빡센 탓에 응용 가능한 것이 사실상 몇 개 없는 관계로($\phi(n)$의 부분합이나 메르텐스 함수 $M(n)$ 정도?) 폭넓게 쓰기는 힘들다. 이러한 단점을 해결한 min_25's sieve라는 알고리즘이 있다. 이에 대해 관심이 있다면 이 글이나(일본어로 작성되어 있다) rkm0959 님의 블로그 글을 참고하라.

  지금까지의 논의를 본 문제에 적용해 보면, ${\rm Id} = \phi \ast 1$이고 $S_{{\rm Id}}(N) = \frac{N(N+1)}{2}$, $S_1(N) = N$이므로 식 (1)은

$$\begin{aligned} S_{\phi}(N) &= \frac{N(N+1)}{2} - \sum_{d=2}^N S_{\phi} {\left(\left\lfloor \frac{N}{d} \right\rfloor \right)} \\ &= \frac{N(N+1)}{2} - \sum_{i=2}^t S_{\phi}(a_i) (e_i - s_i + 1) \end{aligned}$$

가 된다! 「Double Counting과 Harmonic Lemma」 게시글을 참고하여 이를 코드로 작성하면 다음과 같다.

#include <cstdio>
#include <unordered_map>
#define lint long long
#define MAX 500003

lint phi[MAX];
std::unordered_map<int, lint> mp;

lint S(int x) {
	int i, j;
	if (x < MAX) return phi[x];
	if (mp.find(x) != mp.end()) return mp[x];
	lint ret = x * (x + 1LL) / 2;
	for (i = 2; i <= x; i = j+1) {
		j = x / (x / i);
		ret -= (j - i + 1) * S(x / i);
	}
	return mp[x] = ret;
}

int main()
{
	int N, i, j;
	scanf("%d", &N);

	for (i = 1; i < MAX; ++i) {
		phi[i] += i;
		for (j = 2; j*i < MAX; ++j) phi[j*i] -= phi[i];
		phi[i] += phi[i-1];
	}

	printf("%lld", S(N));
	return 0;
}

  위 코드에서 i가 $s_i$, j가 $e_i$, MAX가 $K$ 역할을 하며, 전처리로 오일러 토션트 함수의 값과 그 부분합을 구하는 것을 하나의 반복문으로 처리한 것에 유의하라. 이 코드에 $N$으로 $10^9$를 입력하면 ideone.com 기준 0.1s 정도에 돌아가는 걸 확인할 수 있다. 이로써 우리는, 멀고 먼 길을 거쳐 본 문제를 해결했다!

 

  참고로 위 코드에서 MAX를 $0$으로 바꾸면 전처리 없는 $\mathcal O(N^{3/4})$짜리 xudyh's sieve가 되는데, 이를 ideone.com에서 돌려보면 1.2초 정도의 시간이 걸린다. 전처리 유무만으로 실행 시간이 10배 넘게 차이나는 것이 보이는가? (^^)

 

메르텐스 함수 구하기

  xudyh's sieve의 다른 이름이 Mertens Trick인 만큼, xudyh's sieve로 $M(n)$의 값도 구해 보자. $\mu \ast 1 = \varepsilon$이고 $S_1(N) = N$, $S_{\varepsilon}(N) = 1$이므로 이를 식 (1)에 대입하면

$$\begin{aligned} M(N) &= 1 - \sum_{d=2}^N M {\left(\left\lfloor \frac{N}{d} \right\rfloor \right)} \\ &= 1 - \sum_{i=2}^t M(a_i) (e_i - s_i + 1) \end{aligned}$$

가 되고, 이 식을 바탕으로 코드를 작성하면 다음과 같다.

#include <cstdio>
#include <unordered_map>
#define lint long long
#define MAX 300003

lint mu[MAX] = {0, 1}; // 편의상 mu[0] = 0이라 둔다.
std::unordered_map<int, lint> mp;

lint S(int x) {
	int i, j;
	if (x < MAX) return mu[x];
	if (mp.find(x) != mp.end()) return mp[x];
	lint ret = 1;
	for (i = 2; i <= x; i = j+1) {
		j = x / (x / i);
		ret -= (j - i + 1) * S(x / i);
	}
	return mp[x] = ret;
}

int main()
{
	int N, i, j;
	scanf("%d", &N);

	for (i = 1; i < MAX; ++i) {
		for (j = 2; j*i < MAX; ++j) mu[j*i] -= mu[i];
		mu[i] += mu[i-1];
	}

	printf("%lld", S(N));
	return 0;
}

이 코드 역시 ideone.com에서 입력으로 $10^9$을 넣고 돌려보면 0.11s만에 결과가 나온다!

 

정말 놀랍지 않은가?

 

참고 사이트

1. https://codeforces.com/blog/entry/54150 (코드포스 블로그 글 ─ 영어로 작성되어 있다)

2. https://rkm0959.tistory.com/190 (rkm0959 님 블로그)

3. https://gratus907.com/108 (gratus907 님 블로그)

4. https://www.acmicpc.net/workbook/view/2208 (BOJ Mertens Trick 문제집)