题解-SDOI2015 约数个数和

Problem

题意：设 (d(x)) 为 (x) 的约数个数，给定 (N,M)，求(sum_{i=1}^Nsum_{j=1}^Md(ij))

(1leq N,M,Tleq 5 imes 10^4)

Solution

第一次推出莫反式 ♪(^∇*)

以下部分中小括号代表(gcd)，中括号代表取布尔值

一开始想枚举约数然后看有多少倍数出现过的，发现不好弄，转而想到 (xy) 的因数一定是 (x) 的因数和 (y) 的因数的积，去除冗余后可以得到 (d(xy)=sum_{a|x}sum_{b|y}[(frac xa,frac yb)=1]=sum_{a|x}sum_{b|y}[(a,b)=1])

然后答案为：

[sum_{i=1}^Nsum_{j=1}^Msum_{a|i}sum_{b|j}[(a,b)=1] ]

考虑优先枚举(a,b)

[sum_{a=1}^Nsum_{b=1}^M[(a,b)=1]lfloor frac Na floor lfloor frac Mb floor ]

发现里头那个(gcd)可以套路地搞事了：

构造(f(x)=sum_{a=1}^Nsum_{b=1}^M[(a,b)=x]lfloor frac Na floor lfloor frac Mb floor)

构造(F(x)=sum_{a=1}^Nsum_{b=1}^M[x|(a,b)]lfloor frac Na floor lfloor frac Mb floor)

可以根据定义式得到(F(x)=sum_{x|n}f(n))，经过莫比乌斯反演，得到(f(x)=sum_{x|n}mu(frac nx)F(n))

答案为(f(1)=sum_nmu(n)F(n))

由于(F(n))中可以将(gcd)的限制条件放在前面的(sum)里，所以(F(n)=sum_{n|a}sum_{n|b}lfloor frac Na floor lfloor frac Mb floor)

所以答案为

[f(1)=sum_dmu(d)sum_{d|a}sum_{d|b}lfloor frac Na floor lfloor frac Mb floor ]

发现后面的式子不好求，可以构造(g(x)=sum_{i=1}^xlfloor frac xi floor)

则最终的式子为

[Ans=sum_dmu(d)g(lfloor frac nd floor)g(lfloor frac md floor) ]

分块解决即可

至于如何求(g(x))，可以发现其实(g(n)=sum_{i=1}^nd(i))，利用线性筛完约数函数(d(i))后求个前缀和即可

Code

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

inline void read(int&x){
	char c11=getchar();x=0;while(!isdigit(c11))c11=getchar();
	while(isdigit(c11))x=x*10+c11-'0',c11=getchar();
}

const int N=50101;
int pri[N],is[N];
int u[N],g[N],d[N];
int tp;ll sm;

int main(){
	g[1]=u[1]=1;
	for(int i=2;i<N;++i){
		if(!is[i])pri[++tp]=i,u[i]=-1,g[i]=2,d[i]=1;
		for(int j=1,k;j<=tp and (k=i*pri[j])<N;++j){
			is[k]=1;
			if(i%pri[j])u[k]-=u[i],g[k]=g[i]*2,d[k]=1;
			else {
				u[k]=0,d[k]=d[i]+1;
				g[k]=g[i]/(d[i]+1)*(d[i]+2);
				break;
			}
		}
		g[i]+=g[i-1],u[i]+=u[i-1];
	}
	int n,m,T;read(T);
	while(T--){
		read(n),read(m),sm=0ll;
		if(n>m)swap(n,m);
		for(int i=1,j;i<=n;i=j+1){
			j=min(n/(n/i),m/(m/i));
			sm+=1ll*(u[j]-u[i-1])*g[n/i]*g[m/i];
		}printf("%lld
",sm);
	}return 0;
}