FFT剖析

快速傅里叶变换 (fast Fourier transform)

xn={x0,x1,…xn-1} (num:N)

旋转因子系数：

d=2pik/N

旋转因子

wk(n)=(cos(dn)+isin(dn)) n=[0,N-1]
y(k)= sum(x(n)wk(n),0,N-1)
y(k)={y(0),y(1),…y(N-1)} 傅里叶级数
x(n)wk(n)的级数是：
1.d=2pik/N 这个系数决定转动圈数，不懂看第二个再说
2.y(k)=x(0)wk(0)+x(1)wk(2)+…+x(n-1)wk(n-1)
旋转因子：cos(t)+isin(t)
t=2pikn/N
T=2pi就是一圈，那么可以得出步进量：2pi*k/N, 圈数:k
在一圈中：wk(n)的单位圆上t<=pi/2，是他的展开区 t<=pi是他的对称区

#include <bits/stdc++.h>
#define fu(x , y , z) for(int x = y ; x <= z ; x ++)
#define fd(x , y , z) for(int x = y ; x >= z ; x --)
#define LL long long
using namespace std;
const int N = 4e6 + 5;
const double pi = acos (-1.0);
struct node {
    double x , y;
} a[N] , b[N];
int n , m , len = 1 , r[N] , l;
node operator + (node a, node b) { return (node){a.x + b.x , a.y + b.y};}
node operator - (node a, node b) { return (node){a.x - b.x , a.y - b.y};}
node operator * (node a, node b) { return (node){a.x * b.x - a.y * b.y , a.x * b.y + a.y * b.x};}
int read () {
    int val = 0 , fu = 1;
    char ch = getchar ();
    while (ch < '0' || ch > '9') {
        if (ch == '-') fu = -1;
        ch = getchar ();
    }
    while (ch >= '0' && ch <= '9') {
        val = val * 10 + (ch - '0');
        ch = getchar ();
    }
    return val * fu;
}
void fft (node *A , int inv) {
    for (int i = 0 ; i < len ; i ++)
        if (i < r[i]) 
            swap (A[i] , A[r[i]]);
    for (int mid = 1 ; mid < len ; mid <<= 1) {  //mid={1,2,4,8,16,32,64,...}
        node wn = (node){cos (1.0 * pi / mid) , inv * sin (1.0 * pi / mid)};
        for (int R = mid << 1 , j = 0 ; j < len ; j += R) {
            node w = (node){1 , 0};
            for (int k = 0 ; k < mid ; k ++ , w = w * wn) {
                node x = A[j + k] , y = w * A[j + mid + k];
                A[j + k] = x + y;
                A[j + mid + k] = x - y;
            }
        }
    }
}
int main () {
    n = read () , m = read ();
    fu (i , 0 , n) a[i].x = read ();
    fu (i , 0 , m) b[i].x = read ();
    while (len <= n + m) len <<= 1 , l ++;
    for (int i = 0 ; i < len ; i ++)
        r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
    fft (a , 1);
    fft (b , 1);
    fu (i , 0 , len) 
        a[i] = a[i] * b[i];
    fft (a , -1);
    // for (int i = 0 ; i <= n + m ; i ++) cout << a[i].x << " " << a[i].y << "\n";
    fu (i , 0 , n + m) 
        printf ("%d " , (int)(a[i].x / len + 0.5));
    return 0;
}

计算2^n的逆排长度：2的6次方：64=0b100-0000；逆排长度0b11-1111

while (len <= n + m) len <<= 1 , l ++;
产生逆排：0,32，
for (int i = 0 ; i < len ; i ++)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
逆排交换数据
for (int i = 0 ; i < len ; i ++)
if (i < r[i])
swap (A[i] , A[r[i]]);

