• 《快速傅里叶变换》


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

    快速傅里叶变换,又叫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 ? '
    ' : ' ');
        //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 ? '
    ' : ' ');
       // system("pause");
        return 0;
    }
    View Code
  • 相关阅读:
    让服务器iis支持.apk文件下载的设置方法
    找不到匹配控制器
    Session_Start
    在SQLSERVER中创建聚集索引
    nvm 安装
    web访问命令行
    devmapper: Thin Pool has 162394 free data blocks which is less than minimum required 163840 free data blocks. Create more free space in thin pool or use dm.min_free_space option to change behavior
    服务器使用赛风
    java jpa 报错
    ssh 设置反向代理
  • 原文地址:https://www.cnblogs.com/zwjzwj/p/14762222.html
Copyright © 2020-2023  润新知