fft学习小记

Preface

这东西要打多几遍。

FFT 快速傅里叶变换

核心思想：

利用单位复数根的性质，尝试分治，在$O(nlogn)$的时间内完成点值和插值运算。

单位复数根

从$(1,0)$出发，每次逆时针旋转$\frac{2\pi }{n}$弧得到一个新的复数根。所以单位复数根可以表示为：$$\omega (cos(2\ast i\ast pi/m),sin(2\ast i\ast pi/m))$$

根据复数运算法则，可以知道，两个复数相乘，实际上是幅角相加，模长相乘。

所以有${\omega }_n^{i+j}={\omega }_n^i\ast {\omega }_n^j$

那么单位复数根就有以下性质：

1. 折半引理：$$({\omega }_n^{k+n/2}{)}^2={\omega }_n^{2k+n}={\omega }_n^{2k}{\omega }_n^n={\omega }_n^{2k}=({\omega }_n^k{)}^2$$

这个引理保证了FFT分治策略的时间复杂度。

2. 消去引理：

当$n$不整除$k$时
$$\sum _{j=0}^{n-1}({\omega }_n^k{)}^j=\frac{({\omega }_n^k{)}^n-1}{{\omega }_n^k-1}=\frac{({\omega }_n^n{)}^k-1}{{\omega }_n^k-1}=\frac{(1{)}^k-1}{{\omega }_n^k-1}=0$$

当$n$整除$k$时，就有${\omega }_n^k=1$，所以上式的值显然为$n$.

这个引理推出FFT的插值运算与点值运算的关系，是插值运算的基础。

分治策略

考虑一个多项式$A(x)$，按照系数下标的奇偶来分成两个次数界为$\frac{n}{2}$的多项式$A_0,A_1$，即$$A^{[0]}(x)=a_0+a_{2}x+a_4{x}^2+...+a_{n-2}{x}^{n/2-1}$$$$A^{[1]}(x)=a_1+a_{3}x+a_5{x}^2+...+a_{n-1}{x}^{n/2-1}$$

然后可以写出：$$A(x)=A^{[0]}(x^2)+A^{[1]}(x^2)x$$

然后根据折半引理，可以分成两个子问题求解。

时间复杂度就是：$$T(n)=2T(n/2)+\Theta (n)=\Theta (nlog(n))$$

逆DFT

实际上就是上面的逆运算。

根据矩阵的思想：$$\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ y_3\\ ⋮\\ y_{n-1}\end{array}\right]=\left[\begin{array}{cccccc}1& 1& 1& 1& \dots & 1\\ 1& {\omega }_n& {\omega }_n^2& {\omega }_n^3& \dots & {\omega }_n^{n-1}\\ 1& {\omega }_n^2& {\omega }_n^4& {\omega }_n^6& \dots & {\omega }_n^{2(n-1)}\\ 1& {\omega }_n^3& {\omega }_n^6& {\omega }_n^9& \dots & {\omega }_n^{3(n-1)}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& {\omega }_n^{3(n-1)}& \dots & {\omega }_n^{(n-1)(n-1)}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ ⋮\\ a_{n-1}\end{array}\right]$$来求出新的$\omega$取值。

可以证明就是：$$[V_n^{-1}{V}_n{]}_{j{j}^{\prime }}=\sum _{k=0}^{n-1}({\omega }_n^{-kj}/n)({\omega }_n^{k{j}^{\prime }})=\sum _{k=0}^{n-1}{\omega }_n^{k(j^{\prime }-j)}/n$$
然后做一遍就好了。

Code

#include <bits/stdc++.h>

#define F(i,a,b) for (int i=a;i<=b;i++)
typedef double db;

const int M = 4 * 1e5 + 10;
const db pi = acos(- 1);

using namespace std;

int cnt, n, m, len;

struct Z {
	db x, y;
	Z (db _x=0, db _y=0) { x = _x, y = _y; }
	friend Z operator + (Z a, Z b) { return Z(a.x + b.x, a.y + b.y); }
	friend Z operator - (Z a, Z b) { return Z(a.x - b.x, a.y - b.y); }
	friend Z operator * (Z a, Z b) { return Z(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x); }
} a[M], b[M], t[M];

void Dft(Z *a, int n, int sig) {
	F(i, 0, n - 1) {
		int pos = 0;
		for (int x = i, y = 0; y < cnt; y ++, x >>= 1) pos = (pos << 1) + (x & 1);
		t[pos] = a[i];
	}
	for (int m = 2; m <= n; m <<= 1) {
		int half = m >> 1;
		Z W(cos(2 * sig * pi / m), sin(2 * sig * pi / m));
		for (int i = 0; i < n; i += m) {
			Z o(1, 0);
			for (int j = i; j < i + half; j ++, o = o * W) {
				Z u = t[j + half] * o;
				t[j + half] = t[j] - u;
				t[j] = t[j] + u;
			}
		}
	}
	F(i, 0, n - 1)
		a[i] = t[i];
	if (sig == - 1) {
		F(i, 0, n - 1)
			a[i].x = int(a[i].x / n + 0.5);
	}
}

int main() {
	scanf("%d%d", &n, &m);
	F(i, 0, n) scanf("%lf", &a[i].x);
	F(i, 0, m) scanf("%lf", &b[i].x);
	
	for (len = 1, cnt = 0; len <= n + m + 1; len <<= 1, cnt ++);
	Dft(a, len, 1), Dft(b, len, 1);
	F(i, 0, len - 1) a[i] = a[i] * b[i];
	Dft(a, len, - 1);

	F(i, 0, n + m) printf("%0.lf ", a[i].x);
}