FFT算法学习笔记

写在前边

　　1.辣鸡RRRR_wys之前csdn的博客，千年不更。。。还很水。。。于是开了这个Blog。。。妄图拯救一下自己

　　2.最近接触了一些多项式理论。于是翘掉了愉快的高频自控，通过《算导》稍稍学习了一下

　　3.算法竞赛中，FFT主要解决多项式的乘法等问题

FFT基本概念

　　1.FFT即快速傅里叶变换, 是离散傅里叶变换的加速算法。可以在o(nlogn)的时间内，完成DFT和DFT^-1

　　2.DFT即离散傅里叶变换, 这里主要就是多项式的系数向量转换成点值表示的过程

多项式的表示

　　1.系数表达：对于一个次数上界为n的多项式 A(x)=∑_{(j=0 - n-1)} a_j*x^j 而言，系数表达式便是由系数组成的向量a = (a₀, a₁, a₂, a₃, ... , a_n-1)

         (1)计算A(x₀)的值即：A(x₀) = a₀+x₀(a₁+x₀(a₂+x₀(a₃+x₀(a₄+...+x₀(a_n-2+x₀*a_n-1))...)))  时间为o(n)

　　  (2)通过系数向量计算两个多项式的乘法：设输出向量为c, 而a,b长度n=max(na,nb)那么显然两个c的长度为2n-1，计算方法为：c_j = ∑_{(k=0 ~ j)} a_k*b_j-k，

　　      即向量a,b的卷积

　　2.点值表达式：对于一个次数上界为n的多项式 A(x)=∑_{(j=0 ~ n-1)} a_j*x^j 而言，点值表达式便是由n个点值所构成的集合{ (x₀,y₀), (x₁,y₁), (x₂,y₂), ... , (x_n-1,y_n-1) }

　　    (1)运用之前提到的方法，一个个计算点值时间为o(n²)，之后可以看到如果我们巧妙地选取点x_{k可以使时间变为o(nlogn)}

　　    (2)将点值表达式转换为系数表达式：将已知的点值带入，多项式可已得到一个有n个未知数的n个线性方程

　　　　 a₀ + a₁*x₀ + a₂*x₀² + a₃* x₀³ + ... + a_n-1*x₀^n-1 = y₀ (1)

　　　　 a₀ + a₁*x₁ + a₂*x₁² + a₃* x₁³ + ... + a_n-1*x₁^n-1 = y₁ (2)

　　　　 .......

　　　　 a₀ + a₁*x_n-1 + a₂*x_n-1² + a₃* x_n-1³ + ... + a_n-1*x_n-1^n-1 = y_n-1 (n)

 　　　　 可以写成矩阵形式，既为一个范德蒙矩阵，通过线代知识解这个方程就可以求出系数向量，复杂度o(n³)。另一种方法是采用拉格朗日插值法求解，复杂度o(n²)。因此，n个点的求值运算和插值运算，是定义完备的互逆运算。

　　　 (3)通过点值计算两个多项式的乘积，c = { (x₀,y_A0*y_B0), (x₁,y_A1*y_B1), (x₂,y_A2*y_B2), ... , (x_n-1,y_A(n-1)*y_B(n-1)) }, 可以看到复杂度为o(n)

　　   (4)计算多项式乘法的过程：1）o(nlogn)分别求出两个多项式的点值表达式 2）点值乘法o(n) 3）插值o(nlogn) 4）求出答案多项式的系数表达式

DFT与FFT

　　下面开始讨论，如何在o(nlogn)的时间内，完成求值和插值运算

　　1.单位复数根

    　　n次单位复数根满足ωⁿ=1的复数ω，n此单位复数根的数目恰好有n个，这些根是：ω_n^k = e^2πik/n 这n个根均匀的分布在复平面上，其他n次单位根都是ω_n的幂次。n个n次单位根在乘法意义下形成一个群，该群与模n意义下的加法群，有相同的结构，即ω_nⁱ ω_n^j = ω_n^{(i + j) mod n}

　　　(1) 消去引理：ω_nd^kd =ω_n^k  (n>=0, k>=0, d>0)

