• 【数学】快速傅里叶变换(FFT)


    快速傅里叶变换(FFT)

    FFT 是之前学的,现在过了比较久的时间,终于打算在回顾的时候系统地整理一篇笔记,有写错的部分请指出来啊 qwq。

    卷积

    卷积、旋积或褶积(英语:Convolution)是通过两个函数 (f)(g)​​ 生成第三个函数的一种数学算子

    定义

    (f,g)​ 在 (R1)​ 上可积,那么 (h(x) = int_{-∞}^∞f( au)g(x- au)d au) 称为 (f)(g)​ 的卷积。

    对于整系数多项式域,(n-1) 次多项式 (A,B) 的卷积 (h(x) = A(x)B(x) = sum_{i=0}^{2(n-1)}sum_{j=0}^{min(i,n-1)} a_{j}b_{i-j}x^i)

    系数表示法

    即用多项式各项系数来刻画这个多项式,例如 (n-1) 次多项式就可以写成这样:(A(x) = a_0 + a_1x+...+a_{n-1}x^{n-1})

    点值表示法

    我们知道,(n)​​​​ 个不同的点可以确定一个 (n-1)​​​​ 次的多项式,所以我们可以使用 (n)​​​​ 个(不同)点来刻画一个 (n-1)​​​​ 次多项式。

    这样做会有什么方便呢?

    例如 (f(x) = (x_0, f(x_0)),...(x_n,f(x_n)),g(x) = (x_0, g(x_0)),...(x_n, g(x_n))),那么它们的卷积 (h(x) = (x_0, f(x_0)g(x_0)),...(x_n, f(x_n)g(x_n)))

    这意味着在系数表示法中需要 (O(n^2)) 次的乘法运算在点值表示法中只需要 (O(n)) 次。

    系数表示法转点值表示法(DFT)

    下面考虑如何将 (n-1) 次多项式从系数表示法转为点值表示法。

    因为用普通的方法选取 (n)​​​​​ 个点然后将系数表示法转为点值表示法的复杂度为 (O(n^2))​​​​​(因为需要选 (n)​​​​​ 个点,然后对于每个点 (x) 需要计算共 (n) 项的结果),我们考虑如何优化这一步。

    注意到满足 (w^n=1)​​ 的单位根 (w)​ 有 (n)​ 个,故从这里入手。

    我们记方程 (w^n = 1) 的第 (k) 个单位根为 (w_n^k)​。

    方便起见,设 (n)​ 为 (2)​ 的幂(就算不是也可以看作是高次项的系数为 (0)​)。

    (A(x)) 按照次数的奇偶性分别分成两组 (F(x),G(x)),并表示为 (A(x) = F(x^2) + xG(x^2))

    例如 (A(x) = a_0+a_1x+a_2x^2+a_3x^3)​,那么 (F(x) = a_0+a_2x,G(x)=a_1+a_3x)​。

    (x=w_n^k)​ 代入 (A(x))​,由复数的性质,

    (A(w_n^k) = F(w_{frac{n}{2}}^k) + w_n^k G(w_{frac{n}{2}}^k))​ ,类似地 (A(w_{n}^{k+frac{n}{2}}) = F(w_{frac{n}{2}}^k) - w_n^k G(w_{frac{n}{2}}^k))​。​

    推导:

    (A(w_{n}^{k+frac{n}{2}}) = F(w_n^{2k+n}) + w_n^{k+frac{n}{2}}G(w_n^{2k+n}) \= F(w_{frac{n}{2}}^k) + w_n^{k+frac{n}{2}} G(w_{frac{n}{2}}^k) \= F(w_{frac{n}{2}}^k) - w_n^k G(w_{frac{n}{2}}^k))​​​

    可以发现对于两个相应的单位根 (w_n^k,w_n^{frac{n}{2}+k})​,可以用对应的 (F,G)​ 算出(可以递归地实现这个过程),而且计算的范围折半,所以一共需要计算 (O(logN))​ 层,每一层执行 (O(n))​ 次运算,所以复杂度为 (O(NlogN))​。

    点值表示法转系数表示法(IDFT)

    下面考虑如何将 (n-1)​ 次多项式从点值表示法转为系数表示法。

    因为对于每个点值 (y_i = sum_{k=0}^{n-1}w_n^{ki}),其中 (iin[0, n-1])​​​​,我们可以写出等式:

    [left[ egin{matrix} 1 & 1 & 1 & cdots & 1\ 1 & w_n^1 & w_n^2 & cdots & w_n^{n-1}\ 1 & w_n^2 & w_n^4 & cdots & w_n^{2(n-1)}\ vdots & vdots & vdots & ddots & vdots \ 1 & w_n^{n-1} & w_n^{2(n-1)} & cdots & w_n^{(n-1)(n-1)} end{matrix} ight] left[ egin{matrix} a_0\ a_1\ a_2\ vdots\ a_{n-1} end{matrix} ight] = left[ egin{matrix} y_0\ y_1\ y_2\ vdots\ y_{n-1} end{matrix} ight] ]

    现在我们已经有向量 (y)​​ 了(就是右式),因此,如果要得到向量 (a)​​,只需要两边乘上 (w)​​​ 矩阵的即可。

    这里的 (w)​​ 矩阵正是著名的范德蒙矩阵,它的逆正好是每一项都取倒数,然后除以 (n)​​。

    因此有 (a_i = frac{1}{n}sum_{k=0}^{n-1}w_n^{-ki}),其中 (iin[0, n-1])​。

    有没有发现 (a_i,y_i)​ 的形式非常接近?据此,我们可以在实现的时候在同一个函数中写出逆变换和正变换,然后在得到的结果 res 中除以 (n) 就可以了。(参照下面的代码)

    至此,FFT 的基本原理讲述完毕,下面是优化。

    位逆序置换

    按照上文的讲述,如果不看下面的代码,那么编写出来的是递归版本,但是这个版本的常数太大了,因此运行起来的效果不好,故使用位逆序置换来降低常数。

    我们看看递归过程是什么样的,以 (n=8) 为例:

    [{x_0, x_1, x_2,x_3,x_4,x_5,x_6,x_7}\ {x_0, x_2, x_4,x_6},{x_1,x_3,x_5,x_7}\ {x_0, x_4}, {x_2,x_6},{x_1,x_5},{x_3,x_7}\ {x_0},{x_4},{x_2},{x_6},{x_1},{x_5},{x_3},{x_7}\ ]

    这里就有一个非常神奇的规律:在最后一行中,原下标所对应的二进制数翻转正好是在最后一行的序数。例如 (x_6) 的下标是 (6 = 110_{(2)}),那么它的序数正好是 (011_{(2)} = 3)

    据此,可以处理出 rev 数组,它记录的正是最后一行所有元素对应的下标。

    简单地说,递归形式是自上而下地做 FFT,而利用位逆序置换我们可以自下而上地做 FFT,它们在实际运行中有着常数上的区别。

    模板题及代码

    https://www.luogu.com.cn/problem/P3803

    给定一个 (n) 次多项式 (F(x)),和一个 (m) 次多项式 (G(x))

    请求出 (F(x))(G(x)) 的卷积。

    #include<bits/stdc++.h>
    using namespace std;
    
    const int N=3e5+5;
    const double pi=acos(-1);
    
    int n, m;
    
    // 复数类
    struct Complex{
    	double x, y;
    	Complex operator + (const Complex &o)const { return {x+o.x, y+o.y}; }
    	Complex operator - (const Complex &o)const { return {x-o.x, y-o.y}; }
    	Complex operator * (const Complex &o)const { return {x*o.x-y*o.y, x*o.y+y*o.x}; }
    };
    
    Complex a[N], b[N];
    int res[N];
    
    int rev[N], bit, tot;
    
    void fft(Complex a[], int inv){ // inv 指示正变换、逆变换。
    	for(int i=0; i<tot; i++) if(i<rev[i]) swap(a[i], a[rev[i]]);
    	
    	for(int mid=1; mid<tot; mid<<=1){
    		auto w1=Complex({cos(pi/mid), inv*sin(pi/mid)});
    		for(int i=0; i<tot; i+=mid*2){
    			auto wk=Complex({1, 0});
    			for(int j=0; j<mid; j++, wk=wk*w1){
    				auto x=a[i+j], y=wk*a[i+j+mid];
    				a[i+j]=x+y, a[i+j+mid]=x-y;
    			}
    		}
    	}
    }
    
    int main(){
    	cin>>n>>m;
    	for(int i=0; i<=n; i++) cin>>a[i].x;
    	for(int i=0; i<=m; i++) cin>>b[i].x;
    	
    	while((1<<bit)<n+m+1) bit++; // 结果次数分布在 [0, n+m] 内,一共有 n+m+1 位。
    	
    	tot=1<<bit; // 得到上文所说的 n
    	
    	for(int i=0; i<tot; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    	
    	fft(a, 1), fft(b, 1); // 正变换 DFT
    	
    	for(int i=0; i<tot; i++) a[i]=a[i]*b[i];
    	
    	fft(a, -1); // 逆变换 IDFT
    	
    	for(int i=0; i<=n+m; i++) res[i]=(int)(a[i].x/tot+0.5), printf("%d ", res[i]);
    	
    	return 0;
    }
    
  • 相关阅读:
    Python基础篇【第十一篇】:正则表达式之计算器
    Python基础篇【第十篇】:正则表达式
    Python基础篇【第九篇】:剖析装饰器
    Python基础篇【第八篇】:剖析递归函数
    Python基础篇【第七篇】:文件操作
    Python基础篇【第六篇】:函数补充
    Python基础篇【第五篇】:sed和函数
    TCP/IP详解之:Ping程序、Traceroute程序
    TCP/IP详解之:ICMP协议
    TCP/IP详解之:ARP协议 和 RARP协议
  • 原文地址:https://www.cnblogs.com/Tenshi/p/15434004.html
Copyright © 2020-2023  润新知