快速傅里叶逆变换(IFFT)

本文作者为 JustinRochester。

目录地址

上一篇

下一篇


快速傅里叶逆变换(IFFT)

引入

根据上一篇的内容,我们知道,FFT 可以将多项式转化为点值序列,通过点值序列的乘积来实现多项式的乘积。

而为了保证还原出的多项式就是我们需要求的多项式,我们需要保证 FFT 取点的数量高于多项式的度数。

那么,当我们没保证取点的数量高于多项式数量时,会发生什么呢?

循环卷积

不妨假设一个度数足够大的多项式 \(\displaystyle F(x)=\sum_{i=0}^n a_ix^i\) ,做 FFT 时只取了 \(m\) 个点 \(\omega_m^0, \omega_m^1, \omega_m^2, \cdots, \omega_m^{m-1}\)

那么,当我们代入 \(\omega_m^k\) 时有:

\(\begin{aligned} &F(\omega_m^k) \\=&\sum_{i=0}^n a_i \omega_m^{ik} \\=&\sum_{i=0}^{m-1} a_i\omega_m^{ik}+\sum_{i=m}^{2m-1} a_i\omega_m^{ik} +\cdots +\sum_{i=\lfloor{n\over m}\rfloor\cdot m}^n a_i\omega_m^{ik} \\=&\sum_{i=0}^{m-1} a_i\omega_m^{ik}+\sum_{i=0}^{m-1} a_{i+m}\omega_m^{(i+m)k} +\cdots +\sum_{i=0}^{n\bmod m} a_{i+\lfloor{n\over m}\rfloor\cdot m}\omega_m^{(i+\lfloor{n\over m}\rfloor\cdot m)k} \\=&\sum_{i=0}^{m-1} a_i\omega_m^{ik}+\sum_{i=0}^{m-1} a_{i+m}\omega_m^{ik} +\cdots +\sum_{i=0}^{n\bmod m} a_{i+\lfloor{n\over m}\rfloor\cdot m}\omega_m^{ik} \\=&\sum_{i=0}^{m-1}(\sum_{t=0}^n [t\equiv i\pmod m]a_t) \omega_m^{ik} \end{aligned}\)

这等价于一个只有 \(m\) 次的多项式,其中第 \(i\) 项系数为 \(\displaystyle (\sum_{t=0}^n [t\equiv i\pmod m]a_t)\) 的多项式进行 FFT 的结果。

