快速傅里叶变换

快速傅里叶变换(FFT)，常用于解答多项式乘法相关内容。
吐槽时间：
见了网上太多太多的资料都是物理的QAQ
不给OIers活路这样真的好么……
都是什么频率啊振幅啊什么乱七八糟的东西TAT
欺负我什么都不懂是吧QAQ
还是自己写一个比较好QAQ
背景故事：
在我们平时计算多项式乘法的时候，我们把第一个多项式的每一项都和第二个多项式的每一项相乘，复杂度为O(n ^ 2)，此时我们所使用的表示法就是系数表示法。
现在我们可以有一种比较强大的方式：
点值表示法：
众所周知，假设f(x)的最高次数为$n - 1$，即次数界为$n$，那么我们只要知道了n个不相同的x及f(x)值，就能确定出f(x)的多项式。
有一种算法叫做秦九韶算法，它能得出的结论是：
对于一个n次多项式，至多做n次乘法和n次加法。
所以我们知道了n个点的x值之后，我们在O(n^2)的时间内就能计算出所有的y值，然而如果通过快速傅里叶变换的方法，可以在O(nlogn)的时间内求出所有的y值。
引入新定义：
求值：通过多项式的系数表示法求其点值表示法。
插值：通过多项式的点值表示法求其系数表示法。
显然上面两个定义是互逆的关系，FFT就是用来解决求值的过程的。
引入新定义：
n次单位复根：在复平面内，n次单位复根能把整个平面分成n块。它的严格定义是：满足$\omega^n = 1$的复数$\omega$值，它一共有n个，分别为$\omega^k_n(k = 0...n - 1)$，其数值为$e^{2\pi ik / n}。$
由复数幂的定义，可知：
$(\omega _n^1)^k = \omega _n^k$
它有很多性质：
1.相消引理：

$$\forall n \in N,x\in N,d \in N^*,\omega_{dn}^{dk} = \omega_n^k$$
这是显然的，你可以考虑这个复平面等分成6个平面取第2个根和把这个复平面等分成3个平面取第一个根没什么区别。
2.折半引理：
$$\forall n \in even \bigcap n \in N ^ * ,(\omega_n^{k +n / 2}) ^2 = (\omega _n^k)^2$$
证明:左边 = $\omega_n^{2k + n} = \omega _n^{2k} * \omega _n^n = \omega _n^{2k} = 右边$
3.几何级数：
对于在复数域上的x $\not = 1$,有下列等式成立：
$$\sum ^{n}_{k = 0} x ^ k = 1 + x + x ^ 2 + x ^ 3 + x ^ n = \frac{x ^ {n + 1} - 1}{x - 1}$$
嗯其实就是等比数列求和。
4.求和引理：
对于任意：n >= 1 且 n不能整除非零整数k ,有：
$$\sum _{i = 0}^{n - 1}(\omega _n ^{k}) ^ i = 0$$
证明：
由几何级数的计算过程可知：
$$\sum _{i = 0}^{n - 1}(\omega _n ^{k}) ^ i = \frac {({\omega_n ^ k}) ^ {n} - 1}{\omega_n^k - 1} = \frac {(\omega _n^n) ^ k - 1}{\omega_n^k - 1} = 0$$
证毕。

---

接下来我们分析一下如何得到最终的多项式吧。
1.求A的n个单位根的点值，求B的n个单位根的点值
2.点值相乘，得到C的点值。
3.计算C的多项式。

---

我们先考虑第一步。
我们希望计算多项式：

$$A(x) = \sum _{j = 0}^{n - 1}a_jx^j$$
在n个单位根处的值。
设多项式A以系数形式给出，系数向量 $\vec a = {(a_0,a_1,...,a_{n - 1)}}$,定义其结果
$$y = A(\omega_n^k) = \sum _{j = 0}^{n - 1}a_j\omega_n^{kj}$$
则向量 $\vec y = {y_0,y_1,...,y_{n - 1}}$就是系数向量 $\vec a$的 离散傅里叶变换(DFT),记作：y = $DFT_n(a)$.

