快速傅里叶变换 FFT 学习笔记【多项式】【FFT】

作者: wjyyy 分类: FFT,分治,多项式,学习笔记,数学 发布时间: 2018-12-04 21:58

点击量:58

一、多项式

多项式的定义

    多项式是多个单项式的和。(这是人教版七年级上册的定义)形如\(x^2\),\(6y^5+5y^2+17\)这样的式子就是多项式。

多项式的乘积

    多项式乘积可以类比两个数相乘。对于两个式子中的任意两项,系数相乘,指数相加。

    比如\[\begin{align*}&(x^2+3x-5)(2x^2+6x+1)\\&=2x^2(x^2+3x-5)+6x(x^2+3x-5)+1(x^2+3x-5)\\&=2x^4+6x^3-10x^2+6x^3+18x^2-30x+x^2+3x-5\end{align*}\]

    合并同类项后,原式\(=2x^4+18x^3+9x^2-27x-5\)。

    这一过程可以用竖式来清晰地表示,即\[\begin{align*}
x^2+3x-5 \\
2x^2+6x+1 \\
\overline{\qquad \qquad \qquad \quad x^2+3x-5}\\
6x^3+18x^2-30x\ \quad\\
2x^4+6x^3-10x^2\qquad\qquad\\
\overline{2x^4+12x^3+9x^2-27x-5}
\end{align*}\]

    如果两个多项式的次数(即最高次项的指数)分别为\(n\)和\(m\),那么多项式相乘的时间复杂度就是\(O(nm)\),在上面的竖式中就表现为\(n\)行\(m+n\)列。我们需要一种途径来降低这一过程的时间复杂度。

二、傅里叶变换

    傅里叶变换实际上是一种数学转换,它可以描述一个信号的特征,并转化为三角函数的形式。不过在算法竞赛中我们没有必要掌握它的原理,只需要利用它进行多项式乘法计算。

三、函数的求值与插值

    在二次函数中,有\(f(x)=ax^2+bx+c\),因为其有三个未知数,且函数图像连续,因此我们一旦知道三个图像上的点,就可以联立出三元一次方程组来解它。也就是三个平面直角坐标系上不同的点唯一确定一个二次函数。以此类推,\(n\)次函数为\(f(x)=a_nx^n+a_{n-1}x^{n-1}+\dots a_0x^0\),它有\(n+1\)个未知数,那么我们需要\(n+1\)个点来确定这个函数。

    如果我们能够找出\(n+1\)个点,列出\(n+1\)个方程,利用高斯消元来解出这个方程,时间复杂度是\(O(n^3)\),我们可以利用这一点,并对它进行优化来处理多项式有关的问题。

    对于两个多项式,我们把它们分别视作函数\(g(x)\)和\(h(x)\),他们的积为\(f(x)\)。那么我们在上面取\(n+1\)个点\(\{x_0,x_1,\dots x_n\}\),把它们分别代到\(h(x)\)和\(g(x)\)中去。那么\(f(x_i)=g(x_i)h(x_i)\),得到了一系列点\(\{(x_0,f(x_0)),(x_1,f(x_1)),\dots ,(x_n,f(x_n))\}\),这样就可以唯一确定\(f(x)\)了。

    上面通过\(\{x_0,x_1,\dots x_n\}\)得到\(\{(x_0,f(x_0)),(x_1,f(x_1)),\dots,(x_n,f(x_n))\}\)的过程叫做求值,通过\(\{(x_0,f(x_0)),(x_1,f(x_1)),\dots,(x_n,f(x_n))\}\)解出\(f(x)\)解析式的过程叫做插值

四、复数和单位根

复数

    复数就是既有实部又有虚部的数,形如\(a+b\mathrm{i}\),例如\(3+5\mathrm{i}\),\(9-\mathrm{i}\),就是复数。

    我们可以把复数理解为一个向量,在平面直角坐标系中,横轴代表实部,纵轴代表虚部。任何一个起点为原点的向量都代表着一个复数,和向量一样,复数也可以做加减乘运算,复数还可以做除法。


