快速傅里叶变换,DFT,FFT

发布于 2020-11-07  83 次阅读


快速傅里叶变换FFT

大整数乘法计算的发展历程是缓慢的...在朴素的大整数乘法计算和多项式乘法中,两个多项式f(x), g(x)的任意两项之间都要相乘一次,时间复杂度为O(n^2)。后来出现了用分治法计算大整数乘法,把时间复杂度优化到O(n^{1.585}),然后经历了相当长的时间后才发明出了更快的计算方法--快速傅里叶变换(Fast Fourier Transform),时间复杂度为O(nlogn),为目前计算大数乘法最快的算法。

多项式的表示

一个多项式有两种表示方法,系数表示法和点值表示法。系数表示法就是我们最常用的用来表示一个多项式的方法,即多项式每个项的系数组成的集合。一个最高n次的多项式可以用n+1个系数来表示,也可以用n+1个横坐标不同的点来表示。n+1个不同的点可以唯一表示一个最高n次多项式。

给定两个多项式f(x), g(x),现计算其乘法h(x) = f(x)*g(x),若这两个多项式是用系数表示法表示的,那么需要O(n^2)的时间去计算;若是用点值表示法表示的,假设取到的点的横坐标的都是相同的,即(x1, f(x1)), (x2, f(x2)), (x3, f(x3))....(xn, f(xn))和(x1, g(x1)), (x2, g(x2)), (x3, g(x3))....(xn, g(xn)),那么仅需要O(n)的时间即可计算出h(x)的点值表示法(x1, f(x1)*g(x1)), (x2, f(x2)*g(x2)), (x3, f(x3)*g(x3))....(xn, f(xn)*g(xn))。可以利用这一过程来加速计算多项式乘法。

难点是如何将多项式从系数表示法转化为点值表示法。我们知道计算一个n次多项式的一个f(x)需要O(n)的时间,那转换成点值表示法岂不是也需要O(n^2)的时间。。快速傅里叶变换做的工作就是快速进行转换这个工作。

这里需要引入复数。

复数

如何快速从系数表示法转化为点值表示法呢?首先需要选择n+1个横坐标,这里选择x^n=1在复数意义上的n个单位根来表示,利用单位根的一些性质来快速完成转化!
一个复数a+bi可以像向量一样在二维坐标里表示为坐标(a, b),也可以像极坐标一样用一对(模角,幅角)来表示,模角是复数的长度,幅角是该复数和x正半轴的夹角。两复数相乘时满足:幅角相乘,模角相加。所以x^n=1的n个根其实就是把一个圆心在原点的单位元n等分,得到的n个点。
w_{n}^{k}表示n次复根的第k个单位根。w_{n}^{k}对应的坐标为(cos(\frac{k}{n}2\pi), sin(\frac{k}{n}2\pi)),也就是复数cos(\frac{k}{n}2\pi)+i*sin(\frac{k}{n}2\pi),根据欧拉公式也可表示成e^{i·\frac{k}{n}2\pi}。用欧拉公式可以很容易的解释为什么两个复数相乘结果是模角相乘,幅角相加。
三个重要性质:

w_{n}^{n} = w_{n}^{0} = 1
w_{n}^{k} = w_{2n}^{2k}
w_{n}^{k+\frac{n}{2}} = -w_{n}^{k}

这三个性质都很容易推出来。
第三个性质中,因为w_{n}^{\frac{n}{2}}刚好是-1,所以去掉之后还需要再乘以-1.

蝴蝶变换

如何把复数应用进去应该是FFT最难理解的部分了。

我们现在要做的是:对于一个n-1次多项式f(x),把n个单位复根带入,求得n个复值,如何快速的求出。

我们把f(x)的奇次项的系数和偶此项的系数分别表示成一个多项式,即
A_0(x) = a_0 + a_2x + a_4x^2 + ... + a_{n-2}x^{\frac{n}{2}-1}
A_1(x) = a_1 + a_3x + a_5x^2 + ... + a_{n-1}x^{\frac{n}{2}-1}
那么就有f(x) = A_0(x^2) + x * A_1(x^2).

我们再将n次单位复根带入,设k<=\frac{n}{2},得
f(w_{n}^{k}) = A_0(w_{n}^{2k}) + w_{n}^{k}*A_1(w_{n}^{2k}) = A_0(w_{\frac{n}{2}}^{k}) + w_{n}^{k}*A_1(w_{\frac{n}{2}}^{k})

