前置芝士:FFT 、高斯消元法
众所周知,多项式多点求值可以用转置原理优化。
EI 巨佬 与 rqy 学姐 都写过这方面的 blog,然而读了还是不会 QAQ。
直到读了原论文,才勉强晓得了咋回事。
于是就有了这篇文章。
如有错漏敬请读者在评论区一一指出或喷作者。
初等矩阵
(有基础的同学可以跳过)
要学会转置原理,首先要知道初等矩阵。
通过高斯消元,我们已经知道了矩阵的初等变换:
- 把两行(列)互换。
- 把某行(列)乘一非零常数 k
(其实是零也行,只是不叫初等变换)。
- 把某行(列)加上另一行(列)的 k 倍。
而我们又知道一个特殊的方阵——单位矩阵 I,它长这个样子:
\begin{bmatrix}
1&0&\cdots&0\\
0&1&\cdots&0\\
\vdots&\vdots&\vdots&\vdots\\
0&0&\cdots&1
\end{bmatrix}
它的性质是:任意矩阵乘上单位矩阵等于它本身。
如果作用初等变换于单位矩阵 I,就可得初等矩阵 E:
- 若作初等行变换,E\times A 相当于对 A 作对应的初等行变换。
- 若作初等列变换,A\times E 相当于对 A 作对应的初等列变换。
乘以零的情况与此相同。
因此高斯消元的过程,就是将矩阵分解为初等矩阵的过程。
算法的转置
若一个算法可以看作方阵 A,
输入为向量 \vec{v},输出为 A\vec{v},则称该算法为线性算法。
许多算法都是线性算法,例如 FFT 转化为方阵就是
\begin{bmatrix}
1&1&1&\cdots&1\\
1&w^1_n& w^2_n&\cdots&w^n_n\\
1&w^2_n& w^4_n&\cdots&w^{2n}_n\\
\vdots&\vdots&\vdots&\vdots&\vdots\\
1&w^n_n&w^{2n}_n&\cdots&w^{n^2}_n
\end{bmatrix}
此时,输入为 \vec{v},输出为 A^T\vec{v} 的算法为该算法的转置算法。
转置有几个非常重要的性质(有基础的同学可以跳过),那就是:
- $(AB)^T=B^TA^T$(转置调换矩乘的顺序)
- 对于初等矩阵的转置,
- 前两种的转置就是它本身。
- 最后一种如果是将第 $i$ 行(列)的 $k$ 倍加到第 $j$ 行(列),
其转置就是将第 $j$ 行(列)的 $k$ 倍加到第 $i$ 行(列)。
可以手推慢慢消化一下,这个非常重要。
## 线性算法的判定
然而线性算法的定义并未给我们带来什么别的东西,
这是因为算法的实际流程并未与矩阵建立起对应的关系。
自己有一个比较好的办法,
就是将输入、输出变量与运行过程中的一切辅助变量
(也就是整个内存)拼成一个向量,
也就是已知 $\begin{pmatrix}\underline{\vec{v}}\\\underline{\vec{0}}\\\vec{0}\end{pmatrix}$,求 $\begin{pmatrix}\underline{\vec{0}}\\\underline{\vec{0}}\\A\vec{v}\end{pmatrix}=\begin{bmatrix}O&O&O\\O&O&O\\A&O&O\end{bmatrix}\begin{pmatrix}\underline{\vec{v}}\\\underline{\vec{0}}\\\vec{0}\end{pmatrix}$,
(若 $A$ 为非方阵,也可通过补 $0$ 补成方阵)
这样的话。程序的各种运算与初等变换一一对应起来了。
**为保证行文统一,以下均称向量中的量为变量,矩阵中的量为常量。**
于是具体的,线性算法只包含以下运算:
- 交换两个变量。
- $x\leftarrow m\times x,x$ 为变量,$m$ 为常量。
- tips: $x\gets 0$ 属于该运算。
- $x\leftarrow x+m\times y,x,y$ 为变量,$m$ 为常量。
- tips: $x\gets y$ 为该运算与上方运算的**组合**,
即 $x\leftarrow 0\times x,x\leftarrow x+y
有了这个定义,我们便可以很方便地判定和处理线性算法了。
原算法与转置算法的转换
将矩阵分解为初等矩阵,
设 A=E_1E_2\cdots E_{n-1}E_n,E_1,E_2,\cdots E_n 均为初等矩阵,
则显然有 A^T=E^T_nE^T_{n-1}\cdots E^T_2E^T_1,
再由 E_k^T 性质可知,将算法转置 相当于将运算次序反过来,然后
由此我们得到将算法转置的~~很傻逼的~~机械过程。
# 常见算法的转置
如果上面的“将算法转置”还没有理解,可以看看下面的例子。
## FFT 的转置
上面我们已经知道了,FFT 可写作 $n$ 阶方阵 $F$,其中 $F_{i,j}=w_n^{ij}$ 。
因此 $F^T=F,F^{-1}_{i,j}=2^{-n}w_n^{-ij}$ 。
## 多项式乘法(卷积)的转置
考虑多项式乘法(卷积):
已知 $n+1$ 维向量 $\vec{a}$,$m+1$ 维向量 $\vec{b}$,求 $n+m+1$ 维向量 $\vec{c}$,
使 $c_k=(\vec{a}\times\vec{b})_k=\displaystyle\sum_{\begin{matrix}\scriptsize 0\le i\le n\\\scriptsize 0\le k-i\le m\end{matrix}}a_{i}b_{k-i}\ (0\le k\le n+m)$。
可以很快地写出这个算法:
```cpp
for(int k=0;k<=n+m;++k){
c[k]=0;
for(int i=min(0,k-m);i<=max(k,n);++i)
c[k]+=a[i]*b[k-i];
}
for(int i=0;i<=n;++i) a[i]=0;
```
把 $\vec{b}$ 看作矩阵中的常量,将算法化为线性算法的标准形式:
根据 $\vec{b}$ 构造矩阵 $B$ 使得 $B\begin{pmatrix}\underline{\vec{a}}\\\vec{0}\end{pmatrix}=\begin{pmatrix}\underline{\vec{0}}\\\vec{c}\end{pmatrix}$。
然后写出它的转置:
```cpp
for(int i=n;i>=0;--i) a[i]=0;
for(int k=n+m;k>=0;--k){
for(int i=max(k,n);i>=min(0,k-m);--i)
a[i]+=c[k]*b[k-i];
c[k]=0;
}
```
可以发现,这个算法就是求 $\begin{pmatrix}\underline{\vec{a}'}\\\vec{0}\end{pmatrix}=B^T\begin{pmatrix}\underline{\vec{0}}\\\vec{c}'\end{pmatrix},\vec{c}'$ 为输入向量。
$a'_{i}=\displaystyle\sum_{\begin{matrix}\scriptsize 0\le k\le n+m\\\scriptsize 0\le k-i\le m\end{matrix}}c'_{k}b_{k-i}\ (0\le i\le n)$。
考虑将 $\vec{b}$ 系数翻转,即设 $b^R_i=b_{m-i}$,
则 $a'_i=\displaystyle\sum_{\begin{matrix}\scriptsize 0\le k\le n+m\\\scriptsize i\le k\le m+i\end{matrix}}c'_{k}b^R_{m+i-k}=(c'\times b^R)_{m+i}\ (0\le i\le n)
因此卷积的转置就是已知 n+m+1 维向量 \vec{c}',m 维向量 \vec{b},
求 \vec{c}'\times\vec{b}^R 的第 m 个到第 m+n 个元素。
可以使用长度为 n+m+1 的循环卷积,因为溢出的部分刚好不要。
根据对多项式乘法(卷积)转置的研究,我们得到了以下性质:
- 将向量分段后,某几段分向量的赋值、加法与数乘也是线性的
(其实就是循环赋值、对应位置相加与乘上一个常数),可以当做初等变换看待。
- 转置时输入与输出的位置互换(由改写本身可证明),
因此转置算法时,通常会在最后将输出移动到输入处。
同时要注意以下几点:
- 转置算法的时候要把运算过程(包括内存的调用)一点一点写出来。
- 注意变量与常量的差别对待(变量写在向量里,常量写在矩阵里)
转置原理
由转置的过程本身可知,原算法与转置算法的复杂度相同。
若得到了计算 A\vec{v} 的优化方法,
则必然可以通过如上机械的改写得到 A^T\vec{v} 的优化方法。
反过来,如果要优化 A\vec{v},则可优先考虑 A^T\vec{v},
然后机械改写 A^T\vec{v} 的算法得到 A\vec{v} 的算法。
这便是转置原理,又称特勒根原理。
(也就是说,将算法转置后再优化就是转置原理)
多项式多点求值
回到我们最初的目的:优化多点求值。
多项式多点求值可化为范特蒙德矩阵
y_k=\sum_{i=0}^nf_ix_k^i
\begin{bmatrix}
1&x_0 &x_0^2 &\cdots&x_0^n\\
1&x_1 &x_1^2 &\cdots&x_1^n\\
1&x_2 &x_2^2 &\cdots&x_2^n\\
\vdots&\vdots&\vdots&\vdots&\vdots\\
1&x_m &x_m^2 &\cdots&x_m^n
\end{bmatrix}
\begin{pmatrix}
f_0\\f_1\\f_2\\\vdots\\f_n
\end{pmatrix}=
\begin{pmatrix}
y_0\\y_1\\y_2\\\vdots\\y_m
\end{pmatrix}
首先将该矩阵补全为方阵
(m 不够矩阵添 0,n 不够矩阵与向量一起添 0)
考虑转置该矩阵,
\begin{bmatrix}
1&1&\cdots&1\\
x_0&x_1&\cdots&x_m\\
x_0^2&x_1^2&\cdots&x_m^2\\
\vdots&\vdots&\vdots&\vdots\\
x_0^n&x_1^n&\cdots&x_m^n
\end{bmatrix}
发现转置算法是求 y'_k=\displaystyle\sum_{i=0}^nf_ix_i^k、
写成生成函数,得到 y'(x)=\displaystyle\sum_{i=0}^n\frac{f_i}{1-x_ix}。
我们一般用分治 FFT 解决这个问题。
具体地,设 \displaystyle\sum_{i=l}^{r}\frac{f_i}{1-x_ix}=\frac{f_{l,r}(x)}{g_{l,r}(x)},l,r 均为多项式函数,则答案为 \dfrac{f_{0,n}(x)}{g_{0,n}(x)}
设 l\le mid<r,则有
f_{l,r}(x)=f_{l,mid}(x)\times g_{mid+1,r}(x)+f_{mid+1,r}(x)\times g_{l,mid}(x)
g_{l,r}(x)=g_{l,mid}(x)g_{mid+1,r}(x)
化为线性算法的标准形式:
- 自下而上初始化常量(x_i 在矩阵中,被视为常量):
-
\vec{g}_{l,l}=\dbinom{1}{-x_i}
-
\vec{g}_{l,r}=\vec{g}_{l,mid}\times\vec{g}_{mid+1,r}
- 初始向量为 \begin{pmatrix}\underline{\vec{f}}\\\vec{0}\end{pmatrix},\vec{0} 为初始化的内存。
各个多项式不共用内存(即第一次调用时均为 \vec{0})。
- 赋值 \vec{f}_{l,l}\leftrightarrow(f_i)
- 自下而上分治(\times 为多项式乘法,由上已知它是线性的):
-
\vec{f}_{l,mid}\gets\vec{f}_{l,mid}\times\vec{g}_{mid+1,r}
-
\vec{f}_{mid+1,r}\gets\vec{f}_{mid+1,r}\times\vec{g}_{l,mid}
-
\vec{f}_{l,r}\gets\vec{f}_{l,r}+\vec{f}_{l,mid}
-
\vec{f}_{l,r}\gets\vec{f}_{l,r}+\vec{f}_{mid+1,r}
- 最后计算 \vec{f}\leftarrow\vec{f}_{0,n}\times(\vec{g}^{-1}_{0,n}\bmod x^{n+1})
- 将辅助变量清 0
转置该算法:
- 自下而上初始化常量 \vec{g}_{l,r}
- 初始向量为 \begin{pmatrix}\underline{\vec{f}}\\\vec{0}\end{pmatrix},\vec{0} 为初始化的内存。
- 将辅助变量清 0
- 计算 \vec{f}_{0,n}'\gets \vec{f}\times^T(\vec{g}^{-1}_{0,n}\bmod x^{n+1}),\times^T 为多项式乘法的转置。
- 自上而下分治:
-
\vec{f}'_{mid+1,r}\gets\vec{f}'_{mid+1,r}+\vec{f}'_{l,r}
-
\vec{f}'_{l,mid}\gets\vec{f}'_{l,mid}+\vec{f}'_{l,r}
-
\vec{f}'_{mid+1,r}\gets\vec{f}'_{mid+1,r}\times^T\vec{g}_{l,mid}
-
\vec{f}'_{l,mid}\gets\vec{f}'_{l,mid}\times^T\vec{g}_{mid+1,r}
- 赋值 (f_i)\leftrightarrow\vec{f}'_{l,l}
由此得到了多项式多点求值的 \Theta(n\log^2 n) 快速算法。
P7440「KrOI2021」Feux Follets
题解
转置原理的适用范围
通过 P7440,可总结出转置原理的适用范围。
通常转置原理要结合分治 FFT 和矩阵食用,
因此总是得到一个 \Theta(n\log^2 n) 的解法。
对于二元生成函数的相关题,
转置原理能将计算 [x^k] 的问题转化为计算 [y^k],从而优化递推。
因此遇到二元生成函数时,转置原理是一个不错的选择。
最后,希望更多的人能借我这篇拙劣的文章学会转置原理。QWQ