FFT——快速傅里叶变换

xiejinhao · 2019-07-15 19:27:54 · 个人记录

FFT——快速傅里叶变换

------------ $FFT$全称($Fast \ Fourier \ Transformation$)，中文名：快速傅里叶离散变换。这个名字听起来很高级，实际上也很高级，它是$DFT$和$IDFT$的加速版本，用于加速多项式乘法。 - 接下来我们把某个多项式记为“大写字母+元素”，如多项式： $A(x)=a_nx^n+a_{n-1}x^{n-1}+…+a_2x^2+a_1x^1+a_0

B(x)=b_mx^m+b_{m-1}x^{m-1}+…+b_2x^2+b_1x^1+b_0

你也可以理解为A/B是关于x的函数

现在，我们要求A(x)*B(x)，即两个多项式的乘积(卷积)，那么我们很容易想到的就是O(N^2)的做法：把每一项依次分别相乘再相加。

\sum^{n+m-1}_{i=0}(\sum^{i}_{j=0}a_i*b_{i-j})*x^i

这样的时间复杂度明显不尽人意，那么FFT就是把这种算式计算优化到NlogN的算法。

------------ - **复数** 我们定义虚数 $i$，且 $i^2 = -1$，即 $i = \sqrt {-1}$，那么复数就形如 $a+bi$，$i$ 是一个虚数单位，$a$ 叫做复数的实部， $b$ 叫做复数的虚部，那么 $a+bi$ 在复平面（高中会讲）可以类似于平面直角坐标系，表示为点 $(a, b)$。 ![](https://i.loli.net/2019/07/16/5d2d21a9f0b8789604.png) **复数的运算法则** 1. 加法：**实部和虚部分别相加。** 2. 减法：**实部和虚部分别相减。** 1. 乘法：**像一次多项式一样相乘。** 1. 除法：**分子分母同乘分母的共轭。** **共轭**：两个复数实部相同，虚部相反。上面除法分子分母同乘分母的共轭，把分母有理化成实数。 ```cpp struct complex { //复数 long double a, b; complex() { a = b = 0;} //重载复数四则运算 complex operator + (const complex x) const { complex temp; temp.a = a + x.a; temp.b = b + x.b; return temp; } complex operator - (const complex x) const { complex temp; temp.a = a - x.a; temp.b = b - x.b; return temp; } complex operator * (const complex x) const { complex temp; temp.a = a * x.a - b * x.b; temp.b = a * x.b + b * x.a; return temp; } complex operator / (const complex x) const { double div = x.a * x.a + x.b * x.b; complex temp; temp.a = (a * x.a + b * x.b) / div; temp.b = (b * x.a - a * x.b) / div; return temp; } }; ``` 这边大家直接记一个重要结论： ### **复数相乘，辐角相加，模长相乘。** ![](https://i.loli.net/2019/07/16/5d2cae88d057313942.png) ------------ ### **单位根** 对于复数$w$,若整数 $n$ 满足 $w^n=1$，则称 $w$ 是 $n$ 次单位根。扔个单位圆： ![](https://i.loli.net/2019/07/16/5d2d23a24c3fe95897.png) 为什么扔个单位圆呢？因为单位圆上任意一点到原点的距离，即模长都为$1$。 **那么，对于一个复数，只有模长为$1$的复数才有可能成为$n$次单位根!** 现在，我们在复平面上把复数锁定在了这个单位圆上。设复数$w$模长为$q$，其幅角为$mπ$，$m$为$[0,2]$间实数，那么，其$n$次方幅角为$nmπ$。则有$nmπ=kπ,k∈N$，即$nm=k$。 1. 若$k=0$,因为$n>0$，所以方程有且只有$m=0$一解。 1. 若$k>0$，则$k=1,2,3,…$，易知$m$有$n-1$个不重解综上，复数$w$的$n$次单位根有$n$个。并且，复数的$n$次单位根$n$等分单位圆。 **单位根的书写：** 单位根用符号$ω$表示，从1开始，沿单位圆逆时针把单位根从$0$开始标号。 ![](https://i.loli.net/2019/07/16/5d2d2c8eb2bd792696.png) 然后单位根还有个神奇的东西$ω^{k}_{n}=ω_{n}^{k\mod n}

单位根的性质

对于任意n，ω_n^0=1。
(ω_n^m)^k=(ω_n^k)^m
若 2 \mid n，ω_n^{m+\frac{n}{2}}=-ω_n^m

多项式的表达方法

系数表达法：这个如同我们上面写的多项式即为系数表达法。
点值表达法

定义多项式F(x)：

我们知道，1个n元一次方程最少需要n个式子求解，而我们把系数表达F(x)抽象为一个函数F(x)，那么这个函数至少需要n+1个点确定下来，那么我们任取n+1个不同的数构成集合U={k_1,k_2,…,k_{n+1}}，对F(x)分别求值得到一个新的集合，那么将两个集合合并成点集，有：

T={(k_1,F(k_1)),k_2,F(k_2),…,(k_{n+1},F(n+1))}

这便是F(x)的点值表达式，并且，点值表达法下的乘法运算更加简单：

多项式P,Q分别取点(x,y_1)和(x,y_2)，则P*Q可得到点(x,y_1*y_2)。令S=P*Q，则C(x_i)=y_{1_i}*y_{2_i}。

可见，点值表达法下，多项式乘法是O(N)的。

（补充）乘法本质：卷积

根据我们上面暴力计算A(x)*B(x)，不妨设C=A*B，那么：

C[k]=\sum^{k}_{i=0}A[i]B[k-i]

那么就可以写成：$C[k]=\sum\limits_{i+j=k}A[i]B[j]

那么形如C[k]=\sum\limits_{i⊕j=k}A[i]B[j](⊕为某种运算)的式子称为卷积。

那么，多项式的乘法就是对于加法的卷积。

我们仍然可以对多项式相乘的模板打个暴力：

    #include<bits/stdc++.h>
    using namespace std;

    const int SIZE = 1e6 + 10;

    long long A[SIZE], B[SIZE], C[SIZE];
    int n, m;

    int main() {
        scanf("%d %d", &n, &m);
        n++, m++;
        for (int i = 0; i < n; i++) scanf("%lld", A + i);
        for (int i = 0; i < m; i++) scanf("%lld", B + i);
        for (int k = 0; k < n+m-1; k++)
            for (int i = 0; i <= k; i++)
                C[k] += A[i] * B[k - i];
        for (int i = 0; i< n + m - 1; i++) 
            printf("%lld ", C[i]);
        return 0;
    }

模板题亲测 44 分。

------------ 上面的主要部分还是写了一个$O(N^2)$的暴力，那我们要寻找一个更快的办法： ### **就是它！$FFT$算法** 前面我们已经说过：$FFT$算法是一个$O(NlogN)$的算法，但常数较大(相信你们会卡常)，并且，$FFT$是$DFT$和$IDFT$的快速算法： ------------ ### **思想：** 首先你要会点值表达式，当然我们已经讲过。然后我们来看$A(x)=x^2+3x-2$（红色）与$B(x)=-x^2+6x+4$(绿色)的取点得到点的点值表达式: ![](https://i.loli.net/2019/07/15/5d2c8bf29f89588834.png) 在$A(x)$上取三个点（紫色）$(-4,2)(-2,-4)(1,2)$；在$B(x)$上取三个点$(0,4)(2,12)(6,4)$，然后转为系数表达$O(N^2)$相乘然后再带入$x$得到$C$的点值表达式？这时绝对不行的，我们上面说过，点值表达式相乘可以在$O(n)$内解决，那怎么办呢？我们可以取$3$个$x$相同的点： ![](https://i.loli.net/2019/07/15/5d2c8bf29ed1b15238.png) 这样就行了？那是不可能的。我们知道，两个二次项式相乘得到的是一个四次项式，而$C(x)=y_1*y_2$，只能得到三个点，这肯定是不行的。这很明显我们需要$2n+1$个点，不够怎么办呢？再添两个不就好了？下面添加了$I,J$(横坐标相同)，$G,H$(横坐标相同)，这样就行了。 ![](https://i.loli.net/2019/07/15/5d2c8bf29de0a97307.png) 那么，我们只需要进行$2n+1$次乘法即可，省略常数（都说了常数巨大了），复杂度$O(N)

神仙问题：tql but 这有用吗

当然有用啊：

算法流程：

\text{系数表达式->（DFT）点值表达式->（IDFT）系数表达式}

~~这不就有个用了吗？~~

但这个思路仅仅是DFT+IDFT，而FFT是用来加速DFT和IDFT的。

如果让我们求点值表达式，那我们肯定会带整数进去，因为这样好计算。但是傅里叶，（~~不知道是不是脑抽了~~），把复数带进了多项式中，然后神奇的事情就发生了……

我们把多项式如下拆开：

F(x)=(a_0+a_2x^2+…+a_{n-2}x^{n-2})+(a_1x^1+a_3x^3+…+a_{n-1}x^{n-1})

（保证n=2^k,k∈N^*)

