相信很多人都想到了 $O(n)$ 的解法：线性处理逆元，然后计算。此做法可以通过 $\color{#000000}30$ 分。

观察题目数据发现：本题需要一个低于 $O(n)$ 的解法。

我们先设 $s=\lfloor\sqrt{n}\rfloor$ 且 $h(x)=\sum_{i=1}^s\dfrac{1}{x+i}$，那么

$$\sum_{i=1}^n\dfrac{1}{n}\equiv(\sum_{i=0}^{s-1}h(is))+\sum_{i=s^2+1}^{n}\dfrac{1}{i}\pmod p $$

其中 $\sum_{i=s^2+1}^{n}\dfrac{1}{i}$ 可在 $O(s)$ 内完成。我们希望快速计算 $\sum_{i=0}^{s-1}h(is)$。

然而这个 $h(x)$ 根本就不是整式。我们试着将其通分：

$$\begin{aligned} h(x)&=\sum_{i=1}^s\dfrac{1}{x+i}\\ &=\dfrac{\sum_{i=1}^s\prod_{j\neq i}x+j}{\prod_{i=1}^{s}x+i} \end{aligned} $$

然后我们设分母为 $f(x)$，分子为 $g(x)$。于是我们成功地将非整式转换成了两个整式。

$O(\sqrt n\log^2 n)$ 的做法

我们可以使用多项式连续点值平移在 $O(s\log s)$ 的时间复杂度内得到 $f(x)$ 和 $g(x)$ 的各项系数，但是通过多项式多点求值得到 $f(0),f(s),f(2s),\cdots,f(s^2-s)$ 有 $O(s\log^2 s)$ 的时间复杂度，$g(x)$ 同理。

于是总时间复杂度为 $O(s\log^2 s)=O(\sqrt n\log^2 n)$。

$O(\sqrt n\log n)$ 的做法

令 $f_t(x)=\sum_{i=1}^{t}x+i,g_t(x)=\sum_{i=1}^t\prod_{j\neq i}x+j$。通过倍增即可求出 $f(0),f(s),f(2s),\cdots,f(s^2-s)$ 以及 $g(0),g(s),g(2s),\cdots,g(s^2-s)$。

试试倍增？

• 倍增方程1（乘以 $2$)

$$\begin{aligned} f_{2t}(x)&=f_t(x)\times f_t(t+x)\\ g_{2t}(x)&=\sum_{i=1}^{2t}\prod_{j\neq i}x+j\\ &=\sum_{i=1}^{t}\prod_{j\neq i}(x+j)\times\prod_{i=j+1}^{2j}(x+i)+\sum_{i=t+1}^{2t}\prod_{j\neq i}(x+j)\times\prod_{i=1}^{j}(x+j)\\ &=g_t(x)f_t(t+x)+\sum_{i=1}^{t}\prod_{j\neq i}(x+t+i)\times f_t(x)\\ &=g_t(x)f_t(t+x)+g_t(t+x)f_t(x) \end{aligned} $$

现在看怎么计算。设 $F(x)=f_t(sx),G(x)=g_t(sx)$，则：

1. 现有的点值是 $F(0), F(1), \cdots, F(j)$（即 $f_j(0), f_j(s),\cdots,f_j(js)$）和 $G(0), G(1), \cdots, G(j)$（即 $g_j(0), g_j(s),\cdots,g_j(js)$）
2. 通过多项式连续点值平移得到 $F(j+1), F(j+2), \cdots, F(2j+1)$ 以及 $G(j+1), G(j+2), \cdots, G(2j+1)$（这里由于要求区间不重叠，会多平移出一个 $F(2j+1)$ 和 $G(2j+1)$，否则会因为除以零寄掉）。
3. 将其拼接后得到 $F(0), F(1),\cdots,F(2j)$（即 $f_j(0), f_j(0+s),\cdots,f_j(2js)$）和 $G(0), G(1),\cdots,G(2j)$（即 $g_j(0), g_j(0+s),\cdots,g_j(2js)$）
4. 再次经过平移后得到 $F(0+\dfrac{j}{s}), F(1\dfrac{j}{s}),\cdots,F(2j\dfrac{j}{s})$（即 $f_j(j), f_j(j+s),\cdots,f_j(j+2js)$）和 $G(0+\dfrac{j}{s}), G(1\dfrac{j}{s}),\cdots,G(2j\dfrac{j}{s})$（即 $g_j(j), g_j(j+s),\cdots,g_j(j+2js)$）。
5. 乘起来即得 $f_{2j}(0),f_{2j}(s),\cdots,f_{2j}(2js)$ 和 $g_{2j}(0),g_{2j}(s),\cdots,g_{2j}(2js)$。

