FFT&NTT

发表于 | 分类于 | 阅读次数:

FFT（快速傅里叶变换, Fast Fourier Transformation)

在算法竞赛中的主要应用是加速多项式运算。

以多项式乘法为例：朴素算法需要 $O(N^2)$ 的时间复杂度，而经过 FFT 优化后只需要 $O(Nlog_2N)$ 的时间复杂度。

数学基础

FFT 作为一个数学算法，比起复杂的数据结构，其编程较为简单，但是它对于数学的要求比较高，要想理解此算法必须先有一定的数学基础。

多项式的表示

多项式有两种表示方法：

  • 系数表达法
  • 点值表达法

其实有些像函数的表达，一二次函数为例，可以用解析式表示，也可以用 3 个以上的点来表示。

系数表达

对于次数界为 $n$ 的多项式 $A(x)=\sum_{j=0}^{n-1}a_jx^j$ ，其系数表达是一个由系数组成的向量 $a=(a_0,a_1,\dots,a_{n-1})$ 。

点值表达

一个次数界为 $n$ 的多项式 $A(x)$ 的点值表达就是一个由 $n$ 个点值所组成的集合：

${(x_0,y_0),(x_1,y_1),\dots ,(x_{n-1},y_{n-1})}$

使得集合中无重复元素，且

$y_k=A(x_k),k\in {k\in \mathbb{Z}|0\leq k\leq n-1}$

一个多项式可以有很多不同的点值表达，以为可以采用 n 个不同的点 $x_0,x_1,\dots,x_{n-1}$ 构成的集合作为这种表示方法的基。

插值：求值计算的逆（从一份多项式的点值表达确定其系数表达的形式）

当差值得多项式的次数界等于已知的点值对的数目，插值才是明确的。

FFT计算过程

有了前面的两个概念，我们的运算流程就十分明确了：

按照朴素算法，时间复杂度： $O(N^2)$

而按照 FFT 算法，先进行一次 DFT（求值, Discrete Fourier Transform），时间复杂度： $O(Nlog_2N)$ ，将系数表达法转化成点值表达法；之后进行点值乘法，时间复杂度： $O(N)$ ；之后再进行一次 IDFT（插值, Inverse Discrete Fourier Transform），将点值表达法转换为系数表达法。

单位复数根

至此，大体的框架已经建成了，但在具体实现还要引入一个单位复数根的概念。

首先，介绍一下复数：

  • 我们规定 $i^2=-1$ ，即 $\sqrt{−1}=i$ ，其中 $i$ 就是虚数域中的单位1；
  • 复数集 $\mathbb{C}=\mathbb{R}\cup \mathbb{I}$ ，是现在已知的最大数集；
  • $\forall x \in \mathbb{C}, \exists a,b \in \mathbb{R}$ 使 $x=a+bi$ ，其中 $a$ 叫做实部，$bi$ 叫做虚部；
  • 实数用一维的数轴表示，而复数需要二维的复平面来表示（类似平面直角坐标系），x轴便是正常的实数数轴，y轴是虚数轴. 坐标为 $(a,b)$ 的点P表示复数 $a+b*i$ 。

下面进入正题：

  • n次单位复数根：满足 $\omega^n=1$ 的复数 $\omega$ ；
  • n次单位复数根有恰好有 $n$ 个：对于 $k=0,1,\dots,n-1$ ，这些根是 $e^{2\pi ik/n}$ ；
  • 根据欧拉公式 $e^{iu}=\cos(u)+i\sin(u)$ $n$ 个单位复数根均匀分布在单位圆的圆周上；
  • 主n次单位根： $\omega_n=e^{2\pi i/n}$ ，其他所有n次单位复数根都是 $\omega_n$ 的幂次；

根据以上性质，我们可得到以下两个基本性质：

消去引理

对于任何整数 $n\geq0,k\geq0$ 且 $d>0$ 有 $\omega_{dn}^{dk}=\omega_n^k$

FFT 中需要用到的是 $(\omega_n^k)^2=\omega_{n/2}^k$

折半引理

如果 $n>0$ 为偶数，那么 $n$ 个n次单位复数根的平方的集合就是 $n/2$ 个 $n/2$ 次单位复数根的集合。

即：$n>0$ 为偶数 n次单位复数根集合 ${(\omega_n^k)^2|k\in {k\in \mathbb{Z}|0\leq k\leq n-1}}={\omega_{n/2}^k|k\in{k\in \mathbb{Z}|0\leq k\leq n/2-1}}$ ，即左边集合元素重复出现两次；

折半引理对于用分治策略来对多项式的系数与点值表达法是至关重要的，因为它保证了递归子问题的规模只是递归调用前的一半。

实现思路

其实在数学基础里已经把大体的实现思路有所说明了，这样写是可行的，但是问题在于当数据量很大的时候如果用递归写的话会发生爆栈，所以我们在这里先研究一下递归搜索树，然后进行数组模拟。

dfs树

我们通过研究系数出现的顺序可以的出变换的规律，于是我们考虑进行二进制 DP

  • 当要处理一个 $5$ 位二进制数 $\overline{abcde}$ 时，小于 $\overline{abcde}$ 的数已经处理好二进制反转了；
  • 所以 $\overline{abcde}>>1=\overline{abcd}$ 已有二进制反转 $\overline{0abcd}$；
  • 所以只需处理 $\overline{e}$ 这一位即可；
  • 那么只需 $\overline{0abcd}$ | ( ($\overline{abcde}$&1)<<(bit-1)) 即可；
  • 考虑到 DP 从 $0$ 开始无需反转， $1>>1=0$ 。

int(a[i].real()/n+0.5)

##代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
#include<iostream>
#include<cstdio>
#include<cmath>
#include<complex>
using namespace std;

typedef complex<double> C;

const int MAXN=262145;
const double PI=acos(-1);
int pos[MAXN],bit,n,m;
C a[MAXN],b[MAXN];

void FFT(C * A,int type)
{
	int i,j,k;
	for(i=1;i<n;i++)
		if(i<pos[i])
			swap(A[i],A[pos[i]]);
	for(i=1;i<n;i<<=1)
	{
		C wn(cos(PI/i),type*sin(PI/i));
		for(j=0;j<n;j+=i<<1)
		{
			C w(1,0);
			for(k=0;k<i;k++,w*=wn)
			{
				C x=A[j+k],y=w*A[i+j+k];
				A[j+k]=x+y,A[i+j+k]=x-y;
			}
		}
	}
}
int main()
{
	scanf("%d%d",&n,&m);
	int i,x;
	for(i=0;i<=n;i++)
		scanf("%d",&x),a[i]=x;
	for(i=0;i<=m;i++)
		scanf("%d",&x),b[i]=x;
	m+=n;
	for(n=1;n<m;n<<=1) ++bit;
	for(i=1;i<n;i++)
		pos[i]=(pos[i>>1]>>1)|((i&1)<<(bit-1));
	FFT(a,1),FFT(b,1);
	for(i=0;i<=n;i++)
		a[i]*=b[i];
	FFT(a,-1);
	for(i=0;i<=m;i++)
		printf("%d ",int(a[i].real()/n+0.5));
	return 0;
}