【SDOI2015】约数个数和

题面

求(sum_{i=1}^nsum_{j=1}^md(ij))

(leq 50000)组数据，(1leq n,mleq 50000)。

题目分析

首先，你需要知道一个结论：

[d(ij)=sumlimits_{x|i}sumlimits_{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)，这样一拿一加，自然会在答案中重复。

---

那么，现在我们的问题转化为求

[sumlimits_{i=1}^nsumlimits_{j=1}^msumlimits_{x|i}sumlimits_{y|j}[gcd(x,y)==1] ]

这样还是无法计算，所以我们把枚举因数提前

[sumlimits_{x=1}^nsumlimits_{y=1}^nlfloorfrac n x floorlfloorfrac m y floor[gcd(x,y)==1] ]

现在看起来就可以反演了，设(f(x))表示(gcd(i,j)==x)时的答案，(g(x))表示(gcd(i,j)==kx,kin Z)时的答案，则：

[egin{split} f(x)&=sumlimits_{i=1}^nsumlimits_{j=1}^mlfloorfrac n i floorlfloorfrac m j floor[gcd(i,j)==x]\ g(x)&=sumlimits_{x|d}^nf(d)\ &=sumlimits_{i=1}^nsumlimits_{j=1}^mlfloorfrac n i floorlfloorfrac m j floor[x|gcd(i,j)]\ &=sumlimits_{i=1}^{lfloorfrac n x floor}sumlimits_{j=1}^{lfloorfrac m x floor}lfloorfrac n {ix} floorlfloorfrac m {jx} floor\ &=sumlimits_{i=1}^{lfloorfrac n x floor}lfloorfrac n {ix} floorsumlimits_{j=1}^{lfloorfrac m x floor}lfloorfrac m {jx} floor end{split} ]

---

比较巧的一点是：(sumlimits_{i=1}^nlfloorfrac n i floor)可以表示(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)=sumlimits_{x|d}^nmu(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<<'
';
	}
	return 0;
}