[BZOJ 3994] 约数个数和

题意

求下式的值:

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

其中 (d(x)) 为约数个数函数
(n,mle 5 imes 10 ^ 4, qle 5 imes 10^4)

题解

[egin{align} d(ij)&=sum_{a|i}sum_{b|j}[aperp b] \ ext{Ans}&=sum_isum_jd(ij)\ &=sum_isum_jsum_{a|i}sum_{b|j}[iperp j] \ &=sum_isum_jsum_{a|i}sum_{b|j}sum_{k|a,k|b}mu(k)\ &=sum_ksum_i^{lfloor frac n k floor}sum_j^{lfloor frac m k floor}sum_a^{lfloor frac n {ki} floor}sum_b^{lfloor frac m {kj} floor}mu(k) \ &=sum_kmu(k)sum_i^{lfloor frac n k floor}sum_j^{lfloor frac m k floor}sum_a^{lfloor frac n {ki} floor}sum_b^{lfloor frac m {kj} floor}1 \ &=sum_kmu(k)sum_i^{lfloor frac n k floor}sum_j^{lfloor frac m k floor}leftlfloorfrac n {ki} ight floorleftlfloorfrac m {kj} ight floor \ &=sum_kmu(k)sum_i^{lfloor frac n k floor}sum_j^{lfloor frac m k floor}leftlfloorfrac {leftlfloorfrac n k ight floor} {i} ight floorleftlfloorfrac {leftlfloorfrac m k ight floor} {j} ight floor\ %&=sum_ksum_i^{lfloor frac n k floor}sum_j^{lfloor frac m k floor}leftlfloorfrac {leftlfloorfrac n k ight floor} {i} ight floorleftlfloorfrac {leftlfloorfrac m k ight floor} {j} ight floormu(k) &=sum_kmu(k)left(sum_i^{lfloor frac n k floor}leftlfloorfrac {leftlfloorfrac n k ight floor} {i} ight floor ight)left(sum_j^{lfloor frac m k floor}leftlfloorfrac {leftlfloorfrac m k ight floor} {j} ight floor ight)\ end{align} ]

这时我们可以认为 (g(x)=sumlimits_{i=1}^xleftlfloorfrac x i ight floor), 而由于(n,m)炒鸡小于是可以数论分块+记忆化来求 (g(x)), 然后随便筛一筛 (mu) 的前缀和就行了

代码实现

#include <bits/stdc++.h>

const int MAXN=5e4+10;

int cnt;
int mu[MAXN];
int pr[MAXN];
bool npr[MAXN];
long long g[MAXN];

long long Calc(int);
void EulerSieve(int);

int main(){
	int T;
	scanf("%d",&T);
	EulerSieve(5e4);
	while(T--){
		int n,m;
		scanf("%d%d",&n,&m);
		if(n>m)
			std::swap(n,m);
		long long ans=0;
		for(int i=1,j;i<=n;i=j+1){
			j=std::min(n/(n/i),m/(m/i));
			ans+=(mu[j]-mu[i-1])*Calc(n/i)*Calc(m/i);
		}
		printf("%lld
",ans);
	}
	return 0;
}

long long Calc(int x){
	if(g[x]!=0)
		return g[x];
	else{
		for(int i=1,j;i<=x;i=j+1){
			j=x/(x/i);
			g[x]+=(j-i+1)*(x/i);
		}
		return g[x];
	}
}

void EulerSieve(int n){
	npr[0]=npr[1]=true;
	mu[1]=1;
	for(int i=2;i<=n;i++){


		if(!npr[i]){
			pr[cnt++]=i;
			mu[i]=-1;
		}
		for(int j=0,t;j<cnt&&(t=i*pr[j])<=n;j++){
			npr[t]=true;
			if(i%pr[j])
				mu[t]=-mu[i];
			else{
				mu[t]=0;
				break;
			}
		}
	}
	for(int i=1;i<=n;i++)
		mu[i]+=mu[i-1];
}