然后设两个有n/2次项的多项式FL(x)和FR(x)

FL(x)=a_0+a_2x+…+a_{n-2}x^{\frac{n}{2}-1}

FR(x)=a_1+a_3x+…+a_{n-1}x^{\frac{n}{2}-1}

F(x)=FL(x^2)+xFR(x^2)

设k<\frac{n}{2}，把ω_n^k代入F(x)：F(ω_n^k)=FL(ω_n^{2k})+ω_n^kFR(ω_n^{2k})

接着，就有F(ω_n^k)=FL(ω_{n/2}^{k})+ω_n^kFR(ω_{n/2}^{k})

那么，如果我们知道FL(x),FR(x)分别在ω^0_{n/2},ω^1_{n/2},…,ω^{n/2-1}_{n/2}处的点值表示，那么我们就可以O(N)内求出F(x)在ω^0_{n},ω^1_{n},…,ω^{n/2-1}_{n}处的点值表示了。

接下来我们的目标就是求F(x)在ω^{n/2}_{n},ω^{n/2+1}_{n},…,ω^{n-1}_{n}处的点值表示。

然后……

继续设k<n/2，然后把ω^{k+n/2}_{n}代入F(x)中。

运用一系列单位根的性质，读者可以自行证明出：

F(ω^{k+n/2}_{n})=FL(ω^{k}_{n/2})-ω^{k}_{n}FR(ω^{k}_{n/2})

