自动搬运

查看原文

来自洛谷，原作者为

	Weng_Weijie やー！今日もパンがうまい！

搬运于2025-08-24 22:13:33，当前版本为作者最后更新于2020-02-10 01:16:14，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

---

前言

去我的博客里看可能体验会更好。

如果真的决定要做这个题，建议先去做一下 P5282 快速阶乘算法。

这个算法（包括阶乘算法、调和数算法、组合数前缀和算法）都由 $\text{min\_25}$ 最先提出，具体可以看他的博客。

思路

总体的思路是：先分块，然后快速算出每一块内的倒数和。

首先构造函数：

$$h(x)=\sum_{i=1}^m\dfrac 1{x+i}$$

其中 $m=\mathrm O(\sqrt{n})$。

我们希望求出 $f(0), f(m), \ldots, f(m^2)$，这样就算出了大块内的总和，剩下的一部分直接暴力即可。

然而这个 $h(x)$ 并不是多项式，我们根本没办法表示它。

于是我们将 $h(x)$ 通分。

$$h(x)=\dfrac{\sum_i\prod_{j\neq i}(x+j)}{\prod_i(x+i)} $$

于是令

$$g(x)=\sum_{i=1}^m\prod_{j\neq i}(x+j)$$ $$f(x)=\prod_{i=1}^m(x+i)$$

那么 $h(x)=\dfrac{g(x)}{f(x)}$。

此时便可以直接使用多点求值的做法达到 $\mathrm O(\sqrt n\log^2n)$，但是应该无法通过。

优化

我们设 $f_t(x)=\displaystyle\prod_{i=1}^t(x+i)$， $g_t(x)$ 同理。

发现我们并不需要 $f, g$ 的每一项系数，只需要它在一些位置上的点值。

再设 $F_t=(f_t(0), f_t(m),\dots f_t(tm))$，$G_t$ 同理。

因为 $f_t(x)$ 是一个 $t$ 次多项式，因此 $F_t$ 是 $f_t(x)$ 的一个点值表示。

观察到：

$$f_{2t}(x)=f_t(x)f_t(x+t)$$ $$g_{2t}(x)=f_t(x)g_t(x+t)+f_t(x+t)g_t(x)$$

我们希望在已知 $F_t, G_t$ 的情况下，求出 $F_{2t}, G_{2t}$，这样就可以倍增求出 $F_m, G_m$，从而达成目标。

接下来以 $f$ 为例，$g$ 也是类似的。

在求 $f_{2t}(x)$ 时需要知道 $f_{t}(x), f_{t}(x+t)$，因此求 $F_{2t}$ 需要对每一个 $0\leq x\leq 2t$，求出 $f_t(mx)$ 和 $f_t{(mx+t)}$。

如果令 $p_t(x)=f_t(mx)$，那么 $F_t=(p(0),p(1),\ldots,p(t))$，$f_t(mx)=p_t(x), f_t(mx+t)=p_t\left(x+\dfrac tm\right)$。

如果我们在已知 $p(0),p(1),\ldots,p(t)$ 的情况下求出 $p(k),p(1+k),\ldots,p(t+k)$，那么就可以令 $k=t+1$ 先求出 $p(t+1),\ldots,p(2t)$，再令 $k=\dfrac tm$，求出需要的所有值。

事实上，由 $\text{Lagrange}$ 插值法：

$$p(x+k)=\sum_{i=0}^tp(i)\prod_{j\neq i}\dfrac{x+k-j}{i-j} $$$$=\sum_{i=0}^t\dfrac {p(i)(-1)^{t-i}}{i!(t-i)!}\cdot \prod_{j\neq i} (x+k-j) $$$$=\dfrac{(x+k)!}{(x+k-t-1)!}\sum_{i=0}^t\dfrac{p(i)(-1)^{t-i}}{i!(t-i)!}\cdot\dfrac{1}{x-i+k} $$

可以看到右边是一个卷积的形式，可以用 $\text{FFT}$ 实现。

可以证明不会出现没有逆元的情况（？不知道我有没有伪证）。

因此我们可以由 $F_t, G_t$ 得到 $F_{2t}, G_{2t}$，而得到 $F_{t+1}, G_{t+1}$ 是比较简单的（留给读者思考），这样倍增算出 $F_m, G_m$，就能得到大块内部的倒数和。

因此我们解决了整个问题，总复杂度 $O(\sqrt n\log n)$。

其他

在算阶乘时我们使用了威尔逊定理使 $n$ 的范围，缩小了一半，而在这题里可以发现 $H_n=H_{\mathrm{mod}-n-1}$，也达到了一样的效果。

在算逆元的时候强烈推荐使用离线求逆元（快速幂求逆元 $n$ 次有时会比两次左右的 $\mathrm{FFT}$ 慢）。

代码

#include <bits/stdc++.h>

typedef long long LL;
typedef unsigned long long ULL;

const int N = 131072;

int mod, mod_g, factor[N], ifactor[N];

void reduce(int &x) { x += x >> 31 & mod; }
int pow(int x, int y, int ans = 1) {
	for (; y; y >>= 1, x = (LL) x * x % mod)
		if (y & 1) ans = (LL) ans * x % mod;
	return ans;
}

void init(int n) {
	factor[0] = 1;
	for (int i = 1; i <= n; ++i)
		factor[i] = (LL) factor[i - 1] * i % mod;
	ifactor[n] = pow(factor[n], mod - 2);
	for (int i = n; i; --i)
		ifactor[i - 1] = (LL) ifactor[i] * i % mod;
}

