[国家集训队] Crash的数字表格 / JZPTAB（莫比乌斯反演） | 错题本

题目

[国家集训队] Crash的数字表格 / JZPTAB

分析

不妨设 $n\ge m$：
$$\begin{array}{rl}& \sum _{i=1}^n\sum _{j=1}^m\text{lcm}(i,j)\\ =& \sum _{i=1}^n\sum _{j=1}^m\frac{ij}{\gcd (i,j)}\\ =& \sum _{d=1}^m\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }\frac{ij}{d}[\gcd (i,j)=1]\\ =& \sum _{d=1}^m\frac{1}{d}\sum _{k=1}^{\lfloor \frac{m}{d}\rfloor }\mu (k)\sum _{i=1}^{\lfloor \frac{n}{dk}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{dk}\rfloor }(idk)(jdk)\\ =& \sum _{d=1}^{m}d\sum _{k=1}^{\lfloor \frac{m}{d}\rfloor }\mu (k)k^2\sum _{i=1}^{\lfloor \frac{n}{dk}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{dk}\rfloor }ij\end{array}$$ 令 $f(n,m)=\sum _{k=1}^m\mu (k)k^2\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }ij$，由于 $\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }ij=\left(\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }i\right)\left(\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }j\right)$ 即两个等差数列之和的乘积，所以这个可以数论分块求出。原式即为 $$\sum _{d=1}^{m}df\left(\lfloor \frac{n}{d}\rfloor ,\lfloor \frac{m}{d}\rfloor \right)$$ 亦可数论分块。总时间复杂度 $O(n)$。

代码

#include <bits/stdc++.h>

typedef long long LL;

const int MAXN = 10000000;
const int MOD = 20101009;

int Mu[MAXN + 5];
bool Vis[MAXN + 5];
std::vector<int> Primes;

void Init(int n) {
    Mu[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!Vis[i])
            Mu[i] = -1, Primes.push_back(i);
        for (int j = 0; j < (int)Primes.size() && i * Primes[j] <= n; j++) {
            Vis[i * Primes[j]] = true;
            if (i % Primes[j] == 0) {
                Mu[i * Primes[j]] = 0;
                break;
            }
            Mu[i * Primes[j]] = -Mu[i];
        }
        Mu[i] = ((Mu[i - 1] + (LL)Mu[i] * i * i % MOD) % MOD + MOD) % MOD;
    }
}

int Sum(int n) {
    return (LL)n * (n + 1) / 2 % MOD;
}

int Calc(int n, int m) {
    if (n < m) std::swap(n, m);
    int ret = 0;
    for (int lft = 1, rgt = 1; lft <= m; lft = rgt + 1) {
        rgt = std::min(n / (n / lft), m / (m / lft));
        ret = (ret + (LL)(Mu[rgt] - Mu[lft - 1]) % MOD * Sum(n / lft) % MOD * Sum(m / lft) % MOD) % MOD;
    }
    return (ret + MOD) % MOD;
}

int main() {
    Init(MAXN);
    int N, M; scanf("%d%d", &N, &M);
    if (N < M) std::swap(N, M);
    int Ans = 0;
    for (int lft = 1, rgt = 1; lft <= M; lft = rgt + 1) {
        rgt = std::min(M / (M / lft), N / (N / lft));
        Ans = (Ans + (LL)(lft + rgt) * (rgt - lft + 1) / 2 % MOD * Calc(N / lft, M / lft)) % MOD;
    }
    printf("%d", Ans);
    return 0;
}