FFT入门

FFT入门

背景知识

FFT是快速求解两个多项式乘积的算法,时间复杂度为优秀的(O(nlogn)),主要想法是分治。

FFT = Fast离散傅里叶变换 = Fast DFT

多项式存在系数表示法和点值表示法,二者可以相互转换;
n次多项式的最高次为n,有(n+1)个系数,需要(n+1)个相异的点才能表示

FFT主要思路

  1. 预处理,填充高次为0,使得最终系数数目大于(N+M);(N次多项式和M次多项式)
  2. 系数表示转化成点值表示(代入的x坐标为(x = (w_n^0,w_n^1...w_n^{n-1}));这里的n是需要的点的个数;这个可以通过分治实现(O(nlogn))
  3. 点值表示法相乘,复杂度(O(n))
  4. 点值表示转化成系数表示,该问题可规约成步骤2

下面我们将重点放在2,3两点上

要点

多项式系数表示转化成点值表示

记多项式(A(x)=a_0+a_1x+a_2x^2+...a_{n-1}x^{n-1}), (A(x))是一个(n-1)次多项式,需要(n)个点才能表示。我们选取复数域上(w_n^k = e^{i*frac{2kpi}{n}} = cos(frac{2k pi}{n})+i sin(frac{2k pi}{n}) quad k=0,1,2,..n-1) 这样的i个点。

它们满足两条性质:

  1. (w_n^k = -w_n^{k+frac{n}{2}}quad) 逆时针旋转(pi)
  2. (w_n^k = w_{nd}^{kd}),假设(kd,nd)均为整数

上述两条性质在复数的单位圆上画一下即可知。

假设(n=2^t),实际上步骤0中预处理就会处理成(2^t);我们将(A(x)) 根据系数位置的就分成两部分,即:
(A1(x) = a_0+a_2x^2+..a_{n-2}^x{n-2} = sum_{i=0}^{i=frac{n}{2}-1}a_{2i}x^{2i})

(A2(x) = a_1+a_3x^3+...a_{n-1}x^{n-1} = sum_{i=0}^{i=frac{n}{2}-1}a_{2i+1}x^{2i+1} = xsum_{i=0}^{i=frac{n}{2}-1}a_{2i+1}x^{2i})

(A(x)=sum_{i=0}^{n-1}a_ix^i = A1(x)+A2(x))

注意到(A1(x), A2(x))虽然系数减少了一半,但需要带入的(x)依旧是(n)项;
假设我们将(x=w_n^k)代入(A1(x)) (代入(A2(x))是一样的),那么:

(A1(w_n^k)=sum_{i=0}^{i=frac{n}{2}-1}a_{2i}w_n^{2ki}=sum_{i=0}^{i=frac{n}{2}-1}a_{2i}w_{frac{n}{2}}^{ki} quad k=0,1,2....n-1),(根据性质2)

再一次变化,将k分为前半部分和后半部分处理:

(A1(w_n^k)=sum_{i=0}^{i=frac{n}{2}-1}a_{2i}w_{frac{n}{2}}^{ki} quad k =0,1,2,frac{n}{2}-1), 前半部分

(A1(w_n^{k+frac{n}{2}})=-sum_{i=0}^{i=frac{n}{2}-1}a_{2i}w_{frac{n}{2}}^{ki} quad k =0,1,2,frac{n}{2}-1) 后半部分,(运用性质1)

由此可知我们在求得(A1(w_n^k))的同时已经把(A1(w_n^{k+frac{n}{2}}))求了出来,这样求解(A1(w_n^k))的规模相对于(A(x))将为了一半,多么nice的结果((A2(x))同理)。

所以 (T(n) = 2*(n/2)+O(n)),因为(x)(n)种取值,所以通过(A1(x))(A2(x))求得(A(x))需要(O(n)) 时间

递归写法
void FFT(Complex* a,int len,int opt){
    if(len==1)
        return;
    Complex* a0=new Complex[len/2];
    Complex* a1=new Complex[len/2];
    for(int i=0;i<len;i+=2){
        a0[i/2]=a[i];
        a1[i/2]=a[i+1];
    }
    FFT(a0,len/2);
    FFT(a1,len/2);
    Complex wn(cos(2*Pi/len),opt*sin(2*Pi/len));
    Complex w(1,0);
    for(int i=0;i<(len/2);i++){
        a[i]=a0[i]+w*a1[i];
        a[i+len/2]=a0[i]-w*a1[i];
        w=w*wn;
    }
    delete[] a0;
    delete[] a1;
}

点值转换成系数表示(有人称作IDFT)

思考一下系数表示到点值的矩阵形式:

[egin{bmatrix} w_n^0 & w_n^{0*1} & w_n^{0*2} & ... &w_n^{0*{n-1}} \ w_n^{1*0} & w_n^{1*1} &w_n{1*2} & ... & w_n^{1*{n-1}} \ ... & \ w_n^{{(n-1)}*0} & w_n^{{(n-1)}*1} & w_n^{(n-1)*2} & .. & w_n^{{(n-1)}*{(n-1)} } end{bmatrix} egin{bmatrix} a0 \ a1 \ .. \ a_{n-1} end{bmatrix} = egin{bmatrix} y0 \ y1 \ .. \ y_{n-1} end{bmatrix} ]

我们记作 (W cdot A = Y); 那么从点值到系数的变换即为 (A = W^{-1}cdot Y)
下面不加证明的给出 (W^{-1}=frac{1}{n}[w^{-ij}]), (w^{ij})表示(W) i行j列的元素。即W的逆的W原矩阵的每个复数取逆后的矩阵的(frac{1}{n}). 这样步骤3就规约到了步骤1

FFT 优化

从递归形式到迭代形式,直接给出详细注释的代码,可在luogu3803交题

迭代版本

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int MAXN=1e7+10;
const double Pi=acos(-1.0);
struct complex{
    double x,y;
    complex (double xx=0,double yy=0){x=xx,y=yy;}
}a[MAXN],b[MAXN];
complex operator + (complex a,complex b){ return complex(a.x+b.x , a.y+b.y);}
complex operator - (complex a,complex b){ return complex(a.x-b.x , a.y-b.y);}
complex operator * (complex a,complex b){ return complex(a.x*b.x-a.y*b.y , a.x*b.y+a.y*b.x);}//复数的运算那部分 
int N,M;
int l,r[MAXN];
int limit=1;
void fft(complex *A,int type)
{
    for(int i=0;i<limit;i++) 
        if(i<r[i]) swap(A[i],A[r[i]]);//求出要迭代的序列,即递归到最底层的序列
    for(int mid=1;mid<limit;mid<<=1)//待合并区间的中点,其实(mid<<1)是要求的DFT(某一层)的长度
    {
        complex Wn( cos(Pi/mid) , type*sin(Pi/mid) ); //单位根 本质是 2*PI/(2*mid)
        for(int R=mid<<1,j=0;j<limit;j+=R)//R实际是步长,即同一层每组要处理的起始位置
        {
            complex w(1,0);//幂  单位1
            for(int k=0;k<mid;k++,w=w*Wn)//枚举左半部分 
            {
                 complex x=A[j+k],y=w*A[j+mid+k];//蝴蝶效应(只是装x),实际上什么有没有就是减少一次y的计算(复数乘法)
                A[j+k]=x+y;                                   // A[k]
                A[j+mid+k]=x-y; //mid即位推导中的2/n;此处可以看作 A[k+2/n]
            }
        }
    }
}

int main()
{
    //n次多项式由(n+1)个系数,有n个复数解,最高次为n
    scanf("%d",&N);scanf("%d",&M);
    for(int i=0;i<=N;i++) scanf("%lf",&a[i].x);//复数域的实数部分
    for(int i=0;i<=M;i++) scanf("%lf",&b[i].x);
    while(limit<=N+M) limit<<=1,l++; //最高次为(N+M) 所以至少需要(N+M+1)个点值;为了fft的方便扩展成 2^l(2^l>(N+M)) 高次皆为0
    for(int i=0;i<limit;i++)
        r[i]= ( r[i>>1]>>1 )| ( (i&1)<<(l-1) ) ; //位逆置 r[i>>1]>>1 能反转前面的l-1 位; 然后将末尾的位或到首位即可
    // 在原序列中 i 与 i/2 的关系是 : i可以看做是i/2的二进制上的每一位左移一位得来
    // 那么在反转后的数组中就需要右移一位,同时特殊处理一下复数 
    fft(a,1);// fft保留的是x= w_{limit}^1, w_{limit}^2 .... w_{limit}^{limit}是的取值
    fft(b,1);
    for(int i=0;i<=limit;i++) a[i]=a[i]*b[i]; //点值表示法的优秀乘法
    fft(a,-1);  //点值表示法到系数表示法,会涉及到逆矩阵,所以fft后还有一个1/limit (limit 实际上是矩阵的行列式(n*n的满秩矩阵才有行列式))
    for(int i=0;i<=N+M;i++)
        printf("%d ",(int)(a[i].x/limit+0.1));//直接用实数域
    return 0;
}

参考

blog0 blog1 blog2 blog4

原文地址:https://www.cnblogs.com/fridayfang/p/11094045.html