多项式乘法学习笔记

多项式乘法

给定两个多项式

$A(x) = \sum_{i=0}^{n-1}{a_i x^i}$

$B(x)=\sum^{m-1}_{i=0}{b_i x^i}$

求 $C(x) = A(x) \times B(x) = \sum_{i=0}^{n+m-1}{c_i x^i}$

前置知识

单位根

$n$ 次单位根即为满足 $x^n=1$ 的 x。

由代数基本定理可得 $n$ 次单位根应该有 $n$ 个。

第 $k$ 个记为 $\omega^k_n$。

单位根的一些性质

$$\omega_n^k = \cos {\dfrac{2k\pi}{n}} + \mathcal{i} \sin{\dfrac{2k\pi}{n}} (k = 0, 1, \dots, n - 1) \\ \omega^{2k}_{2n} = \omega^{k}_{n} \\ \omega^{k + \frac{n}{2}}_{n}=-\omega^{k}_{n} \\ \omega^i_n \omega^j_n = \omega^{i+j}_n \\ (\omega^k_n)^j = \omega^{jk}_n \\ $$

原根的性质（$p$ 为质数）

$$\forall i, j \in \Z, g^i \not \equiv g^j \left( mod p \right) \\ g^{p-1} \equiv 1 (modp) \\ $$

多项式的点表示法

众所周知，一个 $n - 1$ 次多项式可以用 $n$ 个互异的点表示出来。

这就叫做多项式的点表示法。

DFT

上面的点表示法启示我们了一种求多项式乘法的方法。

一个 $n$ 次的多项式 $A(x)$，对于取值 $\{x_1, x_2, x_3, \dots, x_n\}$，可以求出$\{y_{A_1}, y_{A_2}, y_{A_3}, \dots, y_{A_n}\}$。

同理一个 $n$ 次的多项式 $B(x)$，也可以求出相对应的 $\{y_{B_1}, y_{B_2}, y_{B_3}, \dots, y_{B_n}\}$。

然后我们对于每一个 $x_i$ 可以求出一个对应的 $y_{C_i}=y_{A_i}\times y_{B_i}$。

得到了 $\{y_{C_1}, y_{C_2}, y_{C_3}, \dots, y_{C_n}\}$，我们就可以用拉格朗日插值法求出多项式$C(x)=A(x) \times B(x)$ 了。

当我们将取值 $\{x_1, x_2, x_3, \dots, x_n\}$ 取为 $\{ \omega^0_n, \omega^1_n, \omega^2_n, \dots, \omega^{n-1}_n \}$ 时我们就可以少算很多的次幂。

利用上述方法求多项式乘法的方法就叫做 DFT。

FFT

上面的 DFT 虽然看起来很妙，但仍然是 $\mathcal{O}(n^2)$ 的，其主要复杂度为求出对应的 $\{y_{A_1}, y_{A_2}, y_{A_3}, \dots, y_{A_n}\}$ 和 $\{y_{B_1}, y_{B_2}, y_{B_3}, \dots, y_{B_n}\}$。

我们考虑优化这个过程。

接下来的过程中我们假设 $n$ 为 $2$ 的整数次幂。

$$A(x) = a_0 + a_1x + a_2x^2 + a_3x^3 + \dots + a_{n-1}x^{n-1}$$

按照 $x$ 次数的奇偶分类。

$$A(x) = (a_0 + a_2x^2 + \dots + a_{n-2}x^{n-2}) + x(a_1 + a_3x^2 + \dots + a_{n-1}x^{n-2})$$

记

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

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

则 $A(x) = A_1(x^2) + xA_2(x^2)$。

显然 $A_1(x)$ 和 $A_2(x)$ 都为 $\dfrac{n}{2}$ 次多项式。

接下来就是 FFT 的核心。

对于所有 $k \leq \frac{n}{2}$，有

$$A(\omega^k_n) = A_1(\omega^{2k}_n ) + \omega^k_n A_2(\omega^{2k}_n) \\ A(\omega^k_n) = A_1(\omega^{k}_{\frac{n}{2}} ) + \omega^{k}_{\frac{n}{2}}A_2(\omega^{k}_{\frac{n}{2}})$$

$$A(\omega^{k + \frac{n}{2}}_n) = A_1(\omega^{2k + n}_n) + \omega^{k + \frac{n}{2}}_n A_2(\omega^{2k + n}_n) \\ A(\omega^{k + \frac{n}{2}}_n) = A_1(\omega^{2k}_n ) - \omega^k_n A_2(\omega^{2k}_n) \\ A(\omega^{k + \frac{n}{2}}_n) = A_1(\omega^{k}_{\frac{n}{2}} ) - \omega^{k}_{\frac{n}{2}}A_2(\omega^{k}_{\frac{n}{2}})$$

