快速傅里叶变换(FFT)与多项式算法学习笔记

参考资料:menci的博客

前言:

最近在学习生成函数,无奈的发现如果我不学习\(O(nlogn)\)的多项式算法的话什么题也做不了qwq

于是就滚来学习FFT了

其实写的很烂,主要是给自己看的

好像整个机房就我不会这玩意了

定义

多项式

形如\(F(x)=\sum\limits_{i=0}^na_ix^i\)的柿子就是一个多项式,这个多项式的次数就是它的最高次数\(n\)

多项式的表示方法

系数表示法

就是用\(\{a_1,a_2,a_i,...,a_n\}\)来表示这个多项式．

点值表示法

就是用n个点\((x_i,y_i)\)来表示这个多项式．

对于任意一个点\(F(x_i)=y_i\)

易知这样能够唯一确定一个多项式．

点值表示法转换成系数表示法可以使用插值.

多项式乘法

定义多项式\(A(x)=\sum\limits_{i=0}^na_ix^i\)与多项式\(B(x)=\sum\limits_{i=0}^nb_ix^i\)的乘积为

\(C(x)=\sum\limits_{k=0}^{2n}(\sum\limits_{i+j=k}{a_ib_j})x^k\)

不足位向高位补零即可．

时间复杂度\(O(n^2)\)

如果使用点值表示法直接把对应的\(y_i\)乘起来就可以了．

时间复杂度\(O(n)\)

那么看到这里你就发现了一种非常优秀的算法，就是使用点值表示法大力相乘，复杂度\(O(n)\)，本篇博客完．

然而我们平常所用到的多项式一般是系数形式的．

非常不幸的告诉你，

将点值/系数表示转化成系数/点值表示时间复杂度是\(O(n^2)\)的．

于是我们就想,有没有一个优秀的算法能够把转化的时间复杂度降下来呢?

当然有了

没错,它就是快速傅里叶变换!

前置知识

复数

令\(i^2=-1\),形如\(a+bi\)的数被称为复数,\(i\)就是虚数单位。

复平面

复平面上的x轴代表实数,y轴代表虚数。

每个复数\(a+bi\)都是复平面上的一个从\((0,0)\)指向\((a,b)\)的向量。

模长为\(\sqrt{a^2+b^2}\),幅角就是从x轴正半轴逆时针转过的有向角度为幅角。

相加遵循平行四边形定则。

相乘时,模长相乘,幅角相加。

单位根

数学上，n次单位根是n次幂为1的复数。它们位于复平面的单位圆上，构成正n边形的顶点，其中一个顶点是1。

以上是百度百科

\(menci\)大佬的解释

在复平面上，以原点为圆心， 为1半径作圆，所得的圆叫做单位圆。以原点为起点，单位圆的\(n\)等分点为终点，作n个向量。

设所得的幅角为正且最小的向量对应的复数为\(\omega_n\)，称为\(n\)次单位根。

由复数乘法的定义（模长相乘，幅角相加）可知，其与的\(n-1\)个向量对应的复数分别为\(\omega^2_n ,\omega^3_n ...\omega^n_n\)，其中 \(\omega^n_n=\omega^0_n=1\)。

欧拉公式

\(\large e^{\pi i}=-1,e^{2\pi i}=1\)

所以\(\large \omega_n= e^{\frac {2\pi i}n},\omega_n^k=(e^\frac{2\pi i}n)^k=\omega_n^k=\cos\frac{2\pi k}n+i\sin\frac{2\pi k}n\)(由向量运算法则可以得到)

单位根的性质

消去引理:\(\omega_{dn}^{dk}=\omega_n^k,k\in N,d,n\in N^+\)

显然成立.

折半引理:若n是偶数,则\(\omega_n^{k+\frac n2}=-\omega_n^k\)

由欧拉公式\(\omega_n^{k+\frac n2 }=(e^\frac{2\pi i}n)^{k+\frac n2}=-(e^{\frac{2\pi i}n})\)

求和引理:\(\sum\limits_{j=0}^{n-1}(\omega_n^k)^j=0,n,k\in N^+,n\nmid k\)

证明:\(\Large \sum\limits_{j=0}^{n-1}(\omega_n^k)^j=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}=0\)

就是等比数列求和公式。

快速傅里叶变换

考虑多项式\(A(x)\)的表示。将\(n\)次单位根的\(0\)到\(n−1\)次幂带入多项式的系数表示,所得点值向量\((A(\omega_n^0),A(\omega_n^1),\ldots,A(\omega_n^{n-1}))\)称为其系数向量\((a_0,a_1,\ldots,a_{n−1})\)的离散傅里叶变换。

利用朴素算法,时间复杂度为\(O(n^2)\)

将多项式按照系数下标的奇偶分为两部分：

\(A(x)=(a_0+a_2x^2+a_4x^4+\cdots+a_{n-2}x^{n-2})+(a_1x+a_3x^3+a_5x^5+\cdots+a_{n-1}x^{n-1})\)