\[x=a+b\mathrm{i},y=c+d\mathrm{i}\]
那么
\[x+y=a+c+(b+d)\mathrm{i}\\
x-y=a-c+(b-d)\mathrm{i}\\
xy=ac-bd+(ad+bc)\mathrm{i}\]
对于除法,我们可以用类似分母有理化的方法来进行化简
\[\frac xy=\frac{a+b\mathrm{i}}{c+d\mathrm{i}}=\frac{(a+b\mathrm{i})(c-d\mathrm{i})}{(c+d\mathrm{i})(c-d\mathrm{i})}=\frac{ac+bd+(bc-ad)\mathrm{i}}{c+d}\]
    不过这些在下面的运用中不是很大,只要把复数看作一个向量就可以了。

单位根

    既然复数可以被看作是向量,那么它就有坐标,对于虚数\(x=a+b\mathrm{i}\),它的横坐标是\(a\),纵坐标是\(b\)。平面直角坐标系中有单位圆,我们可以研究单位圆上复数的性质。

    设\(\omega_n^k\)表示把单位圆分成\(n\)份,取第\(k\)个向量(复数\(1\)(坐标为\((1,0)\))认为是第\(0\)个向量)。或者说,用坐标表示有\(\omega_n^k=\left(\cos\frac{2\pi}{n}k,\sin\frac{2\pi}{n}k\right)\);用数字表示有\(\omega_n^k=\cos\frac{2\pi}{n}k+\mathrm{i}\sin\frac{2\pi}{n}k\)

    单位根有一些性质,比如\(\omega\)右上角的数可以看作指数,相乘时指数相加:\(\omega_n^{k_1}\cdot\omega_n^{k_2}=\omega_n^{k_1+k_2}\),同理幂运算也和数字类似\(\left(\omega_n^p\right)^q=\omega_n^{pq}\)。也因为是在单位圆上,而且也是由三角函数表示的,单位根存在周期,\(\omega_n^k=\omega_n^{k\pmod n}\)。

    还有两个比较重要的性质是\(\omega_{2n}^{2k}=\omega_n^k,\omega_n^k=-\omega_n^{k+\frac2n}\),都可以用三角函数来理解。我们之后要用到的也主要是这两个。

五、离散傅里叶变换

离散傅里叶变换 – DFT

    如果我们把一个\(n-1\)次多项式带入\(n\)个值,用点值来表示这个多项式是可行的。那么这\(n\)个数我们就分别带入\(\omega_n^0,\omega_n^1\dots,\omega_n^{n-1}\)。得到\(f(\omega_n^i)=a_n\omega_n^{in}+a_{n-1}\omega_n^{i(n-1)}+\dots +a_0\)。

    这样用单位根赋值,用点值表示多项式的方法就是离散傅里叶变换。这样做是为了给下面的逆变换做准备和提供方便。

离散傅里叶逆变换 – IDFT

    如果我们已经求出了一个\(n-1\)次多项式的一组点值表示\(\left\{(\omega_n^0,y_0),(\omega_n^1,y_1),\dots,(\omega_n^{n-1},y_{n-1})\right\}\),简写为\(\left\{y_0,y_1,\dots,y_{n-1}\right\}\)。此时要求原多项式的系数\(\left\{a_0,a_1,\dots,a_{n-1}\right\}\),我们可以尝试用逆运算的方法来解决这个问题。

    多项式运算实际上可以和矩阵结合起来,\(\omega_n^k\)是\(x\)的值,\(a_k\)是系数,这组点值\(\left\{y_0,y_1,\dots,y_{n-1}\right\}\)可以认为是矩阵的计算结果,如下。(第一行中的\(1\)实则为\(\left(\omega_n^0\right)^k\))