int wn[N], w[N], lim, s, rev[N];
void fftinit(int len) {
	wn[0] = lim = 1, s = -1; while (lim < len) lim <<= 1, ++s;
	for (int i = 1; i < lim; ++i) rev[i] = rev[i >> 1] >> 1 | (i & 1) << s;
	const int w = pow(mod_g, (mod - 1) / lim);
	for (int i = 1; i < lim; ++i) wn[i] = (LL) wn[i - 1] * w % mod;
}
void fft(int *A, int typ) {
	static ULL tmp[N];
	for (int i = 0; i < lim; ++i) tmp[rev[i]] = A[i];
	for (int i = 1; i < lim; i <<= 1) {
		for (int j = 0, t = lim / i / 2; j < i; ++j) w[j] = wn[j * t];
		for (int j = 0; j < lim; j += i << 1)
			for (int k = 0; k < i; ++k) {
				const ULL x = tmp[k + j + i] * w[k] % mod;
				tmp[k + j + i] = tmp[k + j] + mod - x, tmp[k + j] += x;
			}
	}
	for (int i = 0; i < lim; ++i) A[i] = tmp[i] % mod;
	if (!typ) {
		const int il = pow(lim, mod - 2); std::reverse(A + 1, A + lim);
		for (int i = 0; i < lim; ++i) A[i] = (LL) A[i] * il % mod;
	}
}

void clear(int *a) { std::memset(a, 0, lim << 2); }
void copy(int *a, int *b, int n) { std::memcpy(a, b, n + 1 << 2); }
void multiply(int *a, int *b, int *c) {
	fft(a, 1), fft(b, 1);
	for (int i = 0; i < lim; ++i) c[i] = (LL) a[i] * b[i] % mod;
	fft(c, 0);
}
void poly_shift(int *f, int n, int *g, int k) {
	static int a[N], b[N], q[N]; fftinit(n + n + 1), clear(a), clear(b);
	for (int i = 0; i <= n; ++i) a[i] = (LL) f[i] * ifactor[i] % mod * ifactor[n - i] % mod;
	for (int i = n - 1; i >= 0; i -= 2) a[i] = mod - a[i];
	
	b[0] = k - n;
	for (int i = 1; i <= 2 * n; ++i) b[i] = (LL) b[i - 1] * (k + i - n) % mod;
	q[2 * n] = pow(b[2 * n], mod - 2);
	for (int i = 2 * n; i; --i) q[i - 1] = (LL) q[i] * (k + i - n) % mod;
	for (int i = 2 * n; i; --i) b[i] = (LL) q[i] * b[i - 1] % mod; b[0] = q[0];
	
	multiply(a, b, a);
	static int suf[N], isuf[N]; suf[2 * n + 1] = 1;
	for (int i = 2 * n; ~i; --i)
		suf[i] = (LL) suf[i + 1] * (i + k - n) % mod;
	isuf[0] = pow(suf[0], mod - 2);
	for (int i = 0; i <= 2 * n; ++i)
		isuf[i + 1] = (LL) isuf[i] * (i + k - n) % mod;
	for (int i = 0; i <= n; ++i) {
		if ((i + k) % mod <= n) g[i] = f[(i + k) % mod];
		else g[i] = (LL) isuf[i + n + 1] * a[i + n] % mod * suf[i] % mod;
	}
}

int n, size, f[N], g[N];

void boom(int n) {
	static int a[N], b[N], c[N], d[N];
	poly_shift(f, n, a, n + 1);
	for (int i = 1; i <= n; ++i) f[n + i] = a[i - 1];
	poly_shift(f, 2 * n, b, pow(size, mod - 2, n));
	poly_shift(g, n, c, n + 1);
	for (int i = 1; i <= n; ++i) g[n + i] = c[i - 1];
	poly_shift(g, 2 * n, d, pow(size, mod - 2, n));
	for (int i = 0; i <= 2 * n; ++i) {
		g[i] = ((LL) g[i] * b[i] + (LL) d[i] * f[i]) % mod;
		f[i] = (LL) f[i] * b[i] % mod;
	}
}
void qaq(int n) {
	f[n + 1] = g[n + 1] = 1;
	int tmp = 1, x = (n + 1) * size % mod;
	for (int i = 1; i <= n; ++i)
		f[n + 1] = (LL) f[n + 1] * (x + i) % mod;
	for (int i = 2; i <= n; ++i) {
		tmp = tmp * (x + i - 1LL) % mod;
		g[n + 1] = ((LL) g[n + 1] * (x + i) + tmp) % mod;
	}
	for (int i = 0; i <= n + 1; ++i) {
		x = (LL) i * size % mod;
		g[i] = (g[i] * (x + n + 1LL) + f[i]) % mod;
		f[i] = f[i] * (x + n + 1LL) % mod;
	}
}

void solve(int n) {
	if (n == 1) {f[1] = size + 1, f[0] = g[0] = g[1] = 1; return;}
	solve(n >> 1), boom(n >> 1); if (n & 1) qaq(n - 1);
}

int solve() {
	if (!n) return 0; int ans = 0;
	size = std::sqrt(n), init(size), solve(size);
	int x = 0, y = 1;
	for (int i = 0; i < size; ++i)
		x = ((LL) f[i] * x + (LL) y * g[i]) % mod, y = (LL) y * f[i] % mod;
	for (int i = size * size + 1; i <= n; ++i)
		x = ((LL) i * x + y) % mod, y = (LL) i * y % mod;
	return pow(y, mod - 2, x);
}

void test() {
	std::cin >> n >> mod >> mod_g, n = std::min(n, mod - n - 1);
	std::cout << solve() << '\n';
}

int main() {
	std::ios::sync_with_stdio(0), std::cin.tie(0);
	int tc; std::cin >> tc; while (tc--) test();
	return 0;
}