快速傅里叶变换 - OI version

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

多项式乘法与卷积

​ 快速傅里叶变换(Fast Fourier Transform，FFT)可以用来快速计算多项式乘法的问题。值得一提的是，“快速傅里叶变换”和“傅里叶变换(Fourier Transform)”是有很大区别的。信奥里所讲的“快速傅里叶变换”其实是“离散傅里叶变换(Discrete Fourier Transform, DFT)”，而“傅里叶变换”则是讨论连续函数的。

​ 我们知道多项式乘法的逐项相乘原理，这是基于乘法分配律的。由于

​ $a(b+c)=ab+ac$

​ 所以$(a_0+a_1x+...+a_nx^n)(b_0+b_1x+...+b_mx^m)$首先可以写作

$a_0(b_0+b_1x+...+b_mx^m)+a_1(b_0+b_1x+...+b_mx^m)+...+a_n(b_0+b_1x+...+b_mx^m)$

​ 然后再把每个括号展开，就形成了逐项相乘的形式。合并同类项以后，可以发现结果的系数可以写成这样的表达式：

​ $c_n = \sum\limits_{i=0}^{n}a_ib_{n-i}$

​ 这个形式在以后将会非常常见，如果把多项式的系数$a、b、c$看作关于下标的函数，那么这个运算就好像是通过$a、b$用某种法则生成了函数$c$。我们把$c$叫做$a,b$的卷积(convolution)，因为$c$好像是由$a$和$b$一正一反卷在一起一样。所以我们也可以说，FFT能快速处理求两个函数的卷积的问题。

系数表达式与点值表达式

​ 一个多项式可以看作一个幂函数，一个$n$次函数需要$n+1$个系数来确定。还有没有别的方法来确定一个$n$次函数？我们知道两点确定一个一次函数，三点确定一个二次函数，那么$n+1$个点的坐标就可以确定一个$n$次函数。也就是说，一个多项式不仅可以通过$n+1$个系数来确定（称为系数表达式），也可以通过$n+1$个坐标点以及其对应的多项式的函数值来确定（称为点值表达式）。

​ 多项式系数逐项相乘的复杂度$O(n^2)$，因为我们需要对两个多项式的每一项两两操作。但我们惊奇的发现，点值表达式下的多项式乘法只需要$O(n)$就可以完成——假设两个点值表达式选择的是相同的坐标点，我们要做的只不过是把函数值相乘而已。

​ 所以，想要快速完成多项式乘法，我们要做的就是解决如何快速完成系数表达式与点值表达式之间的转换问题，这已经称为了多项式乘法的复杂度瓶颈。这就是FFT要做的事。把多项式的系数表达式变换成点值表达式，称为“离散傅里叶变换”。反过来称为“离散傅里叶逆变换”。

离散傅里叶变换(DFT)

​ 问题的关键在于我们应该选择怎么样的点来进行代入。如果我们随便选择$n$个点代入，那么每个点计算点值的复杂度是$O(n)$的，总复杂度依然是$O(n^2)$的。而优化就是看我们是否能少做一些点值的计算。也就是去寻找，当我们把不同的点代入后，它们的结果之间是否存在某种联系使得我们可以不用重新计算。

​ 人们发现，用单位根可以解决这个问题。

​ 单位根其实就是$x^n=1$的$n$个复数解，称为$n$次单位根。根据复数乘法“模长相乘，幅角相加”的规则，在复平面内，它们是包括$(1,0)$在内的将单位圆$n$等分的单位向量。记$(1,0)$之后逆时针第一个单位根为$\omega_n$，第$k$个单位根可以记作$\omega_n^k$。

​ 单位根很容易用欧拉公式表示：

​ $\omega_n^k=\cos\frac{2k\pi}{n}+i\sin\frac{2k\pi}{n}$

​ 现在我们来看看离散傅里叶变换是怎么做的。设有多项式的系数表达式

​ $A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}$

