多项式

拉格朗日插值

• 已知最高次数不超过 $n-1$ 的多项式在平面上的 $n$ 个点 $(x_1,y_1),(x_2,y_2),\dots ,(x_n,y_n)$，满足 $\forall 1\le i<j\le n,x_i\ne x_j$，则可还原出多项式：
  $$f(x)=\sum _{i=1}^n{y}_i\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}$$
• 若 $x_i=i$，且 $x>n$，则该式可简化为：
  $$f(x)=\prod _{i=1}^n(x-i)\left(\sum _{i=1}^n\frac{(-1{)}^{n-i}{y}_i}{(x-i)(i-1)!(n-i)!}\right)$$
• 预处理相关逆元和前缀积，即可 $\mathcal{O}(n)$ 计算。
• 若已知二元多项式在空间内的 $nm$ 个点 $(x_{11},y_{11},z_{11}),\dots ,(x_{nm},y_{nm},z_{nm})$，满足任意两点横纵坐标至少一个不相等，且该多项式中 $x$ 的最高次数不超过 $n-1$，$y$ 的最高次数不超过 $m-1$，则可还原出多项式：
  $$f(x,y)=\sum _{i=1}^n\sum _{j=1}^m{z}_{ij}\prod _{k\ne i}\frac{x-x_{kj}}{x_{ij}-x_{kj}}\prod _{l\ne j}\frac{x-x_{il}}{x_{ij}-x_{il}}$$

FFT

• 即快速傅里叶变换。
• 设 $n$ 次单位根 ${\omega }_n^k=\cos \frac{2\pi k}{n}+i\sin \frac{2\pi k}{n}$。
• 将多项式 $A,B$ 的项数补到 2 的整数次幂，设项数为 $n$，采用分治法快速求出 $A,B$ 代入 ${\omega }_n^0,{\omega }_n^1,\dots ,{\omega }_n^{n-1}$ 的点值，将点值相乘再通过类似的过程还原回多项式，即可快速求出多项式 $A\times B$，时间复杂度 $\mathcal{O}(n\log n)$。

DFT

• 即离散傅里叶变换。
• 对于多项式
  $$A(x)=\sum _{k=0}^{n-1}{a}_k{x}^k=\sum _{k=0}^{\frac{n}{2}-1}(a_{2k}{x}^{2k}+a_{2k+1}{x}^{2k+1})=\sum _{k=0}^{\frac{n}{2}-1}{a}_{2k}(x^2{)}^k+x\sum _{k=0}^{\frac{n}{2}-1}{a}_{2k+1}(x^2{)}^k$$
• 设 $A_1(x)=\sum _{k=0}^{\frac{n}{2}-1}{a}_{2k}{x}^k,A_2(x)=\sum _{i=0}^{\frac{n}{2}-1}{a}_{2k+1}{x}^k$，若已求得 $A_1({\omega }_{\frac{n}{2}}^k),A_2({\omega }_{\frac{n}{2}}^k),0\le k\le \frac{n}{2}-1$，则
  $$\begin{gathered} A({\omega }_n^k)=A_1({\omega }_n^{2k})+w_n^k{A}_2(w_n^{2k})=A_1({\omega }_{\frac{n}{2}}^k)+w_n^k{A}_2({\omega }_{\frac{n}{2}}^k) \\ A({\omega }_n^{k+\frac{n}{2}})=A_1({\omega }_n^{2k})-w_n^{k+\frac{n}{2}}{A}_2(w_n^{2k})=A_1({\omega }_{\frac{n}{2}}^k)-w_n^k{A}_2({\omega }_{\frac{n}{2}}^k) \end{gathered}$$
• 不断分治下去，此时我们便求得了 $A({\omega }_n^k),0\le k\le n-1$。

IDFT

