[国家集训队]Crash的数字表格
解法
首先简化题意
[ans=sum_{i=1}^Nsum_{j=1}^Mlcm(i,j)
]
[ans=sum_{i=1}^Nsum_{j=1}^Mfrac{ij}{gcd(i,j)}
]
枚举一下(gcd)
[ans=sum_{k=1}^Nfrac{1}{k}sum_{i=1}^Nsum_{j=1}^M{ij}[gcd(i,j)=k]
]
把(k)提出来
[ans=sum_{k=1}^Nksum_{i=1}^frac{N}{k}sum_{j=1}^frac{M}{k}{ij}[gcd(i,j)=1]
]
后面那一部分莫比乌斯反演走起
定义
[f(x)=sum_{k=1}^Nksum_{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}^Nksum_{i=1}^frac{N}{k}sum_{j=1}^frac{M}{k}{ij}[gcd(i,j)=d]
]
去掉第一个(sum)
[g(x)=sum_{k=1}^Nksum_{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^2sum_{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^2sum_{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^2sum_{i=1}^frac{N}{kd}sum_{j=1}^frac{M}{kd}ij
]
[ans=sum_{d=1}^Nμ(d)sum_{k=1}^Nkd^2sum_{i=1}^frac{N}{kd}isum_{j=1}^frac{M}{kd}j
]
后面两个是等差数列求和
[ans=sum_{d=1}^Nμ(d)sum_{k=1}^Nkd^2frac{(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}^NDsum_{d|D}μ(d)dfrac{(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
",ans);
}
}