　　　(2) ω_n^n/2 =ω₂ =-1

　　　(3) 折半引理：如果n>0为偶数，那么n个n次单位复数根的平方的集合就是n/2个n/2次单位复数根的集合

　　　(4) 求和引理：∑_{(j=0 ~ n-1)} (ω_n^k)^j = 0 (n>=1, 不能被n整除的非负整数k)

　　2.DFT

　　我们希望计算ω_n⁰, ω_n¹, ω_n², .... , ω_n^n-1处的值（在多项式乘法中其实是2n个点值）。假设给定系数a=(a₀, a₁, a₂, ... , a_n-1), 即对于k=0, 1, 2, 3, ... , n-1求出y_k = A(ω_n^k)，记向量y=(y₁, y₂, y₃, ... ,y_n-1)就是系数向量a的离散傅里叶变换DFT，记为 y = DFT_n(a)

　　3. FFT

　　通过快速傅里叶变换FFT，利用单位复数根的性质，我们可以在o(nlogn)的时间内计算DFT_n(a)。首先假设次数n恰好是2的整数幂，不足在高次位置项添0。

　　定义两个新的次数界为n/2的多项式 A_[0] = a₀ + a₂*x + a₄*x² + a₆* x³ + ... + a_n-2*x^n/2-1 , A_[1] = a₁ + a₃*x + a₅*x² + a₇* x³ + ... + a_n-1*x^n/2-1

　　A_[0]包含所有偶数下标的系数，A_[1]包含所有奇数下标的系数，于是有A(x) = A_[0](x2) + x*A_[1](x²) 。

　　所以，求解A(x)的DFN，转换为了求解多项式A_[0], A_[1]在(ω_n⁰)², (ω_n¹)², (ω_n²)² , ... , (ω_n^n-1)² 处的取值。

　　根据折半引理，(ω_n⁰)², (ω_n¹)², (ω_n²)² , ... , (ω_n^n-1)² 并不是由n个不同的值组成，而是由n/2个n/2次单位根复数根所组成，每个根正好出现两次。因此，我们递归的对次数界为n/2的多项式A_[0], A_[1]在n/2个n/2次单位复数根处进行求值。这样我们求出了次数界为n的次多项式在n次单位复数根处的值。

　　注意到除了递归调用的时间外，每次所需时间为o(n)，其中n是输入向量的长度。因此有如下递归式：T(n) = 2T(n/2) + o(n) = o(nlogn)，时间为o(nlogn)。

　　计算在单位复数根处插值，求解系数向量。根据点值表达式，我们可以列出线性方程组：

　　a₀ + a₁ + a₂ + a₃+ ... + a_n-1 = y0 (1)

　　a₀ + a₁* (ω_n)¹+ a₂*(ω_n)² + a₃* (ω_n)³ + ... + a_n-1*(ω_n)^n-1 = y₁ (2)

　　a₀ + a₁* (ω_n)²+ a₂*(ω_n)⁴ + a₃* (ω_n)⁶ + ... + a_n-1*(ω_n)^2(n-1)  =  y₂  (3)

　　......

　　a₀ + a₁* (ω_n)^n-1+ a₂*(ω_n)^2(n-1) + a₃* (ω_n)^3(n-1) + ... + a_n-1*(ω_n)^(n-1)(n-1)  =  y_n-1 (4)　

　　根据这个线性方程组，我们可以令单位复数根组成的系数矩阵为V_n，运用yV_n^-1求解系数。可以证明：V_n^-1的(j, k)处的值为(ω_n)^-kj/n  (j,k = 0, 1, 2, ... , n-1)

　　有了上面的式子我们能推导出DFT^-1：a_j =(1/n)* ∑_{(k=0 ~ n-1)} y_k(ω_n)^-kj  (j = 0, 1, 2, ... , n-1)

　　观察可得，相对于DFT，我们把a与y替换，用ω_n^-1代替ω_n，并将每个元素除以n。这样，我们也可以在o(nlogn)的时间内计算出DFT^-1

