浅谈FFT

多项式的表示

首先学习两种多项式的表达法

系数表达法
点值表达法

系数表达法

比如一个 \(n\) 次多项式 \(A(x)\) 他有 \(n + 1\) 项于是设每一项的系数为 \(a_i\) 则有 \(A(x) = \sum _{i = 0}^ n a_i x^ i\)

利用这种方法计算多项式乘法复杂度为 \(O(n^2)\) 即利用（第一个多项式中每个系数都需要与第二个多项式的每个系数相乘）

用这种方法，计算 \(A(x_0)\) 的值是非常好用的，我们使用霍纳法则即可。

点值表达法

我们知道如果给出 \(n + 1\) 个互不相同的点，他们确定一个 \(n\) 次多项式，这个结论是显然的

我们在下文中对于次数的规定稍微放开，规定 A 的次数界（即最大次数 < 次数界）为 \(K_a\) 而同理，有 \(B_k\) 。通常的，把 \(A(x)*B(x)\) 的次数界规定为 \(K_a + K_b\)

对于许多多项式的操作，点值表达式十分方便，例如有一个点值表达式 A

\[\{\{x_0, y_0\}, \{x_1, y_1\}\cdots\{x_{n - 1}, y_{n-1}\}\} \]

还有一个点值表达式 B

