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

题解（公式推导）

\begin{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)}g\sum_{i=1}^{n}\sum_{j=1}^{m}[\gcd(i,j)=g]\frac{i}{g}\frac{j}{g}\\ &=\sum_{g=1}^{\min(n,m)}g\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{g}\rfloor}[\gcd(i,j)=1]ij\\ &=\sum_{g=1}^{\min(n,m)}g\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}i\sum_{j=1}^{\lfloor\frac{m}{g}\rfloor}j[\gcd(i,j)=1]\\ &=\sum_{g=1}^{\min(n,m)}g\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}i\sum_{j=1}^{\lfloor\frac{m}{g}\rfloor}j\sum_{g|i}\sum_{g|j}\mu(g)\\ \end{aligned}

$$\mathbf{F}(n,m)=\sum\limits_{i=1}^{n}i\sum\limits_{j=1}^{m}j\sum\limits_{g|i}\sum\limits_{g|j}\mu(g)$$，则原式 $$=\sum\limits_{g=1}^{\min(n,m)}g\mathbf{F}(\frac{n}{g},\frac{m}{g})$$

\begin{aligned} \mathbf{F}(n,m) &=\sum_{i=1}^{n}i\sum_{j=1}^{m}j\sum_{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}i\sum_{j=1}^{m}j[g|i][g|j]\\ &=\sum_{g=1}^{\min(n, m)}\mu(g)\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}id\sum_{j=1}^{\lfloor\frac{m}{g}\rfloor}jd\\ &=\sum_{g=1}^{\min(n, m)}\mu(g)g^2\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}i\sum_{j=1}^{\lfloor\frac{m}{g}\rfloor}j\\ \end{aligned}

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

\begin{aligned} \mathbf{F}(n,m) &=\sum_{g=1}^{\min(n, m)}\mu(g)g^2\mathbf{S}(1,\lfloor\frac{n}{g}\rfloor)\mathbf{S}(1,\lfloor\frac{m}{g}\rfloor)\\ \end{aligned}

代码

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;