快速傅里叶变换

快速傅里叶变换(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块。它的严格定义是：满足ωn=1的复数ω值，它一共有n个，分别为ωkn(k=0...n−1)，其数值为e2πik/n。
由复数幂的定义，可知：
(ω1n)k=ωkn
它有很多性质：
1.相消引理：

∀n∈N,x∈N,d∈N∗,ωdkdn=ωkn

这是显然的，你可以考虑这个复平面等分成6个平面取第2个根和把这个复平面等分成3个平面取第一个根没什么区别。
2.折半引理：

∀n∈even⋂n∈N∗,(ωk+n/2n)2=(ωkn)2

证明:左边 = ω2k+nn=ω2kn∗ωnn=ω2kn=右边
3.几何级数：
对于在复数域上的x≠1,有下列等式成立：

∑k=0nxk=1+x+x2+x3+xn=xn+1−1x−1

嗯其实就是等比数列求和。
4.求和引理：
对于任意：n >= 1 且 n不能整除非零整数k ,有：

∑i=0n−1(ωkn)i=0

证明：
由几何级数的计算过程可知：

∑i=0n−1(ωkn)i=(ωkn)n−1ωkn−1=(ωnn)k−1ωkn−1=0

证毕。

---

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

---

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

A(x)=∑j=0n−1ajxj

在n个单位根处的值。
设多项式A以系数形式给出，系数向量a⃗ =(a0,a1,...,an−1),定义其结果

y=A(ωkn)=∑j=0n−1ajωkjn

则向量y⃗ =y0,y1,...,yn−1就是系数向量a⃗ 的离散傅里叶变换(DFT),记作：y = DFTn(a).

---

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

A[0](x)=a0+a2x+a4x2+...+an−2xn/2−1

A[1](x)=a1+a3x+a5x2+...+an−1xn/2−1

则有

A(x)=A[0](x2)+x∗A[1](x2)

我们如果只是这样的话是没办法减少计算量的，因为这只是相当于展开。
我们考虑一个显然的事实：
对于能够计算A(x)的那两个多项式，实际上它们也是可以用FFT求值的，它们的次数界减少了一半。
我们仔细回忆下单位根的性质，发现：

A(k)=A[0](ωkn)2+ωkn∗A[1](ωkn)2;

A(k+n/2)=A[0](ωk+n/2n)2+ωk+n/2n∗A[1](ωk+n/2n)2;

化简一下第二个式子，发现：

A(k+n/2)=A[0](ωkn)2−ωkn∗A[1](ωkn)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;
…
显然，复杂度为nlog2n。
我们终于讲完了如何求出n个单位根的点值，现在来说第二步：

---

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

A(k)∗B(k)=C(k)｛k为定义域上的任意值｝

那么我们把n个点值相乘的复杂度自然是O(n)的。
只需要记录c[i] = a[i] * b[i]即可。

---

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

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢y0y1y2y3...yn−1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢1  111...1 1 ω1nω2nω3nωn−1n  1ω2nω4nω6nω2(n−1)n  1ω3nω6nω9nω3(n−1)n...............1ωn−1nω2(n−1)nω3(n−1)nω(n−1)2n⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢a0a1a2a3...an−1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥

终于挤下来了。
我们发现这个X的矩阵叫做范德蒙德矩阵。
发现我们现在要求的就是a⃗ ，即DFT−1n(y)。
那么怎么求系数矩阵呢？
则：若有矩阵AB = C，已知C,A则有：B = A−1C
证明：A−1∗AB=B=A−1∗C
证毕。
显然的一点是：
范德蒙德矩阵的逆是存在的，我们想求a⃗ 也不是很难，但是，矩阵乘法在这里的复杂度是O(n2)的。
对于中间那个矩阵，我们设其为Vn，则有：
V−1n的(j,k)的元素为ω−kjn/n
证明：
……不想写了。
我们给定逆矩阵V−1n，可以推导出DFT−1n(y):

aj=1n∑k=0n−1ykω−kjn.

然而对比之前的式子……？

y=A(ωkn)=∑j=0n−1ajωkjn

长得一模一样……
用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(nlog2n)，所以我们每次显然不能左移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不开二倍数组长度会跑得更快
//注意：以上内容不保证其正确性与真实性。
//注意：以上内容不保证其正确性与真实性。
//注意：以上内容不保证其正确性与真实性。
我们可以将两个数相乘写成这种形式：

ci=∑j=0iaj∗bi−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_PIusing 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]。
原式 = > ?

Ck=∑i=kn−1ai∗bi−k=∑kn−1ai∗Revbn+k−i−1

让我来冷静一下……
FFT解决的卷积不是：

Fk=∑i=0kai∗bk−i

话说这两个到底有什么相同点啊！
于是我重新冷静了一下：
我们一开始解决的Fk在第k位存储的是能满足把a和b的第i位的元素看作ai∗xi这种形式，实际上根本没有所谓的xi，但是这给了我们一个提示……
从0 - > k：Fk存储的是a和b乘起来之后x的次数等于k的项的系数。
这次的题目k - > n：Ck存储的是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_PIusing 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
算法导论

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