---

也就是说，点值的结果就是DFT，那么我们只需要快速计算DFT的值即可。
如果按照正常的算法，时间复杂度显然是$O(n^2)$的，所以有快速傅里叶变换(FFT)，它采取的是分治的思想。
我们考虑原来的多项式A(x)，定义两个新的次数界为n / 2的多项式：

$$A^{[0]}(x) = a_0 + a_2x + a _4x^2 + ... + a_{n - 2}x^{n / 2 - 1}$$
$$A^{[1]}(x) = a_1 + a_3x + a _5x^2 + ... + a_{n - 1}x^{n / 2 - 1}$$
则有
$$A(x) = A^{[0]}(x^2) + x * A^{[1]}(x^2)$$
我们如果只是这样的话是没办法减少计算量的，因为这只是相当于展开。
我们考虑一个显然的事实：
对于能够计算 $A(x)$的那两个多项式，实际上它们也是可以用FFT求值的， 它们的次数界减少了一半。
我们仔细回忆下单位根的性质，发现：
$$A(k) = A^{[0]}(\omega _n^ k)^2 + \omega _n^ k * A^{[1]}(\omega _n^ k)^2 ;$$
$$A(k + n / 2) = A^{[0]}(\omega _n^ {k + n / 2})^2 + \omega _n^ {k + n / 2} * A^{[1]}(\omega _n ^ {k + n / 2 })^2 ;$$
化简一下第二个式子，发现：
$$A(k + n / 2) = A^{[0]}(\omega _n^ {k})^2 - \omega _n^ {k} * A^{[1]}(\omega _n ^ {k})^2 ;$$
嗯……
$$设A(x) = u + x * v,则A(x + n / 2) = u - x * v.$$
这个玩意叫做 蝴蝶变换。
减少了一半的计算量。
我们来计算下复杂度：
T(n) = 2 * T(n / 2) + n ;
T(n / 2) = 2 * T(n / 4) + n / 2;
…
显然，复杂度为 $nlog_2n$。
我们终于讲完了如何求出n个单位根的点值，现在来说第二步：

---

把点值相乘。
显而易见的一点是，我们如果有两个多项式A(x),B(x),它们相乘得到多项式C(x),则：

$$A(k) * B(k) = C(k) ｛k为定义域上的任意值｝$$
那么我们把n个点值相乘的复杂度自然是O(n)的。
只需要记录c[i] = a[i] * b[i]即可。

---

接下来我们要做的就是插值了。[忘了定义的请自觉往回翻]
因为我不知道拉格朗日插值公式是什么鬼（划掉），所以我决定还是不讲了。
//别打我别打我……我讲…………
嗯我们考虑把DFT写成矩阵的形式：

