FFT快速傅里叶变换原理及实现代码

解决什么问题

求两个多项式乘积或者卷积, 也就是将问题转化为多项式乘积

上述$A(x)\times B(x)$就是求卷积, 直接求时间复杂度$O(n^2)$, $FFT$能优化为$O(n\log n)$

多项式表示法转化为点表示法

性质: 任意一个$n$次多项式$A(x)=a_0+a_{1}x+a_2{x}^2+...+a_n{x}^n$, 在函数图像上任意取$n+1$个点, 唯一能确定一个$n$次多项式

上图就是选取的$n+1$个点

上述就是范德蒙行列式, 结果等于
$$\prod _{1\le i<j\le n+1}(x_i-x_j)$$
系数矩阵是满秩的

$FFT$使用了$n+1$个特殊点, 如果求点的乘积, 对应$y$方向坐标相乘即可, 也就是$A(x_i)B(x_i)$, 非常方便, 因此目标就变为如何快速的将系数表示法转化为点表示法, 在点表示法计算完结果后, 再转化回系数表示法

复数的性质

取复数域上的单位根作为不同的点
复数$a+bi$, $a$是实部, $b$是虚部

复数的运算性质

• 两个复数相加$a+bi+c+di=(a+c)+(b+d)i$, 满足平行四边形
• 两个复数相乘$(a+bi)(c+di)=(ac-bd)+(ad-bc)i$, 相乘得到的复数的模长等于原来两个复数模长相乘, 新的辐角等与原来两个辐角相加

复数域上的单位根

等分成$n$份, 取其中$k$份, ${\omega }_n^k$, 被称为$n$次单位根, 一般来说将$n$变为$2$的整次幂方便运算

$${\omega }_n^i\times {\omega }_n^j={\omega }_n^{i+j}$$

具有上述五条性质

$FFT$正变换

将原多项式转化为点表示法

将原式的奇偶项分开, 使用换元法$x=x^2$, 就有上述$A_1(x)$和$A_2(x)$两个多项式
也就推导出$A(x)=A_1(x^2)+x{A}_2(x^2)$

对于$[\frac{\pi }{2},n-1]$的点, 减去$\frac{\pi }{2}$, 映射到$[0,\frac{\pi }{2}-1]$

想要计算$[0,n-1]$区间内的$\omega$值, 将区间分为左右两份, 第一个区间预处理$A_1({\omega }_{\frac{n}{2}}^k)$, 第二个区间预处理$A_2({\omega }_{\frac{n}{2}}^k)$, 然后就可以根据这两个值计算

本质是基于分治算法, 每一层时间复杂度$O(n)$, 最多递归$\log n$层, 总的时间复杂度$O(n\log n)$

$FFT$逆变换

将点表示法快速的求出原来的多项式
假设已经有点集$({\omega }_n^k,y_k)$, 计算原多项式系数

系数如下

$$C_k=\sum _{i=0}^{n-1}{y}_i({\omega }_n^{-k}{)}^i$$
$A_k=\frac{C_k}{n}$

将$C_k$做上图变换, 求$C_k$等价于再求一次$FFT$

迭代实现

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>

using namespace std;

const int N = 300010;  // 定义数组的最大长度
const double PI = acos(-1);  // 定义π的值

int n, m;  // 两个多项式的最高次数
struct Complex {
	double x, y;  // 复数的实部和虚部
	Complex operator+(const Complex &t) const  // 复数加法
	{
		return {x + t.x, y + t.y};

	}
	Complex operator-(const Complex &t) const  // 复数减法
	{
		return {x - t.x, y - t.y};
	}

	Complex operator*(const Complex &t) const  // 复数乘法
	{
		return {x * t.x - y * t.y, x * t.y + y * t.x};
	}
} a[N], b[N];  // 存储两个多项式的系数
int rev[N], bit, tot;  // rev数组用于存储位逆序置换，bit是二进制位数，tot是总长度

// FFT函数，inv为1时是FFT，inv为-1时是逆FFT
void fft(Complex a[], int inv) {
	// 位逆序置换
	for (int i = 0; i < tot; i++)
		if (i < rev[i])
			swap(a[i], a[rev[i]]);  // 将数组元素按照位逆序排列

	// 迭代实现FFT
	for (int mid = 1; mid < tot; mid <<= 1)  // mid表示当前区间长度的一半
	{
		auto w1 = Complex({cos(PI / mid), inv * sin(PI / mid)});  // 单位根
		for (int i = 0; i < tot; i += mid * 2)  // 遍历每个区间
		{
			auto wk = Complex({1, 0});  // 初始旋转因子
			for (int j = 0; j < mid; j++, wk = wk * w1)  // 遍历区间内的每个元素
			{
				auto x = a[i + j], y = wk * a[i + j + mid];  // 蝶形操作
				a[i + j] = x + y, a[i + j + mid] = x - y;  // 更新结果
			}
		}
	}
}

int main() {
	scanf("%d%d", &n, &m);  // 输入两个多项式的最高次数
	for (int i = 0; i <= n; i++) scanf("%lf", &a[i].x);  // 输入第一个多项式的系数
	for (int i = 0; i <= m; i++) scanf("%lf", &b[i].x);  // 输入第二个多项式的系数

	// 计算tot和bit，tot是大于等于n+m+1的最小的2的幂次
	while ((1 << bit) < n + m + 1) bit++;
	tot = 1 << bit;

	// 初始化rev数组，用于位逆序置换
	for (int i = 0; i < tot; i++)
		rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));

	// 对a和b进行FFT
	fft(a, 1), fft(b, 1);

	// 将a和b的FFT结果相乘
	for (int i = 0; i < tot; i++) a[i] = a[i] * b[i];

	// 对a进行逆FFT
	fft(a, -1);

	// 输出结果，即多项式乘积的系数
	for (int i = 0; i <= n + m; i++)
		printf("%d ", (int) (a[i].x / tot + 0.5));

	return 0;
}