然后我们就可以求出F(x)在ω^0_{n},ω^1_{n},…,ω^{n-1}_{n}处的点值表示，接着，我们就实现了DFT（~~听着好累的样子……~~）

但是我们还不知道FL(x),FR(x)在ω^0_{n},ω^1_{n},…,ω^{n/2-1}_{n}处的点值表示。难道要暴力代入吗？这样反而会使时间复杂度上升。

但是我们想想：FL(x)和FR(x)是由原来问题F(x)转化来的，怎么转化的？看到这里，你应该想到了分治。我们知道，分治的时间复杂度是logN，~~好像完成了一半了？~~

然后，因为FL(x),FR(x)是F(x)的子问题，那么我们可以对它们继续进行分解成一个个小的子问题，然后再合并，类似于归并排序。

### 1. $F(ω^{k}_{n})=FL(ω^{k}_{n/2})+ω^{k}_{n}FR(ω^{k}_{n/2})

2. F(ω^{k+n/2}_{n})=FL(ω^{k}_{n/2})-ω^{k}_{n}FR(ω^{k}_{n/2})

------------ 上文的所有证明都基于$2=2^k,k∈N^*$，但是我们做题的时候并没有这一点，所以我们可以补一些系数为$0$的项，这对计算结果没有影响。 ------------ ### 1. 预处理单位根 **我们知道一个性质**$ω_n^k=(ω_n^1)^k

那么，我们只要知道ω_n^1就可以求出ω_n^0,ω_n^1,ω_n^2,…,ω_n^{n-1}。

怎么求ω_n^1？

解三角形不就好了？已知|OA|=1,幅角是\frac{1}{n}圆周，那么：

(由于c++下三角函数用弧度制表示，以下全部使用弧度制。)

设幅角α=\frac{2π}{n},，则ω_n^1(cos\big(\frac{2π}{n} \big),sin\big(\frac{2π}{n} \big))。