f(w_{n}^{k+\frac{n}{2}}) = A_0(w_{n}^{2k+n}) + w_{n}^{k+\frac{n}{2}}*A_1(w_{n}^{2k+n}) = A_0(w_{\frac{n}{2}}^{k}) - w_{n}^{k}*A_1(w_{\frac{n}{2}}^{k})

f(x)是n-1次多项式,需要计算n次单位复根带入得到的值;A_0(x), A_1(x)\frac{2}{n}-1次多项式,已经计算出将\frac{2}{n}次单位复根带入得到的值,这样就可以在A_0(x), A_1(x)的基础上线性求出f(x)的n个复根!A_0(x), A_1(x)也可以用类似的求法递归求下去,共需要求logn层,所以总的时间复杂度是O(nlogn).
注意到每层循环中n都必须可以被2整除,所以第一层的n必须是2的幂,若所给的多项式的最高次不是2的幂,那么就补0.

快速傅里叶逆变换IFFT

得到目标多项式的点值表示法之后,我们还要把它转化成系数表示法,这一过程叫快速傅里叶逆变换。

把f(x)转化成点值表示法后,再把这些点值作为另一个多项式的系数,并把单位复根的倒数作为根代入,得到的值就是f(x)的系数,所以这里相当于是再做一次FFT即可。

代码:

