【题解】Luogu-P4240 毒瘤之神的考验

P4240 毒瘤之神的考验

Description

• 多测，$t$ 组数据。

• 每次给定两个整数 $n,m$，请求出
  $$\left[\sum _{i=1}^n\sum _{j=1}^m\phi (ij)\right]\mathrm{m}\mathrm{o}\mathrm{d}998244353$$

• $1\le t\le 10^4$，$1\le n,m\le 10^5$。

Solution

引理：
$$\phi (ij)=\frac{\phi (i)\phi (j)\gcd (i,j)}{\phi (\gcd (i,j))}$$

证明：
$$\begin{array}{rl}\phi (i)\phi (j)\gcd (i,j)& =i\prod _{p_1\in \mathbb{P},p_1\mid i}\frac{p_1-1}{p_1}\cdot j\prod _{p_2\in \mathbb{P},p_2\mid j}\frac{p_2-1}{p_2}\cdot \gcd (i,j)\\ & =ij\gcd (i,j)\left(\prod _{p_1\in \mathbb{P},p_1\mid i}\frac{p_1-1}{p_1}\prod _{p_2\in \mathbb{P},p_2\mid j}\frac{p_2-1}{p_{@}}\right)\end{array}$$
你会发现所有 $\gcd (i,j)$ 的质因数都被算了 $2$ 次。
$$\begin{array}{rl}& =\left(ij\prod _{p_1\in \mathbb{P},p_1\mid ij}\frac{p_1-1}{p_1}\right)\cdot \left(\gcd (i,j)\prod _{p_2\in \mathbb{P},p_2\mid \gcd (i,j)}\frac{p_2-1}{p_2}\right)\\ & =\phi (ij)\phi (\gcd (i,j))\end{array}$$
把 $\phi (\gcd (i,j))$ 除过去即可。

不妨设 $n\le m$，根据引理，有
$$\begin{array}{rl}ans& =\sum _{i=1}^n\sum _{j=1}^m\frac{\phi (i)\phi (j)\gcd (i,j)}{\phi (\gcd (i,j))}\\ & =\sum _{d=1}^n\sum _{i=1}^n\sum _{j=1}^m\frac{\phi (i)\phi (j)d}{\phi (d)}[\gcd (i,j)=d]\\ & =\sum _{d=1}^n\frac{d}{\phi (d)}\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }\phi (id)\phi (jd)[\gcd (i,j)=1]\\ & =\sum _{d=1}^n\frac{d}{\phi (d)}\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\phi (id)\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }\phi (jd)[\gcd (i,j)=1]\\ & =\sum _{d=1}^n\frac{d}{\phi (d)}\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\phi (id)\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }\phi (jd)\sum _{k\mid \gcd (i,j)}\mu (k)\\ & =\sum _{d=1}^n\frac{d}{\phi (d)}\sum _{k=1}^{\lfloor \frac{n}{d}\rfloor }\mu (k)\sum _{i=1}^{\lfloor \frac{n}{dk}\rfloor }\phi (idk)\sum _{j=1}^{\lfloor \frac{m}{dk}\rfloor }\phi (jdk)\\ & =\sum _{T=1}^n\sum _{d\mid T}\frac{d}{\phi (d)}\mu \left(\frac{T}{d}\right)\sum _{i=1}^{\lfloor \frac{n}{T}\rfloor }\phi (iT)\sum _{j=1}^{\lfloor \frac{m}{T}\rfloor }\phi (jT)\end{array}$$
其中
$$f(n)=\sum _{d\mid n}\frac{d\mu \left(\frac{n}{d}\right)}{\phi (d)}$$
不是积性函数（它甚至不一定是整数），直接预处理逆元 $\Theta (n\ln n)$ 计算。

对于
$$\sum _{i=1}^{\lfloor \frac{n}{T}\rfloor }\phi (iT)\sum _{j=1}^{\lfloor \frac{m}{T}\rfloor }\phi (jT)$$
发现这东西带两个参数，令
$$g(k,n)=\sum _{i=1}^n\phi (ik)$$
你有递推式
$$g(k,n)=g(k,n-1)+\phi (nk)$$
而此处的 $nk\le n\text{(题目中的 }n\text{)}$，也可以 $\Theta (n\ln n)$ 预处理，$g$ 数组用 vector 存。

代回去
$$ans=\sum _{T=1}^{n}f(T)g\left(T,\lfloor \frac{n}{T}\rfloor \right)g\left(T,\lfloor \frac{m}{T}\rfloor \right)$$
这个 $g$ 带一个 $T$，无法整除分块，只能暴力算。

此时是 $\mathrm{O}(n\ln n)$ 的预处理，$\mathrm{O}(n)$ 的单次回答，过不去。

---

设
$$h(n,m,l)=\sum _{T=1}^{l}f(T)g\left(T,n\right)g\left(T,m\right)$$
这样在整除分块中
$$h\left(\lfloor \frac{n}{l}\rfloor ,\lfloor \frac{m}{l}\rfloor ,r\right)-h\left(\lfloor \frac{n}{l}\rfloor ,\lfloor \frac{m}{l}\rfloor ,l-1\right)$$
又有递推式
$$h(n,m,l)=h(n,m,l-1)+f(l)g(l,n)g(l,m)$$
有 $\max (n,m)l\le n\text{(题目中的 }n\text{)}$。

$h$ 数组同样用 vector 存。

