[SDOI2015][bzoj 3994][Luogu P3327] 约数个数和 (莫比乌斯反演)

题目描述

设$d(x)$为$x$的约数个数，给定$N$、$M$，求 $sum^{N}_{i=1}sum^{M}_{j=1} d(ij)$
$N,M,T<=50000$

题目分析

首先很不显然的有这样一个结论：
$d(ij)=sum_{x|i}sum_{y|j}[(x,y)==1]$

证明如下

将一个数唯一质因数分解为$p_1^{k_1}p_2^{k_2}p_3^{k_3}...p_n^{k_n}$,其约数个数为$(k_1+1)(k_2+1)(k_3+1)...(k_n+1)$

考虑一个质数$p$对$d(ij)$的影响,假设$i$分解质因数后有$k$个$p$,$j$分解质因数后有$q$个$p$,则产生的贡献即为$k+q+1$,接下来是关键的一步(反正我想不到XD)
$k+q+1=sum_{x=0}^ksum_{y=0}^q[(p^x,p^y)==1]$

(此句可自行忽略:只有当$xy=0$时才会有值,这就相当于一个$(k+1)∗(q+1)$的矩形，只取了左下角的$L$字型的一块)

根据乘法原理,有
$d(ij)=sum_{x_1=0}^{k_1}sum_{y_1=0}^{q_1}[(p_1^{x_1},p_1^{y_1})==1] ewline *sum_{x_2=0}^{k_2}sum_{y_2=0}^{q_2}[(p_2^{x_2},p_2^{y_2})==1] ewline ··· ewline *sum_{x_n=0}^{k_n}sum_{y_n=0}^{q_n}[(p_n^{x_n},p_n^{y_n})==1]$
因为必须满足所有$(p_n^{x_n},p_n^{y_n})==1$才会对答案造成贡献,则可以变形为$(prod_{i=1}^np_n^{x_n},prod_{i=1}^np_n^{y_n})==1$

所以$d(ij)=sum_{x|i}sum_{y|j}[(x,y)==1]$

(以上证明摘自:传送门)

证毕

现在就有了$Ans=sum_{i=1}^{N}sum_{j=1}^{M}sum_{x|i}^{}sum_{y|j}^{}[(x,y)==1]$

根据数论中切换枚举次序的套路以及莫比乌斯反演对布尔条件框的替换，我们得到
$Ans=sum_{x=1}^{N}sum_{y=1}^{M}⌊frac{N}{x}⌋⌊frac{M}{y}⌋[(x,y)==1] ewline =sum_{x=1}^{N}sum_{y=1}^{M}⌊frac{N}{x}⌋⌊frac{M}{y}⌋sum_{k|(x,y)}mu(k) ewline =sum_{k=1}^{min(N,M)}mu(k)sum_{k|x}^{}⌊frac{N}{x}⌋sum_{k|y}^{}⌊frac{M}{y}⌋ ewline =sum_{k=1}^{min(N,M)}mu(k)sum_{x=1}^{⌊frac{N}{k}⌋}⌊frac{N}{kx}⌋sum_{y=1}^{⌊frac{M}{k}⌋}⌊frac{M}{ky}⌋$

于是这里使用分块优化,预处理$mu$的前缀和

最后考虑后面的两个Sigma怎么处理
可以观察发现它们都是$sum_{i=1}^{n}⌊frac{n}{i}⌋$的形式,此处可以用分块优化$Theta(nsqrt{n})$预处理

其实还可以$Theta(n)$预处理,可以发现$sum_{i=1}^{n}⌊frac{n}{i}⌋=sum_{i=1}^{n}sum_{i|d}^{n}1$
显然是约数个数和$d(n)$的前缀和函数,于是线性筛出$mu(n)$与$d(n)$,求出前缀和即可

每次询问$Theta(sqrt{n})$, 预处理$Theta(n)$,总时间复杂度为$Theta(nsqrt{n})$

TIPS

线性筛约数个数时存一下最小的质数$p_1$的次方$+1$,即存下$(k_1+1)$就可以方便的线性筛了

线性筛约数和同理,存一下最小的质数$p_1$所对应的首项为$1$,公比为$p_1$,项数为$k_1+1$的等比数列,即存下$sum_{i=1}^{k_1}p_1^i$就可以愉快的线性筛了

AC code

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

const int MAXN = 100001;
namespace Mobius
{
	#define LL long long
	int Prime[MAXN], Cnt, mu[MAXN], d[MAXN], Minseq[MAXN];
	bool IsnotPrime[MAXN];
	LL sum_MU[MAXN], sum_D[MAXN];
	void init()
	{
		mu[1] = d[1] = Minseq[1] = 1;
		for(int i = 2; i < MAXN; ++i)
		{
			if(!IsnotPrime[i])
				Prime[++Cnt] = i, mu[i] = -1, d[i] = Minseq[i] = 2;
			for(int j = 1; j <= Cnt && i * Prime[j] < MAXN; ++j)
			{
				IsnotPrime[i * Prime[j]] = 1;
				if(i % Prime[j] == 0)
				{
					Minseq[i * Prime[j]] = Minseq[i]+1;
					mu[i * Prime[j]] = 0;
					d[i * Prime[j]] = d[i] / Minseq[i] * Minseq[i * Prime[j]];
					break;
				}
				Minseq[i * Prime[j]] = 2;
				mu[i * Prime[j]] = -mu[i];
				d[i * Prime[j]] = d[i]<<1;
			}
		}
		for(int i = 1; i < MAXN; ++i)
			sum_MU[i] = sum_MU[i-1] + mu[i],
			sum_D[i] = sum_D[i-1] + d[i];
	}
	LL calc(int N, int M)
	{
		if(N > M) swap(N, M);
		LL ret = 0;
		for(int i = 1, j; i <= N; i=j+1)
		{
			j = min(N/(N/i), M/(M/i));
			ret += (sum_MU[j]-sum_MU[i-1]) * sum_D[N/i] * sum_D[M/i];
		}
		return ret;
	}
}
using namespace Mobius;
int main ()
{
	init();
	int T, n, m;
	scanf("%d", &T);
	while(T--)
	{
		scanf("%d%d", &n, &m);
		printf("%lld
", calc(n, m));
	}
}