const double PI = acos(-1.0);
void change(Complex y[], int len){
    for (int i=1, j=len/2, k; i<len-1; i++){
        if (i<j)    swap(y[i], y[j]);
        k = len/2;
        while (j>=k){
            j -= k;
            k /= 2;
        }
        if (j<k)    j += k;
    }
}
void FFT(Complex y[], int len, int on){       //快速傅里叶变换,y表示要带入的len个复数,on表示是求正变换还是逆变换。
    change(y, len);                      //直接化成递归的最底一层,即最低一层中各项系数的位置。
    for (int h=2; h<=len; h<<=1){        //h表示该层每个函数有h个元素,开始向上迭代
        Complex wn(cos(PI*2/h), sin(PI*2*on/h));      //wn为第1个单位根/单位根倒数
        for (int j=0; j<len; j+=h){                   //每层的第一个元素
            Complex w(1, 0);                          //w表示第几单位根,求完一次后,都变成下一个单位根。
            for (int k=j; k<j+h/2; k++){
                Complex u = y[k], t = w*y[k+h/2];
                y[k] = u+t, y[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if (on==-1){                        //若为逆变换,还需要再除以n才能得到原多项式的系数。
        for (int i=0; i<len; i++)   y[i].x /= len;
    }
}

复数的使用:

struct Complex{
    double x, y;                //实部和虚部
    Complex(double a = 0.0, double b = 0.0){
        x = a, y = b;
    }
    Complex operator + (const Complex &a) const {
        return Complex(x+a.x, y+a.y);
    }
    Complex operator - (const Complex &a) const {
        return Complex(x-a.x, y-a.y);
    }
    Complex operator * (const Complex &a) const {      //简单推出
        return Complex(x*a.x-y*a.y, x*a.y+y*a.x);
    }
};

题目
P1919 【模板】A*B Problem升级版
大数相乘不同于多项式乘法的地方是在计算完成后需要进位。

#include <bits/stdc++.h>

using namespace std;

const int maxn = 1e6+500000;
const double PI = acos(-1.0);
struct Complex{
    double x, y;
    Complex(double a = 0.0, double b = 0.0){
        x = a, y = b;
    }
    Complex operator + (const Complex &a) const {
        return Complex(x+a.x, y+a.y);
    }
    Complex operator - (const Complex &a) const {
        return Complex(x-a.x, y-a.y);
    }
    Complex operator * (const Complex &a) const {
        return Complex(x*a.x-y*a.y, x*a.y+y*a.x);
    }
};

char s1[maxn], s2[maxn];
struct Complex x1[maxn*2], x2[maxn*2];
int len;
int ans[maxn*2];

void change(Complex y[], int len){
    for (int i=1, j=len/2, k; i<len-1; i++){
        if (i<j)    swap(y[i], y[j]);
        k = len/2;
        while (j>=k){
            j -= k;
            k /= 2;
        }
        if (j<k)    j += k;
    }
}
void FFT(Complex y[], int len, int on){
    change(y, len);                      //直接化成递归的最底一层
    for (int h=2; h<=len; h<<=1){        //h表示该层每个函数有h个元素
        Complex wn(cos(PI*2/h), sin(PI*2*on/h));                    //
        for (int j=0; j<len; j+=h){        //每层的第一个元素
            Complex w(1, 0);                                   //w表示第几单位根
            for (int k=j; k<j+h/2; k++){
                Complex u = y[k], t = w*y[k+h/2];
                y[k] = u+t, y[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if (on==-1){
        for (int i=0; i<len; i++)   y[i].x /= len;
    }
}
int main(){

    while (~scanf ("%s %s", s1, s2)){
        int len1 = strlen(s1), len2 = strlen(s2);
        len = 1;
        while (len < len1+len2) len <<= 1;
        for (int i=0; i<len1; i++)  x1[i] = Complex(s1[len1-i-1]-'0', 0.0);
        for (int i=0; i<len2; i++)  x2[i] = Complex(s2[len2-i-1]-'0', 0.0);
        for (int i=len1; i<len; i++)    x1[i] = Complex(0, 0);
        for (int i=len2; i<len; i++)    x2[i] = Complex(0, 0);
        FFT(x1, len, 1);
        FFT(x2, len, 1);
        for (int i=0; i<len; i++)  x1[i] = x1[i] * x2[i];
        FFT(x1, len, -1);
        for (int i=0; i<len; i++)  ans[i] = int(x1[i].x + 0.5);
        for (int i=0; i<len; i++){
            ans[i+1] += ans[i] / 10;
            ans[i] %= 10;
        }
        while (ans[len-1]==0 && len>1) len--;
        for (int i=len-1; i>=0; i--)   putchar(ans[i]+'0');
        putchar('\n');
    }
    return 0;
}

P3803 【模板】多项式乘法(FFT)

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
const int maxn = 1e6+500000;
const double PI = acos(-1.0);
struct Complex{
    double x, y;
    Complex(double a = 0.0, double b = 0.0){
        x = a, y = b;
    }
    Complex operator + (const Complex &a) const {
        return Complex(x+a.x, y+a.y);
    }
    Complex operator - (const Complex &a) const {
        return Complex(x-a.x, y-a.y);
    }
    Complex operator * (const Complex &a) const {
        return Complex(x*a.x-y*a.y, x*a.y+y*a.x);
    }
};

struct Complex x1[maxn*2], x2[maxn*2];
int len;
ll ans[maxn*2];

void change(Complex y[], int len);
void FFT(Complex y[], int len, int on);
int main(){
    int len1, len2, tmp;
    while (~scanf ("%d %d", &len1, &len2)){
        len = 1, len1++, len2++;
        while (len < len1+len2) len <<= 1;
        for (int i=0; i<len1; i++){
            scanf ("%d", &tmp);
            x1[i] = Complex(tmp, 0.0);
        }
        for (int i=0; i<len2; i++){
            scanf ("%d", &tmp);
            x2[i] = Complex(tmp, 0.0);
        }
        for (int i=len1; i<len; i++)    x1[i] = Complex(0, 0);
        for (int i=len2; i<len; i++)    x2[i] = Complex(0, 0);
        FFT(x1, len, 1);
        FFT(x2, len, 1);
        for (int i=0; i<len; i++)  x1[i] = x1[i] * x2[i];
        FFT(x1, len, -1);
        for (int i=0; i<len; i++)  ans[i] = ll(x1[i].x + 0.5);

        while (ans[len-1]==0 && len>1) len--;
        printf ("%lld", ans[0]);
        len1--, len2--;
        for (int i=1; i<=len1+len2; i++)   printf (" %lld", ans[i]);
    }
    return 0;
}

void change(Complex y[], int len){
    for (int i=1, j=len/2, k; i<len-1; i++){
        if (i<j)    swap(y[i], y[j]);
        k = len/2;
        while (j>=k){
            j -= k;
            k /= 2;
        }
        if (j<k)    j += k;
    }
}
void FFT(Complex y[], int len, int on){
    change(y, len);                      //直接化成递归的最底一层
    for (int h=2; h<=len; h<<=1){        //h表示该层每个函数有h个元素
        Complex wn(cos(PI*2/h), sin(PI*2*on/h));                    //
        for (int j=0; j<len; j+=h){        //每层的第一个元素
            Complex w(1, 0);                                   //w表示第几单位根
            for (int k=j; k<j+h/2; k++){
                Complex u = y[k], t = w*y[k+h/2];
                y[k] = u+t, y[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if (on==-1){
        for (int i=0; i<len; i++)   y[i].x /= len;
    }
}

一沙一世界,一花一天堂。君掌盛无边,刹那成永恒。