这边要牢记，x是实轴，y是虚轴：

    struct complex { //复数 
        long double a, b;
        complex() { a = b = 0;}
            //重载复数四则运算 
        complex operator + (const complex x) const {
            complex temp;
            temp.a = a + x.a;
            temp.b = b + x.b;
            return temp; 
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.a = a - x.a;
            temp.b = b - x.b;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.a = a * x.a - b * x.b;
            temp.b = a * x.b + b * x.a;
            return temp; 
        }

        complex operator *= (const complex &x) {
            *this = *this * x;
            return *this;
        }

        complex operator / (const complex x) const {
            double div = x.a * x.a + x.b * x.b;
            complex temp;
            temp.a = (a * x.a + b * x.b) / div;
            temp.b = (b * x.a - a * x.b) / div;
            return temp; 
        }
    }w[SIZE],temp, buf;

    void unit_roots(int n) {
        const long double pi = 3.14159265358979;
        //complex ;
        temp.a = cos(2*pi/n), temp.b = sin(2*pi/n);
        w[0].a = 1, w[0].b = 0;
        buf.a = 1, buf.b = 0;
        for(int i = 1; i < n ; i++) {
            buf *= temp;
            w[i] = buf; 
        }
        for(int i = 0; i < n; i++)
            printf("w_%d %.3Lf %.3Lf\n", i, w[i].a, w[i].b);    
    }

然后我们就可以实现DFT了

因为还没讲IDFT，先写DFT:

递归DFT实现（复数板子上面有了下面就不写了）

    int n, m;
    Complex f[SIZE];
    void dft(Complex *f, int len) { //当前子问题长度
        if(len == 0) return;
        Complex fl[len + 10], fr[len + 10];
        for(int i = 0 ; i < len; i++)
            fl[i] = f[i << 1], fr[i] = f[i << 1 | 1];
        dft(fl, len >> 1);
        dft(fr, len >> 1);
        len *= 2;
        Complex temp, buf;
        temp.a = cos(2*pi/n), temp.b = sin(2*pi/n);
        buf.a = 1;
        for(int i = 0 ; i < len / 2; i++) {
            f[i] = fl[i] + buf * fr[i];
            f[i + len / 2] = fl[i] - buf * fr[i];
            buf *= temp;
        }
    }
    int main() {

        scanf("%d", &n);
        for(int i = 0; i < n; i++) scanf("%Lf", &f[i].a);
        for(m = 1; m < n; m <<= 1);//补到2的幂次方
        dft(f, m >> 1);
        for(int i = 0; i < m; i++) printf("%.3Lf ", f[i].a);
        return 0;
    }

事实上虚数可以用(小心TLE)：


            complex <double> f;

然后我们得到了一堆点值表达……

这并没有什么用

于是，我们现在来学习IDFT——DFT的逆运算。

刚刚的多项式为F(x)，系数是A数列，点值表示成M，M[k]=\sum \limits_{i=0}^{n-1}(ω_n^{k})^i*F[i]

然后就有F[k]=\frac{1}{n}\sum \limits_{i=0}^{n-1}(ω_n^{-k})^i*M[i]

结论记住就好了，可以自己尝试证明。

我们可以先求出和式，然后再除以n。

接下来我们求出ω_n^0,ω_n^{-1},…ω_n^{-n+1}

并且仍有：ω_n^{-j}=(ω_n^{-1})^j

