2019-03-08

算法笔记

23 分钟读完 (大约 3382 个字)

FFT详解

用于计算一类朴素卷积问题。

笔者学习的是这份博客
内容中可能有很多相同之处，敬请谅解。

现在要计算两个一元$n$次多项式$F(x)$与$G(x)$的乘积，如何计算？
前置知识：多项式的表示方法
一. 系数表示法
对于一个$n$次多项式$F(x)$，它可以被表示成

$F(x) = a_nx^n+a_{n-1}x^{n-1}+...+a_1x^1+a_0x^0.$

更加形式化的来说，它可以表示成

$F(x) = \sum_{i=0}^{n} a_ix^i.$

举个例子，2次多项式，其中$a_0=1,a_1=2,a_2=1$
那么$F(x)=1x^2+2x+1.$这样即可通俗的表示出一个$n$次多项式。
二. 点值表示法
~~众所周知，两个点确定一个一次函数，三个点确定一个二次函数。~~
所以，$n+1$个点确定一个一元$n$次多项式。
所以我们可以通过$n+1$个点来表示它。
那么相乘之后的点值如何计算？
比如说两个2次多项式$F(x)=x^2+2x+1$（红色）与$G(x)=3x^2-4x-2$（蓝色），它们的图像如图所示：

那么观察图中$x=1$时的情况。此时$F(1)=4,G(1)=-3.$
所以，显然，$F(1)\times G(1)=-12$.也就是说在$Z=FG$这一多项式内，带入$1$，得到的结果是$-12$.
等等，好像有哪里不对。如果说$Z=FG$的话，那么Z的次数应该是$2n$.
那$Z$需要$2n+1$个点来确定。但是原来只需要$n+1$个点，咋办？
很简单，在原来的多项式里每个都多加$n$个点即可。反正多项式已知。
这样就可以用点值来进行操作。也就是说先转成点值，再一乘，再转回来，就是计算流程。
但是好像还是很慢。那么如何优化呢?
复数部分
复数，即形如$a+bi$的数，其中$i^2=-1.$ $a$称为实部，$bi$称为虚部。
或者说：在一个数轴上（只有x轴），我们可以表示出任何实数。
那么，多加一维（y轴），也就是类似于平面直角坐标系一样，我们就可以表示出任意一个复数。
所以我们把这个坐标系叫做复平面，其中x轴称为实轴，y轴称为虚轴。
复数运算
复数相加：实部相加，虚部相加，例如

$(a+bi)+(c+di)=(a+c)+(b+d)i.$

复数相减：同理。

$(a+bi)-(c+di)=(a-c)+(b-d)i.$

复数相乘：像一次多项式一样相乘。 注意$i^2=-1$.

$(a+bi)(c+di)=ac+(ad+bc)i-bd=(ac-bd)+(ad+bc)i.$

复数相除：
相信大家都学过共轭根式。同样的，复数也有共轭。
即：$a+bi$的共轭为$a-bi$。
这两个复数卡乘在一起一定是个实数。即

$(a+bi)(a-bi)=a^2-(bi)^2=a^2+b^2.$

所以再除的时候，将分子分母同乘分母的共轭，就可以将分母有理化。
即

$\frac{a+bi}{c+di}=\frac{(a+bi)(c-di)}{c^2+d^2}=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}i.$

复数逆元：

$\frac{1}{a+bi}=\frac{a}{a^2+b^2}-\frac{b}{a^2+b^2}i.$

同样的，复数运算在复平面内也有一定规律可循。
考虑在复平面内的两个复数：（借张图）

表示的是复数$(1+4i)$与复数$(3+2i)$相乘所得的结果：$(-5+14i)$。
设其中$(5,0)$点为位置$P$，则$\angle{POC}=\angle{BOA}.$
还有：$\overline{OB}\times\overline{OC}=\overline{OA}.$
第二个证明：勾股定理。先把$i$消掉。
$\overline{OB}^2=a^2+b^2.$ $\overline{OC}^2=c^2+d^2.$
$\overline{OB}^2\times \overline{OC}^2=(a^2+b^2)(c^2+d^2)=a^2c^2+a^2d^2+b^2c^2+b^2d^2.$
$\overline{OA}^2=(ac-bd)^2+(ad+bc)^2=a^2c^2+a^2d^2+b^2c^2+b^2d^2.$得证。
我们将复数中，复数向量的长度称为模长，向量与x轴正方向的夹角称为幅角。
根据上面的东西：复数相乘时，模长相乘，幅角相加。
单位根
一个n次的单位根即为方程$x^n=1$的复数解。
考虑这样一个图：