也就是说,假设我们将 \(F(x)\) 进行 FFT 后,再进行 IFFT,还原出的不再是 \(F(x)\) 而是上述这个多项式。而是 \(\displaystyle f_k'=\sum_{i\equiv k\pmod {N}}a_ib_j\) ,其中 \(N=2^K\) 是 FFT 分治的上界。

两个多项式乘积时也是类似的:

习惯上,我们称形如 \(\displaystyle c_k=\sum_{i+j=k}a_ib_j\) 形式的卷积为加法卷积,因为它对答案的贡献是加法(同理也有减法、乘法、按位与、按位或、按位异或等卷积)。我们原先以为,FFT 可以通过自己的变换使得加法卷积可以高速的计算。

然而,通过这个例子,我们发现,实际上 FFT 计算的只是 \(\displaystyle c_k=\sum_{i+j\equiv k\pmod {N}}a_ib_j\) 。我们称之为循环加法卷积,在不引起歧义的情况下,可以直接简称为循环卷积。


快速傅里叶逆变换

根据上述的条件,我们知道了,FFT-IFFT 优化多项式 \(A(x), B(x)\) 乘积得到 \(C(x)\) 的结果,其实是得到了 \(\displaystyle c_k=\sum_{i+j\equiv k\pmod {N}}a_ib_j\)

我们假设三者的度数均小于 \(N\) ,代入原式可以进行化简:

\(\begin{aligned} &C(x) \\=&\sum_{k=0}^{N-1} c_kx^k \\=&\sum_{k=0}^{N-1} x^k\sum_{i+j\equiv k\pmod {N}}a_ib_j \\=&\sum_{k=0}^{N-1} x^k\sum_{i=0}^{N-1}\sum_{j=0}^{N-1}a_ib_j[i+j\equiv k\pmod N] \\=&\sum_{k=0}^{N-1} x^k\sum_{i=0}^{N-1}\sum_{j=0}^{N-1}a_ib_j[N\mid (i+j-k)] \end{aligned}\)

这里,我们根据 第三篇 中提到的单位根反演,可以拆解得到:

\(\begin{aligned} &C(x) \\=&\sum_{k=0}^{N-1} x^k\sum_{i=0}^{N-1}\sum_{j=0}^{N-1}a_ib_j[N\mid (i+j-k)] \\=&\sum_{k=0}^{N-1} x^k\sum_{i=0}^{N-1}\sum_{j=0}^{N-1}a_ib_j{1\over N}\sum_{t=0}^{N-1}\omega_N^{(i+j-k)t} \\=&\sum_{k=0}^{N-1} x^k\cdot \sum_{t=0}^{N-1}{1\over N}\omega_N^{-kt}(\sum_{i=0}^{N-1} a_i\omega_N^{it})(\sum_{j=0}^{N-1} b_j\omega_N^{jt}) \end{aligned}\)

这个式子乍一看及其鬼畜,但我们理一下它每部分的含义:

\(\displaystyle \sum_{i=0}^{N-1}a_i\omega_N^{it}\) 其实是 \(A(x)\) 进行 FFT 后的第 \(t\) 项,我们记为 \((\mathcal{F}A)_t\) ,同理 \(\displaystyle \sum_{j=0}^{N-1} b_j\omega_N^{jt}\)\((\mathcal FB)_t\)

根据点值的乘积就是乘积的点值,我们知道 \(\displaystyle (\sum_{i=0}^{N-1} a_i\omega_N^{it})(\sum_{j=0}^{N-1} b_j\omega_N^{jt})=(\mathcal FA)_t(\mathcal FB)_t=(\mathcal FC)_t\) ,即 \(C(x)\) FFT 后的第 \(t\) 项。

于是我们代入就会看到 \(\displaystyle C(x)=\sum_{i=0}^{N-1} x^k\cdot \sum_{t=0}^{N-1}{1\over N}\omega_N^{-kt}(\mathcal FC)_t\)

而我们又知道,\(C(x)\) 的第 \(i\) 项系数 \(\displaystyle c_k=\sum_{t=0}^{N-1}{1\over N}\omega_N^{-kt}(\mathcal FC)_t\)

相比快速傅里叶变换的式子 \(\displaystyle (\mathcal FC)_t=\sum_{k=0}^{N-1}\omega_N^{kt}c_k\) ,我们知道了,上面的式子就是快速傅里叶逆变换的式子,其中 \(\displaystyle {1\over N}\omega_N^{-kt}\) 就是 IFFT 的系数。


快速傅里叶逆变换的实现技巧

当然,我们一般而言,实现的时候并不会直接采用 \({1\over N}\omega_n^{-kt}\) 当作 IFFT 的系数,这样的话复杂度又会退化为 \(O(N^2)\)

考虑到 \(\omega_n^{kt}\) 作为 FFT 的系数可以分治实现,如果 IFFT 的系数是 \(\omega_n^{-kt}\) 的话,也同样可以分治实现了,只是单位根从 \(\omega_n\) 更换为 \(\omega_n^{-1}\) 的问题。

然而,恰好我们可以这样处理:考虑到 \(\displaystyle C(x)=\sum_{i=0}^{N-1} x^k\cdot \sum_{t=0}^{N-1}{1\over N}\omega_N^{-kt}(\mathcal FC)_t={1\over N}\sum_{i=0}^{N-1} x^k\cdot \sum_{t=0}^{N-1}\omega_N^{-kt}(\mathcal FC)_t\)

我们可以先直接按 \(\omega_n^{-kt}\) 进行 IFFT ,求出 \(C(x)\)\(N\) 倍。最后统一 \(O(N)\) 遍历,将所有系数乘上 \({1\over N}\) 还原。

当然,另一个等价的实现方法是直接采用 \({1\over p_{flr}}\omega_n^{-kt}\) 来分治进行 IFFT 。其中 \(p_{flr}\) 表示第 \(flr\) 层的分治数量,在我们这次的篇章中,其恒等于 \(2\)

此外,我们在 FFT 时优先计算出了 \(\omega_N^k\) 数组,通过更新步长的方式,避免了反复计算单位根的计算量。同理,在 IFFT 时,也可以有限计算 \(\omega_N^{-k}\) 数组或 \({1\over 2}\omega_N^{-k}\)

然而,计算 \(\omega_N^{-k}\) 数组时,考虑到 \(\omega_N^{-k}=\omega_N^{N-k}\) 。于是有 \(\omega_N^{-0}=\omega_N^0, \omega_N^{-k}=\omega_N^{N-k}(k>0)\) ,只需要将 \(\omega_N^k\) 数组复制一遍,然后将第 \(2\) 项到第 \(N\) 项首尾翻转即可。(当然,也可以不复制,直接原位上翻转)

当然,如果是计算 \({1\over 2}\omega_N^{-k}\) 就需要翻转后再乘上 \({1\over 2}\)

我们推荐的是采用第一种方法,因为它具有更好的拓展性。


快速傅里叶逆变换是实现

int N;
vir w[MAXN];//单位根
inline void FFT(vir *f, int n) {
	static vir tmp[MAXN];
	if(n==1) return ;
	
	//奇偶分列
	for(int i=0; i<n; ++i) tmp[i]=f[i];
	vir *fl=f, *fr=f+n/2;
	for(int i=0, j=0; i<n; i+=2, ++j) fl[j]=tmp[i];
	for(int i=1, j=0; i<n; i+=2, ++j) fr[j]=tmp[i];
	
	//递归求解
	FFT(fl, n/2); FFT(fr, n/2);
	
	//O(n) 合并
	vir *o=w;
	int pace = N/n;
	for(int i=0; i<n/2; ++i) {
		tmp[i] = fl[i] + *o * fr[i]; 
		tmp[i + n/2] = fl[i] - *o * fr[i];
		o += pace;
	}
	for(int i=0; i<n; ++i) f[i]=tmp[i];
}
inline void IFFT(vir *f, int n) {
	reverse(w+1, w+N);
	FFT(f, n);
	vir invn=!vir(n);//求逆元
	for(int i=0; i<n; ++i) f[i]*=invn;
}
inline void mul(vir *f, vir *g, int n, int m) {
	for(N=1; N<n+m-1; N<<=1);//求最小的 2 的幂次
	vir omega = get_omega(N);
	
	w[0]=1;
	for(int i=1; i<N; ++i) w[i]=w[i-1]*omega;
	FFT(f, N); FFT(g, N);
	for(int i=0; i<N; ++i) f[i]*=g[i];
	IFFT(f, N);
}

例题

洛谷 P3803 【模板】多项式乘法(FFT)

参考代码

其中 \(g=3\) 是质数 \(P=998244353=119\cdot 2^{23}+1\) 的原根。

#include <bits/stdc++.h>
using namespace std;
const int MAXN=1<<21;
const int P=998244353;
typedef unsigned uint;
typedef unsigned long long ull;

inline uint norm(uint v, uint p=P) { return v>=p?v-p:v; }
inline uint mul(uint a, uint b, uint p=P) { return (ull)a*b%P; }
inline uint kpow(uint a, uint x, uint p=P) { uint ans=1; for(; x; x>>=1, a=mul(a, a, p)) if(x&1) ans=mul(ans, a, p); return ans; }
uint exgcd(uint a, uint b, uint &x, uint &y) {
	static uint g;
	return b?(exgcd(b, a%b, y, x), y-=a/b*x, g):(x=1, y=0, g=a);
}
inline uint inv(uint a, uint p=P) {
	uint x, y;
	return (exgcd(a, p, x, y)==1)?norm(x+p):0;
}

struct Z {
	uint v;
	inline Z(uint v_=0):v(norm(v_)) {}
	inline Z& operator += (const Z& x) { v=norm(v+x.v); return *this; }
	inline Z& operator -= (const Z& x) { v=norm(v+P-x.v); return *this; }
	inline Z& operator *= (const Z& x) { v=mul(v, x.v); return *this; }
	
	inline Z operator + (const Z &x) const { Z y=*this; return y+=x; }
	inline Z operator - (const Z &x) const { Z y=*this; return y-=x; }
	inline Z operator * (const Z &x) const { Z y=*this; return y*=x; }
	
	inline Z operator - () const { return Z(P-v); }
	inline Z operator ! () const { return Z(inv(v)); }
	inline operator uint() const { return v; }
};
typedef Z vir;

int n, m;
vir f[MAXN], g[MAXN];
inline vir get_omega(int n) { return kpow(3, (P-1)/n); }
int N;
vir w[MAXN];
inline void FFT(vir *f, int n) {
	static vir tmp[MAXN];
	if(n==1) return ;
	
	for(int i=0; i<n; ++i) tmp[i]=f[i];
	vir *fl=f, *fr=f+n/2;
	for(int i=0, j=0; i<n; i+=2, ++j) fl[j]=tmp[i];
	for(int i=1, j=0; i<n; i+=2, ++j) fr[j]=tmp[i];
	
	FFT(fl, n/2); FFT(fr, n/2);
	
	vir *o=w;
	int pace = N/n;
	for(int i=0; i<n/2; ++i) {
		tmp[i] = fl[i] + *o * fr[i]; 
		tmp[i + n/2] = fl[i] - *o * fr[i];
		o += pace;
	}
	for(int i=0; i<n; ++i) f[i]=tmp[i];
}
inline void IFFT(vir *f, int n) {
	reverse(w+1, w+N);
	FFT(f, n);
	vir invn=!vir(n);
	for(int i=0; i<n; ++i) f[i]*=invn;
}
inline void mul(vir *f, vir *g, int n, int m) {
	for(N=1; N<n+m-1; N<<=1);
	vir omega = get_omega(N);
	
	w[0]=1;
	for(int i=1; i<N; ++i) w[i]=w[i-1]*omega;
	FFT(f, N); FFT(g, N);
	for(int i=0; i<N; ++i) f[i]*=g[i];
	IFFT(f, N);
}
int main() {
	ios::sync_with_stdio(0);
	cin.tie(0); cout.tie(0);
	cin>>n>>m;
	++n; for(uint i=0, v; i<n; ++i) cin>>v, f[i]=v;
	++m; for(uint i=0, v; i<m; ++i) cin>>v, g[i]=v;
	mul(f, g, n, m);
	for(int i=0, I=n+m-1; i<I; ++i) cout<<f[i]<<" \n"[i==I-1];
	cout.flush();
	return 0;
}
posted @ 2023-02-08 09:35  JustinRochester  阅读(615)  评论(0编辑  收藏  举报