洛谷1829：crash的数字表格
洛谷1829：crash的数字表格
题意：
- 求\(\sum_{i=1}^N\sum_{j=1}^Mlcm(i,j)\)。
- 数据范围\(:n,m\leq10^7\)。
思路：
易得：原式
\[\sum_{i=1}^N\sum_{j=1}^M\frac{ij}{gcd(i,j)}
\]
枚举\(gcd(i,j)=d\)。
\[\sum_{d=1}^{min(N,M)}\sum_{i=1}^N\sum_{j=1}^M[gcd(i,j)=d]\frac{ij}{d}
\]
把\(d\)给提出来，其实就是将枚举项\(i,j\)看成是\(di,dj\)。
\[\sum_{d=1}^{min(N,M)}\frac{1}{d}\sum_{i=1}^{\frac{N}{d}}\sum_{j=1}^{\frac{M}{d}}[gcd(i,j)=1]idjd.\\\sum_{d=1}^{min(N,M)}d\sum_{i=1}^{\frac{N}{d}}\sum_{j=1}^{\frac{M}{d}}[gcd(i,j)=1]ij
\]
我们知道莫比乌斯函数的一个重要性质。
\[\sum_{d|n}\mu(d)=[n=1]
\]
代入可得：
\[\sum_{d=1}^{min(N,M)}d\sum_{i=1}^{\frac{N}{d}}\sum_{j=1}^{\frac{M}{d}}\sum_{x|gcd(i,j)}\mu(x)ij
\]
改为枚举\(x\)。
\[\sum_{d=1}^{min(N,M)}d\sum_{x=1}^{min(\frac{N}{d},\frac{M}{d})}\mu(x)\sum_{i=1}^{\frac{N}{d}}\sum_{j=1}^{\frac{M}{d}}ij[x|gcd(i,j)]
\]
因为后面有\(x|gcd(i,j)\)的条件，所以我们可以将枚举\(i,j\)改为枚举\(ix,jx\)。
那么就有:
\[\sum_{d=1}^{min(N,M)}d\sum_{x=1}^{min(\frac{N}{d},\frac{M}{d})}\mu(x)\sum_{i=1}^{\frac{N}{dx}}\sum_{j=1}^{\frac{M}{dx}}ijx^2
\]
接着把\(x^2\)提到前面去：
\[\sum_{d=1}^{min(N,M)}d\sum_{x=1}^{min(\frac{N}{d},\frac{M}{d})}x^2\mu(x)\sum_{i=1}^{\frac{N}{dx}}i\sum_{j=1}^{\frac{M}{dx}}j
\]
后面的两个\(sum\)其实就是两个等差数列求和。
中间的\(sum\)直接\(O(n)\)预处理即可。
我们设\(sum(\frac{N}{d},\frac{M}{d})=\sum_{x=1}^{min(\frac{N}{d},\frac{M}{d})}x^2\mu(x)\sum_{i=1}^{\frac{N}{dx}}i\sum_{j=1}^{\frac{M}{dx}}j\).
那么原式其实就是\(\sum d*sum(\frac{N}{d},\frac{M}{d})\).
两次数论分块。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 1e7+10;
const int mod = 20101009;
const int inv2 = (mod+1)/2;
ll n, m;
int primes[maxn], cnt;
ll mu[maxn];
bool vis[maxn];
void init(int n)
{
mu[1] = 1;
for(int i = 2; i <= n; i++)
{
if(!vis[i])
{
primes[++cnt] = i;
mu[i] = -1;
}
for(int j = 1; primes[j] <= n/i; j++)
{
vis[primes[j]*i] = 1;
if(i % primes[j] == 0) break;
else mu[i*primes[j]] = -mu[i];
}
}
for(ll i = 1; i <= n; i++)
mu[i] = (mu[i-1]+(mu[i]+mod)*i%mod*i%mod)%mod;
}
ll sum(ll x, ll y){
return x%mod*(x+1)%mod*inv2%mod*y%mod*(y+1)%mod*inv2%mod;
}
ll work(ll x, ll y)
{
ll res = 0;
for(ll l = 1, r; l <= min(x, y); l = r+1)
{
r = min(x/(x/l), y/(y/l));
res = (res+((mu[r]-mu[l-1])%mod+mod)%mod*sum(x/l, y/l)%mod)%mod;
}
return res;
}
void solve(ll n, ll m)
{
ll res = 0;
for(ll l = 1, r; l <= min(n, m); l=r+1)
{
r = min(n/(n/l), m/(m/l));
res = (res+(l+r)%mod*(r-l+1)%mod*inv2%mod*work(n/l, m/l)%mod)%mod;
}
cout << res << endl;
}
int main()
{
init(maxn-5);
cin >> n >> m;
solve(n, m);
return 0;
}
