FFT [TPLY]

题目链接
https://www.luogu.org/problemnew/show/1919
https://www.luogu.org/problemnew/show/3803

资料推荐

orz大佬博客
http://www.cnblogs.com/cjoieryl/p/8206721.html (orz YL大佬)
http://blog.csdn.net/iamzky/article/details/22712347 (超级易懂)

知识点
复数:
https://baike.baidu.com/item/%E5%A4%8D%E6%95%B0/254365?fr=aladdin
单位根(单位复数根)
https://baike.baidu.com/item/%E5%8D%95%E4%BD%8D%E6%A0%B9
多项式表示
http://blog.csdn.net/acdreamers/article/details/39005227
卷积
http://blog.csdn.net/bitcarmanlee/article/details/54729807

书籍推荐
ACM/ICPC 算法基础训练教程 7.4快速傅里叶变换
算法导论 30章多项式与快速傅里叶变换

觉得看不懂很正常,知识点很抽象需要逐步理解
而对于本文,本文作用为对两位大佬内容解读并且加入自己看法
读者可以先花时间阅读以上推荐的文字,产生一定理解再读本文.由于作者很弱,可能会产生错误,还请大神帮忙指出.

正文

Part1.初识FFT

作用:求多项式乘法的系数(就是初中学的多项式)
FFT多项式乘法与普通形式有差异

一般多项式乘多项式法则:
多项式与多项式相乘，先用一个多项式的每一项与另一个多项式的每一项相乘，再把所得的积相加。举一个例子由多项式乘多项式法则可以得到(a+b)(c+d)=a(c+d)+b(c+d)=ac+ad+bc+bd

FFT'绕圈'法
FFT处理多项式相乘,先从多项式最普通的表达形式---系数表达形式转化成点值表达形式再进行点值乘法.最后,再将点值乘法的结果转回系数表达形式.
而从-系数表达形式转化成点值表达形式的过程,我们叫做DFT,即离散傅里叶变换.
把点值乘法的结果转回系数表达形式的步骤,我们称之为离散傅里叶变换的逆运算,也叫插值.

下面再对一些可能不理解的概念进行解释
１～对于点值表达式我的理解 :
点值表达形式就是用若干函数上的点的坐标来表示该函数
比如 f(x)=\(x^2\)-x+2={ ( 0 , 2 ) , ( 1 , 2 ) , ( 2 , 4 ) }
但是这些点的数量一定不小于系数的个数n
这是一个常识
就像求解n元方程一定要有n组方程才能得解
２～点值表达式的计算:
假设我们现在有用点乘表示的多项式
A(x)={(0,1),(1,2),(2,7)}
B(x)={(0,2),(1,2),(2,4)}
[P.S.:点乘表达式能运算的条件是两个点要有相同的x值]
C(x)={(0,2×1),(1,2×2),(2,7×4)}={(0,2),(1,4),(2,28)}

FFT流程图如下:(借鉴他人博客)

part2.复数

0.为什么要学习复数?
YL大佬原话:"学习复数,所有这些七七八八的定理,都是为了FFT的奇妙变换!"

所以,为了FFT,好好学习复数吧!(再一次ORZ YL大佬)

1.复数的定义:我们把形如a+bi（a,b均为实数）的数称为复数，其中a称为实部，b称为虚部，i称为虚数单位,i的值为\(\sqrt{-1}\)
2.单位根的定义:数学上，n次单位根是n次幂为1的复数。
单位根的概念不理解很正常,下面举个例子(来自百度百科)
单位的一次根有一个：1
单位的二次根有两个：+1和-1。
单位的三次根是:
{$ 1 $, \(\frac{-1+\sqrt{3}i}{2}\),\(\frac{-1-\sqrt{3}i}{2}\)}
咱们证明一下
\((\frac{-1+\sqrt{3}i}{2})^3\)=\((\frac{-1+\sqrt{3}i}{2})^2\)×\((\frac{-1+\sqrt{3}i}{2})\)
....................=\((\frac{1+3i^2-2\sqrt{3}i}{4})\)×\((\frac{-1+\sqrt{3}i}{2})\)
....................=\((\frac{-2-2\sqrt{3}i}{4})\)×\((\frac{-1+\sqrt{3}i}{2})\)
....................=\((\frac{(-2)×(1+\sqrt{3}i)×(\sqrt{3}i-1)}{8})\)
....................=-\((\frac{1}{4})\)×(\(3i^2-1\))
....................=\(-(\frac{1}{4})\)*(-4)
....................=1
所以\(\frac{-1+\sqrt{3}i}{2}\)是单位的三次根
单位的四次根是{1,+i,-1,-i}