• 即逆离散傅里叶变换，由之前的叙述我们有：
  [ ( ω n 0 ) 0 ( ω n 0 ) 1 ⋯ ( ω n 0 ) n − 1 ( ω n 1 ) 0 ( ω n 1 ) 1 ⋯ ( ω n 1 ) n − 1 ⋮ ⋮ ⋱ ⋮ ( ω n n − 1 ) 0 ( ω n n − 1 ) 1 ⋯ ( ω n n − 1 ) n − 1 ] [ a 0 a 1 ⋮ a n − 1 ] = [ A ( ω n 0 ) A ( ω n 1 ) ⋮ A ( ω n n − 1 ) ] \left[
  (ωn0)0(ωn0)1⋯(ωn0)n−1(ωn1)0(ωn1)1⋯(ωn1)n−1⋮⋮⋱⋮(ωnn−1)0(ωnn−1)1⋯(ωnn−1)n−1" role="presentation" style="position: relative;">(ω0n)0(ω1n)0⋮(ωn−1n)0(ω0n)1(ω1n)1⋮(ωn−1n)1⋯⋯⋱⋯(ω0n)n−1(ω1n)n−1⋮(ωn−1n)n−1$$\begin{array}{cccc}({\omega }_n^0{)}^0& ({\omega }_n^0{)}^1& \cdots & ({\omega }_n^0{)}^{n-1}\\ ({\omega }_n^1{)}^0& ({\omega }_n^1{)}^1& \cdots & ({\omega }_n^1{)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }_n^{n-1}{)}^0& ({\omega }_n^{n-1}{)}^1& \cdots & ({\omega }_n^{n-1}{)}^{n-1}\end{array}$$
  \right] \left[
  a0a1⋮an−1" role="presentation" style="position: relative;">a0a1⋮an−1$$\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{n-1}\end{array}$$
  \right] = \left[
  A(ωn0)A(ωn1)⋮A(ωnn−1)" role="presentation" style="position: relative;">A(ω0n)A(ω1n)⋮A(ωn−1n)$$\begin{array}{c}A({\omega }_n^0)\\ A({\omega }_n^1)\\ ⋮\\ A({\omega }_n^{n-1})\end{array}$$
  \right] ⎣ ⎡​(ωn0​)0(ωn1​)0⋮(ωnn−1​)0​(ωn0​)1(ωn1​)1⋮(ωnn−1​)1​⋯⋯⋱⋯​(ωn0​)n−1(ωn1​)n−1⋮(ωnn−1​)n−1​⎦ ⎤​⎣ ⎡​a0​a1​⋮an−1​​⎦ ⎤​=⎣ ⎡​A(ωn0​)A(ωn1​)⋮A(ωnn−1​)​⎦ ⎤​
• 对矩阵求逆，可得：
  [ a 0 a 1 ⋮ a n − 1 ] = 1 n [ ( ω n − 0 ) 0 ( ω n − 0 ) 1 ⋯ ( ω n − 0 ) n − 1 ( ω n − 1 ) ) 0 ( ω n − 1 ) 1 ⋯ ( ω n − 1 ) n − 1 ⋮ ⋮ ⋱ ⋮ ( ω n − ( n − 1 ) ) 0 ( ω n − ( n − 1 ) ) 1 ⋯ ( ω n − ( n − 1 ) ) n − 1 ] [ A ( ω n 0 ) A ( ω n 1 ) ⋮ A ( ω n n − 1 ) ] \left[
  a0a1⋮an−1" role="presentation" style="position: relative;">a0a1⋮an−1$$\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{n-1}\end{array}$$
  \right] = \frac{1}{n} \left[
  (ωn−0)0(ωn−0)1⋯(ωn−0)n−1(ωn−1))0(ωn−1)1⋯(ωn−1)n−1⋮⋮⋱⋮(ωn−(n−1))0(ωn−(n−1))1⋯(ωn−(n−1))n−1" role="presentation" style="position: relative;">(ω−0n)0(ω−1n))0⋮(ω−(n−1)n)0(ω−0n)1(ω−1n)1⋮(ω−(n−1)n)1⋯⋯⋱⋯(ω−0n)n−1(ω−1n)n−1⋮(ω−(n−1)n)n−1$$\begin{array}{cccc}({\omega }_n^{-0}{)}^0& ({\omega }_n^{-0}{)}^1& \cdots & ({\omega }_n^{-0}{)}^{n-1}\\ ({\omega }_n^{-1}){)}^0& ({\omega }_n^{-1}{)}^1& \cdots & ({\omega }_n^{-1}{)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }_n^{-(n-1)}{)}^0& ({\omega }_n^{-(n-1)}{)}^1& \cdots & ({\omega }_n^{-(n-1)}{)}^{n-1}\end{array}$$
  \right] \left[
  A(ωn0)A(ωn1)⋮A(ωnn−1)" role="presentation" style="position: relative;">A(ω0n)A(ω1n)⋮A(ωn−1n)$$\begin{array}{c}A({\omega }_n^0)\\ A({\omega }_n^1)\\ ⋮\\ A({\omega }_n^{n-1})\end{array}$$
  \right] ⎣ ⎡​a0​a1​⋮an−1​​⎦ ⎤​=n1​⎣ ⎡​(ωn−0​)0(ωn−1​))0⋮(ωn−(n−1)​)0​(ωn−0​)1(ωn−1​)1⋮(ωn−(n−1)​)1​⋯⋯⋱⋯​(ωn−0​)n−1(ωn−1​)n−1⋮(ωn−(n−1)​)n−1​⎦ ⎤​⎣ ⎡​A(ωn0​)A(ωn1​)⋮A(ωnn−1​)​⎦ ⎤​