$$\left [\begin{array} \\ y_0\\ y_1 \\ y_2 \\y_3\\...\\y_{n - 1} \end{array}\right] = \left [\begin{array} \\ 1 \ \ & \ 1\ &\ \ 1 &\ \ 1& ... & 1\\ 1 & \omega_n^1&\omega_n^2&\omega_n^3 & ... & \omega_n^{n - 1}\\ 1& \omega_n^2&\omega_n^4&\omega_n^6 & ... & \omega_n^{2(n - 1)} \\1& \omega_n^3&\omega_n^6&\omega_n^9 & ... & \omega_n^{3(n - 1)}\\...\\1& \omega_n^{n - 1}&\omega_n^{2(n - 1)}&\omega_n^{3(n - 1)} & ... & \omega_n^{(n - 1)^2} \end{array}\right]\left [\begin{array} \\ a_0\\ a_1 \\a_2 \\a_3\\...\\a_{n - 1} \end{array}\right]$$
终于挤下来了。
我们发现这个X的矩阵叫做范德蒙德矩阵。
发现我们现在要求的就是 $\vec a$，即 $DFT^{-1}_n(y)$。
那么怎么求系数矩阵呢？
则：若有矩阵AB = C，已知C,A则有：B = $A^{-1} C$
证明： $A ^{-1}* AB = B = A ^ {-1} * C$
证毕。
显然的一点是：
范德蒙德矩阵的逆是存在的，我们想求 $\vec a$也不是很难，但是，矩阵乘法在这里的复杂度是 $O(n^2)$的。
对于中间那个矩阵，我们设其为 $V_n$，则有：
$V_n^{-1}的(j,k)的元素为\omega_n^{-kj}/n$
证明：
……不想写了。
我们给定逆矩阵 $V_n^{-1}$，可以推导出 $DFT^{-1}_n(y)$:
$$a_j = \frac 1 n \sum _{k = 0}^{n - 1} y_k\omega^{-kj}_{n}.$$
然而对比之前的式子……？
$$y = A(\omega_n^k) = \sum _{j = 0}^{n - 1}a_j\omega_n^{kj}$$
长得一模一样……
用a数组直接存储a数组 * b数组的点值，对当前的a数组做FFT，对其结果除以n就能得到答案。
但是我们现在的东西还停留在递归阶段，我们知道，可以对这个数组进行元素上的调整，这样我们可以先计算那些子节点，省去了递归的过程，这个过程叫做：
雷(二)德(进)算(制)法
大概意思就是把每一个数的二进制头尾翻转过来：
即100 - > 001这样子……
我们可以知道的一件事情就是：0 - > 7这个序列，在递归到叶子节点之后，是：0,4,2,6,1,5,3,7。我们发现把二进制转过来之后按照顺序排序之后，再反转回来就能得到原来的序列了。

---

嗯……发现自己还是介绍一下雷德算法比较好，它有一个O(n)求每个数的二进制的反转的方法，但是请记住，这个O(n)并不如平时我们写的O(nlogn)要快，但是极其优美。
我们考虑一个数x的二进制反转和x / 2的二进制反转，发现如果x是偶数的话，那么[x / 2]的实际二进制和x的实际二进制只差了在二进制中移一次位。用实际二进制数组表达，那么设二进制数组为A[]，则A[x] = A[x >> 1] << 1;
在反转的二进制中,设其数组为Rev[]:

$$Rev[x] = Rev[x >> 1] >> 1;$$
发现，一个数如果是奇数，那么也仅仅是改变了右移之后的最左端的一位。因为：A[x + 1] = A[x] | 1，那么我们就能同样改变反转数组：
$$Rev[x] = Rev[x >> 1] >> 1 | ((i & 1) << (LG - 1))$$
左移的位数大家仔细考虑一下发现，我们只需要使用 < n 的那些数字，也就意味着，在0 - > 3区间内：
Rev[0] = 00,Rev[1] = 10,Rev[2] = 01,Rev[3] = 11。
发现如果我们左移Log次好像也是……
可以的……吗？
显然是不能的，因为我们希望的是：
在将其二进制反转之后，用O(n)的方法对其调整计算的顺序。
我们如果二进制每次左移log位，看样子是可行的，但是我们只用到了0 - > n的，而左移log位就会使得除0以外的二进制反转>=n，因为用快速排序会使得整个程序再次陷入 $O(nlog_2n)$，所以我们每次显然不能左移log位。
我们发现： 左移(log - 1)位之后，O(n)扫一遍整个序列就能排序。

void Rader()
{
    //m 是原来两个多项式中较长的那个的长度的2倍。
    for(n = 1;n <= m;n <<= 1,LG ++);
    for(int i = 1;i < n;i ++)
        Rev[i] = (Rev[i >> 1] >> 1) | ((i & 1) << (LG - 1));
}
void FFT(Virt *F,int ty = 1)
{
    Rep_0(i,n)if(i < Rev[i])swap(F[i],F[Rev[i]]);//每次FFT的过程中根据二进制的反转调整运算顺序。
    ...
}

