《快速傅里叶变换》

前言:写得比较简陋~以后有时间一定来完善。

快速傅里叶变换,又叫FFT。

可以在O(n)内处理多项式乘法的算法。

一个n - 1多项式g(x) 可以用系数表示为$\sum_{i = 0}^{n - 1} ai * x^{i}$

对于我们g(x)多项式中的每一项,我们可以看成有一个函数f(x)。

将每一项的x,x^1 ...都可以得到多项式中的对应的每一项。

那么每一项就可以看成一个点,在f(x)函数的图像上。

那么这n个点就可以确定唯一的一个f(x)函数。

这就是转化为点值法。

假设两个多项式点值法表示为:

F(x) = ((x0,f(x0)),(x1,f(x1)...)
G(x) = ((x0,g(x0)),(x1,g(x1)...)

那么他们的卷积H[x] = (x0,f(x0) * g(x0),....)

所以如果我们知道两个多项式的点值表示,那么就可以O(n)求出他们的卷积的点值表示。

然后引入单位根wn和欧拉定理:这里不赘述,可以自习百度。

然后根据推演:
对于wn^k 的k  < n / 2的情况 = A + wn^k *B

对于wn^k 的k > n / 2的情况 = A - wn^k * B

两者只有一个常数系数的不同。

所以我们在算上面一项的时候可以O(1)算得下面一项的值。

所以我们可以折半分治它,借由蝴蝶操作。

最后我们算得两个多项式卷积的点值表示后,再有推论:ak = ck / n。

可得我们的系数表示 = 点值表示 / n。

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef pair<LL, int> pii;
const int N = 2e6 + 5;
const int M = 2000 + 5;
const LL Mod = 998244353;
#define pi acos(-1)
#define INF 1e9
#define dbg(ax) cout << "now this num is " << ax << endl;
namespace FASTIO {
inline LL read() {
    LL x = 0, f = 1;
    char c = getchar();

    while (c < '0' || c > '9') {
        if (c == '-')
            f = -1;

        c = getchar();
    }

    while (c >= '0' && c <= '9') {
        x = (x << 1) + (x << 3) + (c ^ 48);
        c = getchar();
    }

    return x * f;
}
}
using namespace FASTIO;

complex<double> a[N],b[N];
int limit = 1,n,m;
void FFT(int len,complex<double> *a,int id) {
    if(len == 1) return ;
    complex<double> a1[len >> 1],a2[len >> 1];
    for(int i = 0;i < len;i += 2) {//根据奇偶分类
        a1[i >> 1] = a[i],a2[i >> 1] = a[i + 1];
    }
    FFT(len >> 1,a1,id);
    FFT(len >> 1,a2,id);
    complex<double> wn(cos(2.0 * pi / len),id * sin(2.0 * pi / len)),w(1,0);
    //
    for(int i = 0;i < (len >> 1);++i,w *= wn) {
        a[i] = a1[i] + w * a2[i];
        a[i + (len >> 1)] = a1[i] - w * a2[i];
    }
}
int main() {
    n = read(),m = read();
    for(int i = 0;i <= n;++i) a[i] = read();
    for(int i = 0;i <= m;++i) b[i] = read();
    while(limit <= (n + m)) limit <<= 1;
    FFT(limit,a,1);
    FFT(limit,b,1);
    for(int i = 0;i <= limit;++i) a[i] = a[i] * b[i];
    FFT(limit,a,-1);
    for(int i = 0;i <= n + m;++i) printf("%d%c",(int)(a[i].real() / limit + 0.5),i == n + m ? '\n' : ' ');
    //system("pause");
    return 0;
}
View Code

递归法的缺点很明显,因为要进行奇偶分类,所以申请内存会浪费很多时间。

由此,我们有了迭代法:对分治层操作的序列进行分析后,我们wn^k在最终操作后的序列中的位置就是2进制的反转。

这样我们先处理出反转位置就可以快速优化FFT。

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef pair<LL, int> pii;
const int N = 4e6 + 5;
const int M = 2000 + 5;
const LL Mod = 998244353;
#define pi acos(-1)
#define INF 1e9
#define dbg(ax) cout << "now this num is " << ax << endl;
namespace FASTIO {
inline LL read() {
    LL x = 0, f = 1;
    char c = getchar();

    while (c < '0' || c > '9') {
        if (c == '-')
            f = -1;

        c = getchar();
    }

    while (c >= '0' && c <= '9') {
        x = (x << 1) + (x << 3) + (c ^ 48);
        c = getchar();
    }

    return x * f;
}
}
using namespace FASTIO;

int limit = 1,n,m,r[N];
complex<double> a[N],b[N];
void FFT(complex<double> *a,int id) {
    for(int i = 0;i < limit;++i) {
        if(i < r[i]) swap(a[i],a[r[i]]);
    }
    for(int i = 1;i < limit;i <<= 1) {
        complex<double> x(cos(pi / i),id * sin(pi / i));//这里变化为了pi / i
        for(int j = 0;j < limit;j += (i << 1)) {
            complex<double> y(1,0);
            for(int k = 0;k < i;++k,y *= x) {
                complex<double> p = a[j + k],q = y * a[j + k + i];
                a[j + k] = p + q;
                a[j + k + i] = p - q;
            }
        }
    }
}
int main() {
    n = read(),m = read();
    for(int i = 0;i <= n;++i) a[i] = read();
    for(int i = 0;i <= m;++i) b[i] = read();
    int L = 0;
    while(limit <= (n + m)) limit <<= 1,++L;
    for(int i = 0;i < limit;++i) {
        r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));
    }
    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%c",(int)(a[i].real() / limit + 0.5),i == n + m ? '\n' : ' ');
   // system("pause");
    return 0;
}
View Code

 

posted @ 2021-05-12 22:21  levill  阅读(87)  评论(0编辑  收藏  举报