此时是
$$\begin{array}{rl}\mathrm{O}\left(n\sum _{i=1}^n\frac{n}{i}\right)& =\mathrm{O}\left(n\cdot {\int }_1^n\frac{n}{x}dx\right)\\ & =\mathrm{O}(n^2\ln n)\end{array}$$
的预处理，$\mathrm{O}(\sqrt{n})$ 的单次回答，预处理会 $\text{MLE}+\text{TLE}$。

---

既然两种方法都不行，那么自然想到 用根号分治来平衡时间。

设阈值为 $k$，表示对于 $\le k$ 的部分预处理。

于是当 $\lfloor \frac{m}{l}\rfloor \le k$，即 $l>\lfloor \frac{m}{k}\rfloor$ 的部分用已经预处理过的 $h$ 整除分块算，当 $l\le \lfloor \frac{m}{k}\rfloor$ 时直接暴力算。

• 暴力：$\mathrm{O}(n\ln n)+\mathrm{O}\left(\frac{tn}{k}\right)$
• 分块：$\mathrm{O}(nk\ln k)+\mathrm{O}(t\sqrt{n})$

总时间复杂度为 $\mathrm{O}\left(nk\ln k+t\left(\frac{n}{k}+\sqrt{n}\right)\right)$，提个公因数就是 $k\ln k+\frac{t}{k}$ 最小。

写个程序去算

	int t = 1e4;
	int ans = 0x3f3f3f3f, pos = 0;
	for (int k = 1; k <= t; k++)
	{
		int res = k * log(k) + t / k;
		if (res < ans)
		{
			ans = res, pos = k;
		}
	}
	printf("%d %d\n", ans, pos);

输出

392 47

所以取 $k=47$。

空间复杂度为 $\mathrm{O}(nk\ln k)$。

然而实测 $k$ 取 $60\sim 70$ 更快（

Code

// 18 = 9 + 9 = 18.
#include <iostream>
#include <cstdio>
#include <vector>
#define Debug(x) cout << #x << "=" << x << endl
typedef long long ll;
using namespace std;

const int MAXN = 1e5 + 5;
const int N = 1e5;
const int MAXK = 47 + 5;
const int K = 47;
const int MOD = 998244353;

typedef int arr[MAXN];

int add(int a, int b)
{
	return (a + b) % MOD;
}

int sub(int a, int b)
{
	return (a - b + MOD) % MOD;
}

int mul(int a, int b)
{
	return (ll)a * b % MOD;
}

int qpow(int a, int b)
{
	int base = a, ans = 1;
	while (b)
	{
		if (b & 1)
		{
			ans = (ll)ans * base % MOD;
		}
		base = (ll)base * base % MOD;
		b >>= 1;
	}
	return ans;
}

int inv(int a)
{
	return qpow(a, MOD - 2);
}

arr p, mu, phi, phi_pro, phi_pro_inv, phi_inv, f;
bool vis[MAXN];
vector<int> g[MAXN], h[MAXK][MAXK];

void pre()
{
	mu[1] = phi[1] = 1;
	for (int i = 2; i <= N; i++)
	{
		if (!vis[i])
		{
			p[++p[0]] = i;
			mu[i] = MOD - 1;
			phi[i] = i - 1;
		}
		for (int j = 1; j <= p[0] && i * p[j] <= N; j++)
		{
			vis[i * p[j]] = true;
			if (i % p[j] == 0)
			{
				mu[i * p[j]] = 0;
				phi[i * p[j]] = phi[i] * p[j];
				break;
			}
			mu[i * p[j]] = MOD - mu[i];
			phi[i * p[j]] = phi[i] * phi[p[j]];
		}
	}
	
	phi_pro[0] = phi_inv[1] = 1;
	for (int i = 1; i <= N; i++)
	{
		phi_pro[i] = mul(phi_pro[i - 1], phi[i]);
	}
	phi_pro_inv[N] = inv(phi_pro[N]);
	for (int i = N - 1; i >= 1; i--)
	{
		phi_pro_inv[i] = mul(phi_pro_inv[i + 1], phi[i + 1]);
		phi_inv[i + 1] = mul(phi_pro_inv[i + 1], phi_pro[i]);
	}
	
	for (int i = 1; i <= N; i++)
	{
		for (int j = 1; i * j <= N; j++)
		{
			f[i * j] = add(f[i * j], mul(mul(i, mu[j]), phi_inv[i]));
		}
	}
	
	for (int k = 1; k <= N; k++)
	{
		g[k].resize(N / k + 5);
		for (int n = 1; n * k <= N; n++)
		{
			g[k][n] = add(g[k][n - 1], phi[n * k]);
		}
	}
	
	for (int n = 1; n <= K; n++)
	{
		for (int m = n; m <= K; m++)
		{
			h[n][m].resize(N / m + 5);
			for (int l = 1; m * l <= N; l++)
			{
				h[n][m][l] = add(h[n][m][l - 1], mul(f[l], mul(g[l][n], g[l][m])));
			}
		}
	}
}

int block(int n, int m)
{
	int res = 0;
	for (int i = 1; i <= m / K; i++)
	{
		res = add(res, mul(f[i], mul(g[i][n / i], g[i][m / i])));
	}
	for (int l = m / K + 1, r; l <= n; l = r + 1)
	{
		int k1 = n / l, k2 = m / l;
		r = min(n / k1, m / k2);
		res = add(res, sub(h[k1][k2][r], h[k1][k2][l - 1]));
	}
	return res;
}

int main()
{
	pre();
	int t;
	scanf("%d", &t);
	while (t--)
	{
		int n, m;
		scanf("%d%d", &n, &m);
		if (n > m)
		{
			swap(n, m);
		}
		printf("%d\n", block(n, m));
	}
	return 0;
}