X(k)=sum(x(n)W(kn/N),0,N-1)
X(k)=sum(x(2r)W(k2r/N),0,N/2-1)+sum(x(2r+1)W(k(2r+1)/N),0,N/2-1)
Y1(k)=x1(k)+W(k/N)X2(K)
Y2(k)=x1(k)-W(k/N)X2(K)
例：N=16
偶 0, 4, 2, 6, 奇1, 5, 3, 7 x0,x1,…x7
k=2r
Y(k)=((x0+wkx1)+wk*(x2+wkx3))+wk((x4+wkx5)+wk(x6+wkx7))
k=2r+1
Y(k)=((x0-wkx1)-wk(x2-wkx3))-wk((x4-wkx5)-wk(x6-wkx7))
wk=w0 y(0)
wk=w1 y(1)
wk=w2 y(2)
N=2^L程序怎么写
1.逆排序
2.w(k)循环 0,1，N/2-1
x1={x0+wkx1,x2+wkx3,…xN-2+wkxN-1}
x2={x0-wkx1,x2-wkx3,…xN-2-wkxN-1}
3.L-1层循环
x1={x0+wkx1,x2+wkx3,…xN-2+wkxN-1}
x2={x0-wkx1,x2-wkx3,…xN-2-wk*xN-1}

#include <bits/stdc++.h>
#define fu(x , y , z) for(int x = y ; x <= z ; x ++)
using namespace std;
const int N = 4e6 + 5;
const double pi = acos (-1.0);
int n , m1 , m2 , rev[N];
complex<double> a[N] , b[N];
void fft (complex<double> *a , int type) {
    fu (i , 0 , n - 1) 
        if (i < rev[i]) swap (a[i] , a[rev[i]]);
    for (int j = 1 ; j < n ; j <<= 1) {
        complex<double> W(cos (pi / j) , sin (pi / j) * type);
        for (int k = 0 ; k < n ; k += (j << 1)) {
            complex<double> w(1.0 , 0.0);
            fu (i , 0 , j - 1) {
                complex<double> ye , yo;
                ye = a[i + k] , yo = a[i + j + k] * w;
                a[i + k] = ye + yo;
                a[i + k] = ye + yo;
                a[i + j + k] = ye - yo;
                w *= W;
            }
        }
    }
}
int main () {
    scanf ("%d%d" , &m1 , &m2);
    fu (i , 0 , m1) cin >> a[i];
    fu (i , 0 , m2) cin >> b[i];
    n = m1 + m2;
    int t = 0;
    while (n >= (1 << t)) 
        t ++;
    n = (1 << t);
    fu (i , 0 , n - 1) 
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (n >> 1) : 0);
    fft (a , 1) , fft (b , 1);
    fu (i , 0 , n) 
        a[i] *= b[i];
    fft (a , -1);
    fu (i , 0 , m1 + m2) 
        printf ("%d " , (int)(a[i].real() / (double)n + 0.5));
    return 0;
}

