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

2022-01-01 14:36:32 阅读：184 来源： 互联网

标签：leftarrow 转置 mid textbf 浅尝 vct 求值 quad

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

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

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9

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