多项式与快速傅里叶变换学习笔记
就是抄了一遍书，，，
多项式
是以\(x\)为变量的多项式\(A\)的形式和。
系数表示
这个多项式可以用它的系数\(a\)来表示为\(a=(a_0,a_1,\cdots, a_{n-1})\)，这被称为多项式的系数表示法，\(a\)通常被视为一个向量。其次数为\(n-1\)，次数界为\(n\).
两个用系数表示法表示的多项式\(A(x)\)和\(B(x)\)的乘积\(C(x)\)是一个多项式，其系数向量被称作输入向量\(a\)和\(b\)的卷积，其中\(a,b\)分别是\(A(x)\)和\(B(x)\)的系数向量。
点值表示
一个次数界为\(n\)的多项式\(A(x)\)的点值表示是一个由\(n\)个点值对所组成的集合
使得对\(k=0,1,\cdots,n-1\)，所有的\(x_k\)各不相同。
插值&拉格朗日插值法
由一个多项式的点值表示确定其系数表示被称为插值。
定理(插值多项式的唯一性）：对于任意\(n\)个点值对组成的集合\(\{(x_0,y_0),(x_1,y_1),\cdots,(x_{n-1},y_{n-1})\}\)，其中所有的\(x\)都不同，那么存在唯一的次数界为\(n\)的多项式\(A(x)\),满足\(y_k=A(x_k),k=0,1,\cdots,n-1\)
证明：式\((*)\)等价于矩阵方程
左边的矩阵表示为\(V(x_0,x_1,\cdots,x_{n-1})\)，被称为范德蒙德矩阵，其行列式的值为
这个可以由数学归纳法证得。
因此，若\(\forall 0\le j < k\le n-1,x_k\ne x_j\)，则该矩阵的行列式不为0，为可逆矩阵。因此给定点值表示能够唯一确定系数\(a_j\)
说白话就是把\(n\)对点带回去能解方程，方程有唯一解。解这个方程的过程就是一种插值算法了，，
更快的插值算法是拉格朗日插值法。
套公式复杂度\(O(n^2)\)
用点值表示搞多项式运算非常方便，比如做乘法\(C(x)=A(x)B(x)\),这样直接能搞出n个点对，但\(C(x)\)的度数应该是\(A(x)\)和\(B(x)\)的度数和，如果\(A(x)\)和\(B(x)\)的次数界都为\(n\)，那么\(C\)的次数界为\(2n\)，需要\(2n\)个点值对才能进行插值以确定系数表达。因此，我们需要先对\(A\)和\(B\)的点值表示进行拓展，使得每个多项式都包含\(2n\)个点值对，然后才能求出能用于插值的\(C(x)\)的点值表示。
DFT与FFT
复数
首先回顾一下复数的运算规则。
复数具有一个实部和虚部，你可以把它理解为复平面上的一个向量。其加减运算即为实部虚部分别加减，乘法运算为模长相乘，幅角相加。
C++ STL里给出了复数，但咱也可以自己写一个。
单位复数根
\(n\)次单位复数根是满足\(w^n=1\)的复数\(w\).
\(n\)次单位复数根有恰好\(n\)个：对于\(k=0,1,\cdots, n-1\)，这些根是\(e^{\frac{2\pi ki}{n}}\)。
考虑复数的指数形式的定义
这说明\(n\)个单位复数根均匀分布在以复平面的原点为圆心的单位半径的圆周上。值
被称为主\(n\)次单位根，所有其他的\(n\)次单位根都是\(w_n\)的幂次。\(w_n^n=w_n^0=1\).
\(n\)个\(n\)次单位复数根在乘法意义下形成一个群。即\(w_n^jw_n^k=w_n^{(j+k)\mod n}\)
消去引理：
对任何整数\(n\ge 0, k \ge 0\)，以及\(d > 0\)，
推论：
对任意偶数\(n>0\)，有
折半引理
如果\(n>0\)为偶数，那么\(n\)个\(n\)次单位复数根的平方的集合就是\(n/2\)个\(n/2\)次单位复数根的集合。
证：根据消去引理，对任意非负整数\(k\)，有\((w_n^k)^2=w_{n/2}^k\).如果对所有的\(n\)次单位复数根进行平方，那么获得每个\(n/2\)次单位根正好两次，因为
故\(n\)为偶数时，\(w_n^k\)与\(w_n^{k+n/2}\)平方相同。
求和引理
对任意整数\(n\ge 1\)和不能被\(n\)整除的非负整数\(k\)，有
把复数看成复平面上的向量来理解的话上述引理显然。
DFT
给出次数界为\(n\)的多项式的系数表示，求该多项式在\(n\)个\(n\)次单位复数根处对应的值，记
\(y_k=A(w_n^k)=\sum_{j=0}^{n-1}a_jw_n^{kj}\)
那么向量\(y=(y_0,y_1,\cdots,y_{n-1})\)就是系数向量\(a=(a_0,a_1,\cdots,a_{n-1})\)的离散傅里叶变换。
FFT
通过快速傅里叶变换（\(FFT\)），我们能够在\(O(n\log n)\)的时间内计算出\(DFT\)，前提假设是\(n\)是\(2\)的整数次幂。
\(FFT\)将\(A(x)\)中偶数下标的项和奇数下标的项的系数分成两组，分别分配给两个次数界为\(n/2\)的多项式\(A^{[0]}(x)\)和\(A^{[1]}(x)\)。
于是求\(A(x)\)在\(w_n^0,w_n^1,\cdots,w_n^{n-1}\)处的值的问题转换为：
1.求次数界为\(n/2\)的多项式\(A^{[0]}(x)\)和\(A^{[1]}(x)\)在点
处的取值。
根据消去引理，上面实际上是\(n/2\)个的\(n/2\)次单位复数根，且每个根出现了两次。问题从求多项式在\(n\)次单位复数根处的值转为了求在\(n/2\)次单位复数根处的值，规模减半。因此现在可以递归地对次数界为\(n/2\)的多项式\(A^{[0]}(x)\)和\(A^{[1]}(x)\)在\(n/2\)个\(n/2\)次单位复数根处进行求值，
2.通过上面的式子计算\(A(x)\)
注意\(w_n^{n/2}=-1\)
于是能写出递归版的FFT
void fft(Complex *a, int n) { //n是2的整数次幂
if (n == 1) return;
static Complex buf[maxn];
int m = n / 2;
for (int i = 0; i < m; i++) {
buf[i] = a[i * 2];
buf[i+m] = a[i * 2 + 1];
//分奇偶存
}
for (int i = 0; i < n; i++)
a[i] = buf[i];
fft(a, m);
fft(a + m, m);
Complex wn = Complex(cos(2 * PI / n), sin(2 * PI / n)), w = Complex(1,0);
for (int i = 0; i < m; i++) {
buf[i] = a[i] + w * a[i+m];
buf[i+m] = a[i] - w * a[i+m];
w = w *wn;
}
for (int i = 0; i < n; i++)
a[i] = buf[i];
}
IDFT
考虑把DFT写成矩阵形式\(y=V_na\)，其中\(V_n\)是由\(w_n\)的适当幂次填充而成的范德蒙德矩阵，其第\(i\)行第\(j\)列的元素值为\(w_n^{(i-1)(j-1)}\)。
\(DFT\)的逆变换被记为\(IDFT\)，其可以表示为\(IDFT_y=a=V_n^{-1}y\)，\(V_n^{-1}\)的第\(i\)行第\(j\)列的元素值为\(\frac{w_n^{-(i-1)(j-1)}}{n}\).因此只要对上面的FFT代码进行修改，即可求解\(IDFT\).
void fft(Complex *a, int n, bool f) { //n是2的整数次幂,f用于判断当前求解的是dft还是idft
if (n == 1) return;
static Complex buf[maxn];
int m = n / 2;
for (int i = 0; i < m; i++) {
buf[i] = a[i * 2];
buf[i + m] = a[i * 2 + 1];
}
for (int i = 0; i < n; i++)
a[i] = buf[i];
fft(a, m, f);
fft(a + m, m, f);
Complex wn = Complex(cos(2 * PI / n), sin(2 * PI * f/ n)), w = (Complex){1, 0};
for (int i = 0; i < m; i++) {
buf[i] = a[i] + w * a[i+m];
buf[i+m] = a[i] - w * a[i+m];
w = w * wn;
}
for (int i = 0; i < n; i++) {
a[i] = buf[i];
if (f == -1) a[i].x /= n;
}
}
卷积定理
对任意两个长度为n的向量\(a,b\)，其中\(n\)是2的幂，\(a\)和\(b\)的卷积
\(a*b=DFT_{2n}^{-1}(DFT_{2n}(a)\cdot DFT_{2n}(b))\)
其中向量\(a,b\)用0填充，使其长度到达\(2n\),"\(\cdot\)"表示向量点乘。
迭代版FFT
for (int i = 0; i < m; i++) {
buf[i] = a[i] + w * a[i+m];
buf[i+m] = a[i] - w * a[i+m];
w = w * wn;
}
这里的操作被称为蝴蝶操作，考虑递归版FFT的递归过程：
我们按照元素在叶子节点中出现的顺序排列元素，成对地取出元素，利用蝴蝶变换计算出每一对元素的\(DFT\)，然后再用这\(n/2\)对的DFT代替他们自身，重复上述过程\(n/2\)次（\(n\)是2的幂)，即可得到原输入向量的\(DFT\)
位逆序置换
考虑如何计算元素在叶中出现的顺序：
在最顶层，下标的二进制表示末位为0的元素（偶数下标）被分到了左子树中，下标的二进制表示末位为1的元素（奇数下标）被分到了右子树中。去掉二进制的末位之后，我们就能在子树中重复上述过程。因此，我们能够根据下标的二进制表示计算出其在叶节点中出现的次序。
考虑从底层走向顶层的过程，假设底层下标为\(j\)的元素对应顶层的下标为\(i\)，那么从顶层走向底层的过程，就是逆序遍历\(i\)的二进制表示的每一位的过程。可以发现，\(j=rev(i)\)，其中\(rev(i)\)表示\(i\)的二进制表示各位逆序后得到的整数。
那么可以写出迭代FFT的代码如下(洛谷FFT模板)
#include<bits/stdc++.h>
using namespace std;
const int maxn = (1 << 21) + 10;
const double PI = acos(-1.0);
struct Complex {
double x, y;
Complex() {x = y = 0;}
Complex(double _x, double _y) {
x = _x, y = _y;
}
}a[maxn], b[maxn];
Complex operator + (Complex a, Complex b) {return Complex(a.x+b.x,a.y+b.y);}
Complex operator - (Complex a, Complex b) { return (Complex){a.x - b.x, a.y - b.y};}
Complex operator * (Complex a, Complex b) { return (Complex){a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x};}
Complex conj(Complex a) {return (Complex){a.x, -a.y};}
void fft(Complex *a, int n, int f) {
for (int i = 0, j = 0; i < n; i++) {
//i为原始下标，j为DFT的递归树最底层的下标。
if (i > j) swap(a[i], a[j]); //i比j小的话，i会在某次交换中也被换到它该待的位置。
for (int l = n >> 1; (j ^= l) < l; l >>= 1);
}
//这样循环保证了j的二进制表示是i的二进制表示的逆序，非常nb
for (int i = 2; i <= n; i <<= 1) { //区间长度
int m = i >> 1;
Complex wn = Complex(cos(2*PI/i), sin(2*PI*f/i)), t, w, u;
for (int k = 0; k < n; k += i) {//区间起点
w = Complex(1, 0);
for (int j = 0; j < m; j++) {
t = w * a[j + k + m];
u = a[j + k];
a[j + k] = u + t;
a[j + k + m] = u - t;
w = w * wn;
}
}
}
if (f == -1) {
for (int i = 0; i < n; i++) a[i].x /= n;
}
}
int rd() {
int s = 0, f = 1; char c = getchar();
while (c < '0' || c > '9') {
if (c == '-') f = -1;
c = getchar();
}
while (c >= '0' && c <= '9') {
s = s * 10 + c - '0';
c = getchar();
}
return s * f;
}
int n, m, Lim;
int main() {
n = rd(), m = rd();
for (Lim = 1; Lim < n + m + 1; Lim <<= 1);
for (int i = 0; i <= n; i++) a[i].x = rd();
for (int i = 0; i <= m; i++) b[i].x = rd();
fft(a, Lim, 1); fft(b, Lim, 1);
for (int i = 0; i < Lim; i++) a[i] = a[i] * b[i];
fft(a, Lim, -1);
for (int i = 0; i <= n + m; i++)
printf("%d ", (int)floor(a[i].x + 0.5));
return 0;
}
