【笔记】快速傅里叶变换

0 约定

1. $\exp (x)=e^x$。
2. 有些地方标注有 ?，系本人不太能保证严谨性的部分。

---

1 复数 (Complex)

1.1 三种形式

1. 代数形式：$z=a+bi$，其中 $a,b\in \mathbb{R}$。

2. 三角形式：$z=r(\cos \theta +i\sin \theta )$，其中 $r=\sqrt{a^2+b^2}$（$a,b$ 是在代数形式中定义的）。

3. 指数形式：$z=r\cdot e^{i\theta }$。

   根据 欧拉公式：

   $$e^{ix}=\cos x+i\sin x,x\in \mathbb{C}$$

   令 $x=\theta$，再带入三角形式中，即可得指数形式。

---

1.2 单位根 (unit root)

对于方程 $x^n=1$，在复数意义下，她有 $n$ 个解，我们称这 $n$ 个解都是 单位根 (unit root)。

设 $\boxed{{\displaystyle {\omega }_n=\exp \frac{2\pi i}{n}}}$，则这 $n$ 个单位根可以表示为 ${\omega }_n^k\mid k=0,1\cdots ,n-1$。

在复平面中，$n$ 次单位根把单位圆 $n$ 等分。

1.2.1 性质

1. ${\omega }_n^n=1$

   证：${\omega }_n^n=\exp \left({\displaystyle \frac{2\pi i}{n}}n\right)=(\exp \pi i{)}^2=(-1{)}^2=1$

2. ${\omega }_n^k={\omega }_{pn}^{pk}$

   $\omega$ 的上下标同乘 $p$，相当于 $\exp$ 后面的分子分母同乘 $p$。

3. ${\omega }_{2n}^{k+n}=-{\omega }_{2n}^k$

   几何理解：这实际上相当于将单位园 $2n$ 等分，次数 $+n$ 相当于在复平面中，这两个点关于原点对称，故而她们互为相反数。（？

   证：${\omega }_{2n}^{k+n}=\exp {\displaystyle \frac{2\pi i(k+n)}{2n}}=\exp ({\displaystyle \frac{2\pi ik}{2n}}+{\displaystyle \frac{2\pi in}{2n}})=-\exp \left({\displaystyle \frac{2\pi ik}{2n}}\right)=-{\omega }_{2n}^k$

---

2 快速傅里叶变换 (Fast Fourier Transform)

FFT 可以在 $O(n\log n)$ 的时间复杂度内，计算多项式乘法，优于暴力的 $O(n^2)$。

2.1 定义

对于一个多项式 $f(x)$，定义傅里叶变换为：

$$\boxed{{\displaystyle \mathcal{F}(f(x))=(f(0),f({\omega }_n),\dots ,f({\omega }_n^{n-1}))}}$$

傅里叶变化本质上可以理解为一种「从一个多项式到一个 $n$ 维向量的映射关系」，将多项式的系数表示转化为点值表示。（？

---

2.2 求解

2.2.1 引入（不太重要

考虑如何计算 $\sum _{i=0}^n{a}_i{x}^i$。

一般的方法是，枚举 $i$，然后维护 $x^i$ 的值 $k$，再将 $a_i$ 与 $k$ 相乘并累加。精细统计的话，这样应该需要 $2n$ 次乘法。

另一种办法是：

$$\begin{array}{rl}\sum _{i=0}^n{a}_i{x}^i& =a_0+a_{1}x+\cdots +a_n{x}^n\\ & =a_0+x(a_1+x(a_2+\cdots x{a}_n))& \text{【每次都把后面一部分提一个 }x\text{ 出来】}\end{array}$$

如果是这样计算，总共其实只需要 $n$ 次计算。（似乎没什么用

---

2.2.2 思想

不失一般性的，我们令 $n=2^p(p\in \mathbb{N})$（如果 $n$ 不足 $2^p$，可将更高位的系数视作 $0$。

借用上述思想，对于一个多项式 $f(x)$，我们可以做如下的变换，即将 $f(x)$ 按 $x$ 次数的奇偶分开：

$$\begin{array}{rl}f(x)& =a_0+a_{1}x+\cdots +a_{n-1}{x}^{n-1}\\ & =a_0+a_2{x}^2+\cdots +a_{n-2}{x}^{n-2}+x(a_1+a_3{x}^2+\cdots +a_{n-1}{x}^{n-2})\end{array}$$

