「Bzoj 2154」Crash的数字表格 (莫比乌斯反演)

Bzoj 2154
题意：给定$n, m$，求
$$
\sum_{i=1}^n \sum_{j=1}^m \operatorname{lcm}(i, j)
$$

化为$\gcd$式子
$$
\sum_{d=1}^n d \sum_{i=1}^{n}\sum_{j=1}^{m}\frac{ij}{\gcd(i,j)} \\
\sum_{d=1}^n d \sum_{i=1}^{\lfloor \frac nd \rfloor}\sum_{j=1}^{\lfloor \frac md \rfloor}ij \cdot [\gcd(i,j)=1]
$$
设
$$
f(n,m)=\sum_{i=1}^n \sum_{j=1}^m ij \cdot [\gcd(i,j)=1] \\
=\sum_{d=1}^n \sum_{d|i} \sum_{d|j} \mu(d) \cdot d^2 \sum_{i=1}^{\lfloor \frac nd \rfloor}\sum_{j=1}^{\lfloor \frac md \rfloor}ij
$$
设 (两个等差数列相乘)
$$
g(n,m)=\sum_{i=1}^n \sum_{j=1}^m ij =\frac{n(n+1)}{2} \times \frac{m(m+1)}{2}
$$
则
$$
f(n,m)=\sum_{d=1}^n \mu(d) \cdot d^2 \cdot g(\lfloor \frac nd \rfloor, \lfloor \frac md \rfloor)
$$
即原式为
$$
\sum_{d=1}^m d \cdot f(\lfloor \frac nd \rfloor, \lfloor \frac md \rfloor)
$$

计算$f(n,m)$时整除分块，计算答案时整除分块，时间复杂度$O(n+m)$

#include<cstdio> 
#include<cstring>
#include<algorithm>
#include<iostream>
#include<vector>
#include<queue>
#include<cmath>
#define ms(i, j) memset(i, j, sizeof i)
#define LL long long
#define db long double
#define fir first
#define sec second
#define mp make_pair
using namespace std;

namespace flyinthesky {

    const LL MAXN = 10000000 + 5;
    const LL MO = 20101009;

    LL n, m, mu[MAXN], vis[MAXN], pri[100000 + 5], tot, qzhmu[MAXN];

    LL ksm(LL a, LL b) {
        LL bs = a % MO, ans = 1ll;
        while (b) {
            if (b & 1) ans = (ans * bs) % MO;
            bs = (bs * bs) % MO;
            b >>= 1ll;
        }
        return ans;
    }

    LL g(LL x, LL y) {
        LL tmp1 = x * (x + 1ll) / 2ll % MO;
        LL tmp2 = y * (y + 1ll) / 2ll % MO;
        return tmp1 * tmp2 % MO;
    }
    LL f(LL x, LL y) {
        if (x > y) swap(x, y);
        LL l = 1, ans = 0;
        while (l <= x) {
            LL r = min(x / (x / l), y / (y / l));
            ans = (ans + (qzhmu[r] - qzhmu[l - 1] + MO) % MO * g(x / l, y / l) % MO) % MO;
            l = r + 1;
        }
        return ans;
    }

    void clean() {
        ms(vis, 0), ms(mu, 0), ms(qzhmu, 0), tot = 0;
    }
    int solve() {

        clean();

        scanf("%lld%lld", &n, &m);
        if (n > m) swap(n, m);

        mu[1] = 1, vis[1] = 1;
        for (LL i = 2; i <= m; ++i) {
            if (!vis[i]) vis[i] = 1, pri[++tot] = i, mu[i] = -1;
            for (LL j = 1; j <= tot && i * pri[j] <= m; ++j) {
                vis[i * pri[j]] = 1;
                if (i % pri[j] == 0) {
                    mu[i * pri[j]] = 0;
                    break;
                } else mu[i * pri[j]] = -mu[i];
            }
        }
        for (LL i = 1; i <= m; ++i) qzhmu[i] = (qzhmu[i - 1] + (((mu[i] + MO) % MO * i) % MO) * i % MO) % MO;

        LL l = 1, ans = 0;
        while (l <= n) {
            LL r = min(n / (n / l), m / (m / l));
            ans = (ans + 1ll * (r - l + 1ll) * (l + r) / 2ll % MO * f(n / l, m / l) % MO) % MO;
            l = r + 1;
        }
        printf("%lld\n", ans);
        return 0;
    }
}
int main() {
    flyinthesky::solve();
    return 0;
}