所以我们还是求出ω_0^{-1}即可，并且，ω_0^{-1}与ω_0^{1}关于y轴对称。

    #include<iostream>
    #include<cstdio>
    #include<cmath>
    using namespace std;
    const long double pi = 3.141592653589793;
    const int SIZE = 4e6 + 10;
    struct complex { //复数 
        long double real, imag;
        complex() { real = imag = (double)0;}
            //重载复数四则运算 
        complex operator + (const complex x) const {
            complex temp;
            temp.real = real + x.real;
            temp.imag = imag + x.imag;
            return temp; 
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.real = real - x.real;
            temp.imag = imag - x.imag;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.real = real * x.real - imag * x.imag;
            temp.imag = real * x.imag + imag * x.real;
            return temp; 
        }

        complex operator *= (const complex &x) {
            *this = *this * x;
            return *this;
        }

        complex operator / (const complex x) const {
            double div = x.real * x.real + x.imag * x.imag;
            complex temp;
            temp.real = (real * x.real + imag * x.imag) / div;
            temp.imag = (imag * x.real - real * x.imag) / div;
            return temp; 
        } 
        void conj() { imag = -imag; } 
    }a[SIZE], b[SIZE];

    complex omega(int n, int k) {
        complex temp;
        temp.real = cos(2 * pi * k / n);
        temp.imag = sin(2 * pi * k / n);
        return temp;
    }   
    int n, m, p;        
    void fft(complex *f, int len, int inv) {
        if(len == 1) return;
        int mid = len >>1;
        complex fl[mid + 10], fr[mid + 10];
        for(int i = 0; i < mid; i++)
            fl[i] = f[i << 1], fr[i] = f[i << 1 | 1];
        fft(fl, mid, inv);
        fft(fr, mid, inv);
        complex temp, buf; 
        temp.real = cos(2 * pi / len);
         temp.imag = sin(2 * pi / len);
        buf.real = 1, buf.imag = 0;
        if(inv) temp.conj();
        for(int i = 0; i < mid; i++, buf *= temp) {
            complex s = buf * fr[i];
            f[i] = fl[i] + s;
            f[i + mid] = fl[i] - s;
        }
    }

    int main() {

        scanf("%d %d", &n, &m);
        for(p = 1; p < n + m + 1; p <<= 1);
        for(int i = 0; i <= n; i++) 
            scanf("%Lf", &a[i].real);
        for(int i = 0; i <= m; i++) 
            scanf("%Lf", &b[i].real);
        fft(a, p, 0);
        fft(b, p, 0);
        for(int i = 0; i < p; i++) 
            a[i] *= b[i];
        fft(a, p, 1);
        for(int i = 0; i <= n + m; i++)
            printf("%d ", (int)(a[i].real / p + 0.5));
    return 0;
}

接着TLE了……(~~丝毫没有卡常的欲望~~)

因为递归常数巨大，并且很容易爆空间。(但实际上是因为写了long \ double炸了)

下面是AC的递归代码，~~以后不要写long double 了~~

    #include<iostream>
    #include<cstdio>
    #include<cmath>
    using namespace std;
    const double pi = 3.141592653589793;
    const int SIZE = 4e6 + 10;
    struct complex { //复数 
        double real, imag;
        complex() { real = imag = (double)0;}
            //重载复数四则运算 
        complex operator + (const complex x) const {
            complex temp;
            temp.real = real + x.real;
            temp.imag = imag + x.imag;
            return temp; 
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.real = real - x.real;
            temp.imag = imag - x.imag;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.real = real * x.real - imag * x.imag;
            temp.imag = real * x.imag + imag * x.real;
            return temp; 
        }

        complex operator *= (const complex &x) {
            *this = *this * x;
            return *this;
        }

        complex operator / (const complex x) const {
            double div = x.real * x.real + x.imag * x.imag;
            complex temp;
            temp.real = (real * x.real + imag * x.imag) / div;
            temp.imag = (imag * x.real - real * x.imag) / div;
            return temp; 
        } 

        void conj() { imag = -imag;}
    }a[SIZE], b[SIZE];

    complex omega(int n, int k) {
        complex temp;
        temp.real = cos(2 * pi * k / n);
        temp.imag = sin(2 * pi * k / n);
        return temp;
    }   
    int n, m, p;        
    void fft(complex *f, int len, int inv) {
        if(len == 1) return;
        int mid = len >>1;
        complex fl[mid + 10], fr[mid + 10];
        for(int i = 0; i < mid; i++)
            fl[i] = f[i << 1], fr[i] = f[i << 1 | 1];
        fft(fl, mid, inv);
        fft(fr, mid, inv);
        complex temp, buf; 
        temp.real = cos(2 * pi / len);
         temp.imag = sin(2 * pi / len);
        buf.real = 1, buf.imag = 0;
        if(inv) temp.conj();
        for(int i = 0; i < mid; i++, buf *= temp) {
            complex s = buf * fr[i];
            f[i] = fl[i] + s;
            f[i + mid] = fl[i] - s;
        }
    }

    int main() {

        scanf("%d %d", &n, &m);
        for(p = 1; p < n + m + 1; p <<= 1);
        for(int i = 0; i <= n; i++) 
            scanf("%lf", &a[i].real);
        for(int i = 0; i <= m; i++) 
            scanf("%lf", &b[i].real);
        fft(a, p, 0);
        fft(b, p, 0);
        for(int i = 0; i < p; i++) 
            a[i] *= b[i];
        fft(a, p, 1);
        for(int i = 0; i <= n + m; i++)
            printf("%d ", (int)(a[i].real / p + 0.5));
        return 0;
    }

