[多项式全家桶] 快速傅立叶变换(FFT)

写在之前

因为想要补周五的ECFinal pk赛中2017-2018 ACM-ICPC, Asia Daejeon Regional Contest Rock Paper Scissors，所以打算学习一下FFT的相关知识。然后一不小心看到了很多聚聚们写的其他相关博客，因而又发现了一片有趣的天地所以打算逐渐学习一下多项式全家桶。

---

FFT 快速傅立叶变换

很多博客的算法都是知其然不知其所以然，我认为自为风月马前卒dalao的这一篇博客写的很漂亮
（因为在处理范德蒙德行列式的逆的时候用到了一种蒙人的做法所以很干净简洁）

先大体用我粗糙的语言描述一下我对FFT算法的理解：
一个$n-1$阶的多项式具有$n$个自由变量，也就是说只需要用$n$个自由量就可以进行刻画该多项式。我们所熟悉的表示办法是用$n$个系数进行刻画，但是这种表示办法在计算多项式乘法的时候时间复杂度$O(n^2)$不可接受，所以我们考虑用一种更有利的表示来加速多项式乘法。我们按照一定的规律在多项式上选取$n$个值，用这$n$个值来刻画多项式。问题的要害在于如何选取对应的$x[i]$可以使得对应的$n$个值既好求又具有一些便利的性质可以方便我们的计算。

FFT算法的提出者使用了$n$阶单位根作为$x$的取值，即$x[i]=w_n^i,i=0,1,$$\cdots$$,n-1$.再分离奇偶项将原多项式化归成两个小多项式的和，利用单位根的性质即可在$O(nlogn)$的时间下将系数多项式转化为点值多项式。之后可以将两个多项式进行$O(n)$的乘法获得乘积多项式。
之后的问题就是如何将点值多项式还原到系数多项式，这个过程称为傅立叶逆变换。事实上，我们构造一个n阶矩阵，可以表示出从自变量到因变量之间的变换关系，容易发现这个n阶矩阵是一个范德蒙德矩阵，其中$xi$分别一一被取为$n$次单位根中的一个。

盗来的图：

那么如何从因变量和变换关系求得自变量？那自然是对矩阵取逆。该矩阵有一个优良的性质，就是其逆恰好为每一项取共轭再以$n$，证明方法我不会…因为偷瞄了一眼感觉很复杂的亚子就没有看下去。不过在这个问题的严谨证明上，自为风月马前卒给出了一种办法。

$设(y_1,y_2,...,y_{n-1})$为$(a_1,a_2,...,a_{n-1})$的点值表示，设有另一个向量$(c_1,c_2,...,c_{n-1})$满足
$$c_k=\sum _{i=0}^{n-1}{y}_i(w_n^{-k}{)}^i$$

即是以$y[i]$作为系数、自变量取为n阶单位根的倒数的新多项式的系数。

早晨在还没有思考到逆矩阵的时候看到了这种突如其来的构造实在让人一惊，尽管当时已经有所怀疑它与逆有关，但却不能肯定。当仔细学习完FFT的整个知识内容后才发现，这就是逆矩阵构成的变换关系(但每一项尚未除以$n$)！之后就使用了这个多项式进行推导，得出其与原多项式系数的关系。经过循循善诱的构造与重重艰辛的证明之后，我们得出关系$$a_i=\frac{c_i}{n}$$至此，逆FFT的过程得到了一个较好的解决。

当我们将$n$拓展为2的幂次之后，经过奇偶分离最终得到的其实是一个确定的序列，也就是说，该变换实际上是一个原序列的同构。经过观察可知，底层序列的下标是原序列下标的二进制翻转。因而我们可以利用这个规律更快地进行序列的操作。

关键代码

项数的处理

//分别读入n次方与m次方的系数
int p = 0;
m = n + m; for(n = 1; n <= m; n <<= 1) p++;

注意，如果是两个多项式做乘法，那么应该考虑的是$n+m+1$的最靠近的比他大的二次方数，这也就说明了为什么上面写的是$n<=m$了。这样子等价于$n>m$即$n>=m+1$时跳出循环，这是符合我们的要求的。
同时容易知道，因为要存储的是$n+m+1$最靠近的二的幂次的个数个元素，所以空间(不仅仅是两个多项式A、B，而且inv数组也是！)应该开为$(n+m+1)<<1$，不然会出现RE等问题！