并且我们不妨假定$n$是2的幂——只需要设定多出来的项的系数都是0就好了。现在，我们按照指数的奇偶把它拆成两部分，并构成两个新的多项式。（注意指数也发生了变化）

​ $A_1(x)=a_0+a_2x+...+a_{n-2}x^{\frac{n}{2}-1}$

​ $A_2(x)=a_1+a_3x+...+a_{n-1}x^{\frac{n}{2}-1}$

​ 它们满足

​ $A(x)=A_1(x^2)+xA_2(x^2)$

​ 现在我们将一个$\omega_n^k$代入$A(x)$（其中$k<\dfrac{n}{2}$）得：

​ $A(\omega_n^k)=A_1(\omega_{n}^{2k})+\omega_n^kA_2(\omega_{n}^{2k})$

​ 将相应的$\omega_n^{k+\frac{n}{2}}$代入$A(x)$得：

​ $A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_{n}^{2k+n})+\omega_n^{k+\frac{n}{2}}A_2(\omega_{n}^{2k+n})$

​ 由于我们假定了$n$是偶数，所以$\omega_{n}^{2k}=\omega_{\frac{n}{2}}^{k}$。这就好像“把蛋糕切成$n$份并拿走$2k$份”和“把蛋糕切成$\frac{n}{2}$份并拿走$k$份”是一回事一样。同时，$\omega_n^{\frac{n}{2}}$恰好是走半圈，在复平面上走半圈就是一次关于原点的中心对称，所以又有$\omega_n^{k+\frac{n}{2}}=-\omega_n^k$。

​ 由此上面两个式子可以化简成：

​ $A(\omega_n^k)=A_1(\omega_{\frac{n}{2}}^{k})+\omega_n^kA_2(\omega_{\frac{n}{2}}^k)$

​ $A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_{\frac{n}{2}}^{k})-\omega_n^kA_2(\omega_{\frac{n}{2}}^k)$

​ 两个表达式有了相同的组成部分，只需要完成$A_1(\omega_{\frac{n}{2}}^{k})$和$A_2(\omega_{\frac{n}{2}}^k)$，就可以求得$A(\omega_n^k)$和$A(\omega_n^{k+\frac{n}{2}})$的计算就可以完成$A(\omega_n^k)$和$A(\omega_n^{k+\frac{n}{2}})$的计算。这是一个递归求解的过程，假设求$A(\omega_n^{1..n})$需要计算$f(n)$次，就有$f(n)=2\cdot f(\frac{n}{2})+n$，而$f(1)=1$，所以$f(n)=n \log_2 n + n$，复杂度$O(n \log n)$。

离散傅里叶逆变换(IDFT)

​ 设刚才DFT的结果为点值表达式$(y_0,y_1,...,y_{n-1})$，它对应的系数表达式就是原来的$(a_0,a_1,...,a_{n-1})$

​ 构造以$(y_0,y_1,...,y_{n-1})$为系数的多项式

​ $B(x)=y_0+y_1x+...+y_{n-1}x^{n-1}$

​ 接下来，将单位根的倒数（$\omega_n^0,\omega_n^{-1},...,\omega_n^{1-n}$）代入$B(x)$得到点值表达式$(z_0,z_1,...,z_{n-1})$。（单位根的倒数恰好也是单位根的共轭复数）接着我们来观察$z$和$a$之间的联系。

​ $\begin{aligned} z_k &= \sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i\\ &= \sum_{i=0}^{n-1}[\sum_{j=0}^{n-1}a_j \cdot (\omega_n^{i})^j](\omega_n^{-k})^i\\ &= \sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j \cdot (\omega_n^{j-k})^i\\&= \sum_{j=0}^{n-1}a_j \cdot [\sum_{i=0}^{n-1}(\omega_n^{j-k})^i]\\\end{aligned}$

​ 后半部分为等比数列求和。要对公比是否为1进行分类讨论：

​ i) 当$j \neq k$时：