• 不难发现，将 ${\omega }_n^k$ 替换成 ${\omega }_n^{-k}$ 再做一遍 DFT，最后将结果除以 $n$ 即可实现 IDFT。

位逆序变换

• 为了减小递归实现带来的常数，我们考虑直接得到递归最底层的排列顺序，再逐层向上合并。
• 不难证明，设 $n=2^m$，则系数 $a_i$ 递归到最底层时的下标恰好为 $i$ 在 $m$ 位二进制表示下的对称翻转，我们称这个变换为位逆序变换（蝴蝶变换）。
• 设 $rev[i]$ 表示系数 $a_i$ 递归到最底层时的下标，则显然有递推式

	int m = 0, _n = na + nb; //相乘的两个多项式最高次数分别为 na,nb
	for (n = 1; n <= _n; n <<= 1)
		++m;
	for (int i = 1; i < n; ++i)
		rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (m - 1));
1
2
3
4
5

• 最终我们得到了 FFT 的迭代实现。

typedef long double ld;
typedef complex<ld> com;
const ld pi = acos(-1.0);

inline void DFT(vector<com> &a, int opt)
{
	int n = a.size();
	for (int i = 0; i < n; ++i)
		if (i < rev[i])
			std::swap(a[i], a[rev[i]]);
	for (int k = 1; k < n; k <<= 1)
	{
		com w(cos(pi / k), opt * sin(pi / k));;
		for (int i = 0; i < n; i += k << 1)
		{
			com res(1.0, 0.0);
			for (int j = 0; j < k; ++j)
			{
				com u = a[i + j],
					v = res * a[i + j + k];
				a[i + j] = u + v;
				a[i + j + k] = u - v;
				res = res * w;
			}
		}
	}
	if (opt == -1)
	{
		for (int i = 0; i < n; ++i)
			a[i] /= n;
	}
}
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

NTT

• 即快速数论变换。
• 在模质数 $P$ 意义下，原根 $g$ 具有单位根的性质：
  $$(g^k{)}^{P-1}\equiv 1(\mathrm{mod}P),g^{\frac{P-1}{2}}\equiv -1(\mathrm{mod}P)$$
• 若 $P=2^{x}a+1$，设 $n=2^m$，用 $g^{2^{x-m}a}$ 替换 ${\omega }_n$ 即可实现 NTT。
• 常见的 $P$ 有：
  $$\begin{gathered} P=1004535809=479\times 2^{21}+1,g=3 \\ P=998244353=7\times 17\times 2^{23}+1,g=3 \end{gathered}$$
• 实现时可以预处理 $g^{2^{x-m}a}$ 的幂次，减小常数。

const int mod = 998244353;
const int inv3 = 332748118;

inline void DFT(vector<int> &a, int opt)
{
	int n = a.size(), g = opt == 1 ? 3 : inv3;
	for (int i = 0; i < n; ++i)
		if (i < rev[i])
			std::swap(a[i], a[rev[i]]);
	for (int k = 1; k < n; k <<= 1)
	{
		int w = quick_pow(g, (mod - 1) / (k << 1));
		tw[0] = 1;
		for (int j = 1; j < k; ++j)
			tw[j] = 1ll * tw[j - 1] * w % mod;
		for (int i = 0; i < n; i += k << 1)
		{
			for (int j = 0; j < k; ++j)
			{
				int u = a[i + j],
					v = 1ll * tw[j] * a[i + j + k] % mod;
				add(a[i + j] = u, v);
				dec(a[i + j + k] = u, v);
			}
		}
	}
	if (opt == -1)
	{
		int inv_n = quick_pow(n, mod - 2); 
		for (int i = 0; i < n; ++i)
			a[i] = 1ll * a[i] * inv_n % mod;
	}
}
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