메르텐스 함수 (Mertens Function) 와 그 빠른 계산 (Xudyh Sieve)

메르텐스 함수 (Mertens Function) 은 뫼비우스 함수 (Mobius Function) 으로부터 정의되는, 매우 재밌는 성질을 갖는 함수이다.

1. Mobius Function / Mertens Function의 정의

먼저, 간단히 정의를 살펴보자.

 

정의. Mobius Function
다음과 같은 뫼비우스 함수를 정의한다. $$\mu(n) =\begin{cases} 0, & ^\exists d \in \mathbb{N}, d > 1 \text{ such that } d^2 | n\\ (-1)^k, & n = p_1 p_2 p_3 \dots p_k \end{cases}$$

https://gratus907.com/43?category=869407

한참 예전에 이런 포스팅을 했었는데, 여기에서 뫼비우스 함수를 소개한 적이 있다. Mobius function의 가장 중요한 성질 중 하나는, Dirichlet convolution에서 상수함수 $1(n)$ 의 Inverse라는 특징이다. 저 포스팅에서 이미 이러한 함수가 $1$의 inverse임을 보인 적이 있다 (Dirichlet convolution의 inverse가 유일함은 보이지 않았는데, 대수적인 지식을 활용하면 쉽게 보일 수 있다. 그렇지 않아도 보이기는 어렵지 않은 것 같긴 하다)

 

아무튼, Mertens Function은 이러한 Mobius Function의 부분합으로 정의된 함수이다.

 

정의. Mertens Function
다음과 같은 뫼비우스 함수를 정의한다. $$M(n) = \sum_{i = 1}^{n} \mu(i)$$

2. Asymptotics / Mertens Conjecture

이 함수에 관한 Merten's Conjecture라는 추측이 있다. 엄밀히는, 원래의 mertens conjecture는 이미 틀렸다는 것이 확인되었으나, 조건을 조금 약화한 mertens conjecture가 옳은지는 아직 알려지지 않았다. 

메르텐스 함수의 대략적인 Asymptotic에 대해 알아보자. 먼저, 메르텐스 함수는 1, -1, 0들로 구성된 뫼비우스 함수의 부분합으로 정의되므로, $\abs{M(n)} \leq n$ 임이 자명하다. 실제로 많은 값을 찍어보면 함수가 거의 0 근처에서만 움직임을 알 수 있다. 

출처 : https://en.wikipedia.org/wiki/Mertens_function#/media/File:Mertens.svg

$n$이 1만까지 갈 때 $M(n)$의 절댓값이 50도 벗어나지 않음을 관찰할 수 있다.

이에 대해, 메르텐스의 본래 추측은 다음과 같다.

[Disproved] Merten's Conjecture
$$\abs{M(n)} \leq \sqrt{n}$$

1985년에 이 추측은 (수학적으로) 부정되었으나, 아직 실제로 반례를 찾은 결과는 없다. 이 추측이 오랜 시간 증명되지 않은 것에서 알 수 있듯, 반례를 찾기가 굉장히 빡센데, 현재 알려진 최고의 결과가 반례가 $10^{6.91 \times 10^{64}}$ 이하에 존재한다는 정도만 알려져 있다. (...)

 

이 추측이 특히 중요한 이유 중 하나는, 다음의 (아주 약간 약화된) 추측 때문이다.

[Unsolved] Merten's weaker Conjecture
$$\abs{M(n)} \in O(x^{0.5 + \epsilon})$$ 가 모든 $\epsilon > 0$ 에 대하여 성립하는가?

이 추측이 리만 가설과 동치임이 알려져 있으나, 리만가설이 증명되지 않은 것에서 알 수 있듯이 이것도 전혀 진전이 없다.

 

3. Xudyh's Sieve

참고 : https://codeforces.com/blog/entry/54150 에 설명되어 있고, 이 부분의 몇 문단은 거의 그 글의 번역이다. 참고할만한 자료들의 링크를 맨 끝에 모아 두었다.

본론으로 들어가서, 이러한 함수를 빠르게 계산하는 방법인 Xudyh's Sieve를 알아보자. 대략 이정도면 대회를 위해 공부하는 알고리즘과 그렇지 않은 알고리즘의 경계선 쯤에 있는데, 이 알고리즘이 알려진 대회에 나온 것을 본건 Opencup GP of China에 한번 나온적 있다고 하는데, Opencup + China면 정상적인 대회가 아닌게 분명하므로...

Project Euler Forum 같은데서 Xudyh's Sieve라는 말이 나온거 같은데, 사실 누구 이름인지 뭔지는 잘 모르겠다. Dirichlet convolution trick for partial sums 이라고 부르는 사람도 본 것 같다. 

 

먼저 https://gratus907.com/43?category=869407 에도 작성한 적 있는, Dirichlet Convolution 정의를 적어 놓자.

정의. Dirichlet Convolution
수론적 함수 (대략 수열이라는 뜻이다) $f, g$ 에 대해 다음과 같이 디리클레 합성곱을 정의한다. $$(f*g)(n) = \sum_{d | n}f(d)g(\frac{n}{d})$$

