[SDOI2015]约数个数和
解法
由于(N,M)的顺序对答案没影响
我们规定(N<M)
(ans=sum_{i=1}^Nsum_{j=1}^M d(ij))
给出一个结论
[d(ij)=sum_{x|i}sum_{y|i}[gcd(x,y)=1]
]
所以
[ans=sum_{i=1}^Nsum_{j=1}^Msum_{x|i}sum_{y|i}[gcd(x,y)=1]
]
乱搞一下,改成枚举因子
[ans=sum_{i=1}^Nsum_{j=1}^Mfrac{N}{i}frac{M}{j}[gcd(i,j)=1]
]
显然，来一发莫比乌斯反演
[f(x)=sum_{i=1}^Nsum_{j=1}^Mfrac{N}{i}frac{M}{j}[gcd(i,j)=x]
]
所以(ans=f(1))
[g(x)=sum_{x|d}f(x)
]
[g(x)=sum_{x|d}sum_{i=1}^Nsum_{j=1}^Mfrac{N}{i}frac{M}{j}[gcd(i,j)=d]
]
[g(x)=sum_{i=1}^Nsum_{j=1}^Mfrac{N}{i}frac{M}{j}[gcd(i,j) \% x=0]
]
把(x)提出来我们就可以排除(gcd)的影响
[g(x)=sum_{i=1}^{frac{N}{x}}sum_{j=1}^{frac{M}{x}}frac{N}{ix}frac{M}{jx}[gcd(i,j) \% 1=0]
]
[g(x)=sum_{i=1}^{frac{N}{x}}sum_{j=1}^{frac{M}{x}}frac{N}{ix}frac{M}{jx}
]
根据莫比乌斯反演定理
[f(x)=sum_{x|d}μ(frac{d}{x})g(d)
]
[f(x)=sum_{x|d}μ(frac{d}{x})sum_{i=1}^{frac{N}{d}}sum_{j=1}^{frac{M}{d}}frac{N}{id}frac{M}{jd}
]
而(ans=f(1))所以
[ans=sum_{d=1}^Nμ(d)sum_{i=1}^{frac{N}{d}}sum_{j=1}^{frac{M}{d}}frac{N}{id}frac{M}{jd}
]
[ans=sum_{d=1}^Nμ(d)sum_{i=1}^{frac{N}{d}}frac{N}{id}sum_{j=1}^{frac{M}{d}}frac{M}{jd}
]
先用数论分块(O(nsqrt n))预处理一个$$sum(x)=sum_{i=1}^xfrac{x}{i}$$
于是$$ans=sum_{d=1}^Nμ(d)sum(frac{N}{d})sum(frac{M}{d})$$
显然我们可以把前面一部分预处理前缀和，后面数论分块加上前面的预处理可以(O(sqrt n )的算出每一次的答案)
完整代码
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=1e6,MAX=50000;
int prime[N],isprime[N],miu[N],cnt,sum[N];
void pre(){
miu[1]=1;
for(int i=2;i<=MAX;++i){
if(!isprime[i])miu[i]=-1,prime[++cnt]=i;
for(int j=1;j<=cnt&&i*prime[j]<=MAX;++j){
isprime[i*prime[j]]=1;
if(i%prime[j]==0)break;
miu[i*prime[j]]=-miu[i];
}
}
for(int i=1;i<=MAX;++i)miu[i]+=miu[i-1];
for(int i=1;i<=MAX;++i){
for(int l=1,r;l<=i;l=r+1){
r=i/(i/l);
r=min(r,i);
sum[i]+=(r-l+1)*(i/l);
}
}
}
int main(){
int T=0;
cin>>T;
pre();
while(T--){
int n,m;
ll ans=0;
cin>>n>>m;
if(n>m)swap(n,m);
for(int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans+=1ll*(miu[r]-miu[l-1])*sum[n/l]*sum[m/l];
}
printf("%lld
",ans);
}
return 0;
}