多项式全家桶

学习笔记
7.3k 词

基本概念

对于求和式 ∑anxn\sum a_n x^n，如果是有限项相加，则称为多项式。记错 F(x)=∑i=0naixiF(x) = \sum_{i = 0}^n a_i x^i，其中 aia_i 称为该多项式的 ii 次项系数，记作 [xi]F(x)[x^i]F(x)，有时也用 F[i]F[i] 表示。

定义两个多项式的 ⊕\oplus 运算卷积为：

[xn](F×G)(x)=∑i⊕j=n[xi]F(x)[xj]G(x)[x^n](F \times G)(x) = \sum\limits_{i \oplus j = n} [x^i]F(x)[x^j]G(x)

以下提到卷积一般指加法卷积。

FFT

多项式的点值表达：只需要 n+1n + 1 个点值就可以唯一确定一个最高 nn 次多项式。

两个多项式卷积就对应点值相乘。也就是对于任意 xx，(F×G)(x)=F(x)G(x)(F \times G)(x) = F(x)G(x)。

FFT 的基本思想就是把两个系数表达的多项式转化成点值表达（DFT），点值相乘后再转换为系数表达（IDFT）。

单位根

nn 次单位根（nn 为正整数）是 nn 次幂为 11 的复数。表现在复平面上就是单位圆上幅角 为 2πn\frac{2\pi}{n} 的倍数的点。

显然 nn 次单位根有恰好 nn 个，记作 ωn0,ωn1,ωn2,…,ωnn−1\omega_n^0, \omega_n^1, \omega_n^2, \ldots, \omega_n^{n - 1}。其中 ωni\omega_n^i 的幅角为 in2π\frac{i}{n} 2\pi。

单位根的一些基本性质：

  • ωn0=1\omega_n^0 = 1。
  • ωnk=(ωn1)k\omega_n^k = (\omega_n^1)^k。
  • ωniωnj=ωni+j\omega_n^i \omega_n^j = \omega_n^{i + j}。
  • (ωnk)−1=ωn−k=ωnn−k(\omega_n^k)^{-1} = \omega_n^{-k} = \omega_n^{n - k}。
  • (ωni)j=ωnij(\omega_n^i)^j = \omega_n^{ij}。
  • ω2n2k=ωnk\omega_{2n}^{2k} = \omega_n^k。
  • 若 nn 为偶数，则 ωnk+n/2=−ωnk\omega_n^{k + n / 2} = -\omega_n^k。

DFT

对于 n−1n - 1 次多项式 F(x)F(x)，我们设两个多项式 FL(x),FR(x)FL(x), FR(x)：

FL(x)=F[0]+F[2]x+F[4]x2+…+F[n−2]xn/2−1FR(x)=F[1]+F[3]x+F[5]x2+…+F[n−1]xn/2−1FL(x) = F[0] + F[2]x + F[4]x^2 + \ldots + F[n - 2]x^{n / 2 - 1} \\ FR(x) = F[1] + F[3]x + F[5]x^2 + \ldots + F[n - 1]x^{n / 2 - 1}

这里钦定 nn 是偶数。

则 F(x)=FL(x2)+xFR(x2)F(x) = FL(x^2) + xFR(x^2)。

对于 0≤k<n/20 \le k < n / 2，代入 x=ωnkx = \omega_n^k：

F(ωnk)=FL((ωnk)2)+ωnkFR((ωnk)2)=FL(ωn/2k)+ωnkFR(ωn/2k)\begin{aligned} F(\omega_n^k) &= FL((\omega_n^k)^2) + \omega_n^k FR((\omega_n^k)^2) \\ &= FL(\omega_{n / 2}^k) + \omega_n^k FR(\omega_{n / 2}^k) \\ \end{aligned}

代入 x=ωnk+n/2x = \omega_n^{k + n / 2}：

F(ωnk+n/2)=FL((ωnk+n/2)2)+ωnk+n/2FR((ωnk+n/2)2)=FL(ωn/2k)−ωnkFR(ωn/2k)\begin{aligned} F(\omega_n^{k + n / 2}) &= FL((\omega_n^{k + n / 2})^2) + \omega_n^{k + n / 2} FR((\omega_n^{k + n / 2})^2) \\ &= FL(\omega_{n / 2}^k) - \omega_n^k FR(\omega_{n / 2}^k) \\ \end{aligned}

