【题解】洛谷 P1829 [国家集训队]Crash的数字表格 / JZPTAB

题解（公式推导）

[egin{aligned} sum_{i=1}^{n}sum_{j=1}^{m}operatorname{lcm}(i,j) &=sum_{i=1}^{n}sum_{j=1}^{m}frac{ij}{gcd(i,j)}\ &=sum_{i=1}^{n}sum_{j=1}^{m}sum_{g=1}^{min(i,j)}[gcd(i,j)=g]frac{ij}{g}\ &=sum_{i=1}^{n}sum_{j=1}^{m}sum_{g=1}^{min(i,j)}[gcd(i,j)=g]frac{i}{g}frac{j}{g}g\ &=sum_{g=1}^{min(n,m)}gsum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)=g]frac{i}{g}frac{j}{g}\ &=sum_{g=1}^{min(n,m)}gsum_{i=1}^{lfloorfrac{n}{g} floor}sum_{j=1}^{lfloorfrac{m}{g} floor}[gcd(i,j)=1]ij\ &=sum_{g=1}^{min(n,m)}gsum_{i=1}^{lfloorfrac{n}{g} floor}isum_{j=1}^{lfloorfrac{m}{g} floor}j[gcd(i,j)=1]\ &=sum_{g=1}^{min(n,m)}gsum_{i=1}^{lfloorfrac{n}{g} floor}isum_{j=1}^{lfloorfrac{m}{g} floor}jsum_{g|i}sum_{g|j}mu(g)\ end{aligned} ]

设 (mathbf{F}(n,m)=sumlimits_{i=1}^{n}isumlimits_{j=1}^{m}jsumlimits_{g|i}sumlimits_{g|j}mu(g))，则原式 (=sumlimits_{g=1}^{min(n,m)}gmathbf{F}(frac{n}{g},frac{m}{g}))

考虑化简 (mathbf{F}(n,m))：

[egin{aligned} mathbf{F}(n,m) &=sum_{i=1}^{n}isum_{j=1}^{m}jsum_{g=1}^{min(n, m)}sum_{g|i}sum_{g|j}mu(g)\ &=sum_{g=1}^{min(n, m)}mu(g)sum_{i=1}^{n}isum_{j=1}^{m}j[g|i][g|j]\ &=sum_{g=1}^{min(n, m)}mu(g)sum_{i=1}^{lfloorfrac{n}{g} floor}idsum_{j=1}^{lfloorfrac{m}{g} floor}jd\ &=sum_{g=1}^{min(n, m)}mu(g)g^2sum_{i=1}^{lfloorfrac{n}{g} floor}isum_{j=1}^{lfloorfrac{m}{g} floor}j\ end{aligned} ]

记 (mathbf{S}(l, r)=frac{(l+r)(r-l+1)}{2}) （等差数列求和公式）

那么

[egin{aligned} mathbf{F}(n,m) &=sum_{g=1}^{min(n, m)}mu(g)g^2mathbf{S}(1,lfloorfrac{n}{g} floor)mathbf{S}(1,lfloorfrac{m}{g} floor)\ end{aligned} ]

至此，(mathbf{F}(n,m)) 已经可以用数论分块在 (Theta(sqrt{min(n, m)})) 的时间内求出。

而对原式 (sumlimits_{g=1}^{min(n,m)}gmathbf{F}(frac{n}{g},frac{m}{g})) 也可用用一次数论分块（(Theta(sqrt{min(n, m)}))）

这样的总时间复杂度是 (Theta(min(n,m))=Theta(n)) 的（认为 (n,m) 同阶）

代码

const int N = 1e7 + 7, P = 20101009, inv6 = 16750841;

bool IsntPrime[N];
int primes[N], pcnt;
int mu[N];
int64 sum[N];

void init(int maxx) {
    sum[1] = mu[1] = 1;
    for (int i = 2; i <= maxx; ++i) {
        if (!IsntPrime[i]) primes[++pcnt] = i, mu[i] = -1;
        for (int p = 1; p <= pcnt && i * primes[p] <= maxx; ++p) {
            IsntPrime[i * primes[p]] = true;
            if (i % primes[p] == 0) break;
            mu[i * primes[p]] = -mu[i];
        }
        sum[i] = (sum[i - 1] + (int64)i * i * mu[i]) % P;
    }
}

inline int64 sum_i1(int64 n) { return (n * (n + 1) >> 1) % P; }

#define query_sum(l, r) (sum[r] - sum[(l) - 1])

inline int64 solve(int n, int m) {
    int64 ans = 0;
    for (int l = 1, r; l <= min(n, m); l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        ans += query_sum(l, r) * sum_i1(n / l) % P * sum_i1(m / l) % P;
    }
    return ans % P;
}

signed main() {
    int n, m;
    read >> n >> m;
    init(min(n, m));
    int64 ans = 0;
    for (int l = 1, r; l <= min(n, m); l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        ans += ((int64)(l + r) * (r - l + 1) >> 1) * solve(n / l, m / l) % P;
    }
    write << (ans % P + P) % P << '
';
    return 0;
}