FFT学习笔记

FFT学习笔记

前置知识：

前言：

  • 假设我们现在有一个(n-1)次多项式，通项为(sum_{i=0}^na_i*x^i)。
  • 比如：(A(x)=x^2+2x+1,B(x)=3x^3+2x^2+5x+1)。
  • 我们想将这两个多项式相乘，采用朴素算法只能老老实实的将每一项对应相乘再相加，时间复杂度达到了(O(nm))。
  • 当然可以转换一种思路，将(n)个不同的(x)带入这个(n-1)次多项式，会取得不同的(y)，这(n个)点((x,y))唯一确定了这个多项式。
  • 这里我们可以发现，两个用点值表示的多项式相乘，我们如果(O(n))的时间枚举一下(x_i)，可以得到：(C(x_i)=A(x_i)*B(x_i))。
  • 这样可以在(O(n))的时间内完成多项式乘法。（因为(n)个点可以唯一的确定一个(n-1)次多项式）
  • 但是可惜的是，取(x_i)求(A(x_i))，这样的复杂度是(O(n))的，如果取(n)个(x)，那整体复杂度就是(O(n^2))的，和朴素枚举差不多。
  • 但如果可以快速的将多项式转换为点值表示（以及其反向操作），之后就可以快速地完成多项式乘法。
  • 快速傅里叶变换((FFT))是一种能在(O(nlogn))的时间内将一个多项式转换成点值表示的算法。
  • 所以算法流程是这样的：
    • 将两个式子从系数用(O(nlogn))时间表示转化为点值表示。
    • 然后以(O(n))的时间完成两个式子相乘。
    • 最后以(O(nlogn))时间将点值表示再转换为多项式系数表示。

几个容易混淆的缩写：

  • (DFT:)离散傅里叶变换，(O(n^2))，(FFT)的朴素版。
  • (FFT:)快速傅里叶变换，(O(nlogn)).
  • (FNTT/NTT:FFT)的升级版。
  • (FWT:)快速沃尔什变换。

单位根：

  • 下文中，默认(n)为(2)的正整数次幂。

  • 在二维平面，原点为圆心，(1)为半径做圆，所得圆为单位圆。

  • 那么我们先把圆切成(n)等份，然后对应(n)个向量。

  • 画个图看看吧：

  • 这样就把一个单位圆分为了(8)等份。

  • 根据上一篇博文的知识，其实这就是(z=e^{i})开了(n)次方。

  • (sqrt[n]{z}=e^{ifrac{2kpi}{n}}=cosfrac{2kpi}{n}+isinfrac{2kpi}{n})。

  • 我们设(w_n^k=e^{ifrac{2kpi}{n}}).

  • (w_n^0=e^{0}=cos0+isin0,w_n^1=e^{frac{i2pi}{n}}=cosfrac{2pi}{n}+isinfrac{2pi}{n},...,)

  • (w_n^{n-1}=e^{i2(n-1)pi}=cosfrac{2(n-1)pi}{n}+isinfrac{2(n-1)pi}{n}).

  • 也就对应的是圆上的几个点。

单位根的性质：

  • (1:w_{2n}^{2k}=w_n^k).
    • 证明：(w_{2n}^{2k}=e^{ifrac{2(2k)pi}{2n}}=e^{ifrac{2kpi}{n}}=w_n^k).证毕。
  • (2:w_n^{k+frac{n}{2}}=-w_n^k).
    • 就相当于一个点(+frac{n}{2})变成对面那个点，也就是这个点的反方向。
    • 或者也可以用欧拉公式展开一下用三角函数变换证明。（我太懒了）