也就是只需要知道 FL(x)FL(x) 和 FR(x)FR(x) 在 ωn/20∼ωn/2n/2−1\omega_{n / 2}^0 \sim \omega_{n / 2}^{n / 2- 1} 处的点值表示，就可以知道 F(x)F(x) 在 ωn0∼ωnn−1\omega_n^0 \sim \omega_n^{n - 1} 处的点值表示。

那么对于一个 n−1n - 1 次多项式，先扩展到 2k2^k 次多项式，然后每次分治再合并就可以实现 DFT 了。

IDFT

设 nn 个点值为 gig_i，则：

nF[k]=∑i=0n−1(ωn−k)iG[i]n F[k] = \sum\limits_{i = 0}^{n - 1} (\omega_n^{-k})^i G[i]

证明：

∑i=0n−1(ωn−k)iG[i]=∑i=0n−1(ωn−k)i∑j=0n−1F[j](ωni)j=∑i=0n−1∑j=0n−1ωni(j−k)F[j]=∑j=0n−1F[j]∑i=0n−1ωni(j−k)=∑j=0n−1F[j]×[n∣(j−k)]×n=nF[k]\begin{aligned} \sum\limits_{i = 0}^{n - 1} (\omega_n^{-k})^i G[i] &= \sum\limits_{i = 0}^{n - 1} (\omega_n^{-k})^i \sum\limits_{j = 0}^{n - 1} F[j] (\omega_n^i)^j \\ &= \sum\limits_{i = 0}^{n - 1}\sum\limits_{j = 0}^{n - 1} \omega_n^{i(j - k)} F[j] \\ &= \sum\limits_{j = 0}^{n - 1} F[j] \sum\limits_{i = 0}^{n - 1} \omega_n^{i(j - k)} \\ &= \sum\limits_{j = 0}^{n - 1} F[j] \times [n | (j - k)] \times n \\ &= n F[k] \end{aligned}

第三个等号到第四个等号是单位根反演。

那么就相当于 G[i]G[i] 作为系数再代入一遍。无非这次代入的是 ωn0,ωn−1,ωn−2,…,ωn1−n\omega_n^0, \omega_n^{-1}, \omega_n^{-2}, \ldots, \omega_n^{1 - n}。那么就相当于 DFT 里求单位根后再求个逆元即可。

蝴蝶变换

我们发现 DFT 和 IDFT 差别不大，因此可以把这两个写在一起。

然后考虑优化分治的常数。

考虑过程中系数的变化，这里以 n=8n = 8 为例：

F[0]F[1]F[2]F[3]F[4]F[5]F[6]F[7]F[0]F[2]F[4]F[6]F[1]F[3]F[5]F[7]F[0]F[4]F[2]F[6]F[1]F[5]F[3]F[7]\begin{matrix} F[0] & F[1] & F[2] & F[3] & F[4] & F[5] & F[6] & F[7] \\ F[0] & F[2] & F[4] & F[6] & F[1] & F[3] & F[5] & F[7] \\ F[0] & F[4] & F[2] & F[6] & F[1] & F[5] & F[3] & F[7] \\ \end{matrix}

稍微观察一下可以发现最终序列就是原序列下标二进制翻转。

这个东西显然可以 O(n)O(n) 递推。然后从下往上合并就可以了。

代码实现

复数类：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
struct complex {
	double a, b;

	complex() = default;
	complex(double _a, double _b): a(_a), b(_b) {}
	complex operator+(const complex &x) const {return complex(a + x.a, b + x.b);}
	complex operator-(const complex &x) const {return complex(a - x.a, b - x.b);}
	complex operator*(const complex &x) const {return complex(a * x.a - b * x.b, a * x.b + b * x.a);}
	complex operator/(const complex &x) const {
		double t = b * b + x.b * x.b;
		return complex((a * x.a + b * x.b) / t, (b * x.a - a * x.b) / t);
	}
	complex &operator+=(const complex &x) {return *this = *this + x;}
	complex &operator-=(const complex &x) {return *this = *this - x;}
	complex &operator*=(const complex &x) {return *this = *this * x;}
	complex &operator/=(const complex &x) {return *this = *this / x;}
};

