模板快速傅里叶变换

Posted 2021-05-20 top_secret的小尛博客

P3803

给定一个 \\(n\\) 次多项式 \\(F(x)\\)，和一个 \\(m\\) 次多项式 \\(G(x)\\)。 请求出 \\(F(x)\\) 和 \\(G(x)\\) 的卷积。

\\(n, m \\leq {10}^6\\)

快速傅里叶变换

FFT变换

\\(n=2^k\\)，已知\\(f(x)=\\sum_{i\\in[0,2^k)}a_ix^i\\)，求所有\\(y_j = f(\\omega_{2^k}^{j})\\)。

---

做法：由于

\\[y_j=f(\\omega_{2^k}^{j})=\\sum_{i\\in[0,2^k)}a_i(\\omega_{2^k}^{j})^i\\\\ y_{j+2^{k-1}}=f(\\omega_{2^k}^{j+2^{k-1}})=\\sum_{i\\in[0,2^k)}a_i(\\omega_{2^k}^{j+2^{k-1}})^i=\\sum_{i\\in[0,2^k)}a_i(-\\omega_{2^k}^j)^i \\]

令 \\(f_i(x) = \\sum_{j\\in[0,2^{k-1})}a_{2j+i}x^j\\quad(i\\in\\{0,1\\})\\)：

\\[y_j = f_0((\\omega_{2^k}^j)^2)+\\omega_{2^k}^j f_1((\\omega_{2^k}^j)^2) = f_0(\\omega_{2^{k-1}}^j)+\\omega_{2^k}^j f_1(\\omega_{2^{k-1}}^j) \\\\ y_{j+2^{k-1}} = f_0((-\\omega_{2^k}^j)^2)-\\omega_{2^k}^{j}f_1((-\\omega_{2^k}^j)^2) = f_0(\\omega_{2^{k-1}}^j)-\\omega_{2^k}^j f_1(\\omega_{2^{k-1}}^j) \\]

故只需求出所有\\(y\'_j=f_1(\\omega_{2^{k-1}}^j),y\'\'_j=f_2(\\omega_{2^{k-1}}^j)\\)即可。（成功划分成了两个子问题）

此时\\(y_j=y\'_j+\\omega_{2^k}^jy\'\'_j, y_{j+2^{k-1}}=y\'_j-\\omega_{2^k}^jy\'\'_j\\)。

FFT逆变换

假设已有\\(y_j=\\sum_{i\\in[0,n)}a_i\\omega_n^{ij}\\)，要求\\(a_i\\)。

---

考虑对\\(y_j\\)再做一个反的FFT，得到：

\\[\\begin{aligned} x_i &= \\sum_{j\\in[0,n)}y_j\\omega_n^{-ij} \\\\ &= \\sum_{j\\in[0,n)}\\sum_{k\\in[0,n)}a_k\\omega_n^{kj}\\omega_n^{-ij} \\\\ &= \\sum_{j\\in[0,n)}\\sum_{k\\in[0,n)}a_k\\omega_n^{(k-i)j} \\\\ &= \\sum_{k\\in[0,n)}a_k\\sum_{j\\in[0,n)}(\\omega_n^{k-i})^j \\\\ &= na_{i} \\end{aligned} \\]

即得\\(a_i=x_i/n\\)

递归实现^[1]

void fast_fast_tle(int limit, complex *a, int type) {
	if (limit == 1) return ;
	complex a1[limit >> 1], a2[limit >> 1];
	for (int i = 0; i <= limit; i += 2)
		a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];
	fast_fast_tle(limit >> 1, a1, type);
	fast_fast_tle(limit >> 1, a2, type);
	complex Wn = complex(cos(2.0 * Pi / limit) , type * sin(2.0 * Pi / limit)), w = complex(1, 0);
	for (int i = 0; i < (limit >> 1); i++, w = w * Wn)
		a[i] = a1[i] + w * a2[i],
			   a[i + (limit >> 1)] = a1[i] - w * a2[i];
}

非递归实现^[1:1]

预处理r[i]（i的二进制表示的反转）：

for (int i = 0; i < limit; i++)
	r[i] = ( r[i >> 1] >> 1 ) | ( (i & 1) << (l - 1) ) ;

代码：

void fast_fast_tle(complex *A, int type) {
	for (int i = 0; i < limit; i++)
		if (i < r[i]) swap(A[i], A[r[i]]); //r[r[i]] == i
	for (int mid = 1; mid < limit; mid <<= 1) { //自底向上，mid表示区间长度的一半
		complex Wn( cos(Pi / mid) , type * sin(Pi / mid) ); //mid<<1次单位根。
		for (int R = mid << 1, j = 0; j < limit; j += R) { //R为区间长度，j为当前的最右端点
			complex w(1, 0);
			for (int k = 0; k < mid; k++, w = w * Wn) {
				complex x = A[j + k], y = w * A[j + mid + k];
				A[j + k] = x + y;
				A[j + mid + k] = x - y;
			}
		}
	}
}

有关complex

C++的complex<T>类

手写struct complex^[1:2]

const double Pi = acos(-1.0);
struct complex {
    double x, y;
    complex (double xx = 0, double yy = 0) {x = xx, y = yy;}
} a[MAXN], b[MAXN];
complex operator + (complex a, complex b) { return complex(a.x + b.x , a.y + b.y);}
complex operator - (complex a, complex b) { return complex(a.x - b.x , a.y - b.y);}
complex operator * (complex a, complex b) { return complex(a.x * b.x - a.y * b.y , a.x * b.y + a.y * b.x);}

---

1. 参考文章 ↩︎ ↩︎ ↩︎