Loading...

浅析FFT与NTT

check评论:0 条 remove_red_eye浏览量:385 change_historyTags:编程学习笔记
作者 : deco date_range日期 : 2018-08-07

$FFT$是$DFT$(离散傅里叶变换)的快速算法,将其时间复杂度从$O(n^2)$优化到了$O(n logn)$

前置技能:

  • 系数表示法

多项式$A(x)=\sum_{i=0}^n a_i \times x^i$,这种算法的时间复杂度为$O(n^2)$,但是不好优化,我们需要其他的表示方法

  • 点值表示法

将$n$互不相同的$x$带入多项式,会得到$n$个不同的取值$y$

则该多项式被这$n$个点$(x_1,y_1),(x_2,y_2),\dots,(x_n,y_n)$唯一确定

其中$y_i=\sum_{j=0}^{n-1} a_j\times x_i^j$

利用这种方法计算多项式乘法的时间复杂度仍然为$O(n^2)$,显然的,我们需要优化

如何优化点值表示法

前置技能:

  • 向量:指具有大小和方向的量
  • 复数:目前最大的数域,包含了实数和虚数
  • 向量的复数表示:一个向量的实部和虚部分别对应了复数系点中的横坐标和纵坐标

$Example:\vec a=a+bi$,则它的坐标表示为$(a,b)$

  • 向量的运算
    PsivD0.png

$\overrightarrow {AC}+\overrightarrow {CB}=\overrightarrow {AB}$,形象的理解起来,就是你从$A$走到$C$,再从$C$走到$B$,等于你从$A$走到$B$

$\vec a \cdot \vec b=|a| \times |b| \times cos<\vec a,\vec b>$,$|a|$和$|b|$代表向量的模

  • 单位根

在复平面上,以原点为圆心,$1$为半径作圆,所得的圆叫单位圆。以圆点为起点,圆的$n$等分点为终点,做$n$个向量,设幅角为正且最小的向量对应的复数为$\omega_n$,称为$n$次单位根。

根据复数乘法的运算法则,其余$n-1$个复数为$\omega_n^2,\omega_n^3,\ldots,\omega_n^n$
注意$\omega_n^0=\omega_n^n=1$(对应复平面上以$x$轴为正方向的向量)

那么如何计算它们的值呢?这个问题可以由欧拉公式解决$$\omega_{n}^{k}=\cos\ k \times \frac{2\pi}{n}+i\sin k \times \frac{2\pi}{n}$$

  • 简略证明欧拉公式(这部分无关紧要,可跳过或略看)

$e^x=1+\sum_{i=1}^∞ \frac{x^i}{i!}$

$cos x=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}......$

$sin x=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}......$

若把$e^x$中的$x$转换为$±ix$,$e^{\pm ix}=1 \pm \frac{ix}{1!}-\frac{x^2}{2!} \mp \frac{ix^3}{3!}......$

$e^{\pm ix}=cos x \pm isinx$

则$e^{ix}=cos x + isinx$,证毕

  • 辐角

向量$\vec a$与$x$轴正半轴的夹角称为向量$\vec a$的辐角,记作$Arg\ z$

  • 辐角主值

$Arg\ z$中只有一个$\theta_0$满足$-\pi < \theta_0 \leq \pi$,称作其辐角主值,记作$arg\ z$

  • 单位根性质

$\omega _{n}^{k}=\cos k\frac{2\pi}{n}+i\sin k\frac {2\pi }{n}$(即上面的公式)

$\omega _{2n}^{2k}=\cos 2k\times \frac{2\pi}{2n}+i\sin2k\times \frac{2\pi}{2n}$
$=\omega _{n}^{k}$

$\omega _{n}^{k+\frac{n}{2}}=-\omega _{n}^{k}$

$\omega _{n}^{\frac{n}{2}}=\cos\frac{n}{2}\times \frac{2\pi}{n}+i\sin\frac{n}{2}\times \frac{2\pi}{n}=\cos \pi+i\sin\pi=-1$

$\omega _{n}^{0}=\omega _{n}^{n}=1$

好了,前置技能就这么多了,下面步入正题

$FFT$(快速傅里叶变换)

我们前面提到过,一个$n$次多项式可以被$n$个点唯一确定。
那么我们可以把单位根的$0$到$n-1$次幂带入,这样也可以把这个多项式确定出来。

我们设多项式$A(x)$的系数为$(a_o,a_1,a_2,\ldots,a_{n-1})$

那么$A(x)=\sum_0^{n-1} a_i \times x^i$
将其下标按照奇偶性分类
$A(x)=\sum_{i=2k,k\in N^+}^{n-2}a_i \times x^i+\sum_{i=2k+1,k\in N^+}^{n-1}a_i \times x^i$

设$A_1(x)=a_0+a_2\times{x}+a_4\times{x^2}+\dots+a_{n-2}\times x^{\frac{n}{2}-1}$