FFT：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
void FFT(vector<complex> &f, int flag) const {
	int n = f.size();
	vector<int> swp(n);
	for (int i = 0; i < n; i++) {
		swp[i] = swp[i >> 1] >> 1 | ((i & 1) * (n >> 1));
		if (i < swp[i]) std::swap(f[i], f[swp[i]]);
	}
	for (int mid = 1; mid < n; mid <<= 1) {
		complex w1(cos(pi / mid), flag * sin(pi / mid));
		for (int i = 0; i < n; i += mid << 1) {
			complex w(1, 0);
			for (int j = 0; j < mid; j++, w *= w1) {
				complex x = f[i + j], y = w * f[i + mid + j];
				f[i + j] = x + y, f[i + mid + j] = x - y;
			}
		}
	}
	return ;
}

卷积：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
poly operator*(const poly &b) const {
	const poly &a(*this);
	int n = 1;
	while (n < (int)(a.size() + b.size()) - 1) n <<= 1;
	vector<complex> A(n), B(n);
	for (int i = 0; i < (int)a.size(); i++) A[i] = complex(a[i], 0);
	for (int i = 0; i < (int)b.size(); i++) B[i] = complex(b[i], 0);
	FFT(A, 1), FFT(B, 1);
	vector<complex> C(n);
	for (int i = 0; i < n; i++) C[i] = A[i] * B[i];
	FFT(C, -1);
	poly c;
	for (int i = 0; i < (int)(a.size() + b.size()) - 1; i++) c.push_back(C[i].a / n + 0.5);
	return c;
}

NTT

FFT 依赖于单位根，于是不可避免地会产生精度问题。考虑找一个单位根的平替。

然而数学家已经证明，在复数域下单位根是唯一一类满足要求的数。

好在大部分题目都在模意义下进行。考虑模意义下什么东西能替代单位根。

那么就要思考我们用到了单位根的哪些性质：

  • ωnk=(ωn1)k\omega_n^k = (\omega_n^1)^k。
  • ωn0∼ωnn−1\omega_n^0 \sim \omega_n^{n - 1} 互不相同。
  • ωnk=ωnk mod n\omega_n^k = \omega_n^{k \bmod n}。
  • ω2n2k=ωnk\omega_{2n}^{2k} = \omega_n^k。

我们发现原根能很好地满足要求，更准确地说是令 ωn1=g(p−1)/n\omega_n^1 = g^{(p - 1) / n}。此处 pp 为模数且是质数。

容易发现只要 n∣(p−1)n | (p - 1) 则上述性质都是满足的。

由于 nn 是 22 的幂次，因此我们希望 p−1p - 1 的质因子 22 尽可能多。比如 p=998244353=223×7×17+1p = 998244353 = 2^{23} \times 7 \times 17 + 1。

常见 NTT 模数原根表：https://blog.miskcoo.com/2014/07/fft-prime-table。

代码：

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
void NTT(vector<int> &g, int flag) const {
	int n = g.size();
	vector<unsigned long long> f(g.begin(), g.end());
	vector<int> swp(n);
	for (int i = 0; i < n; i++) {
		swp[i] = swp[i >> 1] >> 1 | ((i & 1) * (n >> 1));
		if (i < swp[i]) std::swap(f[i], f[swp[i]]);
	}
	for (int mid = 1; mid < n; mid <<= 1) {
		int w1 = power(flag ? G : invG, (mod - 1) / mid / 2);
		vector<int> w(mid);
		w[0] = 1;
		for (int i = 1; i < mid; i++) w[i] = (long long)w[i - 1] * w1 % mod;
		for (int i = 0; i < n; i += mid << 1)
			for (int j = 0; j < mid; j++) {
				int t = (long long)w[j] * f[i + mid + j] % mod;
				f[i + mid + j] = f[i + j] - t + mod;
				f[i + j] += t;
			}
		if (mid == 1 << 10)
			for (int i = 0; i < n; i++) f[i] %= mod;
	}
	int inv = flag ? 1 : power(n, mod - 2);
	for (int i = 0; i < n; i++) g[i] = f[i] % mod * inv % mod;
	return ;
}

这里用了一个取模优化：用 unsigned long long 存储结果，只需在中间进行一次取模即可。

留言