【注】參考自算法導論慕爬。
1. 簡介
快速傅里葉變換（FFT）是實現(xiàn)離散傅里葉變換（DFT）和離散傅里葉逆變換（IDFT）的快速算法，其時間復雜度為 隶糕。DFT 在實際生活中有很多應用坝疼，比如通過離散傅里葉變換翠拣，可以將系數(shù)表示的多項式轉(zhuǎn)為點值表示的多項式共耍，從而使得多項式的乘法的復雜度由
降為
虑灰。
2. 實現(xiàn)
FFT 涉及到離散傅里葉變換的知識，F(xiàn)FT 能夠?qū)崿F(xiàn)的基礎在于折半定理痹兜，折半定理提供了求大整數(shù)的 DFT 的一種思路：即通過將大整數(shù)不斷劃分為兩部分穆咐，然后分別求出兩部分的 DFT，最后再合并，本質(zhì)是一種分治的思想对湃。
3. 代碼
3.1 快速傅里葉變換模板
#include <bits/stdc++.h>
using namespace std;
// 復數(shù)
struct Complex{
double re, im;
Complex(): re(0), im(0) {}
Complex(double _re, double _im): re(_re), im(_im) {}
Complex operator + (const Complex cp) const {
return Complex(re+cp.re, im+cp.im);
}
Complex operator - (const Complex cp) const {
return Complex(re-cp.re, im-cp.im);
}
Complex operator * (const Complex cp) const {
return Complex(re*cp.re - im*cp.im, re*cp.im + im*cp.re);
}
Complex operator / (const double div) const {
return Complex(re/div, im/div);
}
void operator = (const Complex cp) {
re = cp.re;
im = cp.im;
}
};
// 快速傅里葉變換
struct FFT{
#ifndef _FFT_
#define _FFT_
#define ll int
#define MAXN 3000005
#endif
// 待計算
ll cnt; // 位數(shù)
ll len; // 數(shù)組大醒陆小（len = 2^cnt）
ll r[MAXN]; // 位逆序置換數(shù)組
const double Pi = acos(-1.0);
FFT(): cnt(0), len(0) {}
// 構(gòu)建 2 的冪次的系數(shù)數(shù)組
void build(Complex a[], ll curlen) {
// 計算 len 所占位數(shù)
cnt = 1;
len = curlen;
while((len >>= 1)) ++cnt;
// 找到最小的大于等于 curlen 的 2 的冪值
len = (ll)1 << cnt;
if(len-curlen == curlen) {
--cnt;
len = curlen;
}
// 超出 curlen 部分的系數(shù)補 0
for(ll i = curlen; i < len; ++i) a[i] = Complex(0,0);
}
// 計算位逆序置換
void bitReverse() {
for(ll i = 0; i < len; ++i) r[i] = i;
// 本質(zhì)是動態(tài)規(guī)劃
// 利用 i>>1 的位逆序推導出 i 的位逆序
for(ll i = 0; i < len; ++i) r[i] = (r[i>>1]>>1) | ((i&1)<<(cnt-1));
}
// 快速傅里葉變換
void DFT(Complex y[]) {
// 位逆序置換
for(ll i = 0; i < len; ++i) {
if(i < r[i]) swap(y[i], y[r[i]]);
}
// 合并 DFT
// 枚舉合并長度
for(ll i = 1; i < len; i<<=1) {
// 主單位復根
Complex wn(cos(Pi/i), sin(Pi/i));
// 枚舉合并組首下標
for(ll j = 0; j < len; j += (i<<1)) {
Complex w(1,0);
// 合并組 DFT
for(ll k = j; k < j+i; ++k, w=w*wn) {
Complex u = y[k];
Complex v = w*y[k+i];
y[k] = u + v;
y[k+i] = u - v;
}
}
}
}
// 快速傅里葉逆變換
void IDFT(Complex a[]) {
// 位逆序置換
for(ll i = 0; i < len; ++i) {
if(i < r[i]) swap(a[i], a[r[i]]);
}
// 合并 DFT
// 枚舉合并長度
for(ll i = 1; i < len; i<<=1) {
// 主單位復根
Complex wn(cos(Pi/i), -sin(Pi/i));
// 枚舉合并組首下標
for(ll j = 0; j < len; j += (i<<1)) {
Complex w(1,0);
// 合并組 DFT
for(ll k = j; k < j+i; ++k, w=w*wn) {
Complex u = a[k];
Complex v = w*a[k+i];
a[k] = u + v;
a[k+i] = u - v;
}
}
}
for(ll i = 0; i < len; ++i) {
a[i] = a[i]/len;
}
}
};
3.2 大整數(shù)乘法模板
#include <bits/stdc++.h>
using namespace std;
// 復數(shù)
struct Complex{
double re, im;
Complex(): re(0), im(0) {}
Complex(double _re, double _im): re(_re), im(_im) {}
Complex operator + (const Complex cp) const {
return Complex(re+cp.re, im+cp.im);
}
Complex operator - (const Complex cp) const {
return Complex(re-cp.re, im-cp.im);
}
Complex operator * (const Complex cp) const {
return Complex(re*cp.re - im*cp.im, re*cp.im + im*cp.re);
}
Complex operator / (const double div) const {
return Complex(re/div, im/div);
}
void operator = (const Complex cp) {
re = cp.re;
im = cp.im;
}
};
// 大整數(shù)乘法（快速傅里葉變換實現(xiàn)）
struct BigIntegersMultiply{
#ifndef _BIGINTEGERSMULTIPLY_
#define _BIGINTEGERSMULTIPLY_
#define ll int
#define MAXN 3000005
#endif
const double Pi = acos(-1.0);
// 待計算
ll cnt; // 位數(shù)
ll len; // 最終乘積數(shù)組大小（len = 2^cnt）
ll r[MAXN]; // 位逆序置換數(shù)組
ll ans[MAXN]; // 最終答案（逆序存儲）
// 待輸入
ll len1, len2; // 初始系數(shù)數(shù)組大小
Complex x1[MAXN], x2[MAXN]; // 兩個大整數(shù)乘數(shù)（逆序儲存）
BigIntegersMultiply(): cnt(0), len(0) {}
// 構(gòu)建 2 的冪次的系數(shù)數(shù)組
void build() {
// 計算兩個乘數(shù)最小的公共 2 的冪次 len 及所占位數(shù) cnt
cnt = 0, len = 1;
while(len < (len1<<1) || len < (len2<<1)) ++cnt, len <<= 1;
// 超出 curlen 部分的系數(shù)補 0
for(ll i = len1; i < len; ++i) x1[i] = Complex(0,0);
for(ll i = len2; i < len; ++i) x2[i] = Complex(0,0);
}
// 計算位逆序置換
void bitReverse() {
for(ll i = 0; i < len; ++i) r[i] = i;
// 本質(zhì)是動態(tài)規(guī)劃
// 利用 i>>1 的位逆序推導出 i 的位逆序
for(ll i = 0; i < len; ++i) r[i] = (r[i>>1]>>1) | ((i&1)<<(cnt-1));
}
// 快速傅里葉變換
void DFT(Complex x[]) {
// 位逆序置換
for(ll i = 0; i < len; ++i) {
if(i < r[i]) swap(x[i], x[r[i]]);
}
// 合并 DFT
// 枚舉合并長度
for(ll i = 1; i < len; i<<=1) {
// 主單位復根
Complex wn(cos(Pi/i), sin(Pi/i));
// 枚舉合并組首下標
for(ll j = 0; j < len; j += (i<<1)) {
Complex w(1,0);
// 合并組 DFT
for(ll k = j; k < j+i; ++k, w=w*wn) {
Complex u = x[k];
Complex v = w*x[k+i];
x[k] = u + v;
x[k+i] = u - v;
}
}
}
}
// 快速傅里葉逆變換
void IDFT(Complex x[]) {
// 位逆序置換
for(ll i = 0; i < len; ++i) {
if(i < r[i]) swap(x[i], x[r[i]]);
}
// 合并 DFT
// 枚舉合并長度
for(ll i = 1; i < len; i<<=1) {
// 主單位復根
Complex wn(cos(Pi/i), -sin(Pi/i));
// 枚舉合并組首下標
for(ll j = 0; j < len; j += (i<<1)) {
Complex w(1,0);
// 合并組 DFT
for(ll k = j; k < j+i; ++k, w=w*wn) {
Complex u = x[k];
Complex v = w*x[k+i];
x[k] = u + v;
x[k+i] = u - v;
}
}
}
for(ll i = 0; i < len; ++i) {
x[i] = x[i]/len;
}
}
// 大整數(shù)乘法
void bigMultiply() {
build();
bitReverse();
DFT(x1), DFT(x2);
for(ll i = 0; i < len; ++i) x1[i] = x1[i]*x2[i];
IDFT(x1);
for(ll i = 0; i < len; ++i) ans[i] = ll(x1[i].re + 0.5);
for(ll i = 0; i < len; ++i) ans[i+1] += ans[i]/10, ans[i] %= 10;
len = len1 + len2 - 1;
while(len && !ans[len]) --len;
++len;
}
};