\[\begin{bmatrix}
y_0\\y_1\\y_2\\\vdots\\y_{n-1}
\end{bmatrix}
=
\begin{bmatrix}
1&1&1&\cdots&1\\
1&\omega_n^1&\omega_n^2&\cdots&\omega_n^{n-1}\\
1&\omega_n^2&\omega_n^4&\cdots&\omega_n^{2(n-1)}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
1&\omega_n^{n-1}&\omega_n^{2(n-1)}&\cdots&\omega_n^{(n-1)^2}
\end{bmatrix}
\begin{bmatrix}
a_0\\a_1\\a_2\\\vdots\\a_n
\end{bmatrix}\]

    此时我们已知上述三个矩阵(分别命名为\(a,b,c\))中的前两个,需要构造出一个矩阵可以代表矩阵\(c\)。既然有了这个等式,我们把它看作一个方程,只要左右两边都把矩阵\(b\)给消掉,等式左边就是矩阵\(c\)了。

    因此我们要构造出矩阵\(b\)的逆矩阵\(d\),即\(b\times d=o\),此时\(o\)代表单位矩阵\(\begin{bmatrix}
1&0&0&\cdots&0\\
0&1&0&\cdots&0\\
0&0&1&\cdots&0\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
0&0&0&\cdots&1
\end{bmatrix}\),即主对角线上都是\(1\),其余位置都是\(0\)的矩阵。此时单独考虑矩阵\(b\),看怎样把它变成矩阵\(o\)。

事实上,矩阵\(b\)的逆矩阵是非常好找的。我们把矩阵\(b\)的每个位置都变成其倒数,得到矩阵\(b’\)。
\[b’=\begin{bmatrix}
1&1&1&\cdots&1\\
1&\omega_n^{-1}&\omega_n^{-2}&\cdots&\omega_n^{-(n-1)}\\
1&\omega_n^{-2}&\omega_n^{-4}&\cdots&\omega_n^{-2(n-1)}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
1&\omega_n^{-(n-1)}&\omega_n^{-2(n-1)}&\cdots&\omega_n^{-(n-1)^2}
\end{bmatrix}\]
此时观察\(b\times b’\)的运算结果,得到矩阵\(f\):
\[f=b\times b’=\begin{bmatrix}
1&1&1&\cdots&1\\
1&\omega_n^1&\omega_n^2&\cdots&\omega_n^{n-1}\\
1&\omega_n^2&\omega_n^4&\cdots&\omega_n^{2(n-1)}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
1&\omega_n^{n-1}&\omega_n^{2(n-1)}&\cdots&\omega_n^{(n-1)^2}
\end{bmatrix}
\begin{bmatrix}
1&1&1&\cdots&1\\
1&\omega_n^{-1}&\omega_n^{-2}&\cdots&\omega_n^{-(n-1)}\\
1&\omega_n^{-2}&\omega_n^{-4}&\cdots&\omega_n^{-2(n-1)}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
1&\omega_n^{-(n-1)}&\omega_n^{-2(n-1)}&\cdots&\omega_n^{-(n-1)^2}
\end{bmatrix}\]

那么对于\(f\)的每一个位置,\(f_{i,j}=\sum_{k=0}^{n-1}b_{i,k}b’_{k,j}\),

通过式子我们可以看出\(b_{i,k}=\left(\omega_n^i\right)^k,b_{k,j}=\left(\omega_n^{-k}\right)^j\),

此时\(f_{i,j}=\sum_{k=0}^{n-1}\omega_n^{k(i-j)}\),

可以看出这是一个首项为\(\omega_n^0=1\),公比为\(\omega_n^{i-j}\),项数为\(n\)的等差数列。

当\(i=j\)时,公比为\(1\),前\(n\)项和为\(S_n=n\);

当\(i\ne j\)时\(S_n=\frac{1-\omega_n^{n(i-j)}}{1-\omega_n^{i-j}}\),在分子上\(\omega_n^{n(i-j)}=\omega_n^0=1\),所以\(S_n=0\)。