其中圆上的所有复数模长都是1，这个圆称为单位圆。
考虑$|x|$的取值范围：
如果$|x|<1$，那么$|x^n|<1$（模长\*模长，越乘越短）如果$|x|>1$，那么$|x^n|>1$（模长*模长，越乘越长）
所以一定有$|x|=1$，即点一定在这个单位圆上。
那么还有一条，就是一个n次的单位根共有n个，并且从(1,0)开始，每次逆时针走$\frac{1}{n}$个圆周。分别称为$\omega_{n}^{0},\omega_{n}^{1},…\omega_{n}^{n-1}.$（逆时针编号！！！）
其实还有$k\geq n$和$k<0$的单位根，本质上就模个n就可以了，这里不再多加考虑。
有关单位根也有很多性质：

$\omega_{n}^{0}=1$
$\omega_{n}^{k}=(\omega_{n}^{1})^k$，很好理解，就是每次转$\frac{1}{n}$圆周。
$\omega_{n}^{j} \times \omega_{n}^{k}=\omega_{n}^{j+k}.$
证明：很简单，代入性质2即可。
根据3，有： $\omega_{n}^{j+k}=\omega_{n}^{j} \times \omega_{n}^{k}.$ $(\omega_{n}^{k})^{-1}=\frac{1}{\omega_{n}^{k}}=\frac{\omega_{n}^{0}}{\omega_{n}^{k}}=\omega_{n}^{-k}=\omega_{n}^{n-k}.$ （根据上面模个n的说法）
$\omega_{pn}^{pk}=\omega_{n}^{k}.$ （将蛋糕分成pn份取pk份=将蛋糕分成n份取k份）
$\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}.$ （相当于绕了半圈$\rightarrow$取了个相反数）
$(\omega_{n}^{k})^j=(\omega_{n}^{j})^k.$ （大小为k，共有j个=大小为j，共有k个）
…差不多就这些了

$FFT$流程：$DFT,IDFT$
$DFT$加速版本：分治
首先，你需要~~控制得当~~将整个n-1次多项式（保证n是2的正整数次幂）按系数奇偶拍成两部分：

$F(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+(a_1x^1+a_3x^3+...+a_{n-1}x^{n-1}).$

再设两个n/2项多项式：

$FL(x)=a_0+a_2x+...+a_{n-2}x^{\frac{n}{2}-1}$ $FR(x)=a_1+a_3x+...+a_{n-1}x^{\frac{n}{2}-1}$

那么有$F(x)=FL(x^2)+x*FR(x^2).$（不懂可以自己设一下代入进去）
下面，我们来代入单位根，看看会有什么神奇的事情发生。
比如说我们代一个$\omega_{n}^{k}(k<\frac{n}{2}):$

$F(\omega_{n}^{k})=FL((\omega_{n}^{k})^2)+\omega_{n}^{k}\times FR((\omega_{n}^{k})^2)$

其中，$(\omega_{n}^{k})^2=\omega_{n}^{2k}=\omega_{\frac{n}{2}}^{k}.$然后代入：

$F(\omega_{n}^{k})=FL(\omega_{\frac{n}{2}}^{k})+\omega_{n}^{k}\times FR(\omega_{\frac{n}{2}}^{k}).$

然后呢？有啥用？
不难发现，如果我们知道了$FL(x),FR(x)$在x取$\omega_{\frac{n}{2}}^{0},\omega_{\frac{n}{2}}^{1}…\omega_{\frac{n}{2}}^{\frac{n}{2}-1}$下的所有值，我们就可以知道$F(x)$在x取$\omega_{n}^{0},\omega_{n}^{1}…\omega_{n}^{\frac{n}{2}-1}$下的所有值（代回原式）。
但是这还不是全部啊，剩下那个部分怎么办？
没关系，考虑$k<\frac{n}{2}$，代入$\omega_{n}^{k+\frac{n}{2}}$时的情况。

$F(\omega_{n}^{k+\frac{n}{2}})=FL((\omega_{n}^{k+\frac{n}{2}})^2)+\omega_{n}^{k+\frac{n}{2}}\times FR((\omega_{n}^{k+\frac{n}{2}})^2).$

其中，$(\omega_{n}^{k+\frac{n}{2}})^2=\omega_{n}^{2k+n}.$

故

$F(\omega_{n}^{k+\frac{n}{2}})=FL(\omega_{n}^{2k+n})+\omega_{n}^{k+\frac{n}{2}}\times FR(\omega_{n}^{2k+n}).$ $F(\omega_{n}^{k+\frac{n}{2}})=FL(\omega_{n}^{2k})+\omega_{n}^{k+\frac{n}{2}}\times FR(\omega_{n}^{2k}).$

（模性质）

$F(\omega_{n}^{k+\frac{n}{2}})=FL(\omega_{\frac{n}{2}}^{k})+\omega_{n}^{k+\frac{n}{2}}\times FR(\omega_{\frac{n}{2}}^{k}).$

（性质3 延展）

$F(\omega_{n}^{k+\frac{n}{2}})=FL(\omega_{\frac{n}{2}}^{k})-\omega_{n}^{k}\times FR(\omega_{\frac{n}{2}}^{k}).$

（取反性质）

然后观察这个式子和前面式子的区别，发现就是把+号改成了-号。也就是说我们可以通过那组数据来算出当$k>\frac{n}{2}$时的情况。这就是加速DFT的原理。
IDFT
至此，我们已经可以将其转化为点值，那么如何还原回来呢？使用IDFT插值。
考虑一个向量$c_0,c_1…c_{n-1}$，满足 $c_k=\sum _{i=0}^{n-1}(\omega_{n}^{-k})^i\times y_i.$
其中$y_i$表示乘出来的点的y坐标。
那么，这个式子也就代表一个以$y$为系数的方程在取$x=\omega_{n}^{-k}$时的复数解。
我们尝试着来化简一波。

$c_k=\sum _{i=0}^{n-1}(\omega_{n}^{-k})^i\times y_i$ $c_k=\sum _{i=0}^{n-1}(\omega_{n}^{-k})^i(\sum_{j=0}^{n-1}a_j(\omega_{n}^{i})^j)$ $c_k=\sum _{i=0}^{n-1}(\omega_{n}^{-k})^i(\sum_{j=0}^{n-1}a_j(\omega_{n}^{j})^i)$ $c_k=\sum _{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_{n}^{j})^i(\omega_{n}^{-k})^i)$ $c_k=\sum _{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_{n}^{j-k})^i)$ $c_k=\sum _{j=0}^{n-1}a_j\color{red}{(\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^i)}$

