rkm0959

무난한 문제. \(f_Y(i,j)\)의 값은 구간의 길이인 \(j-i+1\)에만 영향을 받는다. 

또한, 한 구간의 gcd가 1인 것은 구간에 속한 모든 자연수가 \(p\)의 배수가 되는 소수 \(p|Y\)가 존재하지 않는 것과 동치이다. 

한편, 배열에 들어갈 수 있는 값은 \([1,Y]\) 중 하나이므로, 각 소수를 독립적으로 생각할 수 있다.

여기까지 관찰하면 남는 것은 정말 간단한 단순 계산 문제. 대충 \(\mathcal{O}(N \log NY + \sqrt{Y})\) 정도에 풀린다.

이러한 문제들은 대부분 유클리드 호제법과 비슷한 느낌으로 풀리는 경우가 많다. 

문제를 주어진 \((p,q)\)에 대해 \(\sum_{i=1}^n \lfloor \frac{pi}{q} \rfloor\)를 계산하는 문제로 변환하자.

쉽게 관찰할 수 있는 것은, \((p,q)\)에 대한 문제는 \((p\text{%}q, q)\)에 대한 문제로 간단하게 변환이 가능하다는 것이다.

우리가 \((p,q)\)에 대한 문제를 빠르게 \((q,p)\)에 대한 문제로 변환할 수 있다면, 유클리드 호제법 과정을 그대로 거쳐서 로그 시간에 문제를 해결할 수 있을 것이다. (분모, 분자가 \(10^{18}\) scale인 분수를 double, int128 없이 비교하는 방식도 비슷하다)

이를 위한 대표적인 방법은 문제를 좌표평면에 \(y = \frac{p}{q}x\)를 그리고, \(1 \le x \le n\), \(1 \le y \le \frac{p}{q}x\)를 만족하는 정수점의 개수를 세는 문제로 변환하는 것이다. 이제 역으로 \( x = \frac{q}{p}y\)를 생각하고, \(1 \le y \le \lfloor \frac{pn}{q} \rfloor \), \( 1 \le x \le \frac{q}{p}y\)인 점들을 생각하자. 

두 그룹의 모두 속하는 정점은 \(y = \frac{p}{q}x\) 위에 놓이는 점들로, 이들의 개수는 쉽게 계산된다. 

이 관찰을 모두 하면 \((p,q)\)에서 \((q,p)\)로 문제를 전환하는 것은 어렵지 않다. 여기서부터는 단순 계산.

\(\displaystyle \sum_{1 \le u < v \le n} \frac{uv}{\text{gcd}(u,v)} \)를 계산해야 한다. \(g=\text{gcd}(u,v)\)를 고정하자.

그러면 식은 \(\displaystyle \sum_{g=1}^n \sum_{1 \le u' < v' \le \lfloor \frac{n}{g} \rfloor, \text{gcd}(u',v')=1} gu'v'\)이다.  \(\displaystyle F(x) = \sum_{1 \le u' < v' \le x, \text{gcd}(u',v')=1} u'v' \)이라 하자.

