백준 20940번 시철이가 사랑한 수식

문제 요약

아래 두 수식의 값을 구하면 된다.
$$
\begin{gather}
\sum_{n=1}^N \sum_{i=1}^n \sum_{j=1}^n (gcd(i,j) \times lcm(i,j)) \\
(\sum_{n=1}^N \sum_{i=1}^n \sum_{j=1}^n gcd(i, j)) \times (\sum_{n=1}^N \sum_{i=1}^n \sum_{j=1}^n lcm(i,j))
\end{gather}
$$
N과 K가 주어지는데 N은 최대 100만이고 K는 위의 값들을$\mod K$로 구하라는 뜻이다. K는 소수임이 보장된다.

풀이

먼저 첫번째 식을 보자. $lcm(i,j) = \frac{ij}{gcd(i,j)}$다. 따라서, 시그마 안의 식은 그냥 $ij$다.
$$
\begin{aligned}
\sum_{n=1}^N \sum_{i=1}^n \sum_{j=1}^n ij &= \sum_{n=1}^N \sum_{i=1}^n i \frac{n(n+1)}{2} \\
&= \sum_{n=1}^N(\frac{n(n+1)}{2})^2
\end{aligned}
$$
이 값은 $O(N)$에 구하는 것이 가능하다.

두 번째 식은 각 항을 순서대로 살펴보자. 먼저 이런 값을 생각하자.
$$
sum_n = \sum_{i=1}^{n-1} gcd(i,n)
$$
그러면 첫번째 항은 다음과 같이 바뀐다.
$$
\sum_{n=1}^N 2sum_n+n
$$
$sum_n$만 빠르게 전처리 할 수 있으면 첫번째 항을 구할 수 있을 것 같다.

$gcd(i, n) = g$인 횟수를 $count[g]$라고 하자. 그러면 $sum_n = \sum g \times count[g]$가 된다.

여기서 $gcd(i,n)=g$라면 $gcd(i/g, n/g)=1$이다.

여기서 오일러 파이 함수를 사용할 수 있다. $\phi(N)$은 N보다 작거나 같으면서 $N$과 서로소인 수이다.

이 정의를 이용하면 $i \le n$이면서 $gcd(i,n)=g$인 $i$의 수는 $\phi(n/g)$가 된다. 따라서 $sum_n$은 아래와 같이 식이 변형된다.
$$
sum_n = \sum_{d \vert n} d \times \phi(n/d)
$$
이제 $sum_n$을 하나 구하는 데에 $O(\sqrt{n})$까지는 줄였다. 그럼 $sum_1$부터 $sum_N$까지 구하는데에 $O(N\sqrt{N})$인건데 좀 크다.

이거보다 더 줄이는 것이 가능하다. 에라토스테네스의 체와 비슷한 방식으로 1부터 N까지의 $sum_n$을 구하는 것이 가능하다.

그래서 $\phi$함수의 전처리와 $sum_n$의 전처리에 $O(NloglogN)$이 걸린다.

이제 $lcm$으로 구성된 항을 보자. 여기서도 $gcd$와 비슷하게 아래와 같은 값을 생각하자.
$$
sum_n' = \sum_{i=1}^{n-1} lcm(i,n) = \sum_{i=1}^{n-1} \frac{in}{gcd(i,n)} = n\sum_{i=1}^{n-1} \frac{i}{gcd(i,n)}
$$
저 시그마 안의 값이 또 일일이 구하려면 골치가 아프다. 여기서도 $gcd$에서 적용했던 것과 같은 논리를 적용하자.

$\frac{i}{g}$가 될 수 있는 것들은 $\frac{n}{g}$와 서로소인 것들만 가능하다. 어떤 수 $n$보다 작거나 같으면서 서로소인 수들의 합을 $f(n)$이라 하면 식을 아래처럼 바꿔쓸 수 있는 것이다.
$$
sum_n' = n \sum_{d \vert n} f(d)
$$
놀랍게도 $f(n)$은 $\frac{n\phi(n)}{2}$로 알려져 있다. 오일러 파이 함수의 영문 위키피디아에서 그걸 찾아볼 수 있다.