然后我们就发现了一个很神奇的东西，如果我们知道了 $A_1(\omega^{k}_{\frac{n}{2}})$ 和 $A_2(\omega^{k}_{\frac{n}{2}})$ 我们就可以同时求出 $A(\omega^k_n)$ 和 $A(\omega^{k + \frac{n}{2}}_n)$。

并且由于 $A_1(x)$ 和 $A_2(x)$ 都为 $\dfrac{n}{2}$ 次多项式，我们可以用相同的办法求出他们，并且每次将次数缩小到 $\dfrac{1}{2}$。

当递归到 $1$ 次时直接带入求值即可。

这就是一个类似于分治的复杂度了，不难证明这个过程是 $\mathcal{O}(n \log n)$ 的。

不过我们仍然还有问题需要解决，上面我们只优化了求出对应的点坐标的复杂度，并没有优化求出系数的复杂度，如果我们用拉格朗日插值来求出系数这个算法的复杂度还是 $\mathcal{O}(n^2)$ 的。

但这就是 DFT 的优越之处，它的特殊取值使得它可以同样在 $\mathcal{O}(n \log n)$ 的复杂度内求出对应的系数。

IFFT

IFFT 即将点值变为系数的过程。

有个比较神奇的结论，记我们通过 FFT 算出来的结果的乘积为 $\{y_{C_0}, y_{C_1}, y_{C_2}, \dots, y_{C_{n-1}}\}$。

记多项式 $D(x) = y_{C_0} + y_{C_1}x + y_{C_2}x^2 + \dots + y_{C_{n-1}}x^{n - 1}$ 在 $\{ \omega^{-0}_n, \omega^{-1}_n, \omega^{-2}_n, \dots, \omega^{-(n-1)}_n \}$ 处的取值为 $\{ d_0, d_1, d_2, \dots, d_{n-1} \}$。

有 $c_i = \dfrac{d_i}{n}$。

即 $C(x) = \dfrac{d_0}{n} + \dfrac{d_1}{n}x + \dfrac{d_2}{n}x^2 + \dots + \dfrac{d_{n-1}}{n}x^{n-1}$。

IFFT 证明

$$d_k = \sum^{n-1}_{i=0}{y_{C_i} (\omega^{-k}_n)^i} \\ d_k =\sum^{n-1}_{i=0}{(\sum^{n-1}_{j=0}{c_j (\omega^i_n)^j}) (\omega^{-k}_n)^i} \\ d_k = \sum^{n-1}_{i=0}{(\sum^{n-1}_{j=0}{c_j (\omega^j_n)^i}) (\omega^{-k}_n)^i} \\ d_k =\sum^{n-1}_{i=0}{\sum^{n-1}_{j=0}{c_j (\omega^{j-k}_n)^i}} \\ d_k =\sum^{n-1}_{j=0}{c_j \sum^{n-1}_{i=0}{(\omega^{j-k}_n)^i}} \\ $$

记 $S(k) = \sum_{i=0}^{n-1}{(\omega^k_n)^i}$。

根据等比数列求和公式，首项为 $1$，公比为 $w^k_n$ ，当公比不为 $0$，即 $k \neq 0$ 时，有：

$$S(k) = \dfrac{1 - (\omega^k_n)^n}{1 - \omega^k_n} = \dfrac{1 - (\omega^n_n)^k}{1 - \omega^k_n} = \dfrac{1 - 1^k}{1 - \omega^k_n} = 0$$

当 $k = 0$ 时，有：

$$S(k) = \sum_{i=0}^{n-1}{1^i} = n$$

所以不难得出

$$d_k = \sum^{n-1}_{j=0}S(j-k) c_j \\ d_k = nc_k \\ c_k = \dfrac{d_k}{n}$$

FFT 未优化代码

#include <bits/stdc++.h>
using namespace std;
 
namespace Math {
	#define PI 3.14159265358979323846
 
	struct Complex { double r, i; Complex(double r = 0, double i = 0): r(r), i(i) { } };
	Complex operator + (const Complex &a, const Complex &b) { return { a.r + b.r, a.i + b.i }; }
	Complex operator - (const Complex &a, const Complex &b) { return { a.r - b.r, a.i - b.i }; }
	Complex operator * (const Complex &a, const Complex &b) { return { a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r }; }
 