\(v'\)을 고정하면서 생각하면 \(\displaystyle F(x) = \sum_{i=2}^x \frac{i^2 \phi(i)}{2}\)임을 알 수 있다. 남은 것은 \(\phi(i)\)들을 전처리하고 단순 계산을 하는 것 뿐.

물론, 여기까지 관찰하면 xudyh's sieve를 끼얹어서 문제를 \(\mathcal{O}(n^{2/3})\)에 풀 수 있다. 

BOJ 13358, BOJ 12445, BOJ 12446, BOJ 16214 등 power tower를 다루는 모든 문제에서 사용되는 정리가 필요하다.

예전 블로그에 그 정리에 대한 설명을 써놓았기 때문에, 링크를 걸어놓는 것으로 설명을 대체한다. 마지막 문제다.

자연수가 prime prime power이면, prime prime power로의 표현 방법은 유일하다. 

또한, 지수가 3 이상인 prime prime power는 모두 직접 구할 수 있을 정도로 적으니, 모두 구해놓자.

문제는 지수가 2인 prime prime power이다. 여기서 아이디어는 prime gap이 작다는 것을 이용하는 것이다.

\(n \le 10^9\)이면, \(n\)과 \(n+3 \cdot 10^6\) 사이에는 \(100000\)개 이상의 소수가 있다는 믿음을 가지자.

이러면 남은 문제는 쉽다. \(\lfloor \sqrt{n} \rfloor\)부터 \(\lfloor \sqrt{n} \rfloor + 3 \cdot 10^6\) 사이에 존재하는 수들이 소수인지를 빠르게 확인한다. 

에라토스테네스의 체를 응용하면 어렵지 않게 할 수 있다. 이제 이분탐색을 끼얹으면 빠른 시간안에 문제를 해결할 수 있다. 

먼저 \(x \in [1, n]\), \(y \in [1, m]\)인 \((x,y)\)에 대해서 \(\text{lcm}(x,y)\)를 다 더한 값을 빠르게 구하는 방법을 생각해보자.

명제 \(P\)에 대해서 \([P]\)를 \(P\)가 참이면 1, 아니면 0인 값으로 정의하도록 하자. 

\(\text{gcd}(x,y)=g\)를 고정하면, 계산해야 하는 식은 \(\displaystyle \sum_{g=1}^{\text{min}(n,m)} \sum_{x=1}^{\lfloor \frac{n}{g} \rfloor} \sum_{y=1}^{\lfloor \frac{m}{g} \rfloor} gxy \cdot [\text{gcd}(x,y) =1] \)이다. 

그런데 \([n=1] = \sum_{d|n} \mu(d)\)임은 정수론 문제에서 아주 자주 쓰이는 식이다. 이를 사용하면, $$ \sum_{g=1}^{\text{min}(n,m)} \sum_{x=1}^{\lfloor \frac{n}{g} \rfloor} \sum_{y=1}^{\lfloor \frac{m}{g} \rfloor} gxy \cdot \sum_{d|\text{gcd}(x,y)} \mu(d) $$를 계산하면 충분함을 파악할 수 있다. 식을 \(d\)에 대해서 다시 쓰면 주어진 식은 $$ \sum_{d=1}^{\text{min}(n,m)} \mu(d) \cdot \sum_{g=1}^{\text{min}(n,m)} g \cdot \left( \sum_{x=1}^{\lfloor \frac{n}{gd} \rfloor} d \cdot x \right) \cdot \left( \sum_{y=1}^{\lfloor \frac{m}{gd} \rfloor} d \cdot y \right) = \sum_{d=1}^{\text{min}(n,m)} \mu(d) \cdot \sum_{g=1}^{\text{min}(n,m)} g \cdot \frac{d^2 \cdot \lfloor \frac{n}{gd} \rfloor \left( \lfloor \frac{n}{gd} \rfloor + 1 \right) \cdot \lfloor \frac{m}{gd} \rfloor \left( \lfloor \frac{m}{gd} \rfloor + 1 \right)}{4} $$로 변환된다. 여기서 \(k=gd\)라고 두고, \(k\)에 대해서 식을 정리하면, 주어진 식은 다시 $$ \sum_{k=1}^{\text{min}(n,m)} \frac{\lfloor \frac{n}{k} \rfloor \left(\lfloor \frac{n}{k} \rfloor + 1 \right) \cdot \lfloor \frac{m}{k} \rfloor \left( \lfloor \frac{m}{k} \rfloor + 1 \right)}{4} \cdot \left( \sum_{d|k} k \cdot d \cdot \mu(d) \right) $$와 같음을 알 수 있다. \(f(k) = \sum_{d|k} d \cdot \mu(d)\)라 하면, 간단하게 식을 정리해 $$ \sum_{k=1}^{\text{min}(n,m)} \frac{\lfloor \frac{n}{k} \rfloor \left(\lfloor \frac{n}{k} \rfloor + 1 \right) \cdot \lfloor \frac{m}{k} \rfloor \left( \lfloor \frac{m}{k} \rfloor + 1 \right)}{4} \cdot kf(k)$$를 얻을 수 있다. 이제 준비가 끝났다. \(\lfloor \frac{n}{k} \rfloor\)가 가질 수 있는 값의 개수는 \(\mathcal{O}(\sqrt{n})\)이다. 

그러므로 \(f(k)\)를 먼저 모두 계산한 뒤, \(kf(k)\)의 부분합을 전처리하면 위 식을 \(\mathcal{O}(\sqrt{n}+\sqrt{m})\)에 계산할 수 있다.

이제 \(\text{gcd}(x,y)\)가 squarefree한 경우만을 계산하는 방법을 고민하자. \(g(n,m)\)을 \(\sum_{x=1}^n \sum_{y=1}^m \text{lcm}(x,y)\)라고 하자. 

포함-배제의 원리에서, 우리가 구해야 하는 것은 \(\displaystyle \sum_{t=1}^{\sqrt{\text{min}(n,m)}} \mu(t) \cdot t^2 \cdot g\left(\left\lfloor \frac{n}{t^2} \right\rfloor, \left\lfloor \frac{m}{t^2} \right\rfloor\right)\)이다.

\(g(n,m)\)을 \(\mathcal{O}(\sqrt{n}+\sqrt{m})\) 시간에 구할 수 있으니, 저 식도 대략 \(\mathcal{O}((\sqrt{n}+\sqrt{m}) \log \text{min}(n,m))\)에 구할 수 있다.

전처리 과정이 필요한 부분은 \(f(k)\) 값들과 \(\mu\) 값들이 전부이며, 모두 간단하게 전처리할 수 있다.

전체 시간복잡도는 대략 \(\mathcal{O}(MAX \log MAX + TC \sqrt{MAX} \log MAX)\)이다. 

결국에는 \(\displaystyle \sum_{p \le P, p \in \text{prime}} p^Q\)를 빠르게 계산하는 문제이다. 개인적으로 아직 공부가 덜 된 유형이다 ㅠㅠㅠ

기본적으로 sum of powers를, 즉 \(1^Q+2^Q+ \cdots + N^Q\)를 빠르게 구하는 과정이 필요하다.

개인적으로는 베르누이 수를 많이 사용하는데, 문제는 베르누이 수가 유리수라서 \(\pmod P\) 연산을 하기가 까다롭다.

실제로는 \(Q \le 1000\)이므로, \(P > 5000\)이면 사용하는 베르누이 수들의 분모에 \(P\)가 없게 되어 문제가 없다. 

그러니까 \(P<5000\)이면 naive하게 처리하고, \(P > 5000\)인 경우를 보도록 하자.

베르누이 수를 전처리하면, \(N\)이 주어졌을 때 \(1^Q+2^Q+ \cdots +N^Q\)를 \(\mathcal{O}(Q)\)에 구할 수 있다. 참고

에라토스테네스의 체에서 숫자들을 지워가는 과정을 모방하는 방식으로 문제에 접근할 수 있다.

현재 에라토스테네스의 체에서 현재 소수 \(r\)을 보고 있다고 하자. 

아직 지워지지 않은 수들은 \(r\) 미만의 소수이거나, 아니면 최소 소인수가 \(r\) 이상인 수들이다. 

이제 \(r\)이 아닌 \(r\)의 배수들을 모두 지운다. 새롭게 지워지는 수들은 어떤 형태를 가지고 있을까?

새롭게 지워지는 수들은 \(r \cdot k\) 형태의 수로, \(k\)는 역시 최소 소인수가 \(r\) 이상인 수이다. 

즉, \(r \cdot k\)가 새롭게 지워지는 것은 [\(k\)가 아직 지워지지 않았고, \(r\) 미만의 소수가 아님]과 동치이다.

 \(dp_N\)을 \(N\) 이하의 수 중 아직 지워지지 않은 수들의 \(Q\)제곱의 합이라고 하자.

위 과정을 잘 생각해보면, \(dp_{\lfloor \frac{P}{i} \rfloor}\) 형태의 \(dp\) 값들만 생각해도 모든 계산이 가능함을 알 수 있다. 

물론, \(\lfloor \frac{P}{i} \rfloor\)로 가능한 수의 개수는 \(\mathcal{O}(\sqrt{P})\)이다. \(dp_k\)의 초깃값을 \(2^Q+3^Q+ \cdots +k^Q\)로 초기화한 후, 위에서 언급한 접근을 활용하여 \(dp\)의 값들을 변경해주자. 잘 구현하면, \(\mathcal{O}(P^{3/4} \log Q + \sqrt{P}Q+Q^2)\) 정도로 문제가 풀린다.

굉장히 많은 부분을 (코드 포함) 이 글에서 배웠다. 코드도 사실 저 글에서 거의 복붙한 수준이다.

pps789님의 코드가 굉장히 잘 이해되고 굉장히 빨라서, 그 코드를 사용해서 이 내용을 다시 공부할 생각이다.

핵심 아이디어는 구하고자 하는 식이 \(\displaystyle \sum_{(i,j)=(j,k)=(k,i)=1} \lfloor \frac{a}{i} \rfloor \cdot \lfloor \frac{b}{j} \rfloor \cdot \lfloor \frac{c}{k} \rfloor\)로 변형된다는 사실이다.

앞에서 정의한 것처럼, 명제 \(P\)에 대해서 \([P]\)를 \(P\)가 참이면 1, 아니면 0으로 정의하자.

이제 문제는 굉장히 간단해진다. 먼저 각 \(1 \le N \le ab\)에 대해서 \(\sum_{i=1}^c \lfloor \frac{c}{i} \rfloor \cdot [\text{gcd}(i,N)=1]\)을 전처리하자. 

이 식은 \(\displaystyle \sum_{i=1}^c \lfloor \frac{c}{i} \rfloor \cdot \sum_{d|\text{gcd}(i,N)} \mu(d) = \sum_{d|N} \mu(d) \sum_{i=1}^{\lfloor \frac{c}{d} \rfloor} \lfloor \frac{c}{di} \rfloor\)으로, \(\mu\) 함수를 전처리한 뒤 쉽게 계산할 수 있다.

전처리한 이 식을 \(f(N)\)이라 하자. 구해야 하는 값은 \(\displaystyle \sum_{i=1}^a \sum_{j=1}^b [\text{gcd}(i,j)=1] \cdot \lfloor \frac{a}{i} \rfloor \cdot \lfloor \frac{b}{j} \rfloor \cdot f(i \cdot j)\)다.

\(f\)가 전처리가 되었다면, 이 식은 그냥 naive하게 구할 수 있다. 단순 계산. 

이제 우리가 구하고자 하는 식이 \(\displaystyle \sum_{(i,j)=(j,k)=(k,i)=1} \lfloor \frac{a}{i} \rfloor \cdot \lfloor \frac{b}{j} \rfloor \cdot \lfloor \frac{c}{k} \rfloor\)과 같음을 증명하자.

\(\displaystyle f(a,b,c)=\sum_{i=1}^a \sum_{j=1}^b \sum_{k=1}^c \tau(i \cdot j \cdot k)\)라 하고, \(\displaystyle g(a,b,c)=\sum_{(i,j)=(j,k)=(k,i)=1} \lfloor \frac{a}{i} \rfloor \cdot \lfloor \frac{b}{j} \rfloor \cdot \lfloor \frac{c}{k} \rfloor\)라 하자.

\(f(a,b,c)=g(a,b,c)\)임을 \(a+b+c\)에 대한 강한 귀납법으로 증명한다. Base Case는 생략한다. 우리는 두 식

\(f(a,b,c)-f(a-1,b,c)-f(a,b-1,c)-f(a,b,c-1)+f(a,b-1,c-1)+f(a-1,b,c-1)+f(a-1,b-1,c)-f(a-1,b-1,c-1)\)

\(g(a,b,c)-g(a-1,b,c)-g(a,b-1,c)-g(a,b,c-1)+g(a,b-1,c-1)+g(a-1,b,c-1)+g(a-1,b-1,c)-g(a-1,b-1,c-1)\)이 같음을 보이면 충분하다. 물론, 첫 번째 식은 단순하게 \(\tau(abc)\)와 같다. 

두 번째 식은 \(\displaystyle \sum_{(i,j)=(j,k)=(k,j)=1} \left( \lfloor \frac{a}{i} \rfloor - \lfloor \frac{a-1}{i} \rfloor \right) \left( \lfloor \frac{b}{j} \rfloor - \lfloor \frac{b-1}{j} \rfloor \right) \left(\lfloor \frac{c}{k} \rfloor - \lfloor \frac{c-1}{k} \rfloor \right)\)이다.

그런데 자연수 \(n, k\)에 대해서 \(\lfloor \frac{n}{k} \rfloor - \lfloor \frac{n-1}{k} \rfloor\)은 \(n\)이 \(k\)의 배수일 때만 \(1\)이고 나머지 경우에는 \(0\)이다.

즉, 두 번째 식은 \((i,j)=(j,k)=(k,i)=1\)이고 \(i,j,k\)가 각각 \(a,b,c\)의 약수인 경우의 수와 같다.

첫 번째 식과 두 번째 식의 의미를 생각해보면, 결국 각 소인수에 대해서 독립적으로 생각해도 됨을 알 수 있다.

소수 \(p\)를 고정하자. \(a, b, c\)가 각각 \(p\)를 \(u, v, w\)개 가진다고 하자. 물론 \(u, v, w \ge 0\)이다.

첫 번째 식에서 소수 \(p\)에 의해 곱해지는 값은 \(abc\)가 \(p\)를 \(u+v+w\)개 가지므로 \(u+v+w+1\)과 같다. 

두 번째 식에서 소수 \(p\)에 의해 곱해지는 값도 \(u+v+w+1\)과 같다. \(0 \le i' \le u\), \(0 \le j' \le v\), \(0 \le k' \le w\)인 \(i',j',k'\)를 고르되, (물론, 이 값들은 \(i,j,k\)가 갖는 \(p\)의 개수가 된다) \(i', j', k'\) 중 최대 하나가 양수인 경우의 수를 세면 정확히 \(u+v+w+1\)이기 때문이다. 결국 첫 번째 식과 두 번째 식은 정확히 같다. 증명 끝. 남은 것은 크게 복잡하지는 않은 구현.

주어진 \(n^2\)개의 자연수 중 가장 큰 자연수를 찾으면, 그 수는 \(a_1, a_2, \cdots a_n\) 중 하나여야 함을 쉽게 증명할 수 있다. 이제 그 수를 \(a_1\)이라 하면, \(n^2\)개의 자연수 중 하나는 \(\text{gcd}(a_1, a_1) = a_1\)이어야 한다. 주어진 \(n^2\)개의 자연수에서 \(a_1\) 하나를 제거해주자. 남은 \(n^2-1\)개의 자연수 중에서 가장 큰 자연수를 다시 찾으면, 이 수는 다시 \(a_2, a_3, \cdots a_n\) 중 하나여야 한다. 이 수를 \(a_2\)라 하면, \(n^2-1\)개의 자연수 중에서 하나는 \(a_2\)여야 하고, 두 개는 \(\text{gcd}(a_1, a_2)\)여야 한다. 이 세 개의 수를 목록에서 제거하자. 다시 남은 \(n^2-4\)개의 자연수 중 최대를 찾아 이를 \(a_3\)로 잡는다. \(a_3\) 하나, \(\text{gcd}(a_1, a_3)\) 두 개, \(\text{gcd}(a_2, a_3)\) 두 개를 목록에서 제거한다. 이러한 과정을 반복해주면 \(a_1, a_2, \cdots a_n\)을 모두 계산할 수 있다. std::map, std::set 또는 std::multiset 등을 적절히 활용하면 쉽게 구현이 된다. 시간복잡도는 \(\mathcal{O}(n^2 \log n + n^2 \log MAX)\)이다.

\((a,b,c,t)\)를 \((a,b,c)\)에 시간 \(t\)에 도착하는 입자라고 하자. \((a,b,c,t)\)와 \((a',b',c',t')\)이 같은 입자일 필요충분조건은 \(a:b:c:t = a':b':c':t'\)인 것이다. 결국 우리는 \(1 \le a, b, c, t \le n\)이고 \(\text{gcd}(a,b,c,t)=1\)인 순서쌍의 개수를 세주면 충분하다. 이는 \(\displaystyle \sum_{1 \le a,b,c,t \le n} [\text{gcd}(a,b,c,t)=1] = \sum_{1 \le a,b,c,t \le n} \sum_{g|\text{gcd}(a,b,c,t)} \mu (g)\)과 같다. \(1 \le g \le n\)에 대해서 \(\mu(g)\)는  \(a, b, c, t\)가 모두 \(g\)의 배수인 경우에 추가되므로, 이 식은 \(\displaystyle \sum_{g=1}^n \mu(g) \cdot \lfloor \frac{n}{g} \rfloor^4\)과 같다. 나머지는 단순 전처리 + 계산.  

문제를 적당히 환원하면 결국에는 \(ax+by = \text{gcd}(a,b)\)인 순서쌍을 많이 찾는 문제가 된다. 이것은 전형적인 Extended Euclidean Algorithm 연습문제다. 유클리드 알고리즘을 돌리면서 차례대로 계산하는 것이 정석. Modular Inverse를 Extended Euclidean Algorithm으로 구하는 코드가 템플릿에 있어서 다행. \(a, b\)가 \(0\)인 경우도 적당히 처리해줘야 한다.

\(P\)의 크기에 따라서 풀이가 다르다. 먼저 \(P \le 50\)인 경우에는, \(P\) 이하의 소수 개수가 15개 이하이므로, 포함 배제의 원리를 이용하여 \(P\)를 최소 소인수로 갖는 \(X\) 이하인 자연수의 개수를 \(\mathcal{O}(2^{15} \cdot 15)\) 정도의 연산 횟수로 계산할 수 있다. 여기에 이진탐색을 끼얹으면 문제를 해결할 수 있다. \(P > 50\)인 경우에는, 고려해야 할 수의 개수가 \(\frac{10^9}{P} < 2 \cdot 10^7\)임을 이용한다. 목표는 \(\frac{10^9}{P}\) 이하의 수들 중 \(P\) 미만의 소수를 갖는 것을 빠른 시간 내에 전부 표시하는 것이다. 고려해야 할 소수들은 모두 \(\text{min}\left(P, \frac{10^9}{P}\right) \le \sqrt{10^9}\) 이하이므로, 이 범위에 속한 소수들의 배수를 전부 표시해주면 된다. 이제 답을 구하는 것은 쉽다. 

목표는 \(\text{gcd}(i,j+l) = a_l\)이 모든 \(1 \le l \le k\)에 대해 성립하는 \(1 \le i \le n\)과 \(0 \le j \le m-k\)를 찾는 것이다.

이러한 \((i,j)\)를 문제에 대한 해라고 부르자. \(L=\text{lcm}(a_1, a_2, \cdots a_k)\)라 하자. 

Lemma 1. 해가 존재한다면 있다면, \(i = L\)인 해도 존재한다.

Proof of Lemma 1. 존재하는 해를 \((u,v)\)라 한다면, \(u\)는 자명하게 \(L\)의 배수다.

이제 \((L,v)\) 역시 문제에 대한 해임을 보이자. \(1 \le L \le u \le n\)과 \(0 \le v \le m-k\)는 자명하게 성립한다.

또한, \(a_l|\text{gcd}(L,v+l)|\text{gcd}(u,v+l)=a_l\)이므로, \(\text{gcd}(L,v+l) = a_l\)도 성립한다. 증명 끝.

Lemma 2. (1) 문제의 해가 존재하려면 연립합동식 \(x \equiv -l \pmod{a_l}\), \(1 \le l \le k\)의 해가 존재해야 한다. 

Lemma 2. (2) 또한, 위 합동식의 해가 존재한다면 그 해는 \(\pmod{L}\)에서 유일하다.

Proof of Lemma 2. (1)은 \(j\)가 존재해야 하므로 자명하다. (2)는 해 두 개를 잡으면 그 차이가 \(L\)의 배수이므로 자명.

Lemma 3. 위 합동식의 해가 \(x \equiv T \pmod{L}\)이라 하자. 단, \(0 \le T < L\)이다. 문제의 해가 존재하면, \((L, T)\)도 해다.

Proof of Lemma 3. 문제에 대한 해 \((u,v)\)가 존재한다고 하면, Lemma 1의 증명에서 이미 \((L,v)\)가 해임을 보였다. 

한편, \(v\)는 당연히 \(T \pmod{L}\)이어야 한다. 이제 \(0 \le T \le v \le m-k\)가 성립하고, \(\text{gcd}(L,T+l) = \text{gcd}(L,v+l) = a_l\)이 모든 \(1 \le l \le k\)에 대해서 성립한다. 증명 끝.

이제 문제는 \((L,T)\)가 해인지 여부를 확인하는 것으로 바뀌었다. overflow에 조심하면서 \(L \le n\)인지 여부를 확인해주자. 중국인 나머지 정리를 이용해서 연립합동식을 풀어주자. 연립합동식의 해가 없으면 문제에 대한 해도 없다. \(T\)가 계산되었으니 \(T \le m-k\)인지 확인하고, \(\text{gcd}(L,T+l) = a_l\)인지 여부를 확인하자. 전체 시간복잡도는 \(\mathcal{O}(k \log n)\)이다. 

편의상 \(m =  \text{max} b_i = 10^6\)으로 설정한다. 

\(\text{gcd}(x_1, x_2, \cdots x_n)\)이 \(d\)의 배수인 \((x_1, x_2, \cdots , x_n)\)의 개수를 \(mul_d\)라 하자.

\(\text{gcd}(x_1, x_2, \cdots x_n)\)이 \(d\)인 \((x_1, x_2, \cdots , x_n)\)의 개수를 \(cnt_d\)라 하자.

\(mul\)과 \(cnt\)의 정의에 의해서 \(\displaystyle mul_d = \sum_{k \in \mathbb{N}, kd \le m} cnt_{kd}\)가 당연히 성립한다. 

\(1 \le d \le m\)에 대해 \(mul_d\)를 모두 구해놓으면, \(cnt_d\) 역시 \(\mathcal{O}(m \log m)\)에 모두 구할 수 있다.

한편, \(mul_d\)는 \(\displaystyle \prod_{i=1}^n \left( \lfloor \frac{b_i}{d} \rfloor - \lfloor \frac{a_i -1}{d} \rfloor \right)\) 형태로 단순하게 계산된다. 문제는 이것을 모든 \(d\)에 대해 빠르게 계산하는 것.

그런데 자연수 \(X\)에 대해서 \(\lfloor \frac{X}{i} \rfloor\)의 값은 \(1 \le i \le X\)에서 \(\mathcal{O}(\sqrt{X})\)번 변화한다는 사실이 잘 알려져 있다.

그러므로 \(1 \le d \le m\)에서 \(\lfloor \frac{a_i-1}{d} \rfloor\)이나 \(\lfloor \frac{b_i}{d} \rfloor\)의 값은 (단, \(1 \le i \le n\)) 최대 \(\mathcal{O}(n\sqrt{m})\)번 변화한다.

각 변화가 일어나는 위치를 모두 전처리한 다음 \(d\)를 증가시키면서 순서대로 \(mul_d\)를 구해줄 수 있다.

\(\left( \lfloor \frac{b_i}{d} \rfloor - \lfloor \frac{a_i -1}{d} \rfloor \right)\)의 값이 0이었다가 다시 양수가 되는 경우를 조심하자. 시간복잡도는 \(\mathcal{O}(n\sqrt{m} + m \log m)\)이다. 

먼저 정의를 몇 가지 내리자. \(n=\prod_{i=1}^k p^{e_i}_i\)에 대해서 \(S(n) = \text{gcd}(e_1, e_2, \cdots e_k)\)라 하자. 단, \(n>1\)이다.

또한, \(n\)을 height 2 이상의 power tower로 표현하는 방식을 \(R(n)\)이라고 하자.

\(n\)을 height 2 이상의 power tower로 표현한다고 하자. \(n=a^b\)라고 하면, \(b \neq 1\)이고 \(b|S(n)\)이어야 한다. 

또한, \(b\)는 그대로 둘 수도 있고, 다시 \(b\)를 height 2 이상의 power tower로 표기할 수도 있다.

그러므로 \(R(n) = \sum_{b \neq 1, b|S(n)} \left(1+ R(b)\right)\)가 성립한다. \(b \le S(n) = \mathcal{O}(\log n)\)임에 주목하자.

본 문제로 돌아와서, \(n = a^{b^c}\)라 하자. \(n =u^v\)로 표기한다면, 가능한 \(v\)는 \(t =b^c \cdot S(a)\)의 약수다. 

즉, 우리는 \(\sum_{v|t} R(v)\)를 계산해야 한다. \(t = \prod_{i=1}^k p^{e_i}_i\)라 하고, \(m = \text{max}(e_1, e_2, \cdots , e_k)\)라 하자.

\(t\)를 인수분해하는 것은 \(b, S(a)\)를 인수분해하면 충분하며, 이는 간단하다. 

이제 계산해야 하는 것은 \(\displaystyle \sum_{f_1=0}^{e_1} \sum_{f_2=0}^{e_2} \cdots \sum_{f_k=0}^{e_k} R( \prod_{i=1}^k p^{f_i}_i )\)이다. 위에서 구한 \(R\)에 대한 점화식을 이용해보자.

식은 \(\displaystyle \sum_{f_1=0}^{e_1} \sum_{f_2=0}^{e_2} \cdots \sum_{f_k=0}^{e_k} \sum_{g \neq 1, g|\text{gcd}(f_1, f_2, \cdots f_k)} \left(1+R(g)\right) \)이다. 물론 여기서 \(f_1 = f_2 = \cdots = f_k =0\)인 경우는 무시.

식을 \(g\)에 대한 식으로 변형하면, 이는 \(\displaystyle \sum_{g=2}^{m} (1+R(g)) \cdot \left( \prod_{i=1}^k \left( \lfloor \frac{e_i}{g} \rfloor + 1 \right) - 1\right) \)이다.

한편, \(e_i \le \log_2 t = c \log_2 b + \log_2 S(a) \le 150000\)이므로 \(m \le 150000\)이다. 

즉, \(R(n)\)을 \(n \le 150000\)에 대해서만 전처리해주면 문제가 풀린다. 

\(S(n)\)을 미리 계산해줄 수 있고, 위에서 유도한 점화식을 사용하면 \(R(n)\)까지 전처리할 수 있다. 

\(a=0\)이면 \(n=p\)를 주면 끝난다. 이제 \(a \neq 0\)인 경우를 보자. \(\mathbb{Z}_p\)의 원시근 \(g\)를 하나 잡자.

해가 존재하면 무조건 \(\text{gcd}(n,p)=1\)이다. \(n \equiv k_1 \pmod{p-1}\), \(n \equiv g^{k_2} \pmod {p}\)라 하자. 

우리가 원하는 식은 \(g^{k_1 \cdot k_2} + g^{k_2 \cdot m} \equiv a \pmod {p}\)인 것이다.

이제 \(k_2=1\)이라 해보자. 그러면 식은 \(g^{k_1} \equiv a - g^m \pmod{p}\)가 된다.

만약 \(g^m \equiv a \pmod{p}\)가 아니라면, 원시근의 성질에서 조건을 만족하는 \(k_1\)이 존재한다. 

\(k_1\)을 구해주면, \(k_2=1\)이란 것을 이미 알고 있으니 \(n\)도 계산할 수 있다.

그렇지 않는 경우라면? \(g^m \equiv a \pmod{p}\)라고 해보자. 이때는 \(k_2=1\)을 포기하는 것이 맞는 판단이다.

그런데 \(\text{gcd}(k_2, p-1) \neq 1\)이라면, \(g^{k_2}\)가 \(\mathbb{Z}_p\)의 원시근이 아니게 된다. 

이 경우에는 \(a-g^{k_2 \cdot m} \pmod p\)가 \(0\)이 아니더라도 \(g^{k_1 \cdot k_2}\) 형태로 표현하지 못할 수도 있다.

그러니 \(\text{gcd}(k_2, p-1) = 1\)이도록 하자. 한편, \(g^{k_2 \cdot m} = g^m \equiv a \pmod{p}\)라면 앞선 상황과 똑같은 문제가 발생한다. 그러므로 \(k_2 \cdot m \neq m \pmod{p-1}\)이어야 하고, 이는 \((k_2-1) \cdot m \neq 0 \pmod{p-1}\)임을 의미한다. 

즉, 우리가 원하는 \(k_2\)는 \(\text{gcd}(k_2, p-1)=1\)과 \((k_2-1) \cdot m \neq 0 \pmod{p-1}\)을 동시에 만족한다. 

한편, \(\text{gcd}(k_2, p-1)=1\)을 만족하는 \(1 \le k_2 \le p-1\)의 개수는 \(\phi(p-1)\)이다.

\((k_2-1)m \equiv 0 \pmod{p-1}\)을 만족하는 \(1 \le k_2 \le p-1\)의 개수는 최대 \(m\)개이다. 

그러므로 \(\phi(p-1) > m\)이라면 조건을 만족하는 \(k_2\)가 무조건 존재한다. 

\(m \le 20\)과 \(\phi(n) \ge \frac{n}{4 \log n}\)을 쓰면 \(p>800\)이면 원하는 \(k_2\)가 존재함을 알 수 있다. 이제 \(k_1\)을 구하고, \(n\)을 계산하면 끝. 

\(p \le 800\)이면 \(1 \le n \le p(p-1)\)을 만족하는 모든 \(n\)에 대해서 조건식이 성립하는지 여부를 확인하여 답을 구할 수 있다.

이제 구현에 대한 이야기를 하자. 어떤 \(\text{gcd}(x,p)=1\)인 \(x\)가 주어졌을 때, \(x\)가 \(\mathbb{Z}_p\)의 원시근인지 판별하는 것은 \(p-1\)의 각 소인수 \(p_i\)에 대해서 \(x^{(p-1)/p_i} \not\equiv 1 \pmod{p}\)임을 확인하는 것으로 할 수 있다. 원시근의 수는 \(\phi(p-1)\)개로 충분히 많기 때문에, 랜덤한 자연수를 원시근이 나올 때까지 계속 뽑는 방법을 사용할 수 있다. \(k_2\)를 선택하는 방법도 비슷하게 랜덤으로 할 수 있다. \(k_1\)을 구하는 것은 대표적인 Discrete Logarithm을 구하는 알고리즘인 Baby-Step-Giant-Step 알고리즘으로 할 수 있다. 이제 \(n\)을 구하는 것은 중국인 나머지 정리로 할 수 있다. \(p>800\)인 경우에 시간복잡도는 테케당 \(\mathcal{O}(\sqrt{p} \log p)\) 정도.

\(\text{gcd}(F_i, F_j) = F_{\text{gcd}(i,j)}\)가 성립하므로, 구해야 하는 값은 \(\displaystyle \sum_{1 \le i, j \le n} F_{\text{gcd}(i,j)}\)이다. 

\(1 \le i, j \le n\)이면서 \(\text{gcd}(i,j)=d\)인 \((i,j)\)의 수를 \(cnt_d\)라 하자. 우리가 구해야 하는 값은 \(\displaystyle \sum_{d=1}^n F_d \cdot cnt_d\)가 된다.

그런데 \(cnt_d\)는 \(1 \le i, j \le \lfloor \frac{n}{d} \rfloor\)이면서 \(\text{gcd}(i,j)=1\)인 \((i,j)\)의 수와 같다. 

이를 정리하면 \(\displaystyle cnt_d = 2 \cdot \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \phi(i) - 1\)를 얻는다. 또한, 앞에서 언급했듯이 \(\lfloor \frac{n}{d} \rfloor\)로 가능한 값은 \(\mathcal{O}(\sqrt{n})\)개이다. 

\(\lfloor \frac{n}{d} \rfloor\)로 가능한 값을 \(a_1, a_2 , \cdots a_t\)라 하자. 또한, \(\lfloor \frac{n}{d} \rfloor = a_i\)인 \(d\)의 범위를 \([lef_i, rig_i]\)라 하자.

우리가 구해야 하는 것은 결국 \(\displaystyle \sum_{i=1}^t \left( \sum_{j=lef_i}^{rig_i} F_j \right) \cdot \left( 2 \cdot \sum_{k=1}^{a_i} \phi(k) - 1\right) \)가 된다.

오일러 파이 함수의 부분합은 xudyh's sieve로 빠른 시간에 구할 수 있음이 잘 알려져 있다.

이제 남은 것은 연속한 피보나치 수들의 합을 빠르게 구하는 것이다. 

\(F_1 + F_2 + \cdots + F_n = F_{n+2}-1\)가 성립함은 귀납법으로 쉽게 증명할 수 있고, 이를 활용하면 연속한 피보나치 수들의 합을 빠르게 구할 수 있다. 피보나치 수들은 당연히 로그 시간에 구할 수 있다. 전체 시간복잡도는 \(\mathcal{O}(n^{2/3})\)이다.

위 문제에서 "연속한 피보나치 수의 합"을 BOJ 13716으로 바꾼 버전이다.

오일러 파이의 합, xudyh's sieve 등의 과정은 그대로 남는다. 문제는 BOJ 13716을 빠르게 해결하는 것이다.

\(f(X,t) = \sum_{i=0}^t (X+i)^k \cdot \binom{t}{i} \cdot (-1)^i\)이라고 정의하고, \(h(X) = \sum_{i=0}^k f(X,i) \cdot F_{X+2i+1}\)이라고 정의하자.

문제의 핵심은 \(h(X+1)-h(1) = \sum_{i=1}^X F_i \cdot i^k\)라는 것이다. 증명은 생략한다.

식을 적당히 조작하면, \(h(X) = \sum_{i=0}^k (X+i)^k \cdot (C_i F_X + D_i F_{X-1})\)가 모든 \(X\)에 대해서 성립하는 \(C_i, D_i\)를 \(\mathcal{O}(k^2)\)에 구할 수 있다. 이 과정에서 피보나치 수의 항등식인 \(F_{a+b} = F_{a-1}F_b+F_a F_{b+1}\)이 활용된다.

전처리를 통해서 \(C_i , D_i\)를 계산한 뒤에는, \(\sum_{i=1}^X F_i \cdot i^k\)를 \(\mathcal{O}(k \log k)\)에 구할 수 있다. 

이제 \(\sum_{i=1}^n F_i \cdot i^k\)를 \(n \le 2.5 \cdot 10^6\)까지 전처리하자. 전처리는 선형 시간에 가능하다.

그러면 \(\mathcal{O}(k \log k)\) 시간에 \(\sum_{i=1}^X F_i \cdot i^k\)를 계산해야 하는 횟수는 최대 400번이므로 시간 안에 문제가 풀린다.

UPD: 증명을 생략하는 것도 그렇고, 증명을 다시 생각하기도 그렇고 해서 AOPS에 이 문제를 올렸다.

식 형태에서 꽤 잘 드러나듯이, 차분법을 최대한 활용하면 된다. 링크.