题解：P12395 「RiOI-6」神曲（加强版）

· · 题解

「知らない世界を見せてよ」

初始时值域为 [1,V]，从 1 到 n 依次确定每个区间 [l_i,r_i]，将该区间内的数缩成一个，这和原题意是等价的。设 f_{i,j} 表示确定了前 i 个区间、还有 j 个数的方案数，枚举第 i 个区间缩成的位置及原值域，有 f_{i,j}=j\displaystyle\sum_{k\ge j}f_{i-1,k}，对 j 做后缀和有 s_{i,j}=s_{i,j+1}+js_{i-1,j}，可以做到 \mathcal O(nV)。

考虑优化，把 DP 过程转置（？），有 t_{i,j}=t_{i,j-1}+jt_{i-1,j}，令 t'_{i+j,j}=t_{i,j}，有 t'_{i+j,j}=t'_{i+j-1,j-1}+jt'_{i+j-1,j}，这就是第二类斯特林数的递推式，而转置后答案为 t_{n,V}=t'_{n+V,V}，初值为 t_{0,0}=1，因此所求即为 \displaystyle{n+V\brace V}，现在要做的就是求第二类斯特林数在对角线上的值。

由第二类斯特林数·列的做法，我们知道 \displaystyle{{n\brace m}=\dfrac{1}{m!}\left[\dfrac{x^n}{n!}\right](\text{e}^x-1)^m}，现在目标就是对每个 m 求出 [x^{n+m}](\text{e}^x-1)^m，即 [x^n]\left(\dfrac{\text{e}^x-1}{x}\right)^m。里面这个东西的常数项为 1，只能单独用二项式定理拆出来，变为 \displaystyle\sum_{0\le k\le m}\binom mk[x^n]\left(\dfrac{\text{e}^x-x-1}{x}\right)^k。

这个形式就很另类拉反，设 F=\left(\dfrac{\text{e}^x-x-1}{x}\right)^{\langle-1\rangle}，所求 n 次项系数即为 [x^{n-k}]F'\left(\dfrac xF\right)^{n+1}，原式就化为了卷积形式，现在唯一要做的就是求出 F，如果你常数很小的话，直接 \mathcal O(n\log^2 n) 求复合逆可能就能直接通过原版，但反正我过不了，，

考虑牛顿迭代，有 \text{e}^F=(x+1)F+1，[x^1]F=2，现在由 \bmod\ x^{s+1} 的 F_0 得到 \bmod\ x^{2(s+1)} 的 F=F_0+x^{s+1}R，有：

\begin{aligned} \text{e}^{F_0+x^{s+1}R}&=(x+1)(F_0+x^{s+1}R)+1\\ \text{e}^{F_0}(1+x^{s+1}R)&\equiv(x+1)(F_0+x^{s+1}R)+1&\pmod{x^{2(s+1)}}\\ \text{e}^{F_0}-(x+1)F_0-1&\equiv x^{s+1}(x+1-\text{e}^{F_0})R&\pmod{x^{2(s+1)}} \end{aligned}

注意到右边的最低次项为 x^{s+2}，因此对于左半边 E 也有 [x^{s+1}]E=0；同时 F_0 最高次项为 x^s，因此对于 k>s+1，有 [x^k]E=[x^k]\text{e}^{F_0}。两侧同时除以 x^{s+2}，有：

\dfrac{\text{e}^{F_0}-(\text{e}^{F_0}\bmod x^{s+2})}{x^{s+2}}\equiv-\dfrac{\text{e}^{F_0}-x-1}{x}R\pmod{x^s}

求出 \text{e}^{F_0} 后跑长为 s 的求逆和乘法即可得到 R，这样写看起来常数会比 ri 题解里的写法快一点点。

然后就是卡常，求 exp 最好用 \mathcal O\left(\dfrac{n\log^2 n}{\log\log n}\right) 的多叉分治 NTT，比常规实现的 \mathcal O(n\log n) 要快将近一倍。

放一下核心部分代码。

il poly solve(int n) {
    poly f = {0, 2}; int l = 1; do {
        int s = l; l <<= 1; f.resize(l + 2); poly g = exp(f), u(s), v(s);
        for (int i = 0; i < s; i++) u[i] = g[s + 2 + i], v[i] = g[i + 1]; v[0] += -!v[0] & p, v[0]--;
        g = u * inv(v); for (int i = 0; i < s; i++) f[s + 1 + i] = neg(g[i]);
    } while (l < n); f.resize(n); return f;
}
int main() {
    int n = rd(), m = rd(); poly f = solve(n + 1), tf = f >> 1; // return 0;
    tf.resize(n + 1), f = deri(f) * pow(inv(tf), n + 1); poly g(m + 1, 0), h(m + 1, 0);
    for (int i = 0; i <= m; i++) g[i] = ifac[i], h[i] = i <= n ? (ll)ifac[i] * f[n - i] % p : 0; g = g * h;
    for (int i = 1; i <= m; i++) pt((ll)g[i] * fac[n + i] % p), pc(' ');
}