多项式求除法及求模

问题是这样的：给定一个 $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）

代码实现

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);
}