또한, '곱셈적 함수'(Multiplicative Function) 라는 것을 정의할 필요가 있다. 

정의. 곱셈적 함수 (Multiplicative Function)

수론적 함수 $f$가 서로소인 정수 $a, b$에 대해, $f(ab) = f(a) f(b)$ 를 만족하면 이를 곱셈적 함수라고 정의한다.

Notation으로, $f$의 부분합을 $s_f$ 라 하자. $g$ 같은 다른 함수도 마찬가지.

 

우리는 다음과 같은 목표를 가지고 있다.

- 어떤 함수 $f$의 부분합을 빠르게 계산하고 싶은데, 그 부분합을 빠르게 구하기 어렵다. 구체적으로 $s_f(n)$ 을 Linear보다 빠르게 계산할 방법이 마땅치 않다.

- 다만, 어느 정도 (작은 값 정도) 의 부분합은 미리 계산해 놓을만 하다. 

- 그런데, 이 함수 $f$가 곱셈적 함수임을 안다. 

편하게 생각해서, $s_f$의 linear 계산이 가능할 때, sublinear를 목표로 알고리즘을 발전시킨다고 생각하면 된다.

(실제로는 Quasilinear, 즉 $O(n \log^k{n})$ 정도를 말한다)

 

어떤 함수 $f$에 대해서는, 다음과 같은 함수 $g$를 찾을 수 있다.

(당연히 모든 함수에 대해 되는것일 리 없고, 극히 일부의 함수 $f$에 대해서만 가능하다.)

- $g$의 부분합을 매우 빠르게 계산할 수 있다.

- $f * g$의 부분합을 매우 빠르게 계산할 수 있다.

매우 빠르게라고 함은 $O(1)$ 같은 시간을 의미하는데, 원리적으로는 Complexity가 약간 달라지지만 Poly-Logarithmic 하면 된다. 이 $s_g$와 $s_{f*g}$ 가 polylogarithmic 계산이 가능하면, $s_f$를 Sublinear하게 계산할 수 있음을 보이려고 한다.$f$는 계산하기 힘든데 $f * g$를 계산하는게 쉬운 함수가 있나 싶을 수 있지만, 드물게 찾을 수 있다.

 

아무튼 이런 함수가 잘 있다고 하고, 우리는 다음과 같은 합을 계산하기로 하자.

$$s_{f * g}(n) = \sum_{i = 1}^{n} \sum_{u | i} g(u)f\left(\frac{i}{u}\right)$$

이제 $i = uv$ 로 쓰면, 다음과 같이 간단히 쓸 수 있다.

$$s_{f * g}(n) = \sum_{u = 1}^{n} \sum_{v = 1}^{[n/u]} f(v)g(u)$$

$u = 1$ 인 부분만 따로 계산하기로 하자. 식을 잘라서 정리하면, $f$의 합 부분을 전부 $s_f$ 로 쓸 수 있다. 

$$s_{f*g}(n) = g(1)s_f(n) + \sum_{u = 2}^{n}g(u) s_f\left(\left\lfloor{\frac{n}{u}}\right\rfloor\right) $$

$s_f (n)$ 만 남기면, 다음의 식을 얻게 된다.

$$s_f(n) = \frac{1}{g(1)} (s_{f*g}(n) - \sum_{u = 2}^{n}g(u) s_f\left(\left\lfloor{\frac{n}{u}}\right\rfloor\right) ) $$

나머지 부분들을 쉽게 계산할 수 있다는 것은 이미 가정한 사실이거나 자명하므로, $\sum_{u = 2}^{n}g(u) s_f\left(\left\lfloor{\frac{n}{u}}\right\rfloor\right)$ 를 계산하는 시간이 전체 계산 시간을 지배하게 된다. 이를 어떻게 빨리 계산할지만 고민해 보면 되겠다. 우선, 여기까지 정리한 것으로 계산해 보자.

여기서 이 식을 빠르게 계산하기 위해 필요한 Lemma가 2개 있다. 

Lemma 1
임의의 자연수 $p, q, r$ 에 대하여 다음이 성립한다. $$ \left\lfloor \frac{ \left\lfloor \frac{p}{q}\right\rfloor}{r} \right\rfloor = \left\lfloor \frac{p}{qr}\right\rfloor $$

즉, Integer division을 반복해서 수행할 때 편하게 해도 된다는 의미이다.

Lemma 2
다음과 같은 Harmonic하게 정의된 집합을 생각하자. $$\Setcond{\left\lfloor \frac{n}{x}\right\rfloor}{x \leq n}$$ 이 집합의 원소의 개수는 많아야 $2 \sqrt n$ 개이다.

n/2, n/3, n/4, ... 들을 다 모아봐야 $O(\sqrt n)$ 개밖에 없다는 강력한 명제로, 증명은 대충... $x < \sqrt n$ 인 $x$가 $\sqrt n$ 개밖에 없고, $x > \sqrt n$ 이면 반대로 $\left\lfloor \frac{n}{x} \right\rfloor$ 가 $\sqrt n$ 보다 작으므로 역시 $\sqrt n$ 개밖에 없다.

 