这样的话我们合并相邻的两个，就像线段树一样，这样问题就成功解决了。
请用手写复数类并重载运算符。
请用手写复数类并重载运算符。
请用手写复数类并重载运算符。
题目：
BZOJ2179：给出两个n位10进制整数x和y，你需要计算x*y。
n <= 60000
//据说FFT不开二倍数组长度会跑得更快
//注意：以上内容不保证其正确性与真实性。
//注意：以上内容不保证其正确性与真实性。
//注意：以上内容不保证其正确性与真实性。
我们可以将两个数相乘写成这种形式：

$$c_i = \sum _{j = 0}^{i}a_j * b_{i - j}$$
然后直接输出就可以了。

bzoj2179,FFT模版，据说我这样不加注释会运行的更快。
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#define Rep(i,n) for(int i = 1;i <= n;i ++)
#define RepG(i,x) for(int i = head[x];~ i;i = edge[i].next)
#define RD(i,x,n) for(int i = x;i <= n;i ++)
#define Rep_0(i,n) for(int i = 0;i < n;i ++)
#define Rep0n(i,n) for(int i = 0;i <= n;i ++)
#define PI M_PI
using namespace std;
struct Virt 
{
    double r,i;
    Virt (double r = 0.0,double i = 0.0)
    {
        this -> r = r;
        this -> i = i;
    }
    Virt operator+(const Virt &x)
    {
        return Virt(r + x.r,i + x.i);
    }
    Virt operator-(const Virt &x)
    {
        return Virt(r - x.r,i - x.i);
    }
    Virt operator*(const Virt &x)
    {
        return Virt(r * x.r - i * x.i,x.r * i + r * x.i);
    }
};
const int N = 131072;
int n,m,LG,c[N],Rev[N];
Virt a[N],b[N];
char s[N >> 1];
void Rader()
{
    for(n = 1;n <= m;n <<= 1,LG ++);
    for(int i = 1;i < n;i ++)
        Rev[i] = (Rev[i >> 1] >> 1) | ((i & 1) << (LG - 1));
}
void FFT(Virt *F,int ty = 1)
{
    Rep_0(i,n)if(i < Rev[i])swap(F[i],F[Rev[i]]);
    for(int i = 1;i < n;i <<= 1)
    {
        Virt wn(cos(PI / i),ty * sin(PI / i));
        for(int j = 0;j < n;j += (i << 1))
        {
            Virt w(1,0);
            for(int k = 0;k < i;k ++,w = w * wn)
            {
                Virt u = F[j + k],v = F[j + k + i] * w;
                F[j + k] = u + v,F[j + k + i] = u - v;
            }
        }
    }
    if(ty == -1)Rep_0(i,n)F[i].r = F[i].r / n;
}
int main ()
{
    scanf("%d",&n);
    m = (n - 1) << 1;
    scanf("%s",s);
    Rep_0(i,n)a[i] = s[n - i - 1] - '0';
    scanf("%s",s);
    Rep_0(i,n)b[i] = s[n - i - 1] - '0';
    Rader();
    FFT(a),FFT(b);
    Rep0n(i,n)a[i] = a[i] * b[i];
    FFT(a,-1);
    Rep0n(i,m)c[i] = (int)(a[i].r + 0.5);
    Rep0n(i,m)
    {
        if(c[i] >= 10)
        {
            c[i + 1] += c[i] / 10,c[i] %= 10;
            if(i == m)m ++;
        }
    }
    for(int i = m;i >= 0;i --)printf("%d",c[i]);
    return 0;
}

992ms应该也不算特别慢了吧……
再来看一道模版题:
快速傅立叶之二
请计算C[k]=sigma(a[i]*b[i-k]) 其中 k < = i < n ，并且有 n < = 10 ^ 5。 a,b中的元素均为小于等于100的非负整数。
Input
第一行一个整数N,接下来N行，第i+2..i+N-1行，每行两个数，依次表示a[i],b[i] (0 < = i < N)。
Output
输出N行，每行一个整数，第i行输出C[i-1]。
原式 = > ?