f0=(a0,a1,…an-2)
f1=(a1,a3,…an-1)
f(wnk)=f0(wnk/2)+wnkf1(wkn/2)
f(wn(k+n/2)=f0(wnk/2)-wnkf1(wkn/2)

#include<cstdio>
 2 #include<iostream>
 3 #include<cmath>
 4 #include<cstring>
 5 #include<algorithm>
 6 #include<cstdlib>
 7 using namespace std;
 8 const int mod=1e9+7;
 9 const double pi=acos(-1);
10 struct cn
11 {
12     double x,y;
13     cn (double x=0,double y=0):x(x),y(y) {}
14 }a[300005],b[300005],c[300005];
15 cn operator + (const cn &a,const cn &b) {return cn(a.x+b.x,a.y+b.y);}
16 cn operator - (const cn &a,const cn &b) {return cn(a.x-b.x,a.y-b.y);}
17 cn operator * (const cn &a,const cn &b) {return cn(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
18 void fft(cn a[],int n,int l,int f)                                   
19 {                                                                                
20     int rev[n+5];
21     rev[0]=0;
22     for (int i=1; i<n; i++){
23         rev[i]=(rev[i>>1]>>1)|((i&1)<<l-1);
24         if (i<rev[i]) swap(a[i],a[rev[i]]);
25     }
26     for (int i=1; i<n; i<<=1){ //L层循环
27         cn wi(cos(pi/i),f*sin(pi/i));
28         for (int j=0; j<n; j+=i*2){
29             cn w(1,0);
30             for (int k=0; k<i; k++){
31                 cn x=a[j+k],y=w*a[j+k+i];   //x1=
32                 a[j+k]=x+y;       // x10=a[0]+w0*a[1]  x11=a[2]+w0*a[3]    
33                 a[j+k+i]=x-y;     // x20=a[0]-w0*a[1]  x21=a[2]-w0*a[3]    
34                 w=w*wi;
35             }
36         }
37     }
38     if (f==-1)
39         for (int i=0; i<n; i++){
40             a[i].x/=n; a[i].y/=n;
41         }
42 }
43 int main()
44 {
45     int n,m;
46     scanf("%d%d",&n,&m); n++; m++;
47     for (int i=0; i<n; i++) scanf("%lf",&a[i].x);
48     for (int i=0; i<m; i++) scanf("%lf",&b[i].x);
49     int l=0,N=1;
50     while (N<n+m-1) N<<=1,l++;
51     fft(a,N,l,1);
52     fft(b,N,l,1);
53     for (int i=0; i<N; i++) c[i]=a[i]*b[i];
54     fft(c,N,l,-1);
55     for (int i=0; i<n+m-1; i++) printf("%d ",(int)(c[i].x+0.5));
56     return 0;
57 }

第一层循环控制级数，共M = l o g N M=logNM=logN级，从L = 1 , . . . , M L = 1 , . . . , ML=1,…,M；
第二层循环控制每层的旋转因子，每层有B=2^L-1个不同的旋转因子；
第三层循环控制每个旋转因子对应的蝶形运算。
每个旋转因子会被使用2^M-L次；
每个蝶形运算的两个输入数据相隔B=2^L-1个点，输出也是相隔B个点（蝶形：交叉平行）；
而同一旋转因子对应的两个相邻的蝶形运算相隔D = 2^L 个点；
d=2^i
x(k)=x(k)+w(v/N)*x(k+d) k=(0,N/d); v=(0,i) num=N/d
x(k+d)=x(k)-w(v/N)*x(k+d)

for (l = 1; (f = f / 2) != 1; l++); 	// 蝶形级数
        for (m = 1; m <= l; m++) 		// 控制蝶形结级数
        {
            le = 2 << (m - 1);         	//le蝶形结距离，即第m级蝶形的蝶形结相距le点
            lei = le / 2; 				// 同一蝶形结中参加运算的两点的距离
            u.real = 1.0; 				// u为蝶形结运算系数，初始值为1
            u.imag = 0.0;
            w.real = cos(PI / lei); 	// w为系数商，即当前系数与前一个系数的商
            w.imag = -sin(PI / lei);
            for (j = 0; j <= lei - 1; j++) // 控制计算不同种蝶形结，即计算系数不同的蝶形结
            {
                for (i = j; i <= FFT_N - 1; i = i + le) // 控制同一蝶形结运算，即计算系数相同蝶形结
                {
                    ip = i + lei; 		// i，ip分别表示参加蝶形运算的两个节点
                    t = plural(dat[ip], u); // 蝶形运算，详见公式
                    dat[ip].real = dat[i].real - t.real;
                    dat[ip].imag = dat[i].imag - t.imag;
                    dat[i].real = dat[i].real + t.real;
                    dat[i].imag = dat[i].imag + t.imag;
                }
                u = plural(u, w); 		 // 改变系数，进行下一个蝶形运算
            }
        }

通过变址运算实现F F T FFTFFT算法，变址运算通过雷德算法实现

for (i = 0; i < (FFT_N-1); i++) {
        if (i < j) {    // 如果i<j,即进行变址
            t = dat[j];
            dat[j] = dat[i];
            dat[i] = t;
        }
        k = nv2;
        // 求j的下一个倒位序
        while (k <= j) // 如果k<=j,表示j的最高位为1
        {
            j = j - k; // 把最高位变成0
            k = k / 2; // k/2，比较次高位，依次类推，逐个比较，直到某个位为0
        }
        j = j + k; // 把0改为1
    }