zoukankan      html  css  js  c++  java
  • FFT快速傅里叶变换 [模板]

    • FFTFFT的作用

    O(logN)O(logN) 计算两个多项式的卷积 .

    • FFTFFT的前置

    1. 单位根

    .,.单位根为在单位圆上的点.\ 横坐标为实部, 纵坐标为虚部的复数 .


    如图, 图中落在 xx 正半轴的那个红点称为 wn0w_{n}^0, 被称为 nn 次单位根, 其余的点从 wn0w_n^0 逆时针顺序分别称为 wn1,wn2...wn7w_n^1, w_n^2...w_n^7 .

    从图中可以 感性 得到 单位根22 个性质 downarrow

    • 1:性质1: w2n2k=wnkw_{2n}^{2k}=w_n^k

    • 2:性质2: wnk+n2=wnkw_n^{k+frac{n}{2}} = -w_n^k


    下面给出 理性 证明, 首先要知道

    :wnk=cos(k2πn)+isin(k2πn)欧拉公式:w_n^k=cos(k*frac{2 pi}{n}) + i*sin(k*frac{2 pi}{n})

    其中 ii1sqrt{-1} .


    • 1:性质1: w2n2k=wnkw_{2n}^{2k}=w_n^k
      :证:
      w2n2k=cos(2k2π2n)+isin(2k2π2n)=wnkw_{2n}^{2k} = cos(2k*frac{2 pi}{2n}) + i*sin(2k*frac{2 pi}{2n})=w_n^k

    • 2:性质2: wnk+n2=wnkw_n^{k+frac{n}{2}} = -w_n^k
      :证:
      wnk+n2=cos[(k+n2)2πn]+isin[(k+n2)2πn)]=wnkw_n^{k+frac{n}{2}}\ =cos[(k+frac{n}{2})*frac{2pi}{n}]+i*sin[(k+frac{n}{2})*frac{2pi}{n})]\ =-w_n^k

    • FFTFFT的内容

    • 有多项式 A(x)=a0+a1x+a2x2+​+an1xn1A(x)=a_0 + a_1x+a_2x^2+dotsi+a_{n-1}x^{n-1}

    nn+1.ecause 一个n次多项式可以被n+1个点唯一确定.
    wn0​wnn1. herefore 将 w_n^0 dotsi w_n^{n-1} 代入可以确定该多项式.

    这里默认 nn22 的整数次幂, 所以 n1n-1 为奇数 .

    • 对于上方多项式 A(x)A(x), 将其 系数 按奇偶分开, 得
      A(x)=(a0+a2x2+​+an2xn2)+(a1x+a3x3+​+an1xn1)A(x) = (a_0+a_2*x^2+ dotsi + a_{n-2}*x^{n-2})+(a_1*x+a_3*x^3+ dotsi + a_{n-1}*x^{n-1})

      A1(x)=a0+a2x+a4x2+​+an2xn21 A2(x)=a1+a3x+a5x2+​+an1xn21A_1(x)=a_0+a_2x+a_4x^2+ dotsi+a_{n-2}x^{frac{n}{2}-1}\ \ A_2(x)=a_1+a_3x+a_5x^2+ dotsi+a_{n-1}x^{frac{n}{2}-1}

      A(x)=A1(x2)+xA2(x2)A(x)=A_1(x^2)+x*A_2(x^2)

     wnk    (k<n2) 代入 w_n^k (k < frac{n}{2}) 得
    A(wnk)=A1(wn2k)+wnkA2(wn2k)=A1(wn2k)+wnkA2(wn2k)A(w_n^k) = A_1(w_n^{2k})+w_n^kA_2(w_n^{2k})=A_1(w_{frac{n}{2}}^k)+w_n^kA_2(w_{frac{n}{2}}^k)

     wnk+n2 代入 w_n^{k+frac{n}{2}} 得
    A(wnk+n2) =A1(wn2k+n)+wnk+n2A2(wn2k+n) =A1(wn2kwnn)wnkA2(wn2kwnn) =A1(wn2k)wnkA2(wn2k) =A1(wn2k)wnkA2(wn2k)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}) \ \ =A_1(w_{frac{n}{2}}^k)-w_n^kA_2(w_{frac{n}{2}}^k)

    :结论:
    • A(wnk)=A1(wn2k)+wnkA2(wn2k)A(w_n^k) = A_1(w_{frac{n}{2}}^k)+w_n^kA_2(w_{frac{n}{2}}^k)

    • A(wnk+n2)=A1(wn2k)wnkA2(wn2k)A(w_n^{k+frac{n}{2}})=A_1(w_{frac{n}{2}}^k)-w_n^kA_2(w_{frac{n}{2}}^k)

    可以观察到第二个式子与第一个式子的唯一区别就是正负号 .

    根据上方式子分治处理出每个子项, 再向上合并, 时间复杂度 O(NlogN)O(NlogN) .


    :总结:

    FFTFFT是先将多项式转换为 点值表示, 然后进行分治处理, 再转变为多项式的 O(nlogn)O(nlogn) 算法 .

    推导过程:

    1. 把系数按奇偶分开
    2. 代入 单位根 .
    3. 得到 “分治向上递推式” .

    • FFTFFT的实现

    • 蝴蝶变换

    借用 这位大佬 的一个图,

    在这里插入图片描述

    所谓蝴蝶变换, 就是将原序列进行二进制翻转的变换 .

    rev[i]=(rev[i>>1]>>1)((i&1)<<bit_num1)rev[i] = (rev[i>>1]>>1) | ((i&1) << bit\_num-1)

    详解点击 这里 .


    下面是 模板 代码 .

    #include<bits/stdc++.h>
    
    #define reg register
    
    const int maxn = 1e6 + 5;
    const double Pi = acos(-1);
    
    struct complex{
            double x, y;
            complex(double x = 0, double y = 0):x(x), y(y) {}
    } A[maxn<<2], B[maxn<<2];
    
    complex operator + (complex a, complex b){ return complex(a.x+b.x, a.y+b.y); }
    complex operator - (complex a, complex b){ return complex(a.x-b.x, a.y-b.y); }
    complex operator * (complex a, complex b){ return complex(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x); }
    
    int N, M;
    int rev[maxn<<2];
    
    void FFT(complex *F, short Opt){
            for(int i = 0; i < N; i ++)
                    if(i < rev[i]) std::swap(F[i], F[rev[i]]);
            for(int p = 2; p <= N; p <<= 1){
                    int len = p >> 1;
                    complex tmp(cos(Pi/len), Opt*sin(Pi/len));
                    for(int i = 0; i < N; i += p){
                            complex Buf(1, 0);
                            for(reg int k = i; k < i+len; k ++){
                                    complex Temp = Buf * F[k+len];
                                    F[k+len] = F[k] - Temp;
                                    F[k] = F[k] + Temp;
                                    Buf = Buf * tmp;
                            }
                    }
            }
    }
    
    int main(){
            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();
            int bit_num = 0;
            for(M += N, N = 1; N <= M; N <<= 1) bit_num ++;
            for(int i = 0; i < N; i ++) rev[i] = (rev[i>>1]>>1)|((i&1)?N>>1:0);
            // for(int i = 0; i < N; i ++) rev[i] = (rev[i>>1]>>1) | ((i&1) << bit_num-1);
            FFT(A, 1), FFT(B, 1);
            for(int i = 0; i < N; i ++) B[i] = A[i] * B[i];
            FFT(B, -1);
            for(reg int i = 0; i <= M; i ++) printf("%.0lf ", fabs(B[i].x/N));
            return 0;
    }
    
    • FFTFFT的应用

    万径人踪灭

    礼物


    • FFT &gt;NTTFFT 魔改-&gt; NTT

    1. 将单位根 wnw_n 替换为原根 gmod1leng^{frac{mod-1}{len}}
    2. 最后乘的是 T_lenT\_len 的逆元 .

    其中 gg 在 模数 998244353998244353下取 33 .

    void FFT(complex *F, int opt){
            for(reg int i = 0; i < FFT_Len; i ++)
                    if(i < rev[i]) std::swap(F[i], F[rev[i]]);
            for(reg int p = 2; p <= FFT_Len; p <<= 1){
                    int half = p >> 1;
                    complex t = complex(cos(Pi/half), opt*sin(Pi/half));
                    for(reg int i = 0; i < FFT_Len; i += p){
                            complex buf = complex(1, 0);
                            for(reg int k = i; k < i+half; k ++){
                                    complex Tmp = buf * F[k + half];
                                    F[k + half] = F[k] - Tmp;
                                    F[k] = F[k] + Tmp;
                                    buf = buf * t;
                            }
                    }
            }
    }
    

    void NTT(int *f, int opt){
            for(reg int i = 0; i < T_len; i ++)
                    if(i < rev[i]) std::swap(f[i], f[rev[i]]);
            for(reg int p = 2; p <= T_len; p <<= 1){
                    int half = p >> 1;
                    int wn = Ksm(3, (mod-1)/p);
                    if(opt == -1) wn = Ksm(wn, mod-2);
                    for(reg int i = 0; i < T_len ; i += p){
                            int buf = 1;
                            for(reg int k = i; k < i+half; k ++){
                                    int Tmp_1 = 1ll*buf*f[k+half] % mod;
                                    f[k+half] = (f[k] - Tmp_1 + mod) % mod;
                                    f[k] = (f[k] + Tmp_1) % mod;
                                    buf = 1ll*buf*wn % mod;
                            }
                    }
            }
    }
    
  • 相关阅读:
    hdu 1124 OR toj 1065 简单数论
    this和判断的位置对赋值的要求
    快捷键操作
    常量池和堆的区别
    toString的用法
    使用泛型解决之前的问题
    不使用泛型会运行时会出现的问题
    集合图型
    代码块运行的优先级
    遍历集合的Iterator删除其中的元素
  • 原文地址:https://www.cnblogs.com/zbr162/p/11822533.html
Copyright © 2011-2022 走看看