闲话

总是想封装一个多项式牛顿迭代（
像 fjzzq 的板子里那种的（
感觉下午就能封装完了（

feux follets 是鬼火的意思嘛
我总是想到西奈鬼火（

今日推歌：一人行者 - ilem feat. 心华
挺有感触的……

再再谈排列计数

好我会了（
可能的前置知识：《转置原理的简单介绍》阅读随笔

写完发现不是那么的难？

Feux Follets

承接上文，可以发现我们要计算的就是一个序列 \(\langle g_m\rangle\)，满足

\[g_m = \sum_{i = 0}^k F(i) [x^m y^i]e^{y(-\ln(1 - x)-x)}
\]

最后对每个点值乘 \(m!\) 作为提取 \(x^m\) 时 egf 的贡献即可。

在 \(y\) 上的求和不是很好搞，但是观察后面的 bgf 对 \(x\) 的偏导可以得到一个比较良好的形式，因此考虑转置问题。原问题的变换矩阵 \(A[m, i] = [x^m y^i]e^{y(-\ln(1 - x)-x)}\)，作转置得到 \(A^{\mathsf T} [m, i] = [x^m y^i]e^{y(-\ln(1 - x)-x)}\)；输入向量 \(\bm f\) 就是 \(f_i = F(i)\)。假设输出向量为 \(\bm h\)，转置问题也就是计算序列 \(\langle h_m\rangle\) 满足

\[h_m = \sum_{i\ge 0} F(i) [x^i y^m] e^{y(-\ln(1 - x)-x)}
\]

发现序列 \(h\) 的每一项其实就是 bgf \(G(x, y) = e^{y(-\ln(1 - x)-x)}\) 的各 \(x^i\) 的系数多项式乘上输入向量的加和的各项，则序列 \(h\) 的生成函数 \(H(y)\) 即为

\[H(y) = \sum_{i\ge 0} f_i [x^i] G(x, y)
\]

考虑求解 \(G_i(y) = [x^i]G(x, y)\)。不难发现 \(G\) 对 \(x\) 的偏导有良好形式，也就是

\[\frac{\partial}{\partial x} G(x, y) = \frac{\partial}{\partial x} e^{y(-\ln (1 - x) - x)} = \frac{\partial}{\partial x} \frac{1}{e^{xy} (1 - x)^y} = \frac{y e^{xy}(1 - x)^y - y (1 - x)^{y - 1} e^{xy}}{(e^{xy} (1 - x)^y)^2} = \frac{y - y / (1 - x)}{e^{xy} (1 - x)^y} = \frac{xy}{1 - x}G(x,y)
\]

也就是

\[(1 - x)\frac{\partial}{\partial x} G(x, y) = xyG(x, y)
\]

两边提取 \(x^i\) 项系数可以得到

\[(i + 1)G_{i + 1} - i G_i = y G_{i - 1}
\]

化简得到

\[G_i = \frac{i - 1}{i} G_{i - 1} + \frac{y}{i} G_{i - 2}
\]

容易得到 \(G_0 = 1\)，而假设 \(G_{-1} = 0\) 可以使我们能够用矩阵描述转移，得到

\[\begin{bmatrix} G_{i - 1} \\ G_{i} \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ y / i & (i - 1) / i \end{bmatrix} \begin{bmatrix} G_{i - 2} \\ G_{i - 1} \end{bmatrix}
\]

设 \(\text{M}_i = \begin{bmatrix} 0 & 1 \\ y / i & (i - 1) / i \end{bmatrix}\)，我们可以得到

\[\begin{bmatrix} G_{i - 1} \\ G_{i} \end{bmatrix} = \text{M}_i\times \text{M}_{i - 1}\times \cdots\times \text{M}_1\times \begin{bmatrix} 0 \\ 1 \end{bmatrix}
\]

也就是说

\[\text{M}_i\times \text{M}_{i - 1}\times \cdots\times \text{M}_1 = \begin{bmatrix} \blacksquare & G_{i - 1} \\ \blacksquare & G_{i} \end{bmatrix}
\]

\(\blacksquare\) 代表我们不关心这个位置的值。我们只需让得到的结果矩阵的右下角乘上 \(f_i\) 就是第 \(i\) 项对 \(H(y)\) 的贡献了。

