BZOJ-2154 Crash的数字表格
解法
请注意！下文中默认 (nleq m) ！！！
一步一步推：
[egin{aligned}
f(n,m)=&sumlimits_{i=1}^nsumlimits_{j=1}^m operatorname{lcm}(i,j)
\f(n,m)=&sumlimits_{i=1}^nsumlimits_{j=1}^m frac{icdot j}{gcd(i,j)}
\f(n,m)=&sumlimits_{i=1}^nsumlimits_{j=1}^m gcd(i,j)cdotfrac{i}{gcd(i,j)}cdotfrac{j}{gcd(i,j)}
end{aligned}]
我们调整一下枚举顺序可得（(d) 表示 (gcd)）：
[egin{aligned}
f(n,m)=&sumlimits_{d=1}^n dcdot sum_{i=1}^{lfloorfrac{n}{d}
floor}sum_{j=1}^{lfloorfrac{m}{d}
floor} [gcd(i,j)=1]icdot j
end{aligned}]
我们先不管前面的部分，设 (g(n,m)=sumlimits_{i=1}^{n}sumlimits_{j=1}^{m} [gcd(i,j)=1]icdot j) 我们看到 ([gcd(i,j)=1]) 自然可以想到用 (varepsilon) 表示，再用 (muast 1) 来表示：
[egin{aligned}
g(n,m)=&sumlimits_{i=1}^{n}sumlimits_{j=1}^{m} varepsilon(gcd(i,j))cdot icdot j
\g(n,m)=&sumlimits_{i=1}^{n}sumlimits_{j=1}^{m} (muast 1)(gcd(i,j))cdot icdot j
\g(n,m)=&sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}sumlimits_{dmid gcd(i,j)} mu(d)cdot icdot j
end{aligned}]
故技重施，我们还是先枚举 (d)：
[egin{aligned}
g(n,m)=&sumlimits_{d=1}^nmu(d)cdotsumlimits_{i=1}^{lfloorfrac{n}{d}
floor}sumlimits_{j=1}^{lfloorfrac{m}{d}
floor} (dcdot i)cdot(dcdot j)
\g(n,m)=&sumlimits_{d=1}^nmu(d)cdot d^2cdotsumlimits_{i=1}^{lfloorfrac{n}{d}
floor}sumlimits_{j=1}^{lfloorfrac{m}{d}
floor} icdot j
end{aligned}]
我们再次剥离式子，设 (h(n,m)=sumlimits_{i=1}^{n}sumlimits_{j=1}^{m} icdot j)
[egin{aligned}
h(n,m)=&(sumlimits_{i=1}^{n} i)cdot(sumlimits_{j=1}^{m} j)
\h(n,m)=&frac{n*(n+1)}{2}cdotfrac{m*(m+1)}{2}
end{aligned}]
好的我们的式子已经推到了可以用 (O(1)) 时间可以求解的式子了，我们再来一步一步带回去。
[g(n,m)=sumlimits_{d=1}^nmu(d)cdot d^2cdot h(lfloorfrac{n}{d}
floor,lfloorfrac{m}{d}
floor)
]
可以用数论分块解决，(mu(d)cdot d^2) 可以用前缀和来维护。
[f(n,m)=sumlimits_{d=1}^n dcdot g(lfloorfrac{n}{d}
floor,lfloorfrac{m}{d}
floor)
]
真巧，还是数论分块。总时间复杂度 (O(n+m))
//Don't act like a loser.
//You can only use the code for studying or finding mistakes
//Or,you'll be punished by Sakyamuni!!!
#include<bits/stdc++.h>
#define int long long
using namespace std;
int read() {
char ch=getchar();
int f=1,x=0;
while(ch<'0'||ch>'9') {
if(ch=='-')
f=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9') {
x=x*10+ch-'0';
ch=getchar();
}
return f*x;
}
const int maxn=1e7+10,mod=20101009;
int n,m,mu[maxn],p[maxn/3],cnt,s[maxn];
bool is[maxn];
void sieve(int n) {
mu[1]=1;
for(int i=2;i<=n;i++) {
if(!is[i]) {
p[++cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt;j++) {
if(i*p[j]>n) {
break;
}
is[i*p[j]]=1;
if(i%p[j]==0) {
mu[i*p[j]]=0;
break;
}
mu[p[j]*i]=-mu[i];
}
}
for(int i=1;i<=n;i++) {
s[i]=(s[i-1]+(mu[i]+mod)%mod*i%mod*i%mod)%mod;
}
}
int h(int n,int m) {
return (n*(n+1)/2)%mod*(m*(m+1)/2%mod)%mod;
}
int g(int n,int m) {
int ret=0,j;
for(int i=1;i<=n;i=j+1) {
j=min(n/(n/i),m/(m/i));
ret+=(s[j]-s[i-1]+mod)%mod*h(n/i,m/i)%mod;
ret%=mod;
}
return ret;
}
int f(int n,int m) {
int ret=0,j;
for(int i=1;i<=n;i=j+1) {
j=min(n/(n/i),m/(m/i));
ret+=(j-i+1)*(i+j)/2%mod*g(n/i,m/i)%mod;
ret%=mod;
}
return ret;
}
signed main() {
cin>>n>>m;
if(n>m) {
swap(n,m);//这一步至关重要！！！
}
sieve(n);
printf("%lld
",f(n,m));
return 0;
}