bzoj 3994/luogu P3327 约数个数和

~{→看不见LaTex格式的Dalao们请刷新本页Thanks♪(･ω･)ﾉ←}~

2019寒假集训题

D.约数个数和

题面

设d(x)为x的约数个数，给定N、M，求 $\sum_{i=1}^{N}\sum_{j=1}^{M}d(i×j)$.

Input

输入文件包含多组测试数据。

第一行，一个整数T，表示测试数据的组数。

接下来的T行，每行两个整数N、M。

Output

T行，每行一个整数，表示你所求的答案。

Sample Input

2
7 4
5 6

Sample Output

110

121

Hint

$1\leq N,M \leq 50,000$

$1\leq T \leq 50,000$

哎，今天fhj大佬讲课也是很无奈，莫比乌斯反演虽然是学了，但是确实不太会用哇。

公式恐惧症患者一枚……

第一步

首先，我们证明一条结论： $d(n×m)=\sum_{i|n}\sum_{j|m}[gcd(i,j)==1]$

证明：对于nm每一个质因子$p_i$，我们令$n={p_i}^{a_i}×n_x,m={p_i}^{b_i}×m_x$，则有 $d(n×m)=\Pi_i(a_i+b_i+1)$ ，设 $l$ 和 $t$ 互质且不含质因子$p_i$，那么$(l×{p_i}^{a_i},t),(l×{p_i}^{a_i-1})…(l,t)…(l,t×{p_i}^{b_i-1}),(l,t×{p_i}^{b_i})$都是互质的数对，而 $ l$ 和 $ t $ 的取法刚好是$\Pi_{k\neq i}(a_k+b_k+1)$种。

第二步

下面是公式恐惧症的“福利”？！

$\sum_{i=1}^n\sum_{j=1}^md(i×j)=\sum_{i=1}^n\sum_{j=1}^m\sum_{k|i}\sum_{l|j}[gcd(k,l)==1]$

现在我们枚举每一个数对$(k,l)$,那么它们会被算多少遍呢？每一个$k$的倍数都会算一次$k$，每一个$l$的倍数又会算一次$l$,而用$\frac nk$来表示n以内k的倍数有多少个是显然的，所以：

$\sum_{k=1}^n\sum_{l=1}^m\frac nk×\frac ml[gcd(k,l)==1]$

然后我们能够用莫比乌斯函数的性质：

$\sum_{d|n}\mu(d)=1\ (n==1)$ $\sum_{d|n}\mu(d)=0\ (n>1)$

那我们就可以得到：

$\sum_{k=1}^n\sum_{l=1}^m\frac nk×\frac ml\sum_{d|gcd(k,l)}\mu(d)$

如果$gcd(k,l)==1$，此时后面的这些等于1，否则是0。

我们就能得到：

$\sum_{k=1}^n\sum_{l=1}^m\frac nk× \frac ml\sum_{d|k且d|l}\mu(d)$

我们只要枚举$\mu(d)$，显然，我们只要枚举d的倍数来与$\mu(d)$相乘即可。

$\sum_{d=1}^{min(n,m)}\mu(d)\sum_{i=1}^\frac nd\sum_{j=1}^\frac md\frac n{d_i}×\frac m{d_j}$

第三步

我们由第二步继续变形可得：

$\sum_{d=1}^{min(n,m)}\mu(d)×(\sum_{k=1}^{\frac nd}\frac n{dk})×(\sum_{l=1}^{\frac md}\frac m{dl})$

令$g(x)=\sum_{i=1}^x\frac xi$

$\sum_{d=1}^{min(n,m)}\mu(d)×g(\frac nd)×g(\frac md)$

我们预处理出$g(x)$即可。

观察一下$g(x)$，就会发现其实式子可以这么理解：$x$以内有约数i的数的个数和，也就是可看作原题中$d(x)$的前缀和
要求$d(x)$的前缀和，可以求出每一个$d(x)$,也就是用线性筛求啦，要记录一下$c(x)$:将$x$分解质因数后，最小质因数的指数。
那么参见第一步，x的最小质因数对$d(x)$的贡献为$c(x)+1$。

第四步

代码实现

#include<bits/stdc++.h>
using namespace std;
int read()
{
	int x = 0, w = 0;
	char ch = 0;
	while (!isdigit(ch)) {w |= ch == '-'; ch = getchar();}
	while (isdigit(ch)) {x = (x << 3) + (x << 1) + ch - '0'; ch = getchar();}
	return w ? -x : x;
}
int mu[500010],p[50010],cnt,n,m;
long long a[500010],b[500010];
bool notprime[500010];
void Prime()
{
	mu[1]=a[1]=b[1]=1;
	for(int i=2;i<=50000;++i)
	{
		if(!notprime[i])
		{
			p[++cnt]=i;
			mu[i]=-1,a[i]=1,b[i]=2;
		}
		for(int j=1;j<=cnt&&i*p[j]<=50000;++j)
		{
			notprime[i*p[j]]=1;
			if(i%p[j]==0)
			{
				b[i*p[j]]=b[i]/(a[i]+1)*(a[i]+2);
				a[i*p[j]]=a[i]+1;
				break;
			}
			else
			{
				b[i*p[j]]=b[i]*b[p[j]];
				a[i*p[j]]=1;
				mu[i*p[j]]=-mu[i];
			}
		}
	}
	for(int i=2;i<=50000;++i)
	{
		b[i]+=b[i-1];
		mu[i]+=mu[i-1];
	}
}
int main()
{
	int T=read();
	Prime();
	while(T--)
	{
		long long result=0;
		n=read(),m=read();
		if(n>m) swap(n,m);
		for(int i=1,j;i<=n;i=j+1)
		{
			j=min(n/(n/i),m/(m/i));
			result+=(mu[j]-mu[i-1])*b[n/i]*b[m/i];
		}
		cout<<result<<endl;
	}
	return 0;
}

p.s.

还有很重要的一点，如果初始化用这样a[N]={0,1}，会导致exe文件过大而无法提交，嘤嘤嘤，感谢vhosc巨佬……(/≧▽≦)/