\[\{\{x_0, y'_0\}, \{x_1, y_1'\}\cdots\{x_{n - 1}, y'_{n-1}\}\} \]

那么 \(A + B\) 的表达就是

\[\{\{x_0, y_0 + y'_0\}, \{x_1, y_1 + y'_1\}\cdots\{x_{n - 1}, y_{n-1} + y'_{n-1}\}\} \]

表达法END

复数浅讲

入门

上过高中的请跳过，我是初中生，会讲锅的。

定义一个复数 \(u = a + bi\) 其中 \(i = \sqrt{-1}\) 。你其实可以按向量那么理解。

然后我们复数乘法显得巧妙

比如说一个

复数 \(b\) * 复数 \(c\) 可以理解为就是旋转，比如这个图

看到了吗，旋转，模长相乘。

正经的说

复数包含了实部和虚部。可以表示为一个虚数 \(z\) 可以表示为 \(z = a + bi (a,b \in \R)\) 其中 \(a\) 为 \(z\) 的实部， \(bi\) 为 \(z\) 的虚部.

复数可以表示为平面向量 \((a, b)\) 。类似于共轭根式，我们定义共轭复数 \(z = a - bi\) 有一个复数与其共轭复数的乘积为实数。

定义复数的四则运算如下：

加法 \((a + bi) + (c + di) = (a + c) + (b + d)i\)
减法 \((a + bi) -(c + di) = (a - c) + (b - d)i\)
乘法 \((a + bi) * (c + di) = (ac - bd) + (ad + bc)i\)
除法 \(\frac {a + bi}{c + di} = \frac {(a + bi)* (c - di)}{(c + di) * (c - di)} = \frac {(ac + bd) + (ad - bc)i}{c + d}\)

这几个很基本了。

关于本源单位根

我们还有一类特殊的向量——\(n\) 次本源单位根 \(\omega\)，表示什么，就是方程 \(x^n = i\) 的解，这就很显然了

有几个关于本源单位根的定理

消去引理： \(\omega ^{dk}_{dn} = \omega ^k_n\)
对于任意偶数 \(n > 0\) 有 \(\omega ^{\frac {n}{2}}_n = \omega^1_2 = -1\)
折半引理：如果 \(n\) 为偶数，则有 \((\omega^k_n)^2 = (\omega ^{k + n/2}_n)^2\) 但是 \(\omega^k_n = -\omega ^{k + n/2}_n\)从复数相乘的角度上看这是显然的
求和引理：对于任意整数 \(n\ge1\) 且不能被 \(n\) 整除的 \(k\),有 \(\sum _{j=0}^{n-1} (w^k_n )^j = 0\) 为啥，由于我们可以用等比数列求和公式来计算为 \(\sum ^{n-1}_0 (\omega_n^k)^j = \frac {(w^k_n)^n - 1}{w_n^k -1} = 0\)

显然的，我们可以发现这几个引理是正确的。

对于所有的次数界，不妨把他们看为是 2 的幂次，若不然可以补零解决，~~其实也有不为二的幂次的方法，只是我不会~~

\(DFT\)

~~离散小波变换~~ 离散傅里叶变换。计算次数界为 \(n\) 的的多项式 \(A(x) = \sum _0^{n-1} a_j x^j\) 我们不妨带入每个单位根得到 \(y_k = A(\omega_n^k) = \sum _{j = 0} ^ {n-1}a_j \omega _n ^{kj}\) 。我们发现计算这个还是 \(O(n ^ 2)\) 的，依然没有时间上的优化。这就是离散傅里叶变换，我们由于单位根的一些性质，就可以做到 \(O(n \log n)\) 的时间来计算一个数的点值表达式

\(FFT\)

快速傅里叶变换、利用单位复数根的性质，我们可以在 \(O(n\log n)\) 的时间内计算出 \(DFT_n\)

不妨假设都为 \(n\) 以次数界，且 \(n\) 为 2 的幂次。

不妨计算 \(A(\omega_n^k)\) 然后按下标的奇偶分类为 \(A_1A_2\)

\[A(\omega_n^k) = A_1(w_n^{2k}) + \omega_n^k A_2(\omega_n^{2k}) \]

我们通过观察有

\[A(\omega_n^{k + n/2}) = A_1(w_n^{2k}) - \omega_n^k A_2(\omega_n^{2k}) \]

所以，我们在计算上一个时，可以 \(O(1)\) 得到下一个于是就很明显了

带入 \(\omega ^k_n\) 进入式子 \(A(x) = \sum _0^{n-1} a_j x^j\) 中，我们得到 \(A(\omega ^k_n)= \sum_0^{n - 1} a_j \omega ^{kj}_n\)

考虑展开这个式子，并让下标按奇偶分类有 \(A(\omega ^k_n) = a_o\omega ^0_n + a_2 \omega ^{2k}_n + \ldots + a_{n - 2}\omega ^{k(n - 2)}_n + a_1\omega ^k_n + \ldots + a_{n - 1} \omega ^{k(n - 1)}_n\)

不妨把前边写为 \(A_1(x) = a_ox^0 + a_2 x^1 + \ldots + a_{n - 2}x ^{\frac{n - 2} 2}\) 后边写为 \(A_2(x) = a_1x^0 + \ldots + a_{n - 1} x ^{\frac{n - 2} 2}\) 于是，我们的 \(A(\omega ^k_n) = A_1(\omega ^{2k}_n) + \omega ^{k}_n A2(\omega ^{2k}_n)\)，不妨观察，假如我们计算 \(A(\omega ^{k + \frac{n}{2}}_n)\) 那又如何，很明显 \(A(\omega ^{k + \frac{n}{2}}_n) = A_1(\omega ^{2k}_n) - \omega ^{k}_n A2(\omega ^{2k}_n)\) 由折半引理，我们可以得出来这个式子。

所以我们在计算 \(A(\omega ^k_n)\) 的时候，可以一起计算出 \(A(\omega ^{k + \frac{n}{2}}_n)\) 这就很 nice 了。我们就可以得到一个递归的 \(FFT\) 然后我们完成了把一个系数表达式搞成了点值表达式的工作。

PS: 先不要管 type 的意义。

void SlowfFt(const 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];
	}
	SlowfFt(limit >> 1, a1, type);
	SlowfFt(limit >> 1, a2, type);
  // 做出单位根
	complex w(1, 0), wn(cos(2.0 * Pi / limit), sin(2.0 * Pi / limit) * type);
	for(int i = 0; i < limit >> 1; ++i, w = w * wn)
	{
    // 根据算出的答案直接出
		complex io = w * a2[i];
		a[i] = a1[i] + io, a[(limit >> 1) + i] = a1[i] - io;
	}
	return ;
}

\(DFT^{-1}\)

我们要把他的这个东西变成系数表达法。怎么搞，前方高能！

我们知道，把系数表达式转换为点值表达法的过程叫做求值，而这个逆过程叫做插值，下面引用一个定理：

插值多形式的唯一性：

证明：根据某个矩阵我们等价于矩阵方程

\[\begin{bmatrix}1 & x_0&x_0^2&\cdots&x_0^{n-1}\\1 & x_1&x_0^2&\cdots&x_1^{n-1}\\\vdots&\vdots&\vdots&\ddots&\vdots\\1 & x_{n-1}&x_0^2&\cdots&x_{n-1}^{n-1}\\\end{bmatrix} \begin{bmatrix}a_0\\a_1\\\vdots\\a_{n-1}\\\end{bmatrix}=\begin{bmatrix}y_0\\y_1\\\vdots\\y_{n-1}\\\end{bmatrix} \]

对吧。

所以说，对于 FFT过的向量我们只要把他求一个逆，然后得出的自然就是系数表达法

对于这个矩阵的逆，我们有 \(V_nV_n^{-1} = I_n\) ，如何求出逆矩阵，下面给了你方法

首先我们的矩阵应该是:

\[\begin{bmatrix}y_0\\y_1\\\vdots\\y_{n-1}\\\end{bmatrix}=\begin{bmatrix}1 & 1&1&\cdots&1\\1 & \omega_n^1&\omega_n^{2}&\cdots&\omega_n^{n-1}\\\vdots&\vdots&\vdots&\ddots&\vdots\\1 & \omega_n^{n-1}&\omega_n^{2(n-1)}&\cdots&\omega_n^{(n-1)(n-1)}\\\end{bmatrix} \begin{bmatrix}a_0\\a_1\\\vdots\\a_{n-1}\\\end{bmatrix} \]

证明：对于 \(j, k = 0, 1, 2, \cdots. n - 1\) ，\(V^{-1}_n\) 处的\((j, k)\) 为 \(\omega _n^{-kj}/n\)

我们不难相处这个结论是对的，因为这需要你的手工验证（滑稽）

算了给个证明吧。\([V_n^{-1}V_n]_(j, j') = \sum _0^{n-1} (\omega^{-kj}_n/n)(\omega^{kj'}_n)=\sum _0^{n-1}\omega _n^{k(j'-j)}/n\) 此时由于上边的求和定理除非 \(j = j'\) 这是为 1，否则为0

所以我们可以推导出来 \(DFT^{-1}_n (y)\)：

\[a_j = \sum _0^{n-1}y_k\omega_n^{-kj} \]

我们可以模仿FFT的过程，把 \(a\) 与 \(y\) 互换，用 \(\omega _n^{-kj}\) 代替 \(\omega_n^{kj}\) 这是只要在代码稍加修改即可。

我们就有了通过递归完成的形式。

更高效的实现

接下来是 FFT 的一个常见优化，即不断递归是会降低程序效率的，可以在上面的基础上优化常数。这就是一种叫做蝴蝶变换操作的优化：

我们不断递归的过程中，其实奇偶分类显得不是很优美，这让我们使用了低效的递归实现。事实上，有一种更为高效的迭代实现

图不是我的，我们观察这个过程，在观察他的最后的二进制数的关系，不难发现他是反的，所以有一种实现方式，是把他刚开始时就反着放，即存在一个 \(rev(i)\) 他的二进制与 \(i\) 相反，然后我们只需要枚举区间的半长度，在模拟每一次FFT的过程即可。

关于三次变两次优化

吐槽一下神鱼的题解，~~你这个太敷衍了~~

但是好像就是这样.....

我们关于点值表达式，不妨取他的虚部为答案，我们知道这是对问题没有影响的，因为点值表达式并没有什么定义域的限制。

最后献上完整版代码

#include <bits/stdc++.h>

using namespace std;

template <typename T>
inline T read()
{
	T x = 0;
	char ch = getchar();
	bool f = 0;
	while(ch < '0' || ch > '9')
	{
		f = (ch == '-');
		ch = getchar();
	}
	while(ch <= '9' && ch >= '0')
	{
		x = (x << 1) + (x << 3) + (ch - '0');
		ch = getchar();
	}
	return  f? -x : x;
}

template <typename T>
void put(T x)
{
	if(x < 0)
	{
		x = -x;
		putchar('-');
	}
	if(x < 10) {
		putchar(x + 48);
		return;
	}
	put(x / 10);
	putchar(x % 10 + 48);
	return ;
}

#define reg register
#define rd read <int>
#define pt(i) put <int> (i), putchar(' ')

typedef long long ll;
typedef double db;
typedef long double ldb;
typedef unsigned long long ull;
typedef unsigned int ui;

const int Maxn = 1 << 21;

typedef double db;

int n, m, limit, r[Maxn], l;

namespace myspace
{
	const db Pi = 3.14159265358979323846;
	struct complex
	{
		db x, y;
		complex(db xx = 0, db yy = 0) { x = xx, y = yy; }
		void operator = (const complex &a) { x = a.x , y = a.y; }
		complex operator * (const complex &a) { return complex(x * a.x - y * a.y, x * a.y + y * a.x); }
		complex operator + (const complex &a) { return complex(x + a.x, y + a.y); }
		complex operator - (const complex &a) { return complex(x - a.x, y - a.y); }
		inline db abs() { return sqrt(x * x + y * y); }
	};
	void SlowfFt(const 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];
		}
		SlowfFt(limit >> 1, a1, type);
		SlowfFt(limit >> 1, a2, type);
		complex w(1, 0), wn(cos(2.0 * Pi / limit), sin(2.0 * Pi / limit) * type);
		for(int i = 0; i < limit >> 1; ++i, w = w * wn)
		{
			complex io = w * a2[i];
			a[i] = a1[i] + io, a[(limit >> 1) + i] = a1[i] - io;
		}
		return ;
	}
	void fFt(const int limit, complex *a, int type)
	{
		for(int i = 0; i < limit; ++i) if(i < r[i]) swap(a[i], a[r[i]]);
		for(int mid = 1; mid < limit; mid <<= 1)
		{
      // 这里为啥不乘二，因为 mid 本身就 / 2 了
			complex wn(cos(Pi / mid), sin(Pi / mid) * type);
			for(int R = mid << 1, j = 0; j < limit; j += R)
			{
				complex w(1, 0);
				for(int k = 0; k < mid; ++k, w = w * wn)
				{
					complex io = a[j + k], oi = w * a[j + k + mid];
					a[j + k] = io + oi;
					a[j + k + mid] = io - oi;
				}
			}
		}
	}
}
//注释代码均为不加三次变两次优化
myspace::complex a[Maxn];
// myspace::complex b[Maxn];

int main()
{
#ifdef _DEBUG
	freopen("in.txt", "r", stdin);
#endif
	n = rd();
	m = rd();
	for(reg int i = 0; i <= n; ++i) a[i].x = rd();
	for(reg int j = 0; j <= m; ++j) a[j].y = rd();
	// for(reg int j = 0; j <= m; ++j) b[j].x = rd();
  limit = 1;
	while(limit <= n + m) limit <<= 1, ++l;
	for(reg int i = 0; i < limit; ++i) r[i] = ((r[i >> 1] >> 1) | ((i & 1) << (l - 1)));
	myspace::fFt(limit, a, 1);
	// myspace::fFt(limit, b, 1);
  for(reg int i = 0; i < limit; ++i) a[i] = a[i] * a[i];
	//for(reg int i = 0; i < limit; ++i) a[i] = a[i] * b[i];
  myspace::fFt(limit, a, -1);
	for(reg int i = 0; i <= n + m; ++i) pt((a[i].y / (limit << 1) + 0.5));
	// for(reg int i = 0; i <= n + m; ++i) pt((a[i].x / limit + 0.5));
  return 0;
}

核心思想

运用了点值形式的快速乘法和本源单位根的性质。即有下图所示

FFT 瞎讲