希望你已经理解了什么是单位复数根(单位根),咱们继续

强烈建议学习一下著名的欧拉公式，后面会用到。
单位复数根定义：对于w\(_n\)=\(e^{2πi/n}\)我们称w\(_n\)为n次单位根(前面讲了n次单位跟哦)。

下面有一个定理需要记住(YL大佬的忠告)
PS由于太弱不会证明,YL大佬的博客上由详细证明可以参考
(w上标表示指数)
消去引理
\[\omega^{dk}_{dn}=\omega^{k}_{n}\]
这个引理非常非常重要哦，你后面一定一定会见到！
非常，非常重要！记住！

这些复杂的数学公式可能会使你觉得枯燥无味,但为后面的运算奠定了重要的基础,了解一下下啦!

Part3.DFT与-DFT进阶

在本章内于括号里面的数为什么是\(\omega^{k}_{n}\)而不是x,你要到FFT才能明白.但本章内容中的\(\omega^{k}_{n}\)你可以当成x看(下一章就不行了)

DFT
对于k=0~n-1，定义：
\[y_k=A(\omega^{k}_{n}) = \Sigma^{n-1}_{j=0} a_j(\omega^{k}_{j})^j\]
对于这个式子,你这么看\(y_k=A(x)=\Sigma^{n-1}_{j=0} a_j(x)^j\)
我们知道DFT求的是点乘
点乘就是坐标
平时你是怎么求坐标的?
选定一个x值,把x值带到里面求出y
上面这个式子就是这么做的啦(还不懂就再看看多项式的系数表达式)

所以我们继续分析,对于每一个x,我们都要O(n)地将值带到多项式里求值,我们最少要n个坐标,所以总复杂度是O(n*n)

最后我们将得到的y称为a的离散傅里叶变换，记作\(y=DFT_n(a)\) \((这里的y，a指的是所有的y_k,a_k(k有n种,也就是要取k个))\)

这个定义先记着,还要注意的是,这一步的作用是将系数表达形式变成点乘表达形式\(O(n^2)\).
再用之前的点乘算法算出点乘的结果O(n)(Part1.讲了点乘,复杂度证明很显然).
最后用下面的离散傅里叶逆变换(我叫它-DFT),再换回去\(O(n^2)\)所以总复杂度为\(O(n^2)\)

-DFT
我这里直接摆出结果，有兴趣的可以自己证明一下
\[对于y_k = \Sigma^{n-1}_{i=0}a_i(x)^i\]
\[有a_k = \frac{1}{n}\Sigma^{n-1}_{i=0}y_i(x)^i\]
再对比DFT公式
\[y_k = \Sigma^{n-1}_{j=0} a_j(x)^j\]
可以发现,DFT和-DFT的差别很小,基本上只是y,a调换了一下罢了(也有不同的...)

Part4 FFT水落石出

终于来到了最后,前面所有令人头大的背景知识,都是为了这一章而服务.
所以这一章是至关重要的.
如果前面还没有完全理解,请重新巩固扎实.

网上都说,FFT的复杂度为O(nlgn)
那是怎么来的呢?
答案是:神奇的分治
先把A(x)这个多项式拆分成奇数项与偶数项
\[A^{[0]}(x) = a_0 + a_2x + a_4x^2 + ... + a_{n-2}x^{\frac{n}{2} - 1}\]
\[A^{[1]}(x) = a_1 + a_3x + a_5x^2 + ... + a_{n-1}x^{\frac{n}{2} - 1}\]
\[A(x) = A^{[0]}(x^2) + xA^{[1]}(x^2) \]
这个是正确的,希望读者拿起草稿纸自己演算一遍
咱们把这个多项式拆分成这样子以后
复数闪亮登场
开始用复数的公式了哦
咱们不妨把\(\omega^{k}\)带入到x中去
\[A(\omega^k_n)=A^{[0]}((\omega^k_n)^2) + \omega^k_n A^{[1]}((\omega^k_n)^2)\]
然后,通过之前提到过的单位复数根消去引理我们可以得到（说了消去引理非常重要）
咱们用消去引理把\((\omega^k_n)^2\)化一下
\(（\omega^k_n)^2\)=\(\omega^{2k}_n\)

