快速傅里叶逆变换（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);
}
例题
参考代码
其中 \(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;
}
