# 快速傅里叶逆变换（IFFT）

## 循环卷积

\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}

## 快速傅里叶逆变换

\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$$

## 快速傅里叶逆变换是实现

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);
}


## 例题

### 参考代码

#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编辑  收藏  举报