.............=\(\omega^{\frac{2k}{2}}_{\frac{n}{2}}\)

.............=\(\omega^{k}_{\frac{n}{2}}\)
又因为\[\omega^{\frac{n}{2}}_{n}=\omega_{2}=e^{k\pi i}= cos k\pi + i sin k\pi = -1(有名的欧拉公式)\]
所以
\[\omega^{k+\frac{n}{2}}_n=\omega^{k}_n*\omega^{\frac{n}{2}}_n=-\omega^{k}_n\]
所以
\[A(\omega^{k+\frac{n}{2}}_n)=A^{[0]}((\omega^{k+\frac{n}{2}}_n)^2) + \omega^{k+\frac{n}{2}}_nA^{[1]}((\omega^{k+\frac{n}{2}}_n)^2)\]
\[=A^{[0]}((-\omega^k_{n})^2)-\omega^k_{n}A^{[1]}((-\omega^k_{n})^2)\]
\[=A^{[0]}(\omega^{2k}_{n})-\omega^k_{n}A^{[1]}(\omega^{2k}_{n})\]
\[ 再一次用到咱们的消去引理 \]
\[=A^{[0]}(\omega^k_{\frac{n}{2}})-\omega^k_{n}A^{[1]}(\omega^k_{\frac{n}{2}})。。。（式1）\]
然而我们同时也可以暴力把原式化一下：

\[A(\omega^k_n)=A^{[0]}((\omega^k_n)^2) + \omega^k_n A^{[1]}((\omega^k_n)^2)\]
\[=A^{[0]}(\omega^{2k}_n) + \omega^k_n A^{[1]}(\omega^{2k}_n)\]
\[ 又因为消去引理 \]
\[=A^{[0]}(\omega^{k}_{\frac{n}{2}}) + \omega^k_n A^{[1]}(\omega^k_{\frac{n}{2}})。。。（式2）\]

蝴蝶操作

所以我们把我们推了这么久的（式1）和（式2）拿出来并排看
\[A(\omega^k_n)=A^{[0]}(\omega^{k}_{\frac{n}{2}}) + \omega^k_n A^{[1]}(\omega^k_{\frac{n}{2}}) \]
\[A(\omega^{k+\frac{n}{2}}_n)=A^{[0]}(\omega^k_{\frac{n}{2}})-\omega^k_{n}A^{[1]}(\omega^k_{\frac{n}{2}})\]
然后我们发现
用\[A^{[0]}(\omega^{k}_{\frac{n}{2}}) 和 A^{[1]}(\omega^k_{\frac{n}{2}}) 和\omega^k_{n}\]
可以同时算出\[A(\omega^k_n)和A(\omega^{k+\frac{n}{2}}_n)\]
计算量便减小了一半
我们管这叫蝴蝶操作

蝴蝶操作流程如下

然后我们还可以发现
蝴蝶操作无非就是把奇数项和偶数项合并算出新项再往上继续往上合并直到我们得到最终的式子

也就是这样

是可以递归着做
不断分离奇数项偶数项最后再回溯合并
这是由上而下

那么能不能由下而上呢？
我们看最下层叶子节点数字的编号，通过找规律我们发现每个数和他二进制相反的位置互换
举个例子
4 -> 100
1 -> 001
所以4和1交换
（注意标号从0开始）

接下来。。。

我承认我蝴蝶操作讲的不清楚
由于知识有限只能讲到这里
网上关于蝴蝶操作的好的文章有很多。
大家可以去看一看
下面是一些例题
（链接再最上面）
可以把它当成板子背一背

代码

