[BZOJ 3994] 约数个数和

题意

求下式的值:

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

其中 \(d(x)\) 为约数个数函数
\(n,m\le 5\times 10 ^ 4, q\le 5\times 10^4\)

题解

\[\begin{align} d(ij)&=\sum_{a|i}\sum_{b|j}[a\perp b] \\ \text{Ans}&=\sum_i\sum_jd(ij)\\ &=\sum_i\sum_j\sum_{a|i}\sum_{b|j}[i\perp j] \\ &=\sum_i\sum_j\sum_{a|i}\sum_{b|j}\sum_{k|a,k|b}\mu(k)\\ &=\sum_k\sum_i^{\lfloor \frac n k \rfloor}\sum_j^{\lfloor \frac m k \rfloor}\sum_a^{\lfloor \frac n {ki} \rfloor}\sum_b^{\lfloor \frac m {kj} \rfloor}\mu(k) \\ &=\sum_k\mu(k)\sum_i^{\lfloor \frac n k \rfloor}\sum_j^{\lfloor \frac m k \rfloor}\sum_a^{\lfloor \frac n {ki} \rfloor}\sum_b^{\lfloor \frac m {kj} \rfloor}1 \\ &=\sum_k\mu(k)\sum_i^{\lfloor \frac n k \rfloor}\sum_j^{\lfloor \frac m k \rfloor}\left\lfloor\frac n {ki}\right\rfloor\left\lfloor\frac m {kj}\right\rfloor \\ &=\sum_k\mu(k)\sum_i^{\lfloor \frac n k \rfloor}\sum_j^{\lfloor \frac m k \rfloor}\left\lfloor\frac {\left\lfloor\frac n k\right\rfloor} {i}\right\rfloor\left\lfloor\frac {\left\lfloor\frac m k\right\rfloor} {j}\right\rfloor\\ %&=\sum_k\sum_i^{\lfloor \frac n k \rfloor}\sum_j^{\lfloor \frac m k \rfloor}\left\lfloor\frac {\left\lfloor\frac n k\right\rfloor} {i}\right\rfloor\left\lfloor\frac {\left\lfloor\frac m k\right\rfloor} {j}\right\rfloor\mu(k) &=\sum_k\mu(k)\left(\sum_i^{\lfloor \frac n k \rfloor}\left\lfloor\frac {\left\lfloor\frac n k\right\rfloor} {i}\right\rfloor\right)\left(\sum_j^{\lfloor \frac m k \rfloor}\left\lfloor\frac {\left\lfloor\frac m k\right\rfloor} {j}\right\rfloor\right)\\ \end{align} \]

这时我们可以认为 \(g(x)=\sum\limits_{i=1}^x\left\lfloor\frac x i\right\rfloor\), 而由于\(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\n",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];
}