快速傅里叶变换 FFT
快速求多项式的卷积 (记 $H(x)$ 为结果, 则 $H(x)=F(x)G(x)=\sum_{i=0}^n\sum_{j=0}^i(f_j+g_{i-j})x^i$), 其中多项式最高次为 $n-1$ 次.
例题
建议结合例题 脱发 欣赏.
前置知识
复数, 多项式的点值 / 系数表示法,$n$ 次单位根. 不过不知道也没关系, 因为下面都有讲.
复数
定义
$$ z=a+bi=(\sqrt{a^2+b^2},\arctan\frac ba) $$
其中定义虚数单位 $i=\sqrt{-1}$,$a$ 称为实部 ($\text{real}$),$b$ 称为虚部($\text{imaginary}$, 也记作 $\text{imag}$), 后半段是极坐标表示法(模长和角度, 但是用不上). 如果将其看作向量 $(a,bi)$, 就可以在图上表示出来(实轴上方的那个向量)
共轭复数
对于复数 $z=a+bi=(a,\theta)$, 其共轭复数 $\bar z=a-bi=(a,-\theta)$, 如之前那张图上的两个复数互为共轭复数.
基本运算
记 $z=a+bi=(a,\theta)$, 则有
$$ z_1+z_2=(a_1+a_2)+(b_1+b_2)i\\ z_1z_2=(ac-bd)+(ad+bc)i=(a_1a_2,\theta_1+\theta_2) $$
复数加法满足平行四边形法则.
Code
在<complex>库中提供了复数模板类complex<>, 当然如果没有氧气优化的话还是建议自己写一个以免被卡常
struct comqlex {
double real,imag;
comqlex(double comq_r=0.0,double comq_i=0.0)
:real(comq_r),imag(comq_i) {};
comqlex(comqlex &comq_op) {
real=comq_op.real; imag=comq_op.imag;
}
comqlex &operator =(const comqlex &comq_op) {
if (this!=&comq_op) {
real=comq_op.real; imag=comq_op.imag;
} return *this;
}
bool operator ==(const comqlex &comq_op) const {
return real==comq_op.real&&imag==comq_op.imag;
}
bool operator !=(const comqlex &comq_op) const {
return real!=comq_op.real||imag!=comq_op.imag;
}
comqlex operator +(const comqlex &comq_op) {
comqlex comq_rt=*this; return comq_rt+=comq_op;
}
comqlex operator -(const comqlex &comq_op) {
comqlex comq_rt=*this; return comq_rt-=comq_op;
}
comqlex operator *(const comqlex &comq_op) {
comqlex comq_rt;
comq_rt.real=real*comq_op.real-imag*comq_op.imag;
comq_rt.imag=real*comq_op.imag+imag*comq_op.real;
return comq_rt;
}
comqlex operator /(const comqlex &comq_op) {
comqlex comq_rt;
comq_rt.real=(real*comq_op.real+imag*comq_op.imag)
/(comq_op.real*comq_op.real+comq_op.imag*comq_op.imag);
comq_rt.imag=(real*comq_op.imag+imag*comq_op.real)
/(comq_op.real*comq_op.real+comq_op.imag*comq_op.imag);
return comq_rt;
}
comqlex &operator +=(const comqlex &comq_op) {
real+=comq_op.real; imag+=comq_op.imag;
return *this;
}
comqlex &operator -=(const comqlex &comq_op) {
real-=comq_op.real; imag-=comq_op.imag;
return *this;
}
comqlex &operator *=(const comqlex &comq_op) {
return *this=(*this)*comq_op;
}
comqlex &operator /=(const comqlex &comq_op) {
return *this=(*this)/comq_op;
}
}系数 / 点值表示法
$n-1$ 次多项式可以由 $n$ 对 $(x,f(x))$ 值唯一表示(当然 $x$ 互不相同). 感性理解一下: 这 $n$ 对值给出了 $n$ 个方程,$F(x)$ 总共有 $n$ 个系数(包括常数项), 也就是 $n$ 个方程解 $n$ 个未知数, 显然都可以解出来. 这就是 点值表示法, 记作
$$ F(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),\dots,(x_{n-1},f(x_{n-1}))\} $$
另外, 对应地, 称
$$ F(x)=f_0+f_1x+\dots+f_{n-1}x^{n-1} $$
为 系数表示法 .
$n$ 次单位根
$\omega_n^i$ 表示将单位圆分成 $n$ 等份, 从 $(1,0)$ 开始沿逆时针方向使 $i=0,1,\dots,n-1$. 当然也存在 $i\not\in[0,n)$, 这时候有 $\omega_n^i=\omega_n^{i\bmod n}$(每一个 $n$ 次就是转了一整圈又回到了起点). 注意这个 $i$ 也就是 $\omega_n$ 的 $i$ 次方, 而不是简单的记号. 显然(虽然窝并不清楚为什么) 有
$$ e^{iu}=\cos u+i\sin u,\quad\omega_n^1=e^\frac{2\pi}ni $$
然后可以接着推出以下性质,$\texttt{(1)}$ 用第一条写开即可得证
$$ \omega_n^k=\cos\frac{2k}n\pi+i\sin\frac{2k}n\pi\\ |\omega_n^k|=\cos^2\frac{2k}n\pi+\sin^2\frac{2k}n\pi=1\\ \omega_n^k=\omega_{dn}^{dk},d\in\mathbb N^*\quad\texttt{(1)}\\ \omega_n^{k+\frac n2}=-\omega_n^k\quad\texttt{(2)} $$
最后一条的简要证明如下
$$ \omega_n^\frac n2=\cos\frac{\frac n2}n2\pi+i\sin\frac{\frac n2}n2\pi=\cos\pi+\sin\pi=-1\\ \omega_n^{k+\frac n2}=\omega_n^k\omega_n^\frac n2=-\omega_n^k $$
整体思路
暴力地去求多项式卷积的复杂度显然是 $\Theta(n^2)$ 的.
引自 算法导论
点值乘法
在每一个 $x$ 的取值 $x_0$ 处, 都有卷积后的多项式 $h(x_0)=f(x_0)g(x_0)$(显然对于一个具体的值直接乘起来就好了), 于是卷积可以变为
$$ H(x)=\{(x_0,f(x_0)g(x_0)),(x_1,f(x_1)g(x_1)),\dots,(x_{n-1},f(x_{n-1})g(x_{n-1}))\} $$
这个式子显然可以 $\Theta(n)$ 地算出来, 但是显然得到的多项式和要求的答案都是系数表示法, 所以还需要进行两种表示法之间的转换, 即 快速傅里叶变换 (FFT). 易知 $x$ 的取值对答案并没有影响, 于是现在要做的就是找到合适的 $x$ 简化表示法之间的转换.
DFT
首先是将多项式转换为点值表示法, 有
$$ F(x){=\sum_{i=0}^{n-1}f_ix^i=f_0x^0+f_1x^1+\dots+f_{n-1}^x{n-1}\\ =(f_0x^0+f_2x^2+\dots+f_{n-2}x^{n-2})+x(f_1x^0+f_3x^2+\dots+f_{n-2}x^{n-2}) } $$
从第一行到第二行的变换是奇偶项分开, 发现从后一个括号中提出一个 $x$ 后两个括号形式上已经很接近了, 接下来令 $t=x^2$, 则有
$$ F(x){=(f_0t^0+f_2t^1+\dots+f_{n-2}t^\frac{n-2}2)+x(f_1t^0+f_3t^1+\dots+f_{n-2}t^\frac{n-2}2)\\ =A(t)+xB(t)=A(x^2)+xB(x^2) } $$
其中,$a_0,a_1,\dots$ 依次为 $f_0,f_2,\dots$,$b_0,b_1,\dots$ 依次为 $f_1,f_3,\dots$. 设 $k<\frac n2$, 令 $x=\omega_n^k$, 由 $\texttt{(1)}$ 有第一个式子和第二个式子的第一步, 由 $\texttt{(2)}$, 并 $k+\frac 2n\equiv k\pmod{\frac n2}$, 有第二个式子的最后一步, 即
$$ f(\omega_n^k)=a(\omega_n^{2k})+\omega_n^kb(\omega_n^{2k})=a(\omega_{\frac n2}^k)+\omega_n^kb(\omega_{\frac n2}^k)\\ f(\omega_n^{k+\frac n2}){=a(\omega_n^{2k+n})+\omega_n^{k+\frac n2}b(\omega_n^{2k+n})\\ =a(\omega_{\frac n2}^{k+\frac n2})+\omega_n^{k+\frac n2}b(\omega_{\frac n2}^{k+\frac n2})=a(\omega_{\frac n2}^k)-\omega_n^kb(\omega_{\frac n2}^k) } $$
也就是说, 如果得到了 $a(\omega_{\frac n2}^k)$ 和 $b(\omega_{\frac n2}^k)$, 就可以求出 $F(x)$ 在 $\omega_n^k$ 和 $\omega_n^{k+\frac n2}$ 处的值 $f(\omega_n^k)$ 和 $f(\omega_n^{k+\frac n2})$. 这样就可以递归进行处理. 写出DFT()函数的伪代码
function DFT(*f,n) /*无返回值,依靠修改f[]实现影响结果*/
/*complex f[]表示当前F(x),以指针形式传递(即相当于引用),int n表示项数*/
if n=1 then
return /*如果只有常数项那么函数值就是常数项系数*/
end if
for i←0 to n/2 by 1 do
a[i]←f[i*2]
b[i]←f[i*2+1]
end for
call DFT(a,n/2)
call DFT(b,n/2)
/*在以上两次操作后a[]和b[]中存储的都是a(w^{1...n/2})和b(w^{1...n/2})的值*/
complex wx←complex(1,0) /*complex(real,imag)是复数模板类,wx保存w^x,初值为w^0*/
complex w0←complex(sin(pi/n),cos(pi/n)) /*w0保存w^1,方便每一次乘到wx上,pi为圆周率*/
for i←0 to n/2 by 1 do
f[i]←a[i]+wx*b[i] /*合并a[],b[]*/
f[i+n/2]←a[i]-wx*b[i]
wx←wx*w0 /*计算w^{x+1}*/
end for
return
end functionIDFT
结论是直接将DFT()中complex w0←complex(sin(pi/n),cos(pi/n))改为complex w0←complex(sin(pi/n),-cos(pi/n))( 即cos改为-cos), 然后末尾除一个 $n$ 即可. 至于具体证明过程窝也不是很清楚, 大致看着老师的 PPT 感性理解了一下.
如果将 DFT 写成矩阵乘法形式有
$$ y=V_nf\\ \begin{bmatrix} y_0\\y_1\\y_2\\y_3\\\vdots\\y_{n-1} \end{bmatrix}=\begin{bmatrix} 1&1&1&1&\cdots&1\\ 1&\omega_n&\omega_n^2&\omega_n^3&\cdots&\omega_n^{n-1}\\ 1&\omega_n^2&\omega_n^4&\omega_n^6&\cdots&\omega_n^{2n-2}\\ 1&\omega_n^3&\omega_n^6&\omega_n^9&\cdots&\omega_n^{3n-3}\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\ 1&\omega_n^{n-1}&\omega_n^{2n-2}&\omega_n^{3n-3}&\cdots&\omega_n^{n^2-2n+1}\\ \end{bmatrix}\begin{bmatrix} f_0\\f_1\\f_2\\f_3\\\vdots\\f_{n-1} \end{bmatrix} $$
其中 $V_n$ 是 $\omega_n$ 的范德蒙德矩阵(虽然并不清楚这是什么),$(i,j)$ 处的元素 $v_{i,j}=\omega_n^{ij}$,$1$ 的实质就是 $\omega_n^0$. 试着简单证明, 将矩阵乘法写开, 用 $x$ 换回 $\omega_n^i$ 有
$$ y_i{=f(\omega_n^i)=\sum_{j=0}^{n-1}f_j\omega_n^{ij}\\ =f(x)=\sum_{j=0}^{n-1}f_jx^j } $$
和 $F(x)$ 原本的形式一样, 得证.
那么容易想到
$$ f=V_n^{-1}y $$
即要求 $f_{0,1,\dots}$, 只需要用 $y$ 矩阵乘上 $V_n$ 的逆矩阵 $V_n^{-1}$ 满足 $V_nV_n^{-1}=I_n$ 即可.$I_n$ 表示一个 $n\times n$ 的单位矩阵, 满足 $i_{a,b}=[a=b]$($[x]$ 的值当表达式 $x$ 成立时为 $1$, 否则为 $0$, 即 $[x]=\begin{cases}{1,\quad x\text{ is true}\\0,\quad x\text{is false}}\end{cases}$), 换而言之, 只有对角线上的值为 $1$, 其余部分值为 $0$. 证明以上结论, 有
$$ y=V_nf\\ V_n^{-1}y=V_n^{-1}V_nf=I_nf=f\\ f_i=\sum_{j=0}^{n-1}i_jf_j=\sum_{j=0}^{n-1}[i=j]f_j $$
再 不加证明地 给出一个关于 $V_n^{-1}$ 中 $(i,j)$ 处元素 $v_{i,j}^{-1}$ 的结论
$$ v_{i,j}^{-1}=\frac{\omega_n^{-ij}}n $$
虽然不清楚这个结论是怎么被总结出来的, 但是可以试着推导一下
$$ [V_n^{-1}V_n]_{i,j}=\sum_{k=0}^{n-1}\frac{\omega_n^{kj}\omega_n^{-ki}}n=\sum_{k=0}^{n-1}\frac{\omega_n^{k(j-i)}}n $$
容易看出当 $i=j$ 时, 上式值为 $\sum_{k=0}^{n-1}\frac{\omega_n^0}n=\sum_{k=0}^{n-1}\frac1n=1$; 当 $i\not=j$ 时, 令 $t=j-i$, 有
$$ \sum_{k=0}^{n-1}\frac{\omega_n^{k(j-i)}}n=\frac1n\sum_{k=0}^{n-1}\omega_n^{k(j-i)}=\frac1n\cdot\frac{1-\omega_n^{nt}}{1-\omega_n^t}=0 $$
因为显然 $\omega_n^{nt}=\omega_n^0=1$. 由上可知 $V_n^{-1}V_n$ 确实符合 $I_n$ 的形式, 故得证.
观察发现
$$ v_{i,j}=\omega_n^{ij},\quad v_{i,j}^{-1}=\frac{\omega_n^{-ij}}n $$
两个式子只相差了一个符号和最后的除以 $n$, 所以得到本小节开头给出的结论
结论是直接将DFT中complex w0←complex(sin(pi/n),cos(pi/n))改为complex w0←complex(sin(pi/n),-cos(pi/n))( 即cos改为-cos), 然后末尾除一个 $n$ 即可.
将DFT()修改一下即可变为IDFT(), 将两函数合并写作FFT()
function FFT(*f,n,m) /*int m表示在进行DFT(1)还是IDFT(-1)*/
if n=1 then
return
end if
for i←0 to n/2 by 1 do
a[i]←f[i*2]
b[i]←f[i*2+1]
end for
call FFT(a,n/2,m)
call FFT(b,n/2,m)
complex wx←complex(1,0)
complex w0←complex(sin(pi/n),m*cos(pi/n))
/*上一行有所改动,如果是IDFT那么每次wx乘上的就应该是w^{-1}*/
for i←0 to n/2 by 1 do
f[i]←a[i]+wx*b[i]
f[i+n/2]←a[i]-wx*b[i]
wx←wx*w0
end for
return
end function优化
在FFT()中是递归计算的, 如果改成非递归速度就会快很多. 首先观察FFT()产生的递归树
引自 算法导论
其最终得到的序列与原序列对比有
$$ (000)_2=0,(001)_2=1,(010)_2=2,(011)_2=3,(100)_2=4,(101)_2=5,(110)_2=6,(111)_2=7\\(000)_2=0,(100)_2=4,(010)_2=2,(110)_2=6,(001)_2=1,(101)_2=5,(011)_2=3,(111)_2=7 $$
即将原序列二进制翻转即得到最终序列. 用 $r_i$ 表示 $i$ 翻转后得到的值, 有
for (int i=0;i<n;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(k-1));其中 $k=\lg n$.
Code
这样就可以写出最后的代码
inline void FFT(complex<double> *a,int n,bool m) { /*fast fast TLE*/
int k; complex<double> w0,wx,b,c;
static int r[5000005];
for (k=0;(1<<k)<n;++k);
if (!r[1]) for (int i=0;i<n;++i) /*求过一遍r[]再次调用时就不用再求了*/
r[i]=(r[i>>1]>>1)|((i&1)<<(k-1));
for (int i=0;i<n;++i)
if (i<r[i]) swap(a[i],a[r[i]]); /*如果不加i<r[i]会交换两次,相当于没交换*/
for (int i=1;i<n;i<<=1) { /*i表示当前正在合并的序列的长度*/
w0=complex<double>(cos(pi/i),m?-sin(pi/i):sin(pi/i));
for (int j=0;j<n;j+=(i<<1)) { /*用j来取遍相邻的每两个需要合并的序列*/
wx=complex<double>(1,0);
for (int k=0;k<i;++k) { /*只扫0...i/2-1的部分,得到i/2...i的答案*/
b=a[j+k]; c=wx*a[i+j+k];
a[j+k]=b+c; a[i+j+k]=b-c;
wx*=w0;
}
}
}
if (m) for (int i=0;i<n;++i) a[i]/=(double)n; /*如果是IDFT最后要记得除掉一个n*/
}一个小细节
注意以上的分治全部是基于每一次都完全平分成两份, 这就意味着必须有 $\lg n\in\mathbb Z$, 即 $n$ 恰为 $2$ 的整数次幂. 然而做题不可能这么巧, 所以需要在高位补 $0$ 来使 $n$ 满足要求.