FFT学习笔记
前置知识：
一.复数
1.基础：见高中必修下基础知识
2.复数的表示法
复数有三种表示方法，分别为：
1.标准形式：\(z=a+bi\)，其中 \(i\) 为虚数单位。若 \(a=0\)，则该复数为纯虚数。
2.三角形式：\(z=r(\sin \theta + i \cos \theta )\)。$\theta $ 是复数的幅角（或称为辐角、相位角），它是复数与正实轴之间的夹角。
3.指数形式：\(z=re^{i \theta }\)。这个形式较难理解，但是我们可以通过泰勒展开，然后可以将其和三角形式一一对应。指数形式一般不会用到，在OI中三角形式为最常用的形式。
3.棣莫弗定理
公式：
1.\({(r(\cos \theta+i \sin \theta))}^n = r^n (\cos n \theta + i \sin n\theta )\)
2.\(r_1(\cos \theta_{1}+i\sin \theta_1) r_2(\cos \theta_{2}+i\sin \theta_2)=r_1r_2(\cos (\theta_1+\theta_2) + i \sin (\theta_1+\theta_2))\)
考虑对定理的理解：两复数相乘，即把对应在半径相乘，角度相加即可。当这个半径为 \(1\) 的时候，我们可以看做就是两个角度相加。
此外，我们可以知道：复数乘上它的逆等于半径在平方。
4.单位根（蝶形因子）
以单位圆点为起点，单位圆的 \(n\) 等分点为终点，在单位圆上可以得到 \(n\) 个复数，设幅角为正且最小的复数为 \(\omega\)，称为 \(n\) 次单位根，即 \(\omega_n\)。
则 \(\omega_n^k\) 表示 \(n\)次单位根中第 \(k\) 个单位根。
关于单位根的性质，有以下三条：
1.\(\omega_n^n=\omega_n^0=1\)
2.\(\omega_n^k=\omega_{2n}^{2k}\)
3.\(\omega_{2n}=-\omega_{2n}^k\)
关于这三条性质的证明，我们都可以用欧拉公式将 \(\omega_n^k\) 写成 \(\cos \frac {2 \pi} {n}\) 的形式，化简去证明。在OI方面，通过图像感性理解即可。
5.C++的STL复数容器：complex
C++有自带的复数容器，有以下几个功能：
1.创建复数：创建复数支持三种表示方法，如：
complex<double> c(5,3)
即为创建了一个复数 \(5+3i\)。
我们也可以:
complex<double> st(cos(pi),sin(pi))
这样用三角表示去创建复数，当然，指数表示也是可以的。
2.输入输出：复数可以直接输入与输出，如：
cin>>c;
cout<<c;
输出的结果（以上文的c为例）为：
(5,3)
输入是一定要加上两个圆括号的，否则C++是无法识别出你虚部的内容的。
3.访问实部和虚部：
访问实部： double k_real=c.real()
访问虚部： double k_imag=c.iamg()
4.虚数支持直接加减与乘除，使用四边形法则，除此之外，复数还支持求其共轭(conj)，求其模长(abs)，求其辐角(arg)等操作。
二.多项式表示法
1.两种表示方法
1.系数表示法：对于任何 \(n\) 次多项式 \(f(x)=\sum_{i=0}^{n} a_ix_i\)，我们只需要知道每一项的系数 \(a_i\) 就可以唯一确定这个多项式，故多项式 \(f(x)\) 可以被表示为 \((a_1,a_2,a_3,……,a_n)\) 的形式，即其系数的有序排列，这就是多项式在系数表示。 （其实将其倒过来看，这个函数 \(f(x)\) 可以看作是这个数列的生成函数）
2.点值表示法：我们有结论：\(n\) 次多项式上已知 \(n+1\) 个不同的点可以唯一确定这个多项式。此结论可以通过因式定理与代数基本定理，或者用范德蒙特行列式的性质证明。继而 \(f(x)\) 可以被表示为 \({(x_0,f(x_0)),(x_1,f(x_1)),……,(x_n,f(x_n))}\)。此即为多项式在点值表示法。
3.优劣：系数表示法更直观，而点值表示法在处理多项式相乘时只需要 \(O(n)\) 的时间，因为可以直接将对应的点值相乘，而系数表示则需要用矩阵乘法，时间复杂度为 \(O(n^2)\)。因此，我们在处理多项式乘法时可以多考虑点值表示。
2.两种表示方法转化
1.暴力转化：两者转化的时间复杂度都为 \(O(n^2)\)。将点值表示转化为系数表示可以用 \(O(n^3)\) 的高斯消元，也可以用 \(O(n^2)\) 的矩阵乘法;将系数表示转化为点值表示则可以直接代入，单次代入时间复杂度为 \(O(n)\)，总时间复杂度为 \(O(n^2)\)。
2.优化：快速傅里叶变换：可以将两者的转化时间复杂度都降低到 \(O(n \log n)\)。
快速傅里叶变换
DFT，即离散傅里叶变换。即计算两个多项式的乘积。暴力计算的时间复杂度显然为 \(O(n^2)\) 。而FFT，快速傅里叶变换（Fast Fourier Transform），可以在 \(O(n \log n)\) 的时间内完成这个运算。
我们考虑有多项式 \(f(x)=\sum_{i=0}^{n}a_ix^i\)，令 \(n=2^s\) ，我们就可以把 \(f(x)\) 看作是两个部分 \(F(x)\) 与 \(G(x)\) 相加的结果，这两个部分是按照 \(a_i\) 下标的奇偶来分的。
令：\(F(x)\) 为 \(a_0 x^0+a_2 x^2+a_4 x^4 +…+a_{n} x^{\frac {n} {2} }\)，\(G(x)\) 为 \(a_1 x^0+a_3 x^2+a_5 x^4+……+a_{n-1}^{\frac {n} {2} -1}\)。
故：\(f(x)=F(x^2)+xG(x^2)\)。
我们可以发现，这样只需要求出 \(F(x^2)\) 与 \(G(x^2)\)，我们就可以很简明的求出 \(f(x)\)。但若是直接随便找几个值作为 \(x\) 代入到函数之中，我们会发现求解还是很麻烦。考虑去寻找特殊的 \(x\) 代入。
考虑单位根。代入 $x=\omega_n^{k} $ （\(k < \frac {n} {2}\)），有：\(f(\omega _n^k)=F(\omega_n ^{2k})+\omega _n^k G(\omega _n ^{2k})=F(\omega_{\frac {n} {2}} ^{k})+\omega _n^k G(\omega_{\frac {n} {2}} ^{k}))\)。
考虑代入 \(\omega _n ^{k+\frac {n}{2}}\)，有：\(f(\omega _n ^{k+\frac {n}{2}})=F(\omega _n ^{2k+n})+\omega _n ^{k+\frac {n}{2}}G(\omega _n ^{2k+n})\)。
根据复数相乘的性质，有：\(f(\omega _n ^{k+\frac {n}{2}})=F(\omega _n ^{2k} \cdot \omega _n^n)-\omega _n ^{k} G(\omega_n^{2k} \cdot \omega_n^n)\)。
根据单位根的性质，我们显然可以得出：\(f(\omega _n ^{k+\frac {n}{2}})=F(\omega _{\frac {n}{2}}^k)-\omega _n ^{k} G(\omega _{\frac {n}{2}}^k)\)。
观察第一式与最后一式的结构，我们发现只需要求出 \(F(\omega_{\frac {n} {2}}^k)\) 与 \(G(\omega_{\frac {n} {2}}^k)\) 就可以 \(O(1)\) 求出 \(f(\omega _n^k)\)。显然，这个过程可以递归分解下去，最终的问题变为求 \(F(\omega_1^k)=G(\omega_1^k)=1\)。
递归下去的时间复杂度为 \(O(\log n)\)，故求出所有答案的时间复杂度为 \(O(n \log n)\)。即我们把系数表示转化为点值表示的时间复杂度从 \(O(n^2)\) 降低到了 \(O(n \log n)\)。
快速傅里叶逆变换
考虑将点值表示转化为系数表示，这个过程我们称为快速傅里叶逆变换。
考虑求DFT的过程，用矩阵表示为：
\(\begin{bmatrix} 1& 1& 1& …& 1\\ 1& (\omega _n^1)^1& (\omega _n^1)^2& …& (\omega _n^1)^{n-1}\\ 1& (\omega _n^2)^1& (\omega _n^2)^2& …& (\omega _n^2)^{n-1}\\ …& …& …& … & …\\ 1& (\omega _n^{n-1})^1& (\omega _n^{n-1})^2& …&(\omega _{n-1}^{n-1})^n\end{bmatrix}\begin{bmatrix} a_0\\ a_1\\ a_2\\ …\\a_{n-1} \end{bmatrix}
=\begin{bmatrix} y_0\\ y_1\\ y_2\\ …\\y_{n-1}\end{bmatrix}\)
这其中矩阵 \(\begin{bmatrix} a_0\\ a_1\\ a_2\\ …\\ a_{n-1} \end{bmatrix}\) 为系数表示，矩阵 \(\begin{bmatrix} y_0\\ y_1\\ y_2\\ …\\ y_{n-1} \end{bmatrix}\) 为点值表示。
由于：\(x^n-1=(x-1)(x^{n-1}+x^{n-2}+…+1)\)。
所以当 \(x=1\) 时，有：\(x^{n-1}+x^{n-2}+…+1=n\)。
除此，当 \(x^n-1=0\) 时，有：\(x^{n-1}+x^{n-2}+…+1=0\)。
由前文所述，\(x\) 为 \(n\) 次单位根。
由单位根的性质，有：\(\begin{aligned} {\left[\begin{array}{lllll} 1 & \left(\omega_{n}^{n-i}\right)^{1} & \left(\omega_{n}^{n-i}\right)^{2} & \cdots & \left(\omega_{n}^{n-i}\right)^{n-1} \end{array}\right]\left[\begin{array}{c} 1 \\ \left(\omega_{n}^{j}\right)^{1} \\ \left(\omega_{n}^{j}\right)^{2} \\ \vdots \\ \left(\omega_{n}^{j}\right)^{n-1} \end{array}\right] } & =\sum_{k=0}^{n-1}\left(\omega_{n}^{n-i}\right)^{k}\left(\omega_{n}^{j}\right)^{k} \\ & =\sum_{k=0}^{n-1}\left(\omega_{n}^{n-i+j}\right)^{k} \\ & =\left\{\begin{array}{ll} n & i=j \\ 0 & i \neq j \end{array}\right. \end{aligned}\)
因此，我们就找到了 \(\begin{bmatrix} 1& 1& 1& …& 1\\ 1& (\omega _n^1)^1& (\omega _n^1)^2& …& (\omega _n^1)^{n-1}\\ 1& (\omega _n^2)^1& (\omega _n^2)^2& …& (\omega _n^2)^{n-1}\\ …& …& …& … & …\\ 1& (\omega _n^{n-1})^1& (\omega _n^{n-1})^2& …&(\omega _{n-1}^{n-1})^n \end{bmatrix}\) 的逆矩阵：\(\frac{1}{n} \begin{bmatrix} 1& 1& 1& …& 1\\ 1& (\omega_n^{n-1})^1& (\omega_n^{n-1})^2& …& (\omega_n^{n-1})^{n-1}\\ 1& (\omega_n^{n-2})^1& (\omega_n^{n-2})^2& …& (\omega_n^{n-2})^{n-2}\\ …& …& …& …& … \\ 1& (\omega_n^{1})^1& (\omega_n^{1})^2& …& (\omega_n^{1})^{n-1} \end{bmatrix}\)
该矩阵亦可被写成如下形式：
\(\left[\begin{array}{ccccc} 1 & 1 & 1 & \cdots & 1 \\ 1 & \left(\overline{\omega_{n}^{1}}\right)^{1} & \left(\overline{\omega_{n}^{1}}\right)^{2} & \cdots & \left(\overline{\omega_{n}^{1}}\right)^{n-1} \\ 1 & \left(\overline{\omega_{n}^{2}}\right)^{1} & \left(\overline{\omega_{n}^{2}}\right)^{2} & \cdots & \left(\overline{\omega_{n}^{2}}\right)^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \overline{\left(\omega_{n}^{n-1}\right)^{1}} & \overline{\left(\omega_{n}^{n-1}\right)^{2}} & \cdots & \overline{\left(\omega_{n}^{n-1}\right)^{n-1}} \end{array}\right]\)
其实用汉字来描述这个过程就是将一个多项式在分治的过程中乘的单位根变为其共轭复数，然后分治完的每一项除上 \(n\)。
这样我们就得到了系数表示。
FFT的优化
三步变两步
说起来是很简单一个东西：显然，我们有：
\((a+bi)^2=(a^2-b^2)+2abi\)
所以对于多项式 \(f(x),g(x)\)，我们将 \(g(x)\) 放到 \(f(x)\) 的虚部，求出 \(f(x)^2\) 后把 \(f(x)\) 虚部除以 \(2\) 即可得到 \(f(x) \cdot g(x)\)。
这样我们只需要做两边FFT即可。
蝶形变换
考虑优化FFT过程，发现不断递归的过程很繁琐。考虑优化。通过观察（可以简单证明），我们发现 \(\omega_n^i\) 的最终位置为 \(i\) 二进制翻转下得到的数。
于是我们可以提前处理好这些数，然后不断向上还原即可。
代码
模板题：Luogu P1919
#include<bits/stdc++.h>
#define int long long
#define Comp complex<double>//复数
using namespace std;
const int N=5e6+100;
const double pi=acos(-1);//圆周率pi
int n;
string s1,s2;
Comp a[N],b[N];
int ans[N];
int rev[N];
void init(int k)
{
int len=1<<k;
for(int i=0;i<len;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
}//预处理出需要翻转的数
void FFT(Comp *a,int n,int inv)
{
for(int i=0;i<n;i++)
if(i<rev[i])//保证只翻转一次
swap(a[i],a[rev[i]]);
for(int mid=1;mid<n;mid*=2)
{
Comp st(cos(pi/mid),inv*sin(pi/mid));//单位根
for(int i=0;i<n;i+=mid*2)//合并到了第i位
{
Comp cur(1,0);
for(int j=0;j<mid;j++,cur*=st)
{
Comp x=a[i+j],y=cur*a[i+j+mid];
a[i+j]=x+y,a[i+j+mid]=x-y;//蝶形变换
}
}
}
if(inv==-1)
for(int i=0;i<n;i++)
a[i]/=n;
}
signed main()
{
//cin>>n;
cin>>s1>>s2;
int n1=s1.size(),n2=s2.size();
for(int i=0;i<n1;i++) a[i]=(double)(s1[n1-i-1]-'0');
for(int i=0;i<n2;i++) b[i]=(double)(s2[n2-i-1]-'0');
int k=1,s=2;
while((1<<k)<2*max(n1,n2)-1) k++,s<<=1;
init(k);
FFT(a,s,1),FFT(b,s,1);
//把a,b转化为点值表示法
for(int i=0;i<s;i++) a[i]*=b[i];
FFT(a,s,-1);
//把a转化为系数表示法
for(int i=0;i<s;i++)
ans[i]+=(int)(a[i].real()+0.5),ans[i+1]+=ans[i]/10,ans[i]%=10;
while(s!=-1 && !ans[s]) s--;
if(s==-1) return cout<<0<<endl,0;
for(int i=s;i>=0;i--) cout<<ans[i];
return 0;
}
浙公网安备 33010602011771号