发现后面标红的东西就是一个等比数列，那么我们可以直接计算。即

$\sum_{i=0}^{n-1}X^i=\frac{X^n-1}{X-1}.$

我们设这个东西为公式$Z.$
于是就有

$c_k=\sum _{j=0}^{n-1}a_jZ(\omega_{n}^{j-k}).$

观察函数$Z$的相关性质。
赫然发现：我们带进去的是单位根！单位根啥性质来着？$X^n=1!$
也就是说当$X≠1$，也就是$k≠0$时，原式$=0.$ （分子为0）
那么$X=1$时，也就是$k=0$时呢？这时候这个求和公式就不管用了，带回原式算，发现由于$1^i=1$（i是自然数），那么$Z(x)=n.$
所以说，什么情况下带进去的值$k=0$呢？显然是$j-k$的时候！此时单项$=na_k$.
否则的话会乘一个0，答案固然也是0.
所以有： $c_k=na_k$ （飞跃性的一步！！！）
回顾整个式子，发现我们算$c$的时候连$a$的影子都没有，现在我们竟然推出了一个和$a$有关的关系式！太神奇了！这就是$IDFT$的原理。
至此，$\text{FFT}$的重要部分已经全部介绍完毕，下面进入代码实现与细节处理部分。

代码实现
朴素DFT:

complex < double > omega(const int n , const int k) { // 单位根
	complex < double > ret = {cos(2 * Pi / n * k) , sin(2 * Pi / n * k)};
	if(!inv) return ret;
	else return conj(ret);
}

void DFT(complex < double > a , const int n) {
	if(n == 1) return ;
	static complex < double > buf[MAXN];
	const int m = n / 2;
	for(int i = 0; i < m; i++) {
		buf[i] = a[(i << 1)];
		buf[i + m] = a[(i << 1 | 1)];
	}
	copy(buf , buf + n , a);
	complex < double > *a1 = a , *a2 = a + m;
	DFT(a1 , m); DFT(a2 , m);
	for(int i = 0; i < m; i++) {
		complex < double > w = omega(n , i);
		buf[i] = a1[i] + w * a2[i];
		buf[i + m] = a1[i] - w * a2[i];
	}
	copy(buf , buf + n , a);
}

发现这样的效率显然不行，那么考虑迭代形式。
首先观察分治到最后的二进制位表示：