实际上跑的很快的递归评测记录。

那有没有更优化的办法呢？这是肯定的。

蝴蝶变换

先看一张图

这张图是什么？这不是我们分治的顺序吗？

我在一开始和结束都打上了二进制，你有没有什么发现？

是的，每一个数的二进制都反过来了。

然后我们的任务就转换成求二进制反序：

    rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << l)//l是二进制位数
    //rev[i]是i翻转得到的

结论牢记即可。

还记得DP吗？（~~貌似没什么关系~~）

我们用循环实现FFT，也就像DP一样定义状态、阶段决策：

第一层循环枚举变换区间大小（阶段），第二层循环枚举开头（决策），第三层处理区间。

    #include<iostream>
    #include<cstdio>
    #include<cmath>
    using namespace std;
    const double pi = 3.141592653589793;
    const int SIZE = 4e6 + 10;
    struct complex { //复数 
        double real, imag;
        complex() { real = imag = (double)0;}
            //重载复数四则运算 
        complex operator + (const complex x) const {
            complex temp;
            temp.real = real + x.real;
            temp.imag = imag + x.imag;
            return temp; 
        }

        complex operator += (const complex &x) {
            *this = *this + x;
            return *this;
        }

        complex  operator - (const complex x) const {
            complex temp;
            temp.real = real - x.real;
            temp.imag = imag - x.imag;
            return temp; 
        }

        complex operator * (const complex x) const {
            complex temp;
            temp.real = real * x.real - imag * x.imag;
            temp.imag = real * x.imag + imag * x.real;
            return temp; 
        }

        complex operator *= (const complex &x) {
            *this = *this * x;
            return *this;
        }

        complex operator / (const complex x) const {
            double div = x.real * x.real + x.imag * x.imag;
            complex temp;
            temp.real = (real * x.real + imag * x.imag) / div;
            temp.imag = (imag * x.real - real * x.imag) / div;
            return temp; 
        } 
    }a[SIZE], b[SIZE];
    int n, m, p, L = -1;
    int rev[SIZE];

    void fft(complex *f, int inv) {
        for(int i = 0 ; i < p ; i++) 
            if(i < rev[i]) 
                swap(f[i], f[rev[i]]);
        for(int len = 2; len <= p; len <<= 1) {
            int length = len >> 1;
            complex temp;
            temp.real = cos(pi / length);
            temp.imag = sin(pi / length);
            if(inv) temp.imag = -temp.imag;
            for(int l = 0; l < p; l += len) {
                complex buf; buf.real = 1;
                for(int k = l; k < l + length; k++) {
                    complex t = buf * f[k + length];
                    f[k + length] = f[k] - t;
                    f[k] += t;
                    buf *= temp;
                }
            }
        }
    }

    int main() {
        scanf("%d %d", &n, &m);
        for(int i = 0; i <= n; i++)
            scanf("%lf", &a[i].real);
        for(int i = 0; i <= m; i++)
            scanf("%lf", &b[i].real);
        for(p = 1; p < n + m + 1; p <<= 1, L++);
        //l是位数，因为会多算一位，一开始初值直接给-1
        for(int i = 1; i < p; i++)
            rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L);
        fft(a, 0), fft(b, 0);
        for(int i = 0; i < p; i++)
            a[i] *= b[i];
        fft(a, 1);
        for(int i = 0; i <= n + m; i++)
            printf("%d ", (int)(a[i].real / p + 0.5));
        return 0;
    }

评测记录

END

FFT——快速傅里叶变换

FFT——快速傅里叶变换

多项式的表达方法

神仙问题：tql but 这有用吗

算法流程：

2. F(ω^{k+n/2}_{n})=FL(ω^{k}_{n/2})-ω^{k}_{n}FR(ω^{k}_{n/2})

蝴蝶变换