　　卷积定理：a(卷积)b = DFT_2n^-1( DFT_2n(a)*DFT_2n(b) ) , 其中向量a, b用0填充使其长度达到2n。

[UR#34]多项式乘法

递归实现：参考hzwer的代码

#include <bits/stdc++.h>
#define pi acos(-1.0)
typedef long long ll;
using namespace std;
int n,m;
complex<double> a[300000],b[300000];
void fft(complex<double> *x, int n, int ty){
    if(n==1)return;
    complex<double> l[n>>1],r[n>>1];
    for(int i=0;i<=n;i+=2)
        l[i>>1]=x[i],r[i>>1]=x[i+1];
    fft(l,n>>1,ty);fft(r,n>>1,ty);
    complex<double> wn(cos(2*pi/n),sin(ty*2*pi/n)),w(1,0),t;
    for(int i=0;i<n>>1;i++,w*=wn)
        t=w*r[i],x[i]=l[i]+t,x[i+(n>>1)]=l[i]-t;
}
int main()
{
    scanf("%d%d",&n,&m);
    for(int i=0,x;i<=n;++i)scanf("%d",&x),a[i]=x;
    for(int i=0,x;i<=m;++i)scanf("%d",&x),b[i]=x;
    m=n+m;for(n=1;n<=m;n<<=1);
    fft(a,n,1);fft(b,n,1);
    for(int i=0;i<=n;++i)a[i]=a[i]*b[i];
    fft(a,n,-1);
    for(int i=0;i<=m;++i)
        printf("%d ",(int)(a[i].real()/n+0.5));
    return 0;
}

　　

　　然而，实际应用中，需要算法尽量地快速，因此，我们引出一种小常数FFT的迭代实现。

高效FFT实现

　　1.蝴蝶操作：

    for(int i=0;i<n>>1;i++,w*=wn)
        t=w*r[i],x[i]=l[i]+t,x[i+(n>>1)]=l[i]-t;

　　2.考虑递归调用的过程，将每次输入向量构成一颗树，类似于zkw线段树，如果可以实现自底向上更新，相比递归常数就很优秀了。

　　  (1)首先，我们成对的取出输入向量树最底层的元素，利用蝴蝶操作计算出每对的DFT，这样我们就有了n/2个二元素DFT，同理，下一步取出这n/2个DFT，计算四元素DFT，最终形成n元素DFT

　　  (2)那么如何确定最底层的元素排列呢？尝试一下就能发现，元素出现的顺序就是是一个位逆序置换，即：n = 8时，{0, 4, 2, 6, 1, 5, 3, 7}, 写成二进制形式，即{000, 100, 010, 110, 001, 101, 011, 111}，将所有元素二进制位逆序可得{000, 001, 010, 011, 100, 101, 110, 111} = {0, 1, 2 ,3, 4, 5, 6, 7}

                                                                                        

[UR#34]多项式乘法

迭代实现：参考hzwer的代码

#include <bits/stdc++.h>
#define pi acos(-1.0)
typedef complex<double> E;
using namespace std;
int n,m,L,R[300000];
E a[300000],b[300000];
void fft(E *a,int f){
    for(int i=0;i<n;++i)if(i<R[i])swap(a[i],a[R[i]]);
    for(int i=1;i<n;i<<=1){
        E wn(cos(pi/i),f*sin(pi/i));
        for(int p=(i<<1),j=0;j<n;j+=p){
            E w(1,0);
            for(int k=0;k<i;++k,w*=wn){
                E x=a[j+k],y=w*a[j+k+i];
                a[j+k]=x+y;a[j+k+i]=x-y;
            }
        }
    }
}
int main()
{
    scanf("%d%d",&n,&m);
    for(int i=0,x;i<=n;++i)scanf("%d",&x),a[i]=x;
    for(int i=0,x;i<=m;++i)scanf("%d",&x),b[i]=x;
    m=n+m;for(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]=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(Number Theory Transform)

　　可是，复数运算的常数依然是很巨大的，而且精度不高。当题目给的模数很特殊的时候，我们就可以利用数论变换来加速算法。