快速傅里叶逆变换（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;
}