其他分享
首页 > 其他分享> > Note/Solution - 浅尝转置原理 & 多点求值

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\) 中的指令：

  值得一提的是，实际的算法往往难以表示为如此初等语句的顺序结构。若某个子算法的转置是显然的，我们无需再将其拆开，而可以直接取它的转置。典型的例子是 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)\)，过程如下：

  主过程就三行，不讲啦。

  对于 \(\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 了，新年快乐吖！

标签：leftarrow,转置,mid,textbf,浅尝,vct,求值,quad
来源： https://www.cnblogs.com/rainybunny/p/15755358.html