令 $f_0(x)=a_0+a_2{x}^1+\cdots +a_{n-2}{x}^{n/2-1}$，$f_1(x)=a_1+a_3{x}^1+\cdots +a_{n-1}{x}^{n/2-1}$。

代回原式：

$$f(x)=f_0(x^2)+x{f}_1(x^2)$$

所以我们有：

$$\begin{array}{rl}f({\omega }_n^k)& =f_0({\omega }_n^{2k})+{\omega }_n^k\times f_1({\omega }_n^{2k})\\ & =f_0({\omega }_{n/2}^k)+{\omega }_n^k\times f_1({\omega }_{n/2}^k)& \text{【单位根性质 2】}\end{array}$$

以及：

$$\begin{array}{rlr}f({\omega }_n^{k+n/2})& =f_0((-{\omega }_n^k{)}^2)-{\omega }_n^k\times f_1((-{\omega }_n^k{)}^2)& \text{【单位根性质 3】}\\ & =f_0({\omega }_n^{2k})-{\omega }_n^k\times f_1({\omega }_n^{2k})\\ & =f_0({\omega }_{n/2}^k)-{\omega }_n^k\times f_1({\omega }_{n/2}^k)& \text{【单位根性质 2】}\end{array}$$

然后我们惊奇的发现，这两个式子只有系数差别，而且 $n$ 的规模每次减半。这就意味着我们可以在 $O(n\log n)$ 的时间复杂度内递归求解出 $f({\omega }_n^k)$ 其中 $k=0,1,\dots ,n-1$，即求出 $\mathcal{F}(f(x))$​​。

而上面最后推出的式子被称做：蝴蝶操作 (Butterfly Operation)。下面这张图形象地展示了「蝴蝶」的含义：

---

2.2.3 实现

通过递归我们可以轻易的模拟出上述过程，为减少常数下介绍一种非递归版本。

画出递归树，如下图，并按照下图的方式给每层的 $f$ 编号：

我们发现，最底层的 $f$ 的二进制下标反转过来并转换成十进制，就等于她对应 $a$ 的下标。

所以我们可以直接预处理出最底层的 $f$，然后自下而上递推回来。（具体实现可以参考 4.1 Code 中的 $r$​。

---

3 快速傅里叶逆变换 (Inverse Fast Fourier Transform)

3.1 定义

设 $f^{\prime }$ 为进行了快速傅里叶变化后的点值，即 $\mathcal{F}(f(x))$。

现在考虑还原为原来的系数数列。这样的操作被称之为 快速傅里叶逆变换 (Inverse Fast Fourier Transform)。

推导过程比较复杂，这里直接给出结论：

$$f(x)=\frac{1}{n}{f}^{\prime }(-{\omega }_n^x)$$

---

4 优化多项式乘法

那么傅里叶变换的优势在哪呢？这就要问点值表达了。

由于是用点值表示多项式，两个多项式相乘，显然就是每个点的虚数相乘。

最后再做一次逆变换，相当于插值，就可以还原回最终结果了。

4.1 Code

#include <bits/stdc++.h>
using namespace std;
const int N = 4e6 + 5, logN = 20;
const double PI = acos(-1);
int n, m, lim, len, r[N];
struct Complex {
	double a, b;
	Complex operator + (Complex x) { return {a+x.a, b+x.b}; }
	Complex operator - (Complex x) { return {a-x.a, b-x.b}; }
	Complex operator * (Complex x) { return {a*x.a-b*x.b, a*x.b+b*x.a}; }
} a[N], b[N], c[N];
void fft (Complex *f, int sign) {
	for (int i = 0; i < lim; i++)
		if (i < r[i]) swap(f[i], f[r[i]]);
	for (int mid = 1; mid < lim; mid<<=1) {
		Complex o = {cos(PI/mid), sign*sin(PI/mid)};
		for (int R = mid<<1, j=0; j < lim; j += R) {
			Complex p = {1, 0};
			for (int k = 0; k < mid; k++) {
				Complex f0 = f[j+k], f1 = f[j+k+mid];
				f[j+k] = f0 + p*f1;
				f[j+k+mid] = f0 - p*f1;
				p = p*o;
			}
		}
	}
}
int main () {
	ios::sync_with_stdio(0);
	cin.tie(0), cout.tie(0);
	cin >> n >> m;
	for (int i = 0, x; i <= n; i++) cin >> x, a[i] = {x*1.0, 0};
	for (int i = 0, x; i <= m; i++) cin >> x, b[i] = {x*1.0, 0};
	for (lim=1, len=0; lim <= n+m; lim<<=1, len++);
	for (int i = 0; i <= lim; i++)
		r[i] = (r[i>>1]>>1)|((i&1)<<(len-1));
	fft(a, 1), fft(b, 1);
	for (int i = 0; i <= lim; i++) c[i] = a[i]*b[i];
	fft(c, -1);
	for (int i = 0; i <= n+m; i++)
		cout << ((int)(c[i].a/lim+0.5)) << ' ';
	return 0;
}

