Note/Solution - 转置原理 & 多点求值

$$\newcommand{\vct}[1]{\boldsymbol{#1}} \newcommand{\mat}[1]{\begin{bmatrix}#1\end{bmatrix}} \newcommand{\opn}[1]{\operatorname{#1}} \mathscr{\text{Defining }\LaTeX\text{ Macros...}}$$

  我并没有透彻理解涉及知识点的严谨描述形式，所以本文大量用语是基于让读者理解而非让读者以此为研究依据的，烦请注意。

---

  设现有一个算法 $\mathcal A$，它接受向量 $\vct x$ 为输入，以向量 $\vct y$ 为输出，满足 $\vct y=A\vct x$，其中 $A$ 是常矩阵，则 $\mathcal A$ 为线性算法。转置原理指出，我们可以据此找到一个算法 $\mathcal B$ 用以计算 $\vct x=A^T\vct y$，且所用乘法次数不变，加法次数增加至多 $(m-n)$ 次。即，当已有足够优秀的算法 $\mathcal A$，我们“几乎”得到了同样优秀的算法 $\mathcal B$。（本段并不严谨，若有严重误导性请指正。）

  转置原理基于这样的事实：$(AB)^T=B^TA^T$，我们将 $A$ 分解为若干初等矩阵的乘积，令 $A=A_{k}A_{k-1}\cdots A_1$，那么 $A^T=A_1^TA_2^T\cdots A_k^T$，而初等矩阵的转置是显然的，因此能实现算法间的转化。具体地，从后往前考虑算法 $\mathcal A$ 中的指令：

• $\textbf{output }q_i~r_j$，即将变量 $r_j$ 的值作为输出的 $q_i$。将其改写为 $\textbf{input}~q_i~r_j$，即将变量 $r_j$ 的值作为输入的 $q_i$；
• $\textbf{swap }r_i~r_j$，即交换变量 $r_i,r_j$ 的值，保持不变；
• $\textbf{mul }r_i~k$，即将变量 $r_i$ 的值乘上常数 $k$，保持不变；
• $\textbf{add }r_i~r_j~k$，即将变量 $r_i$ 的值加上 $k$ 倍的 $r_j$。将其改写为 $\textbf{add }r_j~r_i~k$；
• $\textbf{input }p_i~r_j$，将其改写为 $\textbf{output }p_i~r_j$。

  值得一提的是，实际的算法往往难以表示为如此初等语句的顺序结构。若某个子算法的转置是显然的，我们无需再将其拆开，而可以直接取它的转置。典型的例子是 DFT，其插值的矩阵是对称的，所以 DFT 的转置就是 DFT 自身。

---

  举一个经典例子：多点求值。

  对于 $n$ 次多项式 $f(x)$，令

$$\vct f=\mat{ f_0\\ f_1\\ \vdots\\ f_{n-1} },$$

并给定 $\lang a_0,\cdots,a_{m-1}\rang$，令

$$A=\mat{ 1&a_0&\cdots&a_0^{n-1}\\ 1&a_1&\cdots&a_1^{n-1}\\ \vdots&\vdots&\ddots&\vdots\\ 1&a_{m-1}&\cdots&a_{m-1}^{n-1} },$$

求

$$\vct r=A\vct f=\mat{ f(a_0)\\ f(a_1)\\ \vdots\\ f(a_{m-1}) }.$$

  考虑取转置，转置问题形如求 $\vct q=A^T\vct p$。研究 $\vct q$ 的形式

$$q_i=\sum_{j=0}^{m-1}a_j^ip_j.$$

可见

$$q(x)\equiv \sum_{i=0}^{m-1}\frac{p_i}{1-a_ix}\pmod{x^n}.$$

不妨令 $q(x)=\sum_{i=0}^{m-1}\frac{p_i}{1-a_ix}$，整理其形式得

$$q(x)=\frac{[y]\prod_i(p_iy+1-a_ix)}{\prod_i(1-a_ix)}.$$

可以分治 $\mathcal O(n\log^2n)$ 求解。到此，我们对这一算法再次转置即可得到同复杂度，求解原问题的算法。

  设计算法：

$$\text{Algorithm 1: input }\vct p\text{, output }\vct q=A^T\vct p.\\ \begin{array}{} 1& \textbf{function }\opn{solve}(l,r)\\ 2& \quad \textbf{if }l=r\textbf{ then}\\ 3& \quad \quad \textbf{return }p_l\\ 4& \quad \textit{mid}\leftarrow\left\lfloor\frac{l+r}{2}\right\rfloor\\ 5& \quad L\leftarrow\opn{solve}(l,\textit{mid})\\ 6& \quad R\leftarrow\opn{solve}(\textit{mid}+1,r)\\ 7& \quad F_l\leftarrow L\times A_{[\textit{mid}+1,r]}\\ 8& \quad F_r\leftarrow R\times A_{[l,\textit{mid}]}\\ 9& \quad F\leftarrow F_l+F_r\\ 10& \quad \textbf{return }F\\ 11& \textbf{end function}\\ 12\\ 13& \textbf{input }\vct p\\ 14& \vct q\leftarrow \opn{solve}(0,m-1)\times A_{[0,m-1]}^{-1}\\ 15& \textbf{output }\vct q \end{array}$$

其中多项式 $A_{[l,r]}=\prod_{i=l}^r(1-a_ix)$，多项式系数序列与向量不作区分。进行转置得到

$$\text{Algorithm 2: input }\vct q\text{, output }\vct p=A\vct q.\\ \begin{array}{} 1& \textbf{function }\opn{solveT}(l,r,F)\\ 2& \quad \textbf{if }l=r\textbf{ then}\\ 3& \quad \quad p_l\leftarrow F_0\\ 4& \quad \quad \textbf{return null}\\ 5& \quad \textit{mid}\leftarrow\left\lfloor\frac{l+r}{2}\right\rfloor\\ 6& \quad L\leftarrow F\times^T A_{[\textit{mid+1,r}]}\\ 7& \quad R\leftarrow F\times^T A_{[l,\textit{mid}]}\\ 8& \quad \opn{solveT}(l,\textit{mid},L)\\ 9& \quad \opn{solveT}(\textit{mid}+1,r,R)\\ 10& \quad \textbf{return null}\\ 11& \textbf{end function}\\ 12\\ 13& \textbf{input }\vct q\\ 14& \opn{solveT}(0,m-1,\vct q\times^T A_{[0,m-1]}^{-1})\\ 15& \textbf{output }\vct p \end{array}$$

注意到递归的底层仅需要 $F_0$ 的值，所以可以保持 $\deg F=r-l+1$，那么复杂度亦为 $\mathcal O(n\log^2n)$。鉴于我初次理解的困难，这里给出详细的转化步骤。先看 $\text{A1}$ 中的 $\opn{solve}(l,r)$ 函数，我们把它转置成 $\text{A2}$ 中的 $\opn{solveT}(l,r,F)$，过程如下：

• $\begin{array}{}10&\textbf{return }F\end{array}$ 这个“返回值”可以看做函数的“输出”，它应当变为函数的“输入”，所以 $\opn{solveT}$ 额外增加了参数 $F$。
• $\begin{array}{}9&F\leftarrow F_l+F_r\end{array}$ 实际上跳步了。应当是初始 $F\leftarrow0$，后 $F\leftarrow F+1\times F_l$，再 $F\leftarrow F+1\times F_r$，取转置后，得到 $F_l\leftarrow F,F_r\leftarrow F$，所以 $\text{A2}$ 中直接用了 $F$ 而并没有添加变量 $F_l$ 和 $F_r$。
• $\begin{array}{}8&F_r\leftarrow R\times A_{[l,\textit{mid}]}\end{array}$ 注意 $A_{[l,\textit{mid}]}$ 是常量，取转置得到 $R\leftarrow F\times^T A_{[l,\textit{mid}]}$，其中 $\times^T$ 是 $\times$ 的转置。第 $7$ 行同理。
• $\begin{array}{}6&R\leftarrow\opn{solve}(\textit{mid+1},r)\end{array}$ “输出”变“输入”，转置得 $\opn{solveT}(\textit{mid}+1,r,R)$。第 $5$ 行同理。
• $\begin{array}{}3&\textbf{return }p_l\end{array}$ 形象地说，想想你自己写代码的时候，可能在这里才读入 $p_l$ 的值。所以这个实际上是输入 $p_l$（作为 $F_0$），转置为输出 $p_l$（值为 $F_0$）。

  主过程就三行，不讲啦。

  对于 $\times^T$，写出暴力多项式卷积的算法 $A(x)\times B(x)$，将 $B$ 视为常量得到一个关于 $A(x)$ 的线性算法。取转置，发现就是 $A(x)$ 与 $B(x)$ 在做差卷积。所以卷积的转置是差卷积。

---

  先到这里叭，请一定去看看 Tiw 的讲解 w！

  下面是多点求值的代码，挺快的√

/*+Rainybunny+*/

#include <bits/stdc++.h>

#define rep(i, l, r) for (int i = l, rep##i = r; i <= rep##i; ++i)
#define per(i, r, l) for (int i = r, per##i = l; i >= per##i; --i)

typedef std::vector<int> Poly;

const int MAXN = 1 << 17, MOD = 998244353;
int n, m, a[MAXN + 5], ans[MAXN + 5];
Poly A[MAXN << 2], Q;

inline int mul(const int u, const int v) { return 1ll * u * v % MOD; }
inline int sub(int u, const int v) { return (u -= v) < 0 ? u + MOD : u; }
inline int add(int u, const int v) { return (u += v) < MOD ? u : u - MOD; }
inline int mpow(int u, int v) {
    int ret = 1;
    for (; v; u = mul(u, u), v >>= 1) ret = mul(ret, v & 1 ? u : 1);
    return ret;
}

namespace PolyOper {

const int G = 3;
int omega[18][MAXN + 5];

inline void init() {
    rep (i, 1, 17) {
        int* wi = omega[i];
        wi[0] = 1, wi[1] = mpow(G, MOD - 1 >> i);
        rep (j, 2, (1 << i) - 1) wi[j] = mul(wi[j - 1], wi[1]);
    }
}

inline void ntt(Poly& u, const int tp) {
    static int rev[MAXN + 5]; int n = u.size();
    rep (i, 0, n - 1) rev[i] = rev[i >> 1] >> 1 | (i & 1) * n >> 1;
    rep (i, 0, n - 1) if (i < rev[i]) std::swap(u[i], u[rev[i]]);
    for (int i = 1, stp = 1; stp < n; ++i, stp <<= 1) {
        int* wi = omega[i];
        for (int j = 0; j < n; j += stp << 1) {
            rep (k, j, j + stp - 1) {
                int ev = u[k], ov = mul(wi[k - j], u[k + stp]);
                u[k] = add(ev, ov), u[k + stp] = sub(ev, ov);
            }
        }
    }
    if (!~tp) {
        int inv = mpow(n, MOD - 2);
        std::reverse(u.begin() + 1, u.end());
        for (int& a: u) a = mul(a, inv);
    }
}

inline Poly pmul(Poly u, Poly v) {
    int res = u.size() + v.size() - 1, len = 1;
    while (len < res) len <<= 1;
    u.resize(len), v.resize(len);
    ntt(u, 1), ntt(v, 1);
    rep (i, 0, len - 1) u[i] = mul(u[i], v[i]);
    ntt(u, -1);
    return u.resize(res), u;
}

inline Poly pmulT(Poly u, Poly v) {
    int n = u.size(), m = v.size();
    std::reverse(v.begin(), v.end()), v = pmul(u, v);
    rep (i, 0, n - 1) u[i] = v[i + m - 1];
    return u;
}

inline void pinv(const int n, const Poly& u, Poly& r) {
    if (n == 1) return void(r = { { mpow(u[0], MOD - 2) } });
    static Poly tmp; pinv(n >> 1, u, r);
    tmp.resize(n << 1), r.resize(n << 1);
    rep (i, 0, n - 1) tmp[i] = i < u.size() ? u[i] : 0;
    rep (i, n, (n << 1) - 1) tmp[i] = 0;
    ntt(r, 1), ntt(tmp, 1);
    rep (i, 0, (n << 1) - 1) r[i] = mul(r[i], sub(2, mul(tmp[i], r[i])));
    ntt(r, -1), r.resize(n);
}

} // namespace PolyOper.

inline void init(const int u, const int l, const int r) {
    if (l == r) return void(A[u] = { { 1, sub(0, a[l]) } });
    int mid = l + r >> 1;
    init(u << 1, l, mid), init(u << 1 | 1, mid + 1, r);
    A[u] = PolyOper::pmul(A[u << 1], A[u << 1 | 1]);
}

inline void solveT(const int u, const int l, const int r, Poly F) {
    F.resize(r - l + 1);
    if (l == r) return void(ans[l] = F[0]);
    int mid = l + r >> 1;
    solveT(u << 1, l, mid, PolyOper::pmulT(F, A[u << 1 | 1]));
    solveT(u << 1 | 1, mid + 1, r, PolyOper::pmulT(F, A[u << 1]));
}

int main() {
    scanf("%d %d", &n, &m), Q.resize(++n);
    for (int& u: Q) scanf("%d", &u);
    rep (i, 0, m - 1) scanf("%d", &a[i]);

    PolyOper::init();
    init(1, 0, m - 1);
    int len = 1; while (len < n) len <<= 1;
    Poly T; PolyOper::pinv(len, A[1], T);
    
    solveT(1, 0, m - 1, PolyOper::pmulT(Q, T));
    rep (i, 0, m - 1) printf("%d\n", ans[i]);
    return 0;
}

  嘛，2022 了，新年快乐吖！