多项式乘法 FFT

计算 \(f(x) = \sum\limits_{i=0}^{n-1}a_ix^i\) 与 \(g(x) = \sum\limits_{i=0}^{m-1}b_ix^i\) 的乘积 \(h(x)=\sum\limits_{i=0}^{n+m-2}c_ix^i\)

一般可以 \(O(n^2)\) 求出：\(c_i = \sum\limits_{j+k=i}a_jb_k = \sum\limits_{j=0}^{i}a_jb_{i-j}\)

\(\text{FFT}\) 算法可以在 \(O(n\log n)\) 时间内求出，\(O(n) = O(m)\)。

准备工作

多项式

\(f(x) = \sum\limits_{i=0}^{n}a_ix^i = a_0 + a_1x + a_2x^2 + \cdots + a_{n}x^{n}\)

\(a_{n}\neq 0\)

\(f(x)\) 为 \(n\) 次多项式，有 \(\deg f = n\)，\(\deg f\) 为 \(f\) 的最高幂次数。

表示多项式

对于 \(f(x) = \sum\limits_{i=0}^{n-1}a_ix^i\)

系数表示：\((a_0,a_1,a_2,\cdots,a_{n-1})\)

点值表示：\(((x_0,f(x_0)), (x_1,f(x_1)), (x_2,f(x_2)),\cdots, (x_{n},f(x_{n})))\)

\(n - 1\) 次的多项式可以由 \(n\) 个点唯一确定。

单位根

设 \(n\) 为偶数。

取 \(n\) 次单位根 \(w_n = \cos({2\pi \over n}) + i \sin({2\pi \over n})\)

单位根具有以下性质：

• \(w_n^0 = w_n^n = 1\)
• \(w_n^{k + n/2} = -w_n^k\)
• \(w_{2n}^{2k}=w_n^k\)

由欧拉公式即得。

算法

\(\text{FFT}\) 流程可以分为：系数表示 \(\rightarrow\) 点值表示 \(\rightarrow\) 系数表示

系数表示 \(\rightarrow\) 点值表示：快速傅里叶变换

点值表示 \(\rightarrow\) 系数表示：快速傅里叶逆变换

用点值表示可以较快进行乘法：

\(\begin{aligned} &((x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),\cdots,(x_n,f(x_n)))\times\\ &((x_0,g(x_0)),(x_1,g(x_1)),(x_2,g(x_2)),\cdots,(x_n,g(x_n))) =\\ &((x_0,f(x_0)g(x_0)),(x_1,f(x_1)g(x_1)),\cdots,(x_n,f(x_n)g(x_n))) \end{aligned}\)

快速傅里叶变换

设 \(n = 2^k\) 且满足 \(n \geq \deg f + \deg g\) 的最小数，若多项式不满 \(n\) 项，高次幂补 \(0\) 即可。

设 \(f(x) = a_0 + a_1x + a_2x^2 + \cdots + a_{n-1}x^{n-1}\)

将 \(f(x)\) 按 \(x\) 幂次的奇偶性分组：

\(f_0(x) = a_0 + a_2x + a_4x^2 + \cdots + a_{n-2}^{n/2-1}\)

\(f_1(x) = a_1 + a_3x + a_5x^2 + \cdots + a_{n-1}^{n/2-1}\)

则 \(f(x) = f_0(x^2) + xf_1(x^2)\)

取 \(x = w_n^k,0\le k < \frac{n}{2}\)：

\[\begin{aligned} &f(w_n^k) \\ &= f_0(w_n^{2k})+w^k_nf_1(w_n^{2k}) \\ &= f_0(w_{2(n/2)}^{2k})+w^k_nf_1(w_{2(n/2)}^{2k}) \\ &= f_0(w_{n/2}^{k})+w^k_nf_1(w_{n/2}^{k}) \end{aligned}\]

即

\[f(w_n^k) = f_0(w_{n/2}^{k})+w^k_nf_1(w_{n/2}^{k}) \]

取 \(x = w_n^{k + n/2},0\le k < \frac{n}{2}\)：

\[\begin{aligned} &f(w_n^{k+n/2})\\ &=f_0(w_n^{2k+n}) + w_n^{k + n/2}f_1(w_n^{2k+n})\\ &=f_0(w_n^{2k}) - w_n^kf_1(w_n^{2k})\\ &=f_0(w_{n/2}^k) - w_n^kf_1(w_{n/2}^k) \end{aligned}\]

即

\[f(w_n^{k+n/2})=f_0(w_{n/2}^k) - w_n^kf_1(w_{n/2}^k) \]

只要求 \(f_0(w_{n/2}^k)\) 和 \(f_1(w_{n/2}^k)\) 就可求出 \(f(w_n^k)\) 和 \(f(w_n^{k+n/2})\)

要求的 \(k\) 的从 \(0\le k<n\) 变为了 \(0\le k < n/2\)，最多会求 \(O(\log n)\) 层，时间复杂度为 \(O(n \log n)\)。

以上将系数表示

\[(a_0,a_1,a_2,\cdots,a_{n-1}) \]