	void FFT(int len, vector<Complex> &a, int op) {
		if (len == 1) return;
		vector<Complex> a1(len / 2), a2(len / 2);
		for (int i = 0; i < a.size(); i += 2) a1[i / 2] = a[i], a2[i / 2] = a[i + 1];
		FFT(len / 2, a1, op), FFT(len / 2, a2, op);
		Complex step = { cos(PI * 2 / len), op * sin(PI * 2 / len) }, omega = { 1, 0 };
		for (int i = 0; i < len / 2; i ++, omega = omega * step) a[i] = a1[i] + omega * a2[i], a[i + len / 2] = a1[i] - omega * a2[i]; 
	} 
	
	vector<int> polyMul(vector<int> a, vector<int> b) {
		int n = a.size(), m = b.size();
		int len = 1;
		while (len < n + m - 1) len *= 2;
		vector<Complex> ac(len), bc(len);
		Complex init = { 0, 0 };
		for (int i = 0; i < len; i ++) ac[i] = (i < n ? a[i] : init); 
		for (int i = 0; i < len; i ++) bc[i] = (i < m ? b[i] : init); 
		FFT(len, ac, 1), FFT(len, bc, 1); // FFT
		vector<Complex> cc(len);
		for (int i = 0; i < len; i ++) cc[i] = ac[i] * bc[i];
		FFT(len, cc, -1); // IFFT
		vector<int> c(n + m - 1);
		for (int i = 0; i < n + m - 1; i ++) c[i] = cc[i].r / len + 0.5;
		return c;
	}
 
	#undef PI
}
 
int main() {
	int n, m; scanf("%d%d", &n, &m);
	vector<int> a(n + 1), b(m + 1);
	for (int i = 0; i <= n; i ++) scanf("%d", &a[i]);
	for (int i = 0; i <= m; i ++) scanf("%d", &b[i]);
	vector<int> c = Math::polyMul(a, b);
	for (int i = 0; i < c.size(); i ++) printf("%d ", c[i]);
	return 0;
}

虽然看网上博客说分治版本的 FFT 没法通过洛谷的模板题，但我莫名其妙过了……

但是常数确实是比较大的。

FFT 优化

主要是可以优化掉递归的常数。

不难发现可以直接通过下表二进制反转后分组来规避掉递归的常数。（为什么？因为奇偶分组是用最低位来决定的，所以考虑这个序列最终的顺序就应该倒过来）

代码：

#include <bits/stdc++.h>
using namespace std;
 
namespace Math {
	#define PI 3.14159265358979323846
 
	struct Complex { double r, i; Complex(double r = 0, double i = 0): r(r), i(i) { } };
	Complex operator + (const Complex &a, const Complex &b) { return { a.r + b.r, a.i + b.i }; }
	Complex operator - (const Complex &a, const Complex &b) { return { a.r - b.r, a.i - b.i }; }
	Complex operator * (const Complex &a, const Complex &b) { return { a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r }; }
 
	vector<int> rev;
 
	void init_FFT(int len, int cnt) {
		rev.resize(len);
		for (int i = 0; i < len; i ++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (cnt - 1));
	} 
 
	// fast implement
	void FFT(int len, vector<Complex> &a, int op) {
		for (int i = 0; i < len; i ++) if (i < rev[i]) swap(a[i], a[rev[i]]);
		for (int t = 1; t < len; t <<= 1) {
			Complex step = { cos(PI / t), op * sin(PI / t) };
			for (int i = 0; i < len; i += 2 * t) {
				Complex omega = { 1, 0 };
				for (int j = i; j < i + t; j ++, omega = omega * step) {
					Complex x = a[j], y = omega * a[j + t];
					a[j] = x + y, a[j + t] = x - y;
				}
			}
		}
	}
	
	vector<int> polyMul(vector<int> a, vector<int> b) {
		int n = a.size(), m = b.size();
		int len = 1, cnt = 0;
		while (len < n + m - 1) len *= 2, ++ cnt;
		init_FFT(len, cnt);
		vector<Complex> ac(len), bc(len);
		Complex init = { 0, 0 };
		for (int i = 0; i < len; i ++) ac[i] = (i < n ? a[i] : init); 
		for (int i = 0; i < len; i ++) bc[i] = (i < m ? b[i] : init); 
		FFT(len, ac, 1), FFT(len, bc, 1); // FFT
		vector<Complex> cc(len);
		for (int i = 0; i < len; i ++) cc[i] = ac[i] * bc[i];
		FFT(len, cc, -1); // IFFT
		vector<int> c(n + m - 1);
		for (int i = 0; i < n + m - 1; i ++) c[i] = cc[i].r / len + 0.5;
		return c;
	}
 