求解这个仍然可以考虑分治。套路地，我们设

\[P_{[l, r)} = \text{M}_{r - 1} \times \text{M}_{r - 2} \times \cdots \times \text{M}_{l} \qquad Q_{[l, r)} = \sum_{i = l}^{r - 1} f_i\times \text{M}_{i} \times \text{M}_{i - 1} \times \cdots \times \text{M}_{l}
\]

然后转移有

\[P_{[l, r)} = P_{[\text{mid}, r)}\times P_{[l, \text{mid})}
\]

\[Q_{[l, r)} = Q_{[l, \text{mid})} + Q_{[\text{mid}, r)}\times P_{[l, \text{mid})}
\]

边界条件是 \(P_{[l, l + 1)} = \text{M}_l, Q_{[l, l + 1)} = f_l\times \text{M}_l\)。

将该算法转置即可。

总时间复杂度为 \(O(n\log^2 n + k\log^2 k)\)。

Submission.
常数好大！hos_lyric 是咋写的啊（

code

struct mat {
    poly a[2][2];
    poly* operator[] (const int& p) { return a[p]; }
    const poly* operator[] (const int& p) const { return a[p]; }
    mat(int _len = 0) { rep(i,0,1) rep(j,0,1) a[i][j].resize(_len); }
    inline friend mat operator* (const mat& a, const mat& b) {
        mat ret;
        ret[0][0] = a[0][0] * b[0][0] + a[0][1] * b[1][0];
        ret[0][1] = a[0][0] * b[0][1] + a[0][1] * b[1][1];
        ret[1][0] = a[1][0] * b[0][0] + a[1][1] * b[1][0];
        ret[1][1] = a[1][0] * b[0][1] + a[1][1] * b[1][1];
        return ret;
    } 
    inline friend mat operator^ (const mat& a, const mat& b) {
        using namespace polynomial :: __multipoint_operation__;
        mat ret;
        ret[0][0] = _E_Mul(a[0][0], b[0][0]) + _E_Mul(a[0][1], b[0][1]);
        ret[0][1] = _E_Mul(a[0][0], b[1][0]) + _E_Mul(a[0][1], b[1][1]);
        ret[1][0] = _E_Mul(a[1][0], b[0][0]) + _E_Mul(a[1][1], b[0][1]);
        ret[1][1] = _E_Mul(a[1][0], b[1][0]) + _E_Mul(a[1][1], b[1][1]);
        return ret;
    }
} Q[N], S(1);

#define ls (p << 1)
#define rs (p << 1 | 1)

void Init(int p, int l, int r) {
    if (l == r) {
        int invk = qp(l, mod - 2);
        Q[p][0][0] = { 0 };
        Q[p][0][1] = { 1 };
        Q[p][1][0] = { 0, invk };
        Q[p][1][1] = { 1ll * (l - 1) * invk % mod };
        return ;
    } int mid = l + r >> 1;
    Init(ls, l, mid);
    Init(rs, mid + 1, r);
    if (r != n) Q[p] = Q[rs] * Q[ls];
}

void Calc(int p, int l, int r, const mat& P) {
    mat _P(r - l + 2);
    rep(i,0,1) rep(j,0,1) rep(k,0,r-l+1) _P[i][j][k] = (k >= P[i][j].size() ? 0 : P[i][j][k]);
    if (l == r) {
        _P = _P ^ Q[p]; 
        assert(_P[0][0].size() > 0); 
        assert(_P[1][1].size() > 0);
        ans[l] = (_P[0][0][0] + _P[1][1][0]); 
        if (ans[l] >= mod) ans[l] -= mod;
        return;
    } int mid = l + r >> 1;
    Calc(ls, l, mid, _P);
    Calc(rs, mid + 1, r, _P ^ Q[ls]);
}

signed main() {
    cin >> n >> k;
    poly f(k), pt(k); ans.resize(n + 1);
    rep(i,0,k-1) cin >> f[i], pt[i] = i; 
    S[1][1] = f.eval(k, pt);
    Init(1, 1, n);
    Calc(1, 1, n, S);
    rep(i,1,n) cout << 1ll * ans[i] * gfac(i) % mod << ' '; 
}