令
\[ \begin{aligned} A_1(x)&=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac n2-1} \\ A_2(x)&=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\frac n2-1} \end{aligned}\\ \]
则
\[ A(x)=A_1(x^2)+xA_2(x^2) \]
假设\(k<\frac n2\)
\[ A(\omega_n^k)=A1(\omega_n^{2k})+\omega_ n^{k}A2(\omega_n^{2k})\\ A(\omega_n^k)=A1(\omega_{\frac n2}^k)+\omega_n^kA2(\omega_{\frac n2}^k) \]
对于\(\omega_n^{k+\frac n2}\)
\[ A(\omega_n^{k+\frac n2})=A1(\omega_n^{2k+n})+\omega_{n}^{k+\frac n2}A2(\omega_{n}^{2k+n})\\=A1(\omega_{\frac n2}^k)-\omega_n^kA2(\omega_{\frac n2}^k) \]
这样,如果我们知道\(A1,A2\)在\(\omega_{\frac n2}^{0\text{~}\frac n2-1}\)的所有取值,我们就能对于所有的\(k\in [0,n)\)都算出来\(A\)在\(\omega_n^{0\text ~ n-1}\)的所有取值了。

于是我们就可以递归的去求解,记得高位补零。

傅里叶逆变换

将点值表示的多项式转化为系数表示，同样可以使用快速傅里叶变换，这个过程称为傅里叶逆变换。
设\((y_0,y_1,\ldots,y_{n-1})\)为\((a_0,a_1,\dots,a_{n-1})\)的傅里叶变换。

考虑另一个向量\((c_0,c_1,\dots,c_{n-1})\)满足
\[ c_k=\sum\limits_{i=0}^{n-1}y_i(\omega_{n}^{-k})^i \]
展开即得
\[ c_k=\sum\limits_{i=0}^{k-1}(\sum\limits_{j=0}^{k-1}a_j(\omega_{n}^i)^j)(\omega_n^{-k})^i\\ c_k=\sum\limits_{i=0}^{k-1}(\sum\limits_{j=0}^{k-1}a_j(\omega_{n}^j)^i)(\omega_n^{-k})^i\\ c_k=\sum\limits_{i=0}^{k-1}\sum\limits_{j=0}^{k-1}a_j(\omega_{n}^{j-k})^i\\ c_k=\sum\limits_{i=0}^{k-1}a_j(\sum\limits_{j=0}^{k-1}(\omega_{n}^{j-k})^i) \]
根据求和引理,我们知道\(\sum\limits_{j=0}^{k-1}(\omega_n^{j-k})\)当\(j-k\neq 0\)时总为0,当\(j=k\)时为\(n\)。

\(\therefore c_k=na_k\)

\(\therefore a_k=\frac 1n c_k\)

所以，使用单位根的倒数代替单位根，做一次类似快速傅里叶变换的过程，再将结果每个数除以\(n\)，即为傅里叶逆变换的结果。

实现

一般用C++ 自带的complex来实现复数,当然也可以手写。

考虑到单位根的倒数等于其共轭复数，在程序实现中，为了减小误差，通常使用 std::conj() 取得 IDFT 所需的「单位根的倒数」。

优化

有蝴蝶操作和迭代.

这里并不想讲,可以看一下menci聚聚的博客.

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define IL inline
#define RG register
#define gi geti<int>()
#define gl geti<ll>()
#define gc getchar()
#define File(a) freopen(a".in","r",stdin);freopen(a".out","w",stdout)
template<typename T>IL bool chkmax(T &x,const T &y){return x<y?x=y,1:0;}
template<typename T>IL bool chkmin(T &x,const T &y){return x>y?x=y,1:0;}
template<typename T>
IL T geti()
{
    RG T xi=0;
    RG char ch=gc;
    bool f=0;
    while(!isdigit(ch))ch=='-'?f=1:f,ch=gc;
    while(isdigit(ch))xi=xi*10+ch-48,ch=gc;
    return f?-xi:xi;
}
template<typename T>
IL void pi(T k,char ch=0)
{
    if(k<0)k=-k,putchar('-');
    if(k>=10)pi(k/10);
    putchar(k%10+'0');
    if(ch)putchar(ch);
}
const int N=4e6+7;
typedef double db;
class _complex{
    public:
    db x,y;
    _complex(){}
    _complex(db _x,db _y):x(_x),y(_y){}
};
_complex operator + (const _complex&a,const _complex&b){
    return _complex(a.x+b.x,a.y+b.y);
}
_complex operator - (const _complex&a,const _complex&b){
    return _complex(a.x-b.x,a.y-b.y);
}
_complex operator * (const _complex&a,const _complex&b){
    return _complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
_complex& operator *= (_complex &a,const _complex&b)
{
    return a=a*b;
}
const db PI=acos(-1);
int R[N],L,n,m;
inline void FFT(_complex *a,int opt)
{
    for(int i=0;i<n;++i)if(i<R[i])swap(a[i],a[R[i]]);
    for(int j=1;j<n;j<<=1)
    {
        _complex O(cos(PI/j),sin(PI/j)*opt);
        for(int k=0;k<n;k+=(j<<1)){
            _complex o(1,0);
            for(int l=0;l<j;++l,o*=O)
            {
                _complex Nx=a[k+l],Ny=o*a[j+k+l];
                a[k+l]=Nx+Ny;
                a[j+k+l]=Nx-Ny;
            }
        }
    }
}
_complex f[N],g[N];
int main(void)
{
    n=gi,m=gi;
    for(int i=0;i<=n;++i)f[i].x=gi;
    for(int i=0;i<=m;++i)g[i].x=gi;
    for(m+=n,n=1;n<=m;n<<=1,++L);
    for(int i=0;i<=n;++i)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    FFT(f,1),FFT(g,1);
    for(int i=0;i<=n;++i)f[i]*=g[i];
    FFT(f,-1);
    for(int i=0;i<=m;++i)printf("%d ",(int)(f[i].x/n+0.5));
    return 0;
}

转载于:https://www.cnblogs.com/LLCSBlog/p/11622835.html