快速傅里叶变换：

  • 我们知道一个(n-1)次多项式可以被(n)个点确定。
  • 假设(A(x)=(a_0,a_1,...,a_{n-1}))。
  • 那么有(A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1})。
  • 接下来按照下标奇偶性分类。
  • (A(x)=(a_0+a_2x^2+a_4x^4+...+a_{n-2}x^{n-2})+(a_1x+a_3x^3+...+a_{n-1}x^{n-1})).
  • 设
    • (A_1(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{frac{n}{2}-1}).
    • (A_2(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{frac{n}{2}-1}).
  • 那么有：
    • (A(x)=A_1(x^2)+xA_2(x^2)).
  • 将(w_n^k(k<frac{n}{2}))代入得：
    • (A(w_n^k)=A_1(w_n^{2k})+w_n^kA_2(w_n^{2k})).
  • 同理将(w_n^{k+frac{n}{2}})代入得：
    • (A(w_n^{k+frac{n}{2}})=A_1(w_n^{2k+n})+w_n^{k+frac{n}{2}}A_2(w_n^{2k+n})).
    • (=A_1(w_n^{2k}*w_n^n)-w_n^kA_2(w_n^{2k}*w_n^n)).
    • (=A_1(w_n^{2k})-w_n^kA_2(w_n^{2k})).
  • 可以发现，这两个式子只有一个常数项不同。
  • 也就是在枚举第一个式子的时候，可以(O(1))的得到第二个式子的值。
  • 而我们又知第一个式子中(kin[0,frac{n}{2}-1],k+frac{n}{2}in[frac{n}{2},n-1]).
  • 所以也就是将原问题缩小了一半。
  • 满足分治性质，递归后合并求解即可。
  • 于是就做到了在(O(nlogn))时间完成了多项式的系数表示到点值表示。

快速傅里叶逆变换：

  • 接下来我们需要将点值表示法还原回系数表示，这个过程为傅里叶逆变换。

  • 我们先假设((y_0,y_1,...,y_{n-1}))为多项式(A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1})为(A(x))的离散傅里叶变换。

  • 再设(B(x)=y_0+y_1x+y_2x^2+...+y_{n-1}x^{n-1})，现在我们代入单位根的倒数(w_n^0,w_n^{-1},w_n^{-2},...,w_n^{-(n-1)})可以得到一个新的离散傅里叶变换((z_0,z_1,...,z_{n-1}))。

  • 可以得到

    • (z_k=y_0+y_1(w_n^{-1})^1+y_2(w_n^{-2})^2+...+y_{n-1}(w_n^{-(n-1)})^{n-1}=sum_{i=0}^{n-1}y_i(w_n^{-k})^i).
    • (=sum_{i=0}^{n-1}(sum_{j=0}^{n-1}a_j(w_n^i)^j)(w_n^{-k})^i).
    • (=sum_{j=0}^{n-1}a_j(sum_{i=0}^{n-1}(w_n^{j-k})^i)).

  • 可以看一下(sum_{j=0}^{n-1}(w_n^{j-k})^i)。当(j-k=0)时，它等于(n)；其余的时候，可通过等比数列求和得知：

    • (frac{(w_n^{j-k})^n-1}{w_n^{j-k}-1}=frac{(w_n^n)^{j-k}-1}{w_n^{j-k}-1}=frac{1^{j-k}-1}{w_n^{j-k}-1}=0).

  • 那么就可以得知：(z_k=na_k).

    • (a_i=frac{z_i}{n}).

  • 所以我们可以得到一个结论。多项式(A(x))的离散傅里叶变换的另一个多项式(B(x))的系数，取单位根的倒数(w_n^0,w_n^{-1},...,x_n^{-(n-1)})作为(x)代入(B(x))，得到的每个数再除以(n)，得到的是(A(x))的各项系数。这就实现了傅里叶逆变换。

  • 至此(FFT)的理论基础已结束。

代码实现：

  • 根据上述分析可得，一个序列需要划分成两部分后分治递归即可。

  • 但是可以再度优化。

  • 可以发现原序列和后序列分组其实是按照原序列下标的二进制翻转。

  • 因此按照下标进行奇偶性分类是没有必要的，于是我们可以免去递归的过程。

  • 对于二进制翻转要怎么做呢？

  • 这是一个(trick:)蝴蝶定理。

  • 假设即将反转的数字为(i)，在(i)之前的数字都已经翻转好了。

  • 那么对于(i)来说就是右移后翻转后再右移，如果是奇数为在最高位补(1,)也就是(r[i] = (r[i>>1]>>1)|((i&1)<<(l-1)))。

  • 弄不明白可以手摸几个例子。

  • 之后直接枚举子区间后向上合并即可。

  • 枚举分割的中点的时间复杂度为(O(logn))，合并复杂度为(O(n))，总时间复杂度为(O(nlogn))。

  • #include<bits/stdc++.h>
using namespace std;
const int maxn = 1e7 + 10;
const double PI = acos(-1.0);
int n, m, limit, l, r[maxn];
inline int read() {
    char c = getchar(); int x = 0, f = 1;
    while (c < '0' || c > '9') {if (c == '-')f = -1; c = getchar();}
    while (c >= '0' && c <= '9') {x = x * 10 + c - '0'; c = getchar();}
    return x * f;
}

//手写复数类
struct Complex
{
    double x, y;
    Complex (double xx=0, double yy=0){
        x = xx; y = yy;
    }
    Complex operator + (const Complex b) const{
        return Complex(x+b.x, y+b.y);
    }
    Complex operator - (const Complex b) const{
        return Complex(x-b.x, y-b.y);
    }
    Complex operator * (const Complex b) const{
        return Complex(x*b.x-y*b.y, x*b.y+y*b.x);
    }
}a[maxn], b[maxn];

void fft(Complex c[], int type)
{
    for(int i = 0; i < limit; i++)
        if(i < r[i]) swap(c[i], c[r[i]]);
    //枚举待合并区间的中点的长度
    for(int mid = 1; mid < limit; mid <<= 1)
    {
        //设立单位根
        Complex wn(cos(PI/mid), type*sin(PI/mid));
        //R是区间的长度,j表示当前已经到哪个位置了,而且是左端点
        for(int R = mid<<1, j = 0; j < limit; j += R)
        {
            Complex w(1, 0); //初始幂次
            for(int k = 0; k < mid; k++, w = w*wn) //枚举左半部分
            {
                Complex x = c[j+k], y = w*c[j+mid+k];
                c[j+k] = x+y;
                c[j+mid+k] = x-y; //右半部分减去即可
            }
        }
    }
}

int main()
{
    //1e6的读入,需要写快读
    n = read(), m = read();
    for(int i = 0; i <= n; i++) a[i].x = read();
    for(int i = 0; i <= m; i++) b[i].x = read();
    limit = 1;
    //要求limit一定是2的整次幂
    while(limit <= n+m) limit <<= 1, l++;
    
    for(int i = 0; i < limit; i++)
        r[i] = (r[i>>1]>>1)|((i&1)<<(l-1));

    //对a序列和b序列分别处理
    fft(a, 1); fft(b, 1);
    for(int i = 0; i <= limit; i++)
        a[i] = a[i]*b[i];
    fft(a, -1);
    for(int i = 0; i <= n+m; i++) //四舍五入
        printf("%d ", (int)(a[i].x/limit+0.5));
    return 0;
}
【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/zxytxdy/p/12078627.html