题解 P5702 【【模板】调和级数求和】
Weng_Weijie
2020-02-10 01:16:14
## 前言
去我的博客里看可能体验会更好。
如果真的决定要做这个题,建议先去做一下 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;
}
```