【BZOJ2154】Crash的数字表格
【BZOJ2154】Crash的数字表格
题面
题解
不妨设\(n\leq m\)
题目是求：
\[\sum_{i=1}^n\sum_{j=1}^mlcm(i,j)
\]
还是照常推式子qaq
\[\sum_{i=1}^n\sum_{j=1}^m\frac{ij}{gcd(i,j)}\\
\sum_{d=1}^n\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)==d]\frac{ij}{d}\\
\sum_{d=1}^n\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)==1]ijd\\
\sum_{d=1}^nd\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)==1]ij
\]
暂时不管前面的，先搞后面的一坨
设\(x=n/d,y=m/d\)
并且设
\[f(k)=\sum_{i=1}^x\sum_{j=1}^y[gcd(i,j)==k]ij\\
g(k)=\sum_{k|d}f(d)
\]
常见套路
\[g(k)=\sum_{i=1}^x\sum_{j=1}^y[k|gcd(i,j)]ij
\]
提出\(k\)
\[g(k)=k^2\sum_{i=1}^{x/k}\sum_{j=1}^{y/k}[1|gcd(i,j)]ij\\
g(k)=k^2\sum_{i=1}^{x/k}\sum_{j=1}^{y/k}ij
\]
这个就是两个等差数列相乘。
要求的东西变为
\[f(1)=\sum_{i=1}^x\mu(i)g(i)
\]
问题解决很多了，但是还必须要往下
\[ans=\sum_{d=1}^nd\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)==1]ij
\]
看看\(f(1)\)是什么
\[f(1)=\sum_{i=1}^x\mu(i)i^2\sum_{i=1}^{x/d}\sum_{j=1}^{y/d}ij
\]
里外都可以数论分块，总复杂度\(O(n)\)
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
const int Mod = 20101009;
const int MAX_N = 1e7 + 5;
int N = 1e7, M, mu[MAX_N], prime[MAX_N], num;
bool is_prime[MAX_N];
int s[MAX_N];
void sieve() {
for (int i = 1; i <= N; i++) is_prime[i] = 1;
is_prime[1] = 0, mu[1] = 1;
for (int i = 2; i <= N; i++) {
if (is_prime[i]) { prime[++num] = i, mu[i] = -1; }
for (int j = 1; j <= num && i * prime[j] <= N; j++) {
is_prime[i * prime[j]] = 0;
if (i % prime[j]) mu[i * prime[j]] = -mu[i];
else { mu[i * prime[j]] = 0; break; }
}
}
}
int solve(int x, int y) {
int ans = 0;
for (int l = 1, r; l <= x; l = r + 1) {
r = min(x / (x / l), y / (y / l));
int res = 1ll * (1ll * (1 + x / l) * (x / l) / 2 % Mod) * (1ll * (1 + y / l) * (y / l) / 2 % Mod) % Mod;
ans = (1ll * (s[r] - s[l - 1]) % Mod * res % Mod + ans) % Mod;
}
return (ans + Mod) % Mod;
}
int main () {
sieve();
cin >> N >> M;
if (N > M) swap(N, M);
for (int i = 1; i <= N; i++) s[i] = (s[i - 1] + 1ll * i * i % Mod * mu[i] % Mod + Mod) % Mod;
int ans = 0;
for (int l = 1, r; l <= N; l = r + 1) {
r = min(N / (N / l), M / (M / l));
int res = 1ll * (l + r) * (r - l + 1) / 2 % Mod;
ans = (ans + 1ll * solve(N / l, M / l) * res % Mod) % Mod;
}
printf("%d\n", ans);
return 0;
}