转化为点值表示

\[((w_n^0,f(w_n^0)),(w_n^1,f(w_n^1)),\cdots,(w_n^{n-1},f(w_n^{n-1}))) \]

对 \(f, g\) 分别做一遍，将其值相乘得到

\[((w_n^0,f(w_n^0)g(w_n^0)),(w_n^1,f(w_n^1)g(w_n^1)),\cdots,(w_n^{n-1},f(w_n^{n-1})g(w_n^{n-1}))) \]

记为

\[((w_n^0,h(w_n^0)),(w_n^1,h(w_n^1)),\cdots,(w_n^{n-1},h(w_n^{n-1}))) \]

快速傅里叶逆变换

设 \(h(x) = \sum\limits_{i=0}^{n-1}c_ix^i\)

设 \(A(x) = \sum\limits_{i=0}^{n-1}h(w_n^i)x^i\)

将 \((w_n^{0},w_n^{-1},w_n^{-2},\dots,w_n^{-(n-1)})\) 带入 \(A\) 得到点值表示：

\[A(w_n^{-k})\\ =\sum\limits_{i=0}^{n-1}h(w_n^i)w_n^{-ik}\\ =\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}c_jw_n^{i(j-k)}\\ =\sum\limits_{j=0}^{n-1}c_j\sum\limits_{i=0}^{n-1}w_n^{i(j-k)}\\ \]

设 \(S(w_n^t)=\sum\limits_{i=0}^{n-1}w_n^{it}\)

若 \(t \equiv 0 \pmod n\)，则 \(w_n^t = 1\)，则 \(S(w_n^t)=n\)

若 \(t \not\equiv 0 \pmod n\)，由等比数列求和公式可得

\[S(w_n^t) = \frac{w_n^{tn}-1}{w_n^t-1} = 0 \]

则当 \(j = k\) 时，\(\sum\limits_{i=0}^{n-1}w_n^{i(j-k)}\) 才有值，即 \(A(w_n^{-k}) = c_kn\)

则 \(A\) 用点值表示为

\[((w_n^0, c_0n), (w_n^{-1}, c_1n), (w_n^{-2}, c_2n),\cdots, (w_n^{-(n-1)}, c_{n-1}n)) \]

由欧拉公式可得：

\[\begin{aligned} w_n &= \cos\left(\frac{2k\pi}{n}\right) + \text{i}\sin\left(\frac{2k\pi}{n}\right)\\ w_n^{-1} &= \cos\left(\frac{2k\pi}{n}\right) - \text{i}\sin\left(\frac{2k\pi}{n}\right) \end{aligned}\]

做一遍与快速傅里叶变换类似的流程，最后将每一项 \(\times \frac{1}{n}\) 即可得到最终答案。

代码实现

雷德算法

将二进制翻转，例如 \(0010 \rightarrow 0100\)

int n; // 位数, 从 0 开始
int sn = 1 << n;
Rep(i, 0, sn) t[i] = (t[i >> 1] >> 1) | ((i & 1) << 1)

代码

Luogu P3803 【模板】多项式乘法（FFT）

bool _st;

#define N 10000007

int n, m, t, cs, pos[N];
const lf pi = 3.14159265358979323846;
struct cx {
	lf a, b;
	cx() { a = b = 0; }
	cx(lf x, lf y) { a = x, b = y; }
	cx operator + (const cx& rhs) const { return {a + rhs.a, b + rhs.b}; }
	cx operator - (const cx& rhs) const { return {a - rhs.a, b - rhs.b}; }
	cx operator * (const cx& rhs) const { return {a * rhs.a - b * rhs.b, a * rhs.b + b * rhs.a}; }
} a[N], b[N];

void FFT(cx* A, int ty) {
	cx wn, w;
	Rep(i, 0, t) if(i < pos[i]) ::std::swap(A[i], A[pos[i]]);
	for(int mid = 1; mid < t; mid <<= 1)
	{
		wn.a = cos(pi / mid), wn.b = ty * sin(pi / mid);
		for(int R = mid << 1, j = 0; j < t; j += R)
		{
			w.a = 1, w.b = 0;
			for(int k = 0; k < mid; k++, w = w * wn)
			{
				cx x = A[j + k], y = w * A[j + mid + k];
				A[j + k] = x + y;
				A[j + mid + k] = x - y;
			}
		}
	}
}

bool _ed;
int main() { ::std::cerr << (&_ed - &_st) / 1024.0 / 1024 << ::std::endl;
	qr(n), qr(m);
	int x;
	rep(i, 0, n) qr(x), a[i].a = x;
	rep(i, 0, m) qr(x), b[i].a = x;
	t = 1;
	while(t <= n + m) t <<= 1, cs++;
	Rep(i, 0, t) pos[i] = (pos[i >> 1] >> 1) | ((i & 1) << (cs - 1));
	FFT(a, 1);
	FFT(b, 1);
	rep(i, 0, t) a[i] = a[i] * b[i];
	FFT(a, -1);
	rep(i, 0, n + m) qw((int)(a[i].a / t + 0.5)), ps;
	return 0;
}