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 了,新年快乐吖!

上一篇:基本不等式


下一篇:关于『 一月の集训 』