$A_2(x)=a_1+a_3\times{x}+a_5\times{x^2}+ \dots+a_{n-1}\times x^{\frac{n}{2}-1}$

那么不难得到
$A(x)=A_1(x^2)+xA_2(x^2)$

我们将$\omega_n^k (k<\frac{n}{2}) $代入得
$A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})$

$=A_1(\omega_{\frac{n}{2}}^{k})+\omega_n^kA_2(\omega_{\frac{n}{2}}^{k})$

同理,将$\omega_n^{k+\frac{n}{2}}$代入得
$A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}(\omega_n^{2k+n})$

$=A_1(\omega_n^{2k}\times \omega_n^n)-\omega_n^kA_2(\omega_n^{2k}\times \omega_n^n)$

$=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k})$

有没有发现,两个式子只有常数项不同,可以$O(1)$进行转移,于是便满足了分治性质,可以将其计算的复杂度优化到$O(nlogn)$

根据上面推出的公式,可以得出代码:

void FFT(int limit,complex *a,int type)
{
    if(limit==1) 
    {
        return ;
    }
    complex a1[limit>>1],a2[limit>>1];
    for(int i=0;i<=limit;i+=2)
    {
        a1[i>>1]=a[i],a2[i>>1]=a[i+1];
    }
    FFT(limit>>1,a1,type);
    FFT(limit>>1,a2,type);
    complex Wn=complex(cos(2.0*Pi/limit),type*sin(2.0*Pi/limit));
    complex w=complex(1,0);
    for(int i=0;i<(limit>>1);i++,w=w*Wn)
    {
        a[i]=a1[i]+w*a2[i],
        a[i+(limit>>1)]=a1[i]-w*a2[i];
    }
}

很明显,这样常数太大,我们需要优化,而基本上使用的就是蝴蝶操作来优化,会在最后讲

快速傅里叶逆变换

当然,我们还需要由点值表示法转回系数表示法,这是需要用到快速傅里叶逆变换

在单位复数根处插值,用$\omega_i^{-1}$来替代$\omega_i$,将每个数值除以$n$即可得到原答案,我们只需将上面代码的$type$改成$-1$即可

蝴蝶操作

原序列后序列
$a_{0,1,2,3,4,5,6,7}$$a_{0,4,2,6,1,5,3,7}$

我们观察一下其二进制变化:

$000$ $001$ $010$ $011$ $100$ $101$ $110$ $111$

$000$ $100$ $010$ $110$ $001$ $101$ $011$ $111$

不难发现,新序列二进制是原数列二进制下的反转,那么我们就可以$O(n)$求得新序列了

我们以洛谷$P3803$为例,代码如下

#include <bits/stdc++.h>
using namespace std;
#define maxn 5000010
#define pi acos(-1)
#define complex_t complex<double> 
int n,m,l,r[maxn];
complex_t a[maxn],b[maxn];
void FFT(complex_t *c, int type) 
{
    for(int i=0;i<n;i++) 
    {
        if(i<r[i]) 
        {
            swap(c[i],c[r[i]]);
        }
    }
    for(int i=1;i<n;i<<=1) 
    {
        complex_t x(cos(pi/i),type*sin(pi/i));
        for(int j=0;j<n;j+=(i<<1)) 
        {
            complex_t y(1,0);
            for(int k=0;k<i;k++,y*=x) 
            {
                complex_t p=c[j+k],q=y*c[i+j+k];
                c[j+k]=p+q,c[i+j+k]=p-q;
            }
        }
    }
}
int main()
{
    scanf("%d%d",&n,&m);
    for (int i = 0; i <= n; i++) 
    {
        int x;
        scanf("%d",&x);
        a[i]=x;
    }
    for(int i=0;i<=m;i++) 
    {
        int y;
        scanf("%d",&y);
        b[i]=y;
    }
    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(a,1);
    FFT(b,1);
    for(int i=0;i<=n;i++) 
    {
        a[i]*=b[i];
    }
    FFT(a,-1);//逆变换
    for(int i=0;i<=m;i++) 
    {
        printf("%d ", (int)(a[i].real()/n+0.5));
    }
    return 0;
}

$NTT$(快速数论变换)

对于质数$p=qn+1$,$n=2^m$,原根$\omega$,则$\omega^{qn}\equiv1\pmod p$

我们可发现一些性质,如$\omega_n^n\equiv 1\pmod p$等

需要注意的事,这里的$p=k \times 2^n +1$($p$为质数)

任意模数$NTT$

我们只需要将一个大于模数的数分解为三个满足$NTT$的数相乘,最后用中国剩余定理合并即可

由于博主非常菜,就讲这么多吧

暂无评论

正在回复给  
去登陆?

   点击刷新验证码

标签云

文章留名