​ $\begin{aligned} \sum_{i=0}^{n-1}(\omega_n^{j-k})^i &= \dfrac{(\omega_n^{j-k})^n-1}{\omega_n^{j-k}-1}\\&= \dfrac{(\omega_n^{n})^{j-k}-1}{\omega_n^{j-k}-1}\\&= \dfrac{1-1}{\omega_n^{j-k}-1}\\&= 0\\ \end{aligned}$

​ ii) 当$j=k$时：

​ $\begin{aligned} \sum_{i=0}^{n-1}(\omega_n^{j-k})^i &= \sum_{i=0}^{n-1}1^i\\&=n\\ \end{aligned}$

​ 综上

​ $z_k=n \cdot a_k$

​ 所以

​ $a_k=\dfrac{z_k}{n}$

​ 奇妙之处就在于，我们通过代入单位根的倒数恰好消掉了大部分复杂的式子，这使得我们只需要在DFT的基础上略作修改就可以完成IDFT：我们只要把点值表达式当作系数表达式，把单位根取倒数后再做一次$DFT$，最后除以$n$，就得到了系数表达式。

数论变换(Number Theoretic Transforms, NTT)

​ FFT需要用复数运算，存在精度问题。NTT是模意义下的FFT。它的思想核心就是模数的原根来代替单位根。

​ 在FFT中，我们选择单位根$\omega_n^k$是由于单位根具有如下几个我们所需要的性质：

​ 1） $(\omega_n^k)^n=1$
​
​ 2） $\omega_{2n}^{2k}=\omega_n^k$
​
​ 3） $\omega_n^{k+\frac{n}{2}}=-\omega_n^k$

​ 4）$n$个单位根互不相同

​ 在模意义下，我们有阶的概念。满足$a^n \equiv 1 \pmod P$的最小的$n$被称为$a$在模$P$意义下的阶（暂时假定$P$为质数）。我们知道费马小定理：

​ $a^{P-1} \equiv 1 \pmod P$

​ 当$a$在模$P$意义下的阶为$P-1$时，我们便称$a$为模数$P$的原根，记作$g$。常见地，模数998244353的原根是3。

​ 由此，我们可以用$g_n^k$来代替原来的$\omega_n^k$。令：

​ $g_n^k=g^{\frac{P-1}{n} \cdot k}$

​ （注意，这同时要求$n$是$P-1$的约数。我们之前还假定了$n$是2的幂。恰好998244353是这样一个满足条件的约数。）

​ 由此可知性质一：

​ $(g_n^k)^n \equiv g^{\frac{P-1}{n}\cdot k n}\equiv (g^{P-1})^k \equiv 1 \pmod P$

​ 性质二：

​ $g_{2n}^{2k} \equiv g^{\frac{P-1}{2n} \cdot 2k} \equiv g^{\frac{P-1}{n} \cdot k} \equiv g_n^k \pmod P$

​ 性质三：

​ $g_n^{k+\frac{n}{2}} \equiv g_n^k \cdot g_n^{\frac{n}{2}}\equiv g_n^k \cdot g^{\frac{P-1}{n}\cdot \frac{n}{2}} \equiv g_n^k \cdot g^{\frac{P-1}{2}}$

​ 由于$(g^\frac{P-1}{2})^2 \equiv 1 \pmod P$，把它当作一个二次剩余方程可以观察到

​ $g^{\frac{P-1}{2}}=-1\pmod P$

​ 也就有

​ $g_n^{k+\frac{n}{2}}=-g_n^k$

​ 性质四：

​ 如果存在$r<s<P-1$使得$g^r \equiv g^s \equiv 1 \pmod P$，则$P \mid g^{s-r}$。于是$g^{(P-1)-(s-r)} \equiv 1 \pmod P$，这就与“原根的阶是$P-1$”矛盾了。

​ 综上，我们发现原根具有我们需要的单位根的性质。在FFT中用原根代替单位根，也就是NTT了。注意在IDFT中我们用到了单位根的倒数，这就对应着原根的乘法逆元。