$$C_k = \sum _{i = k}^{n - 1}a_i*b_{i - k} = \sum_{k}^{n - 1}a_i * Revb_{n + k - i - 1}$$
让我来冷静一下……
FFT解决的卷积不是：
$$F_k = \sum _{i = 0}^{k}{a_i * b_{k - i}}$$
话说这两个到底有什么相同点啊！
于是我重新冷静了一下：
我们一开始解决的 $F_k$在第k位存储的是能满足把a和b的第i位的元素看作 $a_i * x^i$这种形式，实际上根本没有所谓的 $x^i$，但是这给了我们一个提示……
从0 - > k：$F_k$存储的是a和b乘起来之后x的次数等于k的项的系数。
这次的题目k - > n：$C_k$存储的是a和b乘起来之后x的次数等于n + k的系数。
那就直接输出n - 1 - > 2 * (n - 1)位的结果即可。
其实这个题也给我提了个醒：
不要关心卷积的过程是什么，你只需要知道：
这个东西是一个卷积。
所以你只要求一下对应的次数的那些卷积就可以了。

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#define Rep(i,n) for(int i = 1;i <= n;i ++)
#define RepG(i,x) for(int i = head[x];~ i;i = edge[i].next)
#define RD(i,x,n) for(int i = x;i <= n;i ++)
#define Rep_0(i,n) for(int i = 0;i < n;i ++)
#define Rep0n(i,n) for(int i = 0;i <= n;i ++)
#define PI M_PI
using namespace std;
struct Virt 
{
    double r,i;
    Virt (double r = 0.0,double i = 0.0)
    {
        this -> r = r;
        this -> i = i;
    }
    Virt operator+(const Virt &x)
    {
        return Virt(r + x.r,i + x.i);
    }
    Virt operator-(const Virt &x)
    {
        return Virt(r - x.r,i - x.i);
    }
    Virt operator*(const Virt &x)
    {
        return Virt(r * x.r - i * x.i,x.r * i + r * x.i);
    }
};
const int N = 131072 * 2;
int n,m,LG,c[N],Rev[N];
Virt a[N],b[N];
char s[N >> 1];
void Rader()
{
    for(n = 1;n <= m;n <<= 1,LG ++);
    for(int i = 1;i < n;i ++)
        Rev[i] = (Rev[i >> 1] >> 1) | ((i & 1) << (LG - 1));
}
void FFT(Virt *F,int ty = 1)
{
    Rep_0(i,n)if(i < Rev[i])swap(F[i],F[Rev[i]]);
    for(int i = 1;i < n;i <<= 1)
    {
        Virt wn(cos(PI / i),ty * sin(PI / i));
        for(int j = 0;j < n;j += (i << 1))
        {
            Virt w(1,0);
            for(int k = 0;k < i;k ++,w = w * wn)
            {
                Virt u = F[j + k],v = F[j + k + i] * w;
                F[j + k] = u + v,F[j + k + i] = u - v;
            }
        }
    }
    if(ty == -1)Rep_0(i,n)F[i].r = F[i].r / n;
}
int main ()
{
    scanf("%d",&n);
    m = (n - 1) << 1;
    Rep_0(i,n)
        scanf("%lf%lf",&a[i].r,&b[n - i - 1].r);
    Rader();
    FFT(a),FFT(b);
    Rep0n(i,n)a[i] = a[i] * b[i];
    FFT(a,-1);
    for(int i = m / 2;i <= m;i ++)printf("%d\n",c[i] = (int)(a[i].r + 0.5));
    return 0;
}

参考资料：
ACdreamers的博客
http://blog.csdn.net/acdreamers/article/details/39005227
zky的博客
http://blog.csdn.net/iamzky/article/details/22712347
算法导论

//请扫描以下二维码关注本人微信。