결국 우리가 위 $s_f$ 의 $n/x$ 값들을 계산하는 데 있어서, Lemma 1과 2를 적용하면 계산해야 할 서로 다른 값은 $2 \sqrt n$ 개밖에 없고, 이 값들을 모두 알 때 합치는 시간이 $O(\sqrt n)$ 만큼 들기 때문에, 다음과 같은 식으로 전체 Complexity를 계산할 수 있다. 이 Complexity를 확실히 하기 위해서는, 반드시 Lemma 2의 Harmonic set을 미리 전부 찾은 다음, 작은 값부터 순서대로 메모이제이션하면서 계산해야 한다. 

$$T(n) = O(\sqrt {n}) + \sum_{i = 1}^{\sqrt{n}} \sqrt{i} + \sum_{i = 1}^{\sqrt{n}} \sqrt{\left\lfloor \frac{n}{i}\right\rfloor}$$

이 식을 적분으로 열심히 전개해 보면, $T(n) = O(n^{\frac{3}{4}})$ 를 얻는다.

 

이미 Sub-linear 알고리즘을 얻긴 했지만, 이를 더 개선할 방법이 없는지 찾아 보자. $s_f$를 계산하기가 어렵다고 하긴 했지만, 일단 Linear 비슷 하게는 어떻게든 계산할 수 있다. 그러면, 미리 $K$를 하나 잡아서 거기까지는 Precomputation을 해놓기로 하자. $K \geq \sqrt{n}$ 을 잡는다면, 다음과 같이 복잡도 식을 정리할 수 있겠다.

$$T(n) = O(\sqrt {n}) + O(k) + \sum_{i = 1}^{\frac{n}{k}} \sqrt{\frac{n}{i}}$$ (어차피 Asymptotic Behavior에 대한 논의이므로 floor는 이제 별로 중요하지 않다. 앞에서 floor를 계속 명시한 이유는 그 값이 몇개 없다는 중요한 관찰이 있기 때문이다)

이는 다시,

$$O(k) + O\left(\frac{n}{\sqrt{k}} \right)$$ 이러한 복잡도를 얻게 된다. 이 식이 Asymptotically 최소화되는 상수 $K$를 찾으면 $K = n^{\frac{2}{3}}$ 을 택하여, 전체 알고리즘을 $O(n^{\frac{2}{3}})$ 으로 줄일 수 있다.

 

* 참고

- $K$를 편하게 그렇게 잡았지만, 실제로는 Quasilinear 하게 $s_f$를 계산해야 하므로 약간 달라질 수도 있다. 메르텐스 함수의 경우, $K = n^{\frac{2}{3}} (\log n) ^ {\frac{1}{3}}$ 같은 값을 택하는 것이 최적이고, 등등... 이와 같이, 실제로 $K$를 어떻게 잡을지는 $s_f$를 계산하는 시간에 따라 달려 있다. 다만 Quasilinear $s_f$ 계산만 가능했던 것이, $O(n^{\frac{2}{3}})$ 에 $\log$ 항이 몇개 붙은, Quasi-$n^{2/3}$ 시간에 계산된다는 점은 유지된다.

- $g$를 찾는 일반적인 방법론은 딱히 없는 것 같다. Dirichlet convolution에 관한 몇가지 특수함수들은 기억하면 좋고, 아니면 어차피 $s_g$를 매우 빨리 계산할 수 있는 함수가 몇개 없음을 기억하자. 기껏해야 $Id$, $\epsilon$, $1$ 정도나 될 것 같다. $id$의 partial sum은 $\frac{n(n+1)}{2}$, $1$은 $n$, $\epsilon$은 $1$ 임이 자명.

- 잘은 모르겠지만, 메르텐스 함수처럼 매우 특수한 함수들은 Analytic Number Theory의 여러 이론들을 이용해서 복잡도를 더 ($O(\sqrt n)$ 까지 내릴 수 있다고 한다. 그런 논문에서는 위 접근을 'Combinatorial' 하다고 써놓고 있다.

* Reference

[1] https://codeforces.com/blog/entry/54150 

[2] Deléglise, Marc; Rivat, Joöl. Computing the summation of the Möbius function. (1996)

 

3-2. Mertens Function에의 적용 / 구현

위 트릭을 잘 이용해서 실제로 메르텐스 함수를 구하는 것은 매우 쉽다. 메르텐스 함수는 $\mu$의 부분합이므로 $f = \mu$ 를 택하고, $g$로는 constant function $1$ 을 택하자. 뫼비우스 함수를 $1$의 역함수로 정의했으므로 $\mu * 1 = \epsilon$ 이고, 둘 모두 부분합을 $O(1)$ 에 구할 수 있는 좋은 함수이다. 이제 이를 실제로 구현하면 된다.

------

구현 코드 추가 예정