000 001 010 011 100 101 110 111
 0   1   2   3   4   5   6   7
 0   2   4   6 - 1   3   5   7
 0   4 - 2   6 - 1   5 - 3   7
 0 - 4 - 2 - 6 - 1 - 5 - 3 - 7
000 100 010 110 001 101 011 111

然后…我们惊奇的发现最后的二进制位就是之前的二进制位倒过来！TQL!
接着考虑空间优化。
不难发现，由于DFT的两个式子十分相近，所以我们可以只算一次，同时通过技巧性的操作直接覆盖原数组，而不再新开一个。这个操作被称为蝴蝶操作。
蝴蝶操作的代码：

for(int l = 2; l <= n; l <<= 1) {
			int m = l / 2;
			for(complex < double > *p = a; p != a + n; p += l) {
				for(int i = 0; i < m; i++) {
					complex < double > t = omega[n / l * i] * p[m + i];
					p[m + i] = p[i] - t;
					p[i] += t;
				}
			}
		}

然后穿起来，与以前的FFT模板结合即可。

#include <bits/stdc++.h>
#define dbg(x) cerr << #x " = " << x << "\n"
#define INF 0x3f3f3f3f
/* do not give up ! */
/* try your best! */
/* Read the meaning clearly! */
typedef long long LL;
typedef long double ld;
typedef unsigned long long ULL;

using namespace std;

template < typename T > inline void inp(T& t) {
	char c = getchar(); T x = 1; t = 0; while(!isdigit(c)) {if(c == '-') x = -1; c = getchar();}
	while(isdigit(c)) t = t * 10 + c - '0' , c = getchar();t *= x;
}
template < typename T , typename... Args > inline void inp(T& t , Args&... args) {inp(t); inp(args...);}
template < typename T > inline void outp(T t) {
	if(t < 0) putchar('-') , t = -t; T y = 10 , len = 1;
	while(y <= t) y *= 10 , len++; while(len--) y /= 10 , putchar(t / y + 48) , t %= y;
}
template < typename T > inline T mul(T x , T y , T MOD) {x=x%MOD,y=y%MOD;return ((x*y-(T)(((ld)x*y+0.5)/MOD)*MOD)%MOD+MOD)%MOD;}

const int MAXN = 2097152 + 10;
const double Pi = acos(-1.0);
int ans[MAXN];
struct FastFourierTransform {
	complex < double > omega[MAXN] , omegaInverse[MAXN];

	void init(const int n) {
		for(int i = 0; i < n; i++) {
			omega[i] = complex < double > (cos(2 * Pi / n * i) , sin(2 * Pi / n * i));
			omegaInverse[i] = conj(omega[i]);
		}
	}
	void Transform(complex < double > *a , const int n , const complex < double > *omega) {
		int k = 0;
		while((1 << k) < n) ++k;
		for(int i = 0; i < n; i++) {
			int t = 0;
			for(int j = 0; j < k; j++) if(i & (1 << j)) t |= (1 << (k - j - 1));
			if(i < t) swap(a[i] , a[t]);
		}
		for(int l = 2; l <= n; l <<= 1) {
			int m = l / 2;
			for(complex < double > *p = a; p != a + n; p += l) {
				for(int i = 0; i < m; i++) {
					complex < double > t = omega[n / l * i] * p[m + i];
					p[m + i] = p[i] - t;
					p[i] += t;
				}
			}
		}
	}

	void DFT(complex < double > *a , const int n) {
		Transform(a , n , omega);
	}
	void IDFT(complex < double > *a , const int n) {
		Transform(a , n , omegaInverse);
		for(int i = 0; i < n; i++) a[i] /= n;
	}
}FFT;
int main() {
	int nn , mm , x; inp(nn , mm);
	nn++ , mm++;
	int n = 1; while(n < nn + mm) n <<= 1;
	complex < double > a[MAXN] , b[MAXN];
	for(int i = 0; i < nn; i++) {
		inp(x); a[i].real(x);
	}
	for(int i = 0; i < mm; i++) {
		inp(x); b[i].real(x);
	}
	FFT.init(n);
	FFT.DFT(a , n);
	FFT.DFT(b , n);

	for(int i = 0; i < n; i++) a[i] *= b[i];

	FFT.IDFT(a , n);
	for(int i = 0; i < nn + mm - 1; i++) {
		ans[i] = static_cast < int > (floor(a[i].real() + 0.5));
		cout << ans[i] << " ";
	}
	return 0;
}

至此，该代码已经可以效率很高地通过P3803一题！完结撒花！
FFT例题将会另外放在一个帖子~

# FFT, 数学

FFT详解

评论

分类

Your browser is out-of-date!