[国家集训队]Crash的数字表格
[国家集训队]Crash的数字表格
解法
首先简化题意
\[ans=\sum_{i=1}^N\sum_{j=1}^Mlcm(i,j)
\]
\[ans=\sum_{i=1}^N\sum_{j=1}^M\frac{ij}{gcd(i,j)}
\]
枚举一下\(gcd\)
\[ans=\sum_{k=1}^N\frac{1}{k}\sum_{i=1}^N\sum_{j=1}^M{ij}[gcd(i,j)=k]
\]
把\(k\)提出来
\[ans=\sum_{k=1}^Nk\sum_{i=1}^\frac{N}{k}\sum_{j=1}^\frac{M}{k}{ij}[gcd(i,j)=1]
\]
后面那一部分莫比乌斯反演走起
定义
\[f(x)=\sum_{k=1}^Nk\sum_{i=1}^\frac{N}{k}\sum_{j=1}^\frac{M}{k}{ij}[gcd(i,j)=x]
\]
\[g(x)=\sum_{x|d}f(d)
\]
所以\(ans=f(1)\)
\[g(x)=\sum_{x|d}\sum_{k=1}^Nk\sum_{i=1}^\frac{N}{k}\sum_{j=1}^\frac{M}{k}{ij}[gcd(i,j)=d]
\]
去掉第一个\(\sum\)
\[g(x)=\sum_{k=1}^Nk\sum_{i=1}^\frac{N}{k}\sum_{j=1}^\frac{M}{k}{ij}[gcd(i,j)\%x=0]
\]
提出\(x\)，排除\(gcd\)的影响
\[g(x)=\sum_{k=1}^Nkx^2\sum_{i=1}^\frac{N}{kx}\sum_{j=1}^\frac{M}{kx}{ij}
\]
根据莫比乌斯定理
\[f(x)=\sum_{x|d}μ(\frac{d}{x})g(d)
\]
\[f(x)=\sum_{x|d}μ(\frac{d}{x})\sum_{k=1}^Nkd^2\sum_{i=1}^\frac{N}{kd}\sum_{j=1}^\frac{M}{kd}{ij}
\]
而\(ans=f(1)\)
\[ans=\sum_{d=1}^Nμ(d)\sum_{k=1}^Nkd^2\sum_{i=1}^\frac{N}{kd}\sum_{j=1}^\frac{M}{kd}ij
\]
\[ans=\sum_{d=1}^Nμ(d)\sum_{k=1}^Nkd^2\sum_{i=1}^\frac{N}{kd}i\sum_{j=1}^\frac{M}{kd}j
\]
后面两个是等差数列求和
\[ans=\sum_{d=1}^Nμ(d)\sum_{k=1}^Nkd^2\frac{(\frac{N}{kd})(\frac{N}{kd}+1)(\frac{M}{kd})(\frac{M}{kd}+1)}{4}
\]
(我一开始做的时候是在这里就结束了，两次数论分块时间复杂度变\(O(n)\),但后来做了一道一样的但是有多组询问的题目，所以就有了以下的优化)
设\(D=kd\)，枚举\(D\),考虑到若\(D>N\)，对答案没有贡献,所以\(D\)只要枚举到\(N\)
\[ans=\sum_{D=1}^ND\sum_{d|D}μ(d)d\frac{(\frac{N}{D})(\frac{N}{D}+1)(\frac{M}{D})(\frac{M}{D}+1)}{4}
\]
中间那个线性筛一下
然后前缀和加数论分块是\(O(\sqrt n)\)的。
完整代码(未优化的)(下面有优化过的)
#include<bits/stdc++.h>
#define ll long long
#define il inline
#define rg register
using namespace std;
const ll N=1e7+20,p=20101009;
__int128 ljq=1;
ll pre[N],prime[N],n,m,s;
bool isprime[N];
void work(){
cin>>n>>m;if(n>m)swap(n,m);
for(rg ll i=2;i<=n;++i){
if(!isprime[i])pre[i]=-1,prime[++s]=i;
for(rg ll j=1;j<=s&&i*prime[j]<=n;++j){
isprime[i*prime[j]]=1;
if(i%prime[j])pre[i*prime[j]]=-pre[i];
else break;
}
}
pre[1]=1;
for(int i=1;i<=n;++i)
pre[i]=(ljq*pre[i]*i*i+pre[i-1])%p;
}
ll js(rg int x,rg int y){
ll ans=0;
for(rg int l=1,r;l<=x;l=r+1){
r=min(x/(x/l),y/(y/l));
ans+=(ljq*(pre[r]-pre[l-1])*(x/l+1)*(x/l)/2*(y/l)*(y/l+1)/2)%p;
ans%=p;
}
return ans;
}
int main(){
work();
ll ans=0;
for(int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans+=(ljq*(l+r)*(r-l+1)/2*js(n/l,m/l))%p;
ans%=p;
}
cout<<((ans%p)+p)%p<<endl;
return 0;
}
优化代码(注意模数不一样)
#include<bits/stdc++.h>
#define ll long long
#define rg register
using namespace std;
const int N=1e7+20,yl=1e8+9,M=1e7+10;
ll dis[N],isprime[N],cnt,prime[N];
void eular(){
for(rg int i=2;i<=M;++i){
if(!isprime[i])prime[++cnt]=i,dis[i]=1-i;
for(int j=1;j<=cnt&&i*prime[j]<=M;++j){
isprime[i*prime[j]]=1;
if(i%prime[j])dis[i*prime[j]]=dis[i]*dis[prime[j]]%yl;
else {dis[i*prime[j]]=dis[i];break;}
}
}
dis[1]=1;
for(int i=1;i<=M;++i){
dis[i]=(dis[i]*i+dis[i-1])%yl;
}
}
int main(){
eular();
int T;
cin>>T;
while(T--){
int n,m;
scanf("%d%d",&n,&m);
if(n>m)swap(n,m);
ll ans=0;
for(ll l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans=(ans+(dis[r]-dis[l-1])*((n/l)*(n/l+1)/2%yl)%yl*((m/l)*(m/l+1)/2%yl))%yl;
}
if(ans<0)ans+=yl;
printf("%lld\n",ans);
}
}
