[SDOI2015]约数个数和（莫比乌斯反演）

problem

• 设 \(d(x)\) 为 \(x\) 的约数个数，求：$$ans=\sum_{i=1}^{{n}\sum_{j=1}}{m}d(ij)$$
• 每个读入文件有 \(T\) 组测试数据，\(T,n,m≤50000\)。

Solution

• 众所周知：

• \[d(ij)=\sum_{x|i}\sum_{y|j}[gcd(x,y)=1] \]

• 证明：
• 考虑每个数对 \((x,y)\) 对应的约数是什么。
• 因为 \(gcd(x,y)=1\)，所以对于任意一个质因子 \(p\)，要么 \(x\) 含有 \(p\)，要么 \(y\) 含有 \(p\)，不可能 \(x,y\) 同时含有 \(p\)。
• 如果包含 \(p\) 的是 \(x\) 且 \(x\)含有 \(a\) 个质因子，那么对应约数中含有 \(a\) 个质因子 \(p\)。
• 否则设 \(i\) 含有 \(b\) 个质因子 \(p\)，\(y\) 含有 \(c\) 个质因子 \(p\)，那么对应约数中含有 \(b+c\) 个质因子 \(p\)。
• 因此数对 \((x,y)\) 一一对应了 \(ij\) 的每个约数。
• 那么下面开始推式子：

• \[ans=\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{x|i}\sum_{y|j}[gcd(x,y)=1] \]

• 我们把 \([gcd(x,y)=1]\) 换成含 \(μ\) 的形式，然后我们枚举 \(x,y\) 分别是 \(d\) 的几倍：

• \[ans=\sum_{d=1}^{min(n,m)}μ(d)\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{\lfloor\frac{n}{d}\rfloor}{x}\rfloor\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\lfloor\frac{\lfloor\frac{n}{d}\rfloor}{y}\rfloor \]

• 显然不同的 \(\lfloor\frac{n}{d}\rfloor\) 只有 \(O(\sqrt n)\) 种（把 \(n\) 换成 \(m\) 也一样）。
• 那么我们用 \(o(n\sqrt n)\) 的时间预处理出所有的 \(\sum_{x=1}^{n}\lfloor\frac{\lfloor\frac{n}{d}\rfloor}{x}\rfloor\)，然后对 \(\lfloor\frac{n}{d}\rfloor\) 整除分块，就可以 \(O(T\sqrt n)\) 回答所有询问。

Code

#include <bits/stdc++.h>

using namespace std;

#define ll long long

const int e = 1e5 + 5, lim = 50000;
int n, m;
ll ans, f[e], sum[e], miu[e];
bool b[e];

inline void init()
{
   int i, j;
   miu[1] = 1;
   for (i = 2; i <= lim; i++) miu[i] = 1;
   for (i = 2; i <= lim; i++)
   if (!b[i])
   {
   		miu[i] = -1;
   		for (j = 2 * i; j <= lim; j += i)
   		{
   			if ((j / i) % i == 0) miu[j] = 0;
   			else miu[j] *= -1;
   			b[j] = 1;
   		}
   }
   for (i = 1; i <= lim; i++)
   for (j = i; j <= lim; j += i)
   f[j]++;
   for (i = 1; i <= lim; i++)
   {
   		sum[i] = sum[i - 1] + miu[i];
   		f[i] += f[i - 1];
   }
}

int main()
{
   init();
   int i, tst, j;
   scanf("%d", &tst);
   while (tst--)
   {
   	scanf("%d%d", &n, &m);
   	ans = 0;
   	int tmp = min(n, m);
   	for (i = 1; i <= tmp; i = j + 1)
   	{
   		j = min(n / (n / i), m / (m / i));
   		ans += (sum[j] - sum[i - 1]) * f[n / i] * f[m / i];
   	}
   	printf("%lld\n", ans);
   }
   return 0;
}