翻转序列的获取：

for(int i = 0; i <= n; i++) inv[i] = (inv[i >> 1] >> 1) | ((i & 1) << (p - 1));

时间复杂度优秀，利用之前获得的结果进行计算。

FFT

void fft(Complex * a,int f, int n){
    for(int i = 0; i < n; i++) if(i < inv[i]) std::swap(a[i], a[inv[i]]);
    for(int i = 1; i < n; i <<= 1){
        Complex wn = {cos(pi / i), f * sin(pi / i)};
        for(int j = 0, r = i << 1; j < n; j += r){
            Complex w{1, 0};
            for(int k = 0; k < i; k++, w = wn * w){
                Complex x = a[k + j], y = w * a[k + j + i];
                a[k + j] = x + y; a[k + j + i] = x - y;
            }
        }
    }
}

有一些值得注意的地方⚠️

1. 在三层循环中，$i$是某一层的半长度，这一点的作用在第三层循环内的蝴蝶操作中可以很清晰地看出。$j$跑遍该层对应的区间所有起点，因为$i$是半长，所以应该再开一个$r=i<<1$来作为累加的全周期长度；在每一个二层循环内，都应该开设一个0度的n阶单位根作为初始单位根。第三层循环中的$k$起到一个计数器的功能，用于遍历半长区间并进行修改操作；随着k的加$1$，w也在不断的逆时针旋转，当$k$取到$i$时，w也转遍了所有的无重角度；在三层循环之内，因为在修改之前$a[i]$数组中存储的是$w_{\frac{n}{2}}^m$即$n-1$阶单位根，所以在当前的$n$阶单位根时，只需要直接获取数组中的原元素即可迭代更新。
2. 板子、原理要熟悉。$x$直接取为更小的情况，$y$要再乘上一个$w$.
3. 因为$i$是半长，所以$i$对应的$n$阶单位根的角度实际上应该是$$\frac{2\pi }{2i}=\frac{\pi }{i}$$这一个地方要牢记，不然很容易出现问题。
   
   

点值求积操作

for(int i = 0; i < n; i++) a[i] = a[i] * b[i];

答案输出

for(int i = 0; i <= m; i++) printf("%d ", (int)(a[i].x / n + 0.5));

范围int和longlong也要注意，这里的范围采取的是int。四舍五入是为了保证精度以免发生一些截断误差etc。

洛谷P3803 AC代码

参考博客：
倾心讲解 FFT 多项式与快速傅里叶变换 && 迭代模板带注释 （模版学习自此）
快速傅里叶变换(FFT)详解 (逆矩阵的严谨证明学习于此)
十分简明易懂的FFT（快速傅里叶变换） (单位根的性质说明较为清晰)

---

补题：Rock Paper Scissors

AC代码

写的比较惨。。花了很久的时间de了无数的小bug总算是过去了。。
总结一下踩到的坑。。这样子以后尽量避免再去犯错误：

1. 在进行多次FFT时要进行初始化，而该初始化是对复数的初始化，所以不仅仅要对x初始化也要考虑y的初始化！
2. 如果想要将一个长度为$n$的字符串转化为一个$n-1$次多项式而进行操作n--，那么前后一定要注意好当前的n与原始的n所推导出的公式之间的关系，仔细而妥善的处理好其中的等价变换。
3. 原板子中的n、m变换实在是不堪阅读，我们还是应该再开设一个意义较为明确的变量(比如N)来存储最接近n + m的幂次。与此同时我们也要再一次明确N代表的含义——N代表的是拓展后的元素个数！也就是说我们在进行与N有关的循环时应该遍历0到N-1，虽然遍历到N并不会引起过多的错误(前提是数组大小不能开的恰恰好)，但是这样子写不利于展现我们对于该算法的理解，所以为此还是应该写的专业一些。
4. 针对于本题而言，允许m区间超过n区间的末端，也就是说可以仅匹配前若干项并且取得较大值。因而，我们应该将第一个字符串延长后再进行匹配。因为n、m均在$1e5$以内，所以有关FFT数组大小至少应该开到$6e5$！同时延长的字符串大小应该开到$2e5$！