所以当且仅当\(i=j\),有\(f_{i,j}=n\),得到矩阵\(d’​\)如下:
\[\begin{bmatrix}
n&0&0&\cdots&0\\
0&n&0&\cdots&0\\
0&0&n&\cdots&0\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
0&0&0&\cdots&n
\end{bmatrix}\]

这时原等式为
\[\begin{bmatrix}
1&1&1&\cdots&1\\
1&\omega_n^{-1}&\omega_n^{-2}&\cdots&\omega_n^{-(n-1)}\\
1&\omega_n^{-2}&\omega_n^{-4}&\cdots&\omega_n^{-2(n-1)}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
1&\omega_n^{-(n-1)}&\omega_n^{-2(n-1)}&\cdots&\omega_n^{-(n-1)^2}
\end{bmatrix}
\begin{bmatrix}
y_0\\y_1\\y_2\\\vdots\\y_{n-1}
\end{bmatrix}
=
\begin{bmatrix}
n&0&0&\cdots&0\\
0&n&0&\cdots&0\\
0&0&n&\cdots&0\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
0&0&0&\cdots&n
\end{bmatrix}
\begin{bmatrix}
a_0\\a_1\\a_2\\\vdots\\a_n
\end{bmatrix}\]

右边为
\[\begin{bmatrix}
na_0&0&0&\cdots&0\\
0&na_1&0&\cdots&0\\
0&0&na_2&\cdots&0\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
0&0&0&\cdots&na_{n-1}
\end{bmatrix}\]

此时将得到的多项式每一项除以\(n\),得到的就是原多项式的系数了。

六、快速傅里叶变换 – FFT

操作:

上面的计算如果全部用多项式进行,复杂度当然是\(O(n^2)\)的,当然矩阵运算可以达到O(n³)。然而FFT可以把这一过程优化到\(O(n\log n)\)。

对于\(f(x)=a_0+a_1x+a_2x^2+a_3x^3+\dots+a_{n-1}x^{n-1}​\),考虑分治。

它可以被化为\(f(x)=(a_0+a_2x^2+\dots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+\dots+a_{n-1}x^{n-2})\)。

可以发现,在奇数项统一提出一个\(x\)后,两个括号中的指数变得很相似。我们把\(x=\omega_n^k\)带入式子得
\[\begin{align*}
f(\omega_n^k)&=\left(a_0+a_2\left(\omega_n^k\right)^2+\dots+a_{n-2}\left(\omega_n^k\right)^{n-2}\right)+\omega_n^k\left(a_1+a_3\left(\omega_n^k\right)^2+\dots+a_{n-1}\left(\omega_n^k\right)^{n-2}\right)\\
&=\left(a_0+a_2\omega_n^{2k}+\dots+a_{n-2}\omega_n^{k\left(n-2\right)}\right)+\omega_n^k\left(a_1+a_3\omega_n^{2k}+\dots+a_{n-1}\omega_n^{k\left(n-2\right)}\right)
\end{align*}\]

此时所有的指数都是偶数,那么根据上面单位根中讲到的性质,\(\omega_{2n}^{2k}=\omega_n^k\)、\(\omega_n^k=-\omega_n^{k+\frac n2}\)。

则当\(k<\frac n2\)时,定义两个新函数\(f_1,f_2\)分别表示两个系数不同的子多项式
\[\begin{align*}
f(\omega_n^k)&=\left(a_0+a_2\omega_{\frac n2}^k+\dots+a_{n-2}\omega_{\frac n2}^{\frac{n-2}2k}\right)+\omega_n^k\left(a_1+a_3\omega_{\frac n2}^k+\dots+a_{n-1}\omega_{\frac n2}^{\frac{n-2}2k}\right)\\
&=f_1(\omega_\frac n2^k)+\omega_n^kf_2(\omega_\frac n2^k)
\end{align*}\]

