题解 P5702 【【模板】调和级数求和】

Weng_Weijie

2020-02-10 01:16:14

Solution

## 前言 去我的博客里看可能体验会更好。 如果真的决定要做这个题,建议先去做一下 P5282 快速阶乘算法。 这个算法(包括阶乘算法、调和数算法、组合数前缀和算法)都由 $\text{min\_25}$ 最先提出,具体可以看他的[博客](https://min-25.hatenablog.com)。 ## 思路 总体的思路是:先分块,然后快速算出每一块内的倒数和。 首先构造函数: $$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}$ 慢)。 ## 代码 ```cpp #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; } ```