引入

第一篇博文，写一个FFT的简介。过几天肯定用得到，也符合我的爱好嘛……
首先，FFT是快速傅里叶变换的简称。傅里叶变换是一个数学中使用的对函数的变换，而离散傅里叶变换顾名思义就是对于离散点的傅里叶变换。
那么，在OI中，傅里叶变换的最大意义，就在于加速向量卷积，即多项式乘法。

多项式乘法

假如我们有两个多项式：

$$A(x)=\sum_{i=0}^{n}a_i*x^i B(x)=\sum_{i=0}^{n}b_i*x^i$$

我们想求得：

$$C(x)=\sum_{i=0}^{2n-1} (\sum_{j=0}^{i}a_j*b_{i-j}) x^i$$

使用朴素的做法肯定需要 $O(n^2)$ 的复杂度。
我们将上述的 $A(x)$ 的系数表示为向量 $\boldsymbol{A}$ ，则称这种表示法为系数表示法。
现在来考虑另一种表示法，即点值表示法。
设向量 $\boldsymbol{A} = (A(p_1),A(p_2),A(p_3),...,A(p_n))$ ，则 $\boldsymbol{A}$ 称作 $A(x)$ 的点值表示。
点值表示有什么优点？首先，n维点值表示与n维系数表示一一对应（证明似乎有点麻烦……），其次，点值表示下的多项式运算得到了化简。
以乘法为例：

$$\begin{eqnarray} C(x)&=&\sum_{i=0}^{2n-1} (\sum_{j=0}^{i}a_j*b_{i-j}) x^i \\ &=&\sum_{i=0}^{2n-1} \sum_{j=0}^{i}(a_j*x^j)*(b_{i-j}*x^{i-j}) \\ &=&\sum_{j=0}^{2n-1}(a_j*x^j)*\sum_{i=0}^{2n-1}(b_{i}*x^{i}) \end{eqnarray}$$

这样便转换为了点值表示的两个式子。而且我们可以发现，完成点值表示下的乘法复杂度是 $O(n)$ 的。
那么加速多项式乘法的过程，变为了：
系数表示 $\rightarrow$ 点值表示 $\rightarrow$ 乘法 $\rightarrow$ 点值表示 $\rightarrow$ 系数表示
而快速傅里叶变换，就是从系数表示转换到点值表示的桥梁。
我们观察到，点值表示下所取的点 $p_1,p_2,...,p_n$ 是可以任取的，那么我们利用一些特殊值，是否能加速运算？

快速傅里叶变换

此处我们引入单位根的概念。即满足：

$$x^n=1$$

的复数，记为 $\omega_n$。它可以看做是复平面上x轴正方向绕逆时针方向旋转 $\frac{2\pi}{n}$ 的复数。
单位根有一些很好的性质，这些性质都可以利用复平面上的向量比较直观地感受。

$$\begin{eqnarray} \omega_{2n}^{2m}&=&\omega_{n}^{m}\\ \omega_{2n}^{m}&=&-\omega_{2n}^{m+n} \end{eqnarray}$$

假如我们取单位根的幂进行转换，会有什么效果？设 $A_0(x)$ 是 $A(x)$ 的偶次项的和， $A_1(x)$ 是奇次项的和。那么：

$$\begin{eqnarray} A(\omega_{n}^{m})&=&A_0((\omega_{n}^{m})^{2})+\omega_{n}^{m}A_1((\omega_{n}^{m})^{2}) \\ &=& A_0(\omega_{\frac{n}{2}}^{m})+\omega_{n}^{m}A_1(\omega_{\frac{n}{2}}^{m}) \\ A(\omega_{n}^{m+\frac{n}{2}})&=&A_0((\omega_{n}^{m})^{2})+\omega_{n}^{m+\frac{n}{2}}A_1((\omega_{n}^{m})^{2}) \\ &=& A_0(\omega_{\frac{n}{2}}^{m})-\omega_{n}^{m}A_1(\omega_{\frac{n}{2}}^{m}) \end{eqnarray}$$

即我们只要有了 $A_0(x),A_1(x)$ 的点值表示，就能在 $O(n)$ 时间内算出 $A(x)$ 的点值表示。
每一层对于一个确定的 $m$，我们用对应的两个值来更新出当前的值，称这个操作为“蝴蝶变换”。
分析一下复杂度，$A(x)$ 的规模为 $T(n)$ ，而他的奇、偶次项分别的和的规模都为 $T(\frac{n}{2})$ ，而合并的时间为 $O(n)$ ，即 $T(n) = 2T(\frac{n}{2}) + O(n)$ ，根据主定理，这个问题的复杂度为 $O(n \log n)$.

