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

前言

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

如果真的决定要做这个题，建议先去做一下 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;
}