题面
求\(\sum_{i=1}^n\sum_{j=1}^md(ij)\)
\(\leq 50000\)组数据,\(1\leq n,m\leq 50000\)。
题目分析
首先,你需要知道一个结论:
\[ d(ij)=\sum\limits_{x|i}\sum\limits_{y|j}[gcd(x,y)==1] \]
你可以认为\(x,y\)表示你选择的因数为\(\frac i x \cdot y\),即:\(x\)表示\(i\)中不要的部分,\(y\)表示\(j\)中要的部分。
如果\(gcd(x,y)==p_i\),那么\(\frac i x\)表示在约数中拿掉\(p_i\),\(y\)表示在约数中加入\(p_i\),这样一拿一加,自然会在答案中重复。
那么,现在我们的问题转化为求
\[ \sum\limits_{i=1}^n\sum\limits_{j=1}^m\sum\limits_{x|i}\sum\limits_{y|j}[gcd(x,y)==1] \]
这样还是无法计算,所以我们把枚举因数提前
\[ \sum\limits_{x=1}^n\sum\limits_{y=1}^n\lfloor\frac n x\rfloor\lfloor\frac m y\rfloor[gcd(x,y)==1] \]
现在看起来就可以反演了,设\(f(x)\)表示\(gcd(i,j)==x\)时的答案,\(g(x)\)表示\(gcd(i,j)==kx,k\in Z\)时的答案,则:
\[ \begin{split} f(x)&=\sum\limits_{i=1}^n\sum\limits_{j=1}^m\lfloor\frac n i\rfloor\lfloor\frac m j\rfloor[gcd(i,j)==x]\\ g(x)&=\sum\limits_{x|d}^nf(d)\\ &=\sum\limits_{i=1}^n\sum\limits_{j=1}^m\lfloor\frac n i\rfloor\lfloor\frac m j\rfloor[x|gcd(i,j)]\\ &=\sum\limits_{i=1}^{\lfloor\frac n x\rfloor}\sum\limits_{j=1}^{\lfloor\frac m x\rfloor}\lfloor\frac n {ix}\rfloor\lfloor\frac m {jx}\rfloor\\ &=\sum\limits_{i=1}^{\lfloor\frac n x\rfloor}\lfloor\frac n {ix}\rfloor\sum\limits_{j=1}^{\lfloor\frac m x\rfloor}\lfloor\frac m {jx}\rfloor \end{split} \]
比较巧的一点是:\(\sum\limits_{i=1}^n\lfloor\frac n i\rfloor\)可以表示\(1 \sim n\)的约数个数的前缀和。
约数个数可以在线性筛中预处理,原理如下:
对于\(x=p_1^{a_1}p_2^{a_2}p_3^{a_3}...p_n^{a_n}\),\(x\)的约数个数为\((a_1+1)\cdot(a_2+1)\cdot(a_3+1)\cdot...\cdot(a_n+1)\)
由于在线性筛中,每个数只会被它最小的质因子更新,所以:
如果\(i\%prime[j]==0\),说明\(i\)中含有\(prime[j]\),此时\(x\)中\(prime[j]\)的个数为\(i\)中\(prime[j]\)的个数\(+1\),\(x\)的约数个数=\(i\)的约数个数/(\(i\)中\(prime[j]\)的个数)*(\(i\)中\(prime[j]\)的个数\(+1\));
如果\(i\%prime[j]!=0\),说明\(prime[j]\)在\(x\)中只有\(1\)个,\(x\)的约数个数=\(i\)的约数个数*\(2\)。
这样一来\(g(x)\)可以进行预处理,然后\(O(1)\)计算。
反演得\(f(x)=\sum\limits_{x|d}^n\mu(\frac dx)g(d)\),为了针对多组数据,整除分块即可。
代码实现
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdio>
#include<iomanip>
#include<cstdlib>
#define MAXN 0x7fffffff
typedef long long LL;
const int N=50005;
using namespace std;
inline int Getint(){register int x=0,f=1;register char ch=getchar();while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}return x*f;}
int mu[N],prime[N];
bool vis[N];
int ys[N],lw[N],g[N];
int main(){
mu[1]=g[1]=1;
for(int i=2;i<=50000;i++){
if(!vis[i])prime[++prime[0]]=i,mu[i]=-1,ys[i]=2,lw[i]=1;
for(int j=1;j<=prime[0]&&i*prime[j]<=50000;j++){
vis[i*prime[j]]=1;
if(i%prime[j]==0){
ys[i*prime[j]]=ys[i]/(lw[i]+1)*(lw[i]+2);
lw[i*prime[j]]=lw[i]+1;
break;
}
mu[i*prime[j]]=-mu[i];
ys[i*prime[j]]=ys[i]*2,lw[i*prime[j]]=1;
}
mu[i]+=mu[i-1],g[i]=g[i-1]+ys[i];
}
int T=Getint();
while(T--){
int n=Getint(),m=Getint();
if(n>m)swap(n,m);
LL ans=0;
for(int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans+=1ll*(mu[r]-mu[l-1])*g[n/l]*g[m/l];
}
cout<<ans<<'\n';
}
return 0;
}