• 倍增方程2（加上 $1$）

$$\begin{aligned} f_{t+1}(x)&=f_t(x)\times (x+t+1)\\ g_{t+1}(x)&=\sum_{i=1}^{t+1}\prod_{j\neq i}x+j\\ &=\sum_{i=1}^{t}\prod_{j\neq i}(x+j)\times\prod_{i=t+1}^{t+1}(x+i)+\sum_{i=t+1}^{t+1}\prod_{j\neq i}(x+j)\times\prod_{i=1}^{t}(x+i)\\ &=g_t(x)\times(x+t+1)+f_t(x) \end{aligned} $$

然后就可以倍增了。

我们只需要计算 $s$ 个点值，所以时间复杂度为 $O(s\log s)=O(\sqrt n\log n)$。

优化

1. 由于我们在平移时已经计算出了 $g_t(2t+1)$ 的值，所以如果要加上 $1$，可以直接使用这个值。
2. 根据威尔逊定理，我们得到 $h(x)=h(p-x-1)$。这意味着我们可以将常数减半。

上代码：

inline int harmonic(int n, int p) {
	(n > p - n - 1) && (n = p - n - 1);
	vector<int> mc[2], md[2];
	int pos = 0, s = sqrt(n) + 1e-6;
	mc[0].reserve(s);
	mc[1].reserve(s);
	md[0].reserve(s);
	md[1].reserve(s);
	int invs = fpow(s, p - 2, p);
	static vector<int> st;
	st.resize((int)log2(s) + 5);
	for (int i = s; i > 1; i >>= 1) {
		st[++pos] = i;
	}
	mc[0].resize(2);
	mc[1].resize(2);
	mc[0][0] = mc[1][0] = mc[1][1] = 1;
	mc[0][1] = s + 1;
	for (int l = st[pos]; pos; l = st[--pos]) {
        LagrangeInterpolation_ex(l >> 1, (l >> 1) + 1, p, mc[0], md[0]);
        mc[0].resize(mc[0].size() << 1);
        std::copy(md[0].begin(), md[0].end(), mc[0].begin() + md[0].size());
        LagrangeInterpolation_ex(mc[0].size() - 1, bt(1ll * invs * (l >> 1), p), p, mc[0], md[0]);
        LagrangeInterpolation_ex(l >> 1, (l >> 1) + 1, p, mc[1], md[1]);
        mc[1].resize(mc[1].size() << 1);
        std::copy(md[1].begin(), md[1].end(), mc[1].begin() + md[1].size());
        LagrangeInterpolation_ex(mc[1].size() - 1, bt(1ll * invs * (l >> 1), p), p, mc[1], md[1]);
		for (int i = (int)mc[0].size() - 1; ~i; i--) {
			mc[1][i] = bt(bt(1ll * mc[1][i] * md[0][i], p) + bt(1ll * mc[0][i] * md[1][i], p), p);
			mc[0][i] = bt(1ll * mc[0][i] * md[0][i], p);
		}
		if (l & 1) {
			for (int i = 0, x; i <= l; i++) {
				x = bt(bt(1ll * i * s, p) + l, p);
				mc[1][i] = bt(bt(1ll * mc[1][i] * x, p) + mc[0][i], p);
				mc[0][i] = bt(1ll * mc[0][i] * x, p);
			}
		}
		else {
			mc[0].resize(l + 1);
			mc[1].resize(l + 1);
		}
	}
	int res1 = 1, res2 = 0;
	for (int i = 0; i < s; i++) {
		res2 = bt(bt(1ll * res1 * mc[1][i], p) + bt(1ll * res2 * mc[0][i], p), p);
		res1 = bt(1ll * res1 * mc[0][i], p);
	}
	for (int i = s * s + 1; i <= n; i++) {
		int x = i;
		res2 = bt(res1 + bt(1ll * x * res2, p), p);
		res1 = bt(1ll * x * res1, p);
	}
	return bt(1ll * fpow(res1, p - 2, p) * res2, p);
}
int main() {
	int T, n, p;
	std::cin.tie(NULL)->sync_with_stdio(false);
	std::cin >> T;
	while (T--) {
		std::cin >> n >> p >> NTT::G;
		brt = ((unsigned __int128)1 << 64) / p; // 这题不用 barrett 取模约减 NTT 过不了
		NTT::invG = fpow(NTT::G, p - 2, p);
		NTT::mod = p;
		std::cout << harmonic(n, p) << std::endl;
	}
}