FFT

多项式求除法及求模

Posted by miskcoo on May 22, 2015

问题是这样的:给定一个 $n$ 次多项式 $A(x)$ 和一个 $m(m \leq n)$ 次多项式 $B(x)$,要求求出两个多项式 $D(x), R(x)$,满足

\[\begin{equation} \label{div0} A(x) = D(x)B(x) + R(x) \end{equation}\]

并且 $deg D \leq degA - degB = n - m$,$degR < m$

在解决这个问题之前你需要掌握多项式求逆,利用快速傅里叶变换可以在 $\mathcal O(n\log n)$ 的复杂度内求出这个问题的解

多项式除法在很多问题中都有应用,例如多项式的扩展欧几里得算法、线性递推的优化等等

多项式除法

求解过程

首先看见 $(\ref{div0})$ 会觉得麻烦的是 $R(x)$,现在要想办法消除 $R(x)$ 的影响,如果能做到这一点式子也会变得顺眼很多嗯!

现在给出了一个 $n$ 次多项式 $A(x)$,我们来看看下面这个运算是在干什么,记

\[A^R(x) = x^n A(\frac{1}{x})\]

来拿一个多项式来试试看,比如说

\[A(x) = x^3 + 2x^2 + 4x + 1\]

那么将刚刚的操作应用上来就会变成

\[A^R(x) = x^3 A(\frac{1}{x}) = 1 + 2x + 4x^2 + x^3\]

这相当于将 $A(x)$ 的系数反转

好,现在来干点有趣的事情,把 $(\ref{div0})$ 的 $x$ 全部用 $\frac{1}{x}$ 来替代,然后等式两边同时乘上 $x^n$,现在等式变成了(下面将 $R(x)$ 看成 $m - 1$ 次多项式,$D(x)$ 看成是 $n - m$ 次多项式,高次项不存在认为是 $0$)

\[\begin{eqnarray*} x^n A(\frac{1}{x}) &=& x^{n - m}D(\frac{1}{x}) x^mB(\frac{1}{x}) + x^{n - m + 1}x^{m - 1}R(\frac{1}{x}) \\ A^R(x) &=& D^R(x)B^R(x) + x^{n - m + 1}R^R(x) \end{eqnarray*}\]

现在情况变得很有趣了,$D(x)$ 是要求的元素之一,反转后次数仍然是不会高于 $n - m$,而 $x^{n - m + 1}R(x)$ 这个部分的最低次项高于 $n - m$,我们把上面的式子放到 $\bmod x^{n - m + 1}$ 意义下,$R(x)$ 的影响被消除了,并且不会影响到 $D(x)$,又由于 $A(x), B(x)$ 是已知的元素所以不会有任何问题,现在式子变成

\[A^R(x) = D^R(x)B^R(x) \pmod {x^{n - m + 1}}\]

这样你只需要求一个 $B^R(x)$ 的逆元(过程看开头的那篇文章),就可以得到 $D_R(x)$,然后反转回来回代,便可以得到 $R(x)$,最后的复杂度就是 $\mathcal O(n\log n)$(然后膜拜 Picks

代码实现

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
34
35
36
void polynomial_division(int n, int m, long long *A, long long *B, long long *D, long long *R)
{
	static long long A0[MaxN], B0[MaxN];

	int p = 1, t = n - m + 1;
	while(p < t << 1) p <<= 1;

	fill(A0, A0 + p, 0);
	reverse_copy(B, B + m, A0);
	polynomial_inverse(t, A0, B0);
	fill(B0 + t, B0 + p, 0);
	transform(p, B0);

	reverse_copy(A, A + n, A0);
	fill(A0 + t, A0 + p, 0);
	transform(p, A0);

	for(int i = 0; i != p; ++i)
		A0[i] = A0[i] * B0[i] % mod_v;
	inverse_transform(p, A0);
	reverse(A0, A0 + t);
	copy(A0, A0 + t, D);

	for(p = 1; p < n; p <<= 1);
	fill(A0 + t, A0 + p, 0);
	transform(p, A0);
	copy(B, B + m, B0);
	fill(B0 + m, B0 + p, 0);
	transform(p, B0);
	for(int i = 0; i != p; ++i)
		A0[i] = A0[i] * B0[i] % mod_v;
	inverse_transform(p, A0);
	for(int i = 0; i != m; ++i)
		R[i] = (A[i] - A0[i]) % mod_v;
	fill(R + m, R + p, 0);
}