代码实现

​ 根据刚才的理论，可以直接写出递归实现的代码。注意单位根可以在计算的过程中递推，而不用每次重新求。

void DFT(int* a, int n, bool IDFT){
	if(n == 1) return;
	for(int i = 0; i < n; ++i){
		if(i & 1){
			b[(n>>1)+(i>>1)] = a[i];
		}else{
			b[i>>1] = a[i];
		} 
	}
	for(int i = 0; i < n; ++i) a[i] = b[i];
	DFT(a, n>>1, IDFT);
	DFT(a+(n>>1), n>>1, IDFT);
	gn = ((IDFT) ? qpow(qpow(3,(P-1)/n),P-2) :qpow(3,(P-1)/n));
	w0 = 1;
	for(int k = 0; k < (n>>1); ++k){
		a1 = a[k], a2 = a[k+(n>>1)];
		tmp = mul(w0, a2);
		a[k] = add(a1, tmp);
		a[k+(n>>1)] = sub(a1, tmp);
		w0 = mul(w0, gn);
	}
}

​ 能不能迭代实现？这也就要求我们从叶节点出发向上完成计算。首先需要确定叶节点的值。在上述过程中，我们每一轮都会对节点奇偶分类。在最顶层的分类中，偶数被分到前面，奇数被分到后面。（下标从0开始）对于被分到前面的偶数，此时的第偶数个被分到前面，第奇数个被分到后面。如果说第一次看的是奇偶，那第二次看的就是偶数的偶数和偶数的奇数、奇数的偶数和奇数的奇数。仔细思考一下可以发现，这就是二进制的第一位和第二位。最终所有节点是按照二进制的倒序排序的。

​ 在计算二进制倒序的时候，不需要暴力计算，而是可以利用历史信息：$r_i$通过$r_{i/2}$来转移。

	for(int i = 1; i < N; ++i){
		r[i] = (r[i>>1]>>1) | ((i&1) ? (N>>1) : 0);
	}

​ 迭代实现的代码如下：

inline void DFT(int* a, int n, bool IDFT){
	for(int i = 1; i < n; ++i){
		if(r[i] < i) swap(a[r[i]],a[i]);
	}
	for(int i = 2; i <= n; i = i<<1){
		gn = ((IDFT) ? qpow(qpow(3,(P-1)/i),P-2) :qpow(3,(P-1)/i));
		for(int j = 0; j < n; j = j+i){
			w0 = 1;
			for(int k = 0; k < (i>>1); ++k){
				a1 = a[j+k], a2 = a[j+k+(i>>1)];
				tmp = mul(w0, a2);
				a[j+k] = add(a1, tmp);
				a[j+k+(i>>1)] = sub(a1, tmp);
				w0 = mul(w0, gn);
			}
		}
	}
}

任意模数NTT

​ 上面假定了$P$是一个非常特殊的模数，它首先是指数，$P-1$还是2的幂次的倍数，它还有一个比较简单的原根。

​ 但假设要求模数任意，怎么解决？在解决这个问题时，人们发现，我们并不需要对我们的算法进行某种修改，而是对我们的结果进行一些变换就可以解决这个问题。对于最终输出的某个答案$x$，假设我们知道$x$在其他“特殊模数”下的答案，我们可以通过解同余方程组来求出$x$。

​ 我们选定满足普通NTT要求的三个模数$P_1,P_2,P_3$，分别求得$x$对应的结果$x_1,x_2,x_3$，于是我们就有

​ $\left \{ \begin{matrix} x ≡ x_1 \pmod {P_1} \\ x ≡ x_2 \pmod {P_2} \\ x ≡ x_3 \pmod {P_3} \\ \end{matrix}\right.$

我们通过中国剩余定理可以得到$x \equiv ans \pmod{ P_1P_2P_3}$，假若$P_1P_2P_3$已经大于给定的数据的最大范围，那么我们得到的答案就可以被认为是正确的了。

（还没有实践过）