	#undef PI
}
 
int main() {
	int n, m; scanf("%d%d", &n, &m);
	vector<int> a(n + 1), b(m + 1);
	for (int i = 0; i <= n; i ++) scanf("%d", &a[i]);
	for (int i = 0; i <= m; i ++) scanf("%d", &b[i]);
	vector<int> c = Math::polyMul(a, b);
	for (int i = 0; i < c.size(); i ++) printf("%d ", c[i]);
	return 0;
}

优化掉了一秒，可喜可贺。

NTT

FFT 中由于要使用单位根，精度是一个很大的问题，所以就有了 NTT。

在下文中，$n$ 为 $2$ 的整数次幂，$p$ 为质数且 $n \mid p - 1$，$g$ 为 $p$ 的一个原根。

我们设 $h^k_n = (g^{\frac{p-1}{n}})^k$。

那么 $\{h^0_n, h^1_n, \dots, h^{n-1}_n\}$ 就有一些和单位根一样优秀的性质。

$$h^0_n = 1 \\ h^n_n = g^{p-1} \equiv 1 (定义) \\ h^{\frac{n}{2}}_n = g^{\frac{p-1}{2}} \equiv 1 (因为互异性所以不是 1 ) \\ h^{ak}_{an} = g^{\frac{p-1}{an}ak} = (g^{\frac{p-1}{n}})^k = h^k_n$$

可以发现在 FFT 中使用到的单位根的性质原根全部满足。

所以我们就可以使用原根来代替单位根了。

质数我们一般选择 $998244353$ 因为 $998244352$ 有因子 $2^{22}$。

原根我们一般使用 $3$，好记。

NTT 代码

#include <bits/stdc++.h>
using namespace std;
 
namespace Math {
	typedef long long ll;
	const ll MOD = 998244353;
	const ll G = 3;
	const ll GINV = 332748118;
 
 
	vector<int> rev;
	void init_NTT(int len, int cnt) {
		rev.resize(len);
		for (int i = 0; i < len; i ++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (cnt - 1));
	} 
 
	ll qpow(ll a, ll b) {
		ll ret = 1;
		for (; b; b >>= 1) {
			if (b & 1) ret = ret * a % MOD;
			a = a * a % MOD;
		}
		return ret;
	}
 
	void NTT(int len, vector<int> &a, int op) {
		for (int i = 0; i < len; i ++) if (i < rev[i]) swap(a[i], a[rev[i]]);
		for (int t = 1; t < len; t <<= 1) {
			ll h = qpow(op ? G : GINV,  (MOD - 1) / (t * 2));
			for (int i = 0; i < len; i += 2 * t) {
				ll w = 1;
				for (int j = i; j < i + t; j ++, w = w * h % MOD) {
					ll x = a[j], y = w * a[j + t] % MOD;
					a[j] = (x + y) % MOD, a[j + t] = (x - y + MOD) % MOD;
				}
			}
		}
	}
	
	vector<int> polyMul(vector<int> a, vector<int> b) {
		int n = a.size(), m = b.size();
		int len = 1, cnt = 0;
		while (len < n + m - 1) len *= 2, ++ cnt;
		init_NTT(len, cnt);
		a.resize(len), b.resize(len);
		NTT(len, a, 1), NTT(len, b, 1);
		vector<int> c(len);
		for (int i = 0; i < len; i ++) c[i] = 1ll * a[i] * b[i] % MOD;
		NTT(len, c, 0);
		ll inv = qpow(len, MOD - 2);
		c.resize(n + m - 1);
		for (int i = 0; i < n + m - 1; i ++) c[i] = inv * c[i] % MOD;
		return c;
	}
}
 
int main() {
	int n, m; scanf("%d%d", &n, &m);
	vector<int> a(n + 1), b(m + 1);
	for (int i = 0; i <= n; i ++) scanf("%d", &a[i]);
	for (int i = 0; i <= m; i ++) scanf("%d", &b[i]);
	vector<int> c = Math::polyMul(a, b);
	for (int i = 0; i < c.size(); i ++) printf("%d ", c[i]);
	return 0;
}

确实快了不止一点。

完结撒花！

CREDIT

快速傅里叶变换(FFT)详解 - 自为风月马前卒 - 博客园 (cnblogs.com)

十分简明易懂的FFT（快速傅里叶变换）_路人黑的纸巾的博客-CSDN博客_fft