//P3803 【模板】多项式乘法（FFT）
#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <queue>
#include <cmath>#define rg register int
#define ll long long
#define db double
#define il inline#define INF 2147483647
#define SZ 10000001using namespace std;const double pi=acos(-1);struct Complex
{db r,i;Complex(){}Complex(db a,db b):r(a),i(b) {}Complex operator+(const Complex B)const{return Complex(r+B.r,i+B.i);}Complex operator-(const Complex B)const{return Complex(r-B.r,i-B.i);}Complex operator*(const Complex B)const{return Complex(r*B.r-i*B.i,r*B.i+i*B.r);}
};
Complex X,Y,a[SZ],b[SZ];int r[SZ],n,m,l;
il void FFT(Complex *a,int x)
{for(rg i=0;i<n;++i) if(i<r[i]) swap(a[i],a[r[i]]);for(rg i=1;i<n;i<<=1){Complex wn(cos(pi/i),x*sin(pi/i));for(rg j=0;j<n;j+=(i<<1)){Complex w(1,0);for(rg k=0;k<i;++k,w=w*wn){X=a[j+k],Y=a[i+j+k]*w;a[j+k]=X+Y,a[i+j+k]=X-Y;}}}if(x==-1) for(rg i=0;i<n;++i) a[i].r=a[i].r/n;
}int main()
{scanf("%d%d",&n,&m);for(rg i=0;i<=n;++i) scanf("%lf",&a[i].r);for(rg i=0;i<=m;++i) scanf("%lf",&b[i].r);m+=n;for(n=1;n<=m;n<<=1) ++l;for(rg i=0;i<n;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));FFT(a,1);FFT(b,1);for(rg i=0;i<=n;++i) a[i]=a[i]*b[i];FFT(a,-1);for(rg i=0;i<=m;++i) printf("%d ",(int)(a[i].r+0.5));return 0;
}

P1919 【模板】A*B Problem升级版（FFT快速傅里叶）
#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <queue>
#include <cmath>#define rg register int
#define ll long long
#define il inline
#define ldb long double#define INF 2147483647
#define SZ 1000001using namespace std;
const ldb pi=acos(-1);struct Complex
{ldb r,i;Complex(){}Complex(ldb a,ldb b):r(a),i(b){}il Complex operator+(const Complex B)const{return Complex(r+B.r,i+B.i);}il Complex operator-(const Complex B)const{return Complex(r-B.r,i-B.i);}il Complex operator*(const Complex B)const{return Complex(r*B.r-i*B.i,r*B.i+i*B.r);}
};
Complex X,Y,a[SZ],b[SZ];
int n,m,l,r[SZ];
il void FFT(Complex *a,int x)
{for(rg i=0;i<n;++i) if(i<r[i]) swap(a[i],a[r[i]]);for(rg i=1;i<n;i<<=1){Complex wn(cos(pi/i),x*sin(pi/i));for(rg j=0;j<n;j+=(i<<1)){Complex w(1,0);for(rg k=0;k<i;++k,w=w*wn){X=a[j+k],Y=a[i+j+k]*w;a[j+k]=X+Y,a[i+j+k]=X-Y;}}}if(x==-1) for(rg i=0;i<n;++i) a[i].r=a[i].r/n;
}char ch1[SZ],ch2[SZ];
int res[SZ];
int main()
{scanf("%d",&n);cin>>ch1;cin>>ch2;for(rg i=0;i<n;++i){a[i].r=ch1[n-i-1]-'0';b[i].r=ch2[n-i-1]-'0';}m=n+n-2;for(n=1;n<=m;n<<=1) ++l;for(rg i=0;i<n;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));FFT(a,1);FFT(b,1);for(rg i=0;i<=n;++i) a[i]=a[i]*b[i];FFT(a,-1);for(rg i=0;i<=m;++i) res[i]=(int)(a[i].r+0.5);for(rg i=0;i<=m;++i) res[i+1]+=res[i]/10,res[i]%=10;rg top=m+1;while(res[top]>9){res[top+1]=res[top]/10;res[top]%=10;++top;}while(!res[top]) --top;for(rg i=top;i>=0;--i)printf("%d",res[i]);return 0;
}

几点注意
pi 最好不要手写（当然如果你能背到小数点后十几位我也不栏你）
不然很可能出精度问题
最好手写Complex速度快很多
注意一些从0开始的地方
加油哦!

END