逆离散傅里叶变换

现在我们已经可以快速从系数表示转换到点值表示了。接下来的问题就是如何从点值表示转换到系数表示。
我们转移到点值表示的过程可以看做是一个矩阵乘法：

$$\left( \begin{array}{ccccc} \omega_{n}^{0} & \omega_{n}^{0} & \omega_{n}^{0} & ... & \omega_{n}^{0} \\ \omega_{n}^{0} & \omega_{n}^{1} & \omega_{n}^{2} & ... & \omega_{n}^{n} \\ \omega_{n}^{0} & \omega_{n}^{2} & \omega_{n}^{4} & ... & \omega_{n}^{2n} \\ ... & ... & ... & ... & ... \\ \omega_{n}^{0} & \omega_{n}^{n-1} & \omega_{n}^{2n-2} & ... & \omega_{n}^{n^2-n} \end{array} \right) * \left( \begin{array}{c} a_0 \\ a_1 \\ a_2 \\ ... \\ a_{n-1} \end{array} \right) = \left( \begin{array}{c} b_0 \\ b_1 \\ b_2 \\ ... \\ b_{n-1} \\ \end{array} \right)$$

那么关键就是找到左边那个矩阵的逆矩阵。这个逆矩阵是存在的……算导上也给出了这个矩阵，即对他的每一个元素取共轭复数，然后除n。为什么会这么想？我们需要逆矩阵的一行乘原矩阵的一列恰为 $[i=j]$ ，由于是单位根所以有：

$$\omega_{n}^{m} * \omega_{n}^{-m} = |\omega_{n}^{m}|=1 \\ \sum_{i=0}^{n-1}\omega_{n}^{i}=[n=1] \\ \omega_{n}^{m}=\omega_{n}^{m\ mod\ n}$$

代入后发现当 $i = j$ 时，由上性质，正好约去。而 $i \neq j$ 时，我们将所有得到的结果的指数都化到模 $n$ 意义下，由数论知识，恰为一组 $0 \sim n-1$ 的自然数，由于 $n \neq 1$，值变为0.
这样，我们只需要把带进去的那个 $\omega_{n}^{m}$ 的向量翻转一下，然后进行FFT即可。

迭代版FFT

依上述内容便可实现递归版的FFT了，但是这样常数很大。
观察到我们的分治是将奇偶分类，即对于每一层我取这层对应位数的0放在一起，1放在一起。
那么最后元素 $i$ 所处的位置，就是将 $i$ 的二进制表示翻转的位置。
这里我觉得引入 Cooley-Tukey 算法的图比较好：

可以看到，我们每次取对应的两个元素，进行一次“蝴蝶变换”，即完成快速傅里叶变换章节中的对两个元素进行的操作。

Number Theory Transform

复数运算的FFT，常数巨大，精度不高。
当题目所给的模数非常带感时：即 $2^k|P-1$，其中 $2^k > n$ 时，我们便可以使用数论变换来加速运算。数论变换中的关键为一个等价：

$$g^{\frac{P-1}{m}}\equiv \omega_{m} \pmod{P}$$

其中 $g$ 是 $P$ 的原根。
你可以手动验证一下上述等价的正确性。
原根的求法不赘述了，枚举即可……

Code

说了这么多，理解起来可能也有点困难，倒是代码比较简短。

void FFT_Init() {
    for ( K = 1; K < n << 1; K <<= 1 ); inv_K = inver( K );
    w[ 0 ][ 0 ] = w[ 0 ][ K ] = w[ 1 ][ 0 ] = w[ 1 ][ K ] = 1;
    int G = fpm( g, ( P - 1 ) / K );
    FOR( i, 1, K - 1 ) {
        w[ 0 ][ i ] = (ll)w[ 0 ][ i - 1 ] * G % P;
    }
    FOR( i, 0, K ) {
        w[ 1 ][ i ] = w[ 0 ][ K - i ];
    }
}

void FFT( int X[], int k, int v ) {
    int i, j, l;
    for ( i = j = 0; i < k; i++ ) {
        if ( i > j ) swap( X[ i ], X[ j ] );
        for ( l = k >> 1; ( j ^= l ) < l; l >>= 1 );
    }
    for ( i = 2; i <= k; i <<= 1 )
        for ( j = 0; j < k; j += i )
            for ( l = 0; l < i >> 1; l++ ) {
                int t = (ll)X[ j + l + ( i >> 1 ) ] * w[ v ][ ( K / i ) * l ] % P;
                X[ j + l + ( i >> 1 ) ] = ( (ll)X[ j + l ] - t + P ) % P;
                X[ j + l ] = ( (ll)X[ j + l ] + t ) % P;
            }
}