此时如果计\(f(\omega_n^{k+\frac n2}\),由上面两个性质和单位根的周期性\(\omega_n^k=\omega_n^{k+n}\)可得
\[\begin{align*}
f(\omega_n^{k+\frac n2})&=\left(a_0+a_2\omega_{\frac n2}^{k+\frac n2}+\dots+a_{n-2}\omega_{\frac n2}^{\frac{n-2}2\left(k+\frac n2\right)}\right)+\omega_n^{k+\frac n2}\left(a_1+a_3\omega_{\frac n2}^{k+\frac n2}+\dots+a_{n-1}\omega_{\frac n2}^{\frac{n-2}2\left(k+\frac n2\right)}\right)\\
&=\left(a_0+a_2\omega_{\frac n2}^k+\dots+a_{n-2}\omega_{\frac n2}^{\frac{n-2}2k}\right)-\omega_n^k\left(a_1+a_3\omega_{\frac n2}^k+\dots+a_{n-1}\omega_{\frac n2}^{\frac{n-2}2k}\right)\\
&=f_1(\omega_\frac n2^k)-\omega_n^kf_2(\omega_\frac n2^k)\\
\end{align*}\]

对比后看出两式子化简后只有后一项的系数不同。那么当我们要求\(f(\omega_n^0),f(\omega_n^1),\dots,f(\omega_n^{n-1})\)时,我们只需要求出\(f(\omega_n^i),i<\frac n2\),那么\(i\ge \frac n2\)的部分就直接求出来了。

不过因为我们分治过程中要保证每一次的\(n\)均为偶数,所以一开始的\(n\)一定要是\(2\)的整数幂。我们对于两个项数分别为\(i,j\)的多项式,需要开\(2^{\left\lceil\log_2(i+j+1)\right\rceil}\)个位置才能顺利地进行分治。

    一次FFT的过程如下图所示:黄色箭头表示右边相应位置的值可以通过左边来计算;橙色箭头表示\(f_1,f_2\)的来源。

    那么我们分别对两个多项式进行FFT操作,之后把点值的纵坐标相乘,再通过FFT的过程代入\(\omega\)的倒数,即可求得两式相乘后的系数。

Code:

const double pi=acos(-1);
struct cpl//复数类
{
    double x,y;
    cpl(double x,double y)
    {
        this->x=x;
        this->y=y;
    }
    cpl()
    {
        x=0,y=0;
    }
    friend cpl operator +(cpl a,cpl b)
    {
        return cpl(a.x+b.x,a.y+b.y);
    }
    friend cpl operator -(cpl a,cpl b)
    {
        return cpl(a.x-b.x,a.y-b.y);
    }
    friend cpl operator *(cpl a,cpl b)
    {
        return cpl(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
    }
}a[1<<21],b[1<<21];
cpl omega(int n,int k)//生成复数ω
{
    double ang=2*pi*k/n;
    return cpl(cos(ang),sin(ang));
}
void fft(cpl *f,int len,int flag)//flag=1表示DFT,-1表示IDFT
{
    if(len==1)
        return;
    int tot=(len>>1);
    cpl a1[len>>1],a2[len>>1];
    for(int i=0;i<tot;++i)
        a1[i]=f[i<<1],a2[i]=f[i<<1|1];
    fft(a1,tot,flag);
    fft(a2,tot,flag);
    for(int i=0;i<tot;++i)
    {
        cpl omg=omega(len,i*flag);
        f[i]=a1[i]+omg*a2[i];
        f[i+tot]=a1[i]-omg*a2[i];
    }
    return;
}

七、优化

1.非递归写法

    递归耗时比较大,可以改为非递归。在每个递归过程中,处于奇数位置的数和偶数位置的数分别要归于自己的位置,在上面的代码中,直到进入最后一层前都一直在进行交换操作。因此递归的作用就是进行有层次的交换。

    有人总结出来一个规律,一个数一开始如果位于\(p\),那么它最终就会位于\(p\)在二进制下翻转一次后表示的数。例如在长为\(8\)的变换中,二进制长度为\(3\)。\(110_{(2)}=6\)变换后成为了\(011_{(2)}=3\),事实证明确实如此。我们实际上可以\(O(n)\)递推预处理出每个位置变换后的位置,pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1))。然后直接循环处理,理解起来比较简单,实现时要注意细节。

