BZOJ 3994 约数个数和

Description

设(d(x))为(x)的约数个数，给定(N,M)，求$$sum_{i=1}^{N}sum_{j=1}{M}d(ij)$$。

Input

输入文件包含多组测试数据。
第一行，一个整数(T)，表示测试数据的组数。
接下来的(T)行，每行两个整数(N,M)。

Output

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

Sample Input

2
7 4
5 6

Sample Output

110
121

HINT

(1 le N, M le 50000)
(1 le T le 50000)

这题有个很屌的结论：$$sum_{i=1}^{N}sum_{j=1}{M}d(ij)=sum_{i=1}^{N}sum_{j=1}{M}lfloorfrac{N}{i} floorlfloorfrac{M}{j} floorlbrack gcd(i,j)=1 brack$$
根据PoPoQQQ博客所说的，我们可以先证明这个式子的成立：$$d(nm)=sum_{i mid n}sum_{j mid m}lbrack gcd(i,j)=1 brack$$
我们可以证明一下：我们对每个质数(p)单独算贡献，设(n=n’ imes p^{k_{1}})，(m=m’ imes p^{k_{2}})。那么，该质数(p)对答案的贡献显然为(k_{1}+k_{2}+1)。于是我们考虑$$d(nm)=sum_{i mid n}sum_{j mid m}lbrack gcd(i,j)=1 brack$$这个式子，发现(p)对之有贡献的数对((i,j))仍然是$$(p^{{k_{1}},1),(p}{k_{1}-1},1) cdots (1,1) cdots (1,p^{{k_{2}-1}),(1,p}{k_{2}})$$这(k_{1}+k_{2}+1)个，因此得证。
代入得$$ sum_{n = 1}^{N}sum_{m = 1}^{M}d(nm) = sum_{n = 1}^{N}sum_{m = 1}^{M}sum_{i mid n} sum_{j mid m}[gcd(i,j)=1]$$
我们转变枚举量，先枚举(i,j)就有

[sum_{i=1}^{N}sum_{j=1}^{M}lfloorfrac{N}{i} floorlfloorfrac{M}{j} floorlbrack gcd(i,j)=1 brack ]

于是$$sum_{i=1}^{N}sum_{j=1}{M}lfloorfrac{N}{i} floorlfloorfrac{M}{j} floorlbrack gcd(i,j)=1 brack$$这个式子我们可以上反演了。
反演化为$$sum_{i=1}^{N}sum_{j=1}{M}lfloorfrac{N}{i} floorlfloorfrac{M}{j} floor sum_{g mid i;g mid j} mu(g)$$
转而枚举(g)，于是就可得到$$sum_{g=1}^{{N}mu(g)sum_{i=1}}{lfloor frac{N}{g} floor}sum_{j=1}^{lfloor frac{M}{g} floor}lfloor frac{N}{ig} floorlfloor frac{M}{jg} floor$$
再化一下就可得到$$sum_{g=1}^{{N}mu(g)sum_{i=1}}{lfloor frac{N}{g} floor}lfloor frac{N}{ig} floorsum_{j=1}^{lfloor frac{M}{g} floor}lfloor frac{M}{jg} floor$$
又有$$lfloor frac{N}{ab} floor=lfloor frac{lfloor frac{N}{a} floor}{b} floor$$
于是我们发现(sum_{i=1}^{lfloor frac{N}{g} floor}lfloor frac{N}{ig} floor)只与(lfloor frac{N}{g} floor)有关，我们可以(O(n sqrt{n}))预处理$$f_{x}=sum_{i=1}^{x}lfloor frac{x}{i} floor$$
有了这个后再化简$$sum_{g=1}^{N}mu(g)f_{lfloor frac{N}{g} floor}f_{lfloor frac{M}{g} floor}$$就可在(O(sqrt{n}))分段求了。皆大欢喜。

#include<iostream>
#include<cstdio>
#include<cstdlib>
using namespace std;

typedef long long ll;
#define maxn (50010)
int f[maxn],mu[maxn],prime[maxn],n,m,tot; bool exist[maxn];

inline int calc(int x)
{
	int ret = 0;
	for (int i = 1,last;i <= x;i = last+1)
	{
		last = min(x,x/(x/i));
		ret += (x/i)*(last-i+1);
	}
	return ret;
}
inline void ready()
{
	mu[1] = 1;
	for (int i = 2;i <= 50000;++i)
	{
		if (!exist[i]) { prime[++tot] = i; mu[i] = -1; }
		for (int j = 1;j <= tot&&prime[j]*i <= 50000;++j)
		{
			exist[i*prime[j]] = true;
			if (i % prime[j] == 0) { mu[i*prime[j]] = 0; break; }
			mu[i*prime[j]] = -mu[i];
		}
	}
	for (int i = 1;i <= 50000;++i) mu[i] += mu[i-1],f[i] = calc(i);
}

inline ll work()
{
	if (n > m) swap(n,m);
	ll ret = 0;
	for (int i = 1,last;i <= n;i = last+1)
	{
		last = min(n,min(n/(n/i),m/(m/i)));
		ret += (ll)(mu[last]-mu[i-1])*((ll)f[n/i]*f[m/i]);
	}
	return ret;
}

int main()
{
	freopen("3994.in","r",stdin);
	freopen("3994.out","w",stdout);
	ready();
	int T; scanf("%d",&T);
	while (T--) scanf("%d %d",&n,&m),printf("%lld
",work());
	fclose(stdin); fclose(stdout);
	return 0;
}