FFT 瞎讲

浅谈FFT

多项式的表示

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

1. 系数表达法
2. 点值表达法

系数表达法

比如一个 (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)从复数相乘的角度上看这是显然的
• 求和引理：对于任意整数 (nge1) 且 不能被 (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(nlog 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_oomega ^0_n + a_2 omega ^{2k}_n + ldots + a_{n - 2}omega ^{k(n - 2)}_n + a_1omega ^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})

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

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

插值多形式的唯一性：

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

[egin{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} egin{bmatrix}a_0\a_1\vdots\a_{n-1}\end{bmatrix}=egin{bmatrix}y_0\y_1\vdots\y_{n-1}\end{bmatrix} ]

对吧。

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

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

首先我们的矩阵应该是:

[egin{bmatrix}y_0\y_1\vdots\y_{n-1}\end{bmatrix}=egin{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} egin{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_komega_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;
}

核心思想

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