【BZOJ2154】Crash的数字表格
题面
题解
不妨设(nleq m)
题目是求：
[sum_{i=1}^nsum_{j=1}^mlcm(i,j)
]
还是照常推式子qaq
[sum_{i=1}^nsum_{j=1}^mfrac{ij}{gcd(i,j)}\
sum_{d=1}^nsum_{i=1}^nsum_{j=1}^m[gcd(i,j)==d]frac{ij}{d}\
sum_{d=1}^nsum_{i=1}^{n/d}sum_{j=1}^{m/d}[gcd(i,j)==1]ijd\
sum_{d=1}^ndsum_{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}^xsum_{j=1}^y[gcd(i,j)==k]ij\
g(k)=sum_{k|d}f(d)
]
常见套路
[g(k)=sum_{i=1}^xsum_{j=1}^y[k|gcd(i,j)]ij
]
提出(k)
[g(k)=k^2sum_{i=1}^{x/k}sum_{j=1}^{y/k}[1|gcd(i,j)]ij\
g(k)=k^2sum_{i=1}^{x/k}sum_{j=1}^{y/k}ij
]
这个就是两个等差数列相乘。
要求的东西变为
[f(1)=sum_{i=1}^xmu(i)g(i)
]
问题解决很多了，但是还必须要往下
[ans=sum_{d=1}^ndsum_{i=1}^{n/d}sum_{j=1}^{m/d}[gcd(i,j)==1]ij
]
看看(f(1))是什么
[f(1)=sum_{i=1}^xmu(i)i^2sum_{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
", ans);
return 0;
}