그래서 $sum_n'$또한 $O(NloglogN)$으로 구하는 것이 가능하고 문제를 풀 수 있다.

참고로, 주어진 식에서 제일 바깥 시그마를 떼어낸 식은 GCD extreme이나 LCM extreme과 같은 이름으로 부르기도 한다.

코드에서 $sum_n$은 gcd_ex[n]으로 1부터 n까지 구한거다. FastMod는 Barett Reduction을 적용한것으로 모듈로 연산을 빠르게 해주는 역할을 한다. 물론 이거 없이도 AC를 받을 수 있으니 그냥 모듈로 연산이라고 생각해도 무방하다.

#include<bits/stdc++.h>
using namespace std;
using ll = unsigned long long;
const int MAXN = 1e6;
typedef unsigned long long ull;
typedef __uint128_t L;
struct FastMod {
    ull b, m;
    FastMod(ull b) : b(b), m(ull((L(1) << 64) / b)) {}
    ull reduce(ull a) {
        ull q = (ull)((L(m) * a) >> 64);
        ull r = a - q * b; // can be proven that 0 <= r < 2*b
        return r >= b ? r - b : r;
    }
};
FastMod F(2);
ll mod, N;
ll phi[MAXN + 5];
ll gcd_ex[MAXN + 5], lcm_ex[MAXN + 5];

ll ipow(ll a, ll b) {
    ll ret = 1;
    while (b) {
        if (b & 1) ret = F.reduce(ret * a);
        a = F.reduce(a*a); b >>= 1;
    }
    return ret;
}

void calc_phi() {
    phi[0] = 0;
    phi[1] = 1;
    for (int i = 2; i <= N; ++i) phi[i] = i;
    for (int i = 2; i <= N; ++i) {
        if (phi[i] == i) {
            for (int j = i; j <= N; j += i) phi[j] -= phi[j] / i;
        }
    }
}

void calc_gcd_ex() {
    for (int i = 1; i <= N; ++i) {
        for (int j = i; j <= N; j += i) {
            gcd_ex[j] += F.reduce(i*phi[j/i]);
            gcd_ex[j] = F.reduce(gcd_ex[j]);
        }
    }
    for (int i = 1; i <= N; ++i) {
        gcd_ex[i] = F.reduce(gcd_ex[i] + mod - i);
        gcd_ex[i] = F.reduce(2*gcd_ex[i] + i);
    }
}

void calc_lcm_ex() {
    ll two_inv = ipow(2, mod - 2);
    for (ll i = 2; i <= N; ++i) {
        for (ll j = i; j <= N; j += i) {
            ll t = F.reduce(i*phi[i]);
            t = F.reduce(t * two_inv);
            lcm_ex[j] += F.reduce(t*j);
            lcm_ex[j] = F.reduce(lcm_ex[j]);
        }
    }
    ll cnt = 0;
    for (int i = 1; i <= N; ++i) {
        lcm_ex[i] = F.reduce(lcm_ex[i] * 2 + i);
    }
}

ll out1, out2, out3;
int main() {
    cin >> N >> mod;
    F = FastMod(mod);
    calc_phi();
    calc_gcd_ex();
    calc_lcm_ex();
    ll cnt = 0;
    for (int i = 1; i <= N; ++i) {
        cnt = F.reduce(cnt + i);
        out1 += F.reduce(cnt * cnt);
        out1 = F.reduce(out1);
    }
    cout << out1 << '\n';
    for (int i = 1; i <= N; ++i) {
        gcd_ex[i] = F.reduce(gcd_ex[i] + gcd_ex[i-1]);
        out2 = F.reduce(out2 + gcd_ex[i]);
        lcm_ex[i] = F.reduce(lcm_ex[i] + lcm_ex[i-1]);
        out3 = F.reduce(out3 + lcm_ex[i]);
    }
    cout << F.reduce(out2 * out3) << '\n';
    return 0;
}