2.蝴蝶优化

    这是一个常数上的优化,在非递归写法中,如果直接转移,可能会带来后效性的影响,一般会额外开一个数组进行copy。有了蝴蝶操作,我们就要先更新\(f\left[k+\frac n2\right]\),再更新\(f[k]\),这样在本次操作中没有紊乱,在之后的操作中也不会调用自己,所以不需要用一个额外的数组了。

    事实上我们把\(f_1\)和\(\omega_n^kf_2\)分别拿临时变量存一下也能达到同样的效果。

3.输出优化

    输出时使用"%.0lf"会比"%d"慢很多,最好先转为整型。

八、全代码

#include<cstdio>
#include<cstring>
#include<cmath>
const double pi=acos(-1);
struct cpl
{
    double x,y;
    cpl(double x,double y)
    {
        this->x=x;
        this->y=y;
    }
    cpl()
    {
        x=0,y=0;
    }
    cpl friend operator +(cpl a,cpl b)
    {
        return cpl(a.x+b.x,a.y+b.y);
    }
    cpl friend operator -(cpl a,cpl b)
    {
        return cpl(a.x-b.x,a.y-b.y);
    }
    cpl friend operator *(cpl a,cpl b)
    {
        return cpl(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
    }
}a[1<<21],b[1<<21];
cpl omega(int n,int k)
{
    double ang=2*pi*k/n;
    return cpl(cos(ang),sin(ang));
}
cpl f[1<<21],t1,t2;
int pos[1<<21],cnt;
void fft(cpl *f,int len,int flag)
{
    int tmp=(1<<len);
    for(int i=0;i<tmp;++i)
        if(i<pos[i])
        {
            cpl t=f[i];
            f[i]=f[pos[i]];
            f[pos[i]]=t;
        }
    for(int bs=1;bs<=len;++bs)
    {
        int L=(1<<bs);
        int g=(L>>1);
        cpl step=omega(L,flag),ori=omega(L,0);
        for(int i=0;i<tmp;i+=L)
        {
            cpl omg=ori;
            for(int j=0;j<g;++j,omg=omg*step)
            {
                t1=f[i+j],t2=omg*f[i+j+g];
                f[i+j]=t1+t2;
                f[i+j+g]=t1-t2;
            }
        }
    }
    return;
}
int main()
{
    int n,m,u;
    scanf("%d%d",&n,&m);
    int len=1;
    while((1<<len)<=m+n)
        ++len;
    for(int i=1;i<(1<<len);++i)
        pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
    for(int i=0;i<=n;++i)
    {
        scanf("%d",&u);
        a[i]=cpl(u,0);
    }
    for(int i=0;i<=m;++i)
    {
        scanf("%d",&u);
        b[i]=cpl(u,0);
    }
    fft(a,len,1);
    fft(b,len,1);
    for(int i=0;i<(1<<len);++i)
        a[i]=a[i]*b[i];
    fft(a,len,-1);
    for(int i=0;i<=n+m;++i)
        printf("%d ",int(a[i].x/(1<<len)+.5));
    return 0;
}

九、例题

1.A*B Problem

2.[ZJOI2014]力

2
说点什么

avatar
1 Comment threads
1 Thread replies
0 Followers
 
Most reacted comment
Hottest comment thread
2 Comment authors
wjyyy脑瘫OIer cby Recent comment authors
  Subscribe  
最新 最旧 得票最多
提醒
脑瘫OIer cby
游客
脑瘫OIer cby

%%%

/* */