---

5 例题

「BZOJ4503」两个串

$n=|S|,m=|T|$

首先将 $\{?,a,b,\dots ,z\}\to \{0,1,2,\dots ,26\}$​。然后将 $T$ reverse 下。

在 $s$ 的第 $i$ 匹配成功，可以转换为：

$$c[i]=\sum _{j=0}^{min(m-1,i)}(s[i-j]-t[j]{)}^{2}t[j]=0$$

平方是为了防止负数可能的抵消，如果是绝对值不能化简。

然后暴力拆开看下：

$$c[i]=\sum _{j=0}^{min(m-1,i)}(s[i-j{]}^{2}t[j])+(t[j{]}^3)-(2s[i-j]t[j{]}^2)$$

然后发现，这就是两个卷积和一个普通累加相加，fft 加速算即可。

---

「ZJOI2014」力

Warning: 不同于题面，下面的数组均从 $0$ 开始标号。

令：

$$\begin{array}{rl}{A}_i& =\sum _{j=0}^{i-1}\frac{q_j}{(i-j{)}^2}\\ B_i& =\sum _{j=i+1}^{n-1}\frac{q_j}{(i-j{)}^2}\end{array}$$

于是：

$$E_i=A_i-B_i$$

考虑分开计算 $A,B$。

再令：

$$p_i={\displaystyle \frac{1}{i^2}},\text{特别地}{p}_0=0$$

由于 $p_0=0$，于是 $i=j$​ 可以被纳入运算，而不造成影响。

再把 $q$ reverse 下，令作 $q^{\prime }$。

于是：

$$\begin{array}{rl}A[i]& =\sum _{j=0}^{i}q[j]p[i-j]\\ B[i]& =\sum _{j=i}^{n-1}q[j]p[j-i]\\ & =\sum _{j=0}^{n-i-1}q[i+j]p[j]\\ & =\sum _{j=0}^{n-i-1}{q}^{\prime }[n-1-(i+j)]p[j]\end{array}$$

又令：$t=n-i-1$。

有：

$$B[i]=\sum _{j=0}^t{q}^{\prime }[t-j]p[j]$$

fft 加速算即可。

---

万径人踪灭

假定 $a_i=0(i\ge n)$。

令：$t0=2i,t1=2i+1$。

设：

• $f_i$ ：以 $i$ 为对称中心，两边最多能匹配的对数（包括自己）；
• 而 $g$ 是枚举分割线。

则有

$$\begin{array}{rl}\{a,b\}& \to \{-1,1\}\\ f_i& =\sum _{j=0}^{min(i,n-i-1)}{\displaystyle \frac{a[i-j]\times a[i+j]+1}{2}}\\ & ={\displaystyle \frac{1}{4}}((\sum _{j=0}^{t0}a[t0-j]\times a[j])+(a[i]\times a[i]+1)+min(t0,n-1)-max(0,t0-n+1)+1)\\ g_i& =\sum _{j=1}^{min(i,n-i)}{\displaystyle \frac{a[i-j]\times a[i+j-1]+1}{2}}\\ & ={\displaystyle \frac{1}{4}}((\sum _{j=0}^{t1}a[t1-j]\times a[j])+min(t1n-1)-max(0,t0-n+1)+1)\end{array}$$

然后答案就是：

$$Ans=\sum _{i=0}^{n-1}(2^{f[i]}-1)+\sum _{i=1}^{n-1}(2^{g[i]}-1)-res$$

其中 $res$ 为字符串中回文串个数，这个可以用 Manacher 算。

$f,g$ 卷积算完之后加上后面的那部分即可。

---

「BZOJ 3771」Triple

---

「Luogu5488」差分与前缀和

算一次的前缀：

$$\begin{array}{rl}a1[i]& =\sum _{j=0}^{i}a[j]\\ a2[i]& =\sum _{j=0}^{i}a1[j]=\sum _{j=0}^i\sum _{k=0}^{j}a[j]\end{array}$$

考虑每个 $a[i]$ 对 $Ans$ 的贡献。