[SDOI2015]约数个数和

题目

解法

由于\(N,M\)的顺序对答案没影响
我们规定\(N<M\)
\(ans=\sum_{i=1}^N\sum_{j=1}^M d(ij)\)

给出一个结论
\[d(ij)=\sum_{x|i}\sum_{y|i}[gcd(x,y)=1]\]
所以
\[ans=\sum_{i=1}^N\sum_{j=1}^M\sum_{x|i}\sum_{y|i}[gcd(x,y)=1]\]
乱搞一下,改成枚举因子
\[ans=\sum_{i=1}^N\sum_{j=1}^M\frac{N}{i}\frac{M}{j}[gcd(i,j)=1]\]
显然，来一发莫比乌斯反演
\[f(x)=\sum_{i=1}^N\sum_{j=1}^M\frac{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}^N\sum_{j=1}^M\frac{N}{i}\frac{M}{j}[gcd(i,j)=d]\]
\[g(x)=\sum_{i=1}^N\sum_{j=1}^M\frac{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(n\sqrt n)\)预处理一个\[sum(x)=\sum_{i=1}^x\frac{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\n",ans);
    }
    return 0;
    
}