多项式除法

问题

给出\(n\)次多项式\(A(x)\)，\(m\)次多项式\(B(x)\)，求多项式\(D(x)\)，\(R(x)\)使得$$A(x)=B(x)D(x)+R(x)$$，满足\(deg\le n-m,deg\ R<m\)。

即求多项式\(A(x)\)对\(B(x)\)的带余除法。

做法

首先，对于任意多项式$$P(x)=\sum _{i=0}^na_ixi$$，显然有$$P^R(x)=xnP(\frac{1}{x})=\sum _{i=0}^na_{n-i}xi$$。

这就是说，通过这个操作，我们可以把多项式的系数翻转。例如：

\[P(x)=x^2+2x+3 \\ x^2P(\frac{1}{x})=3x^2+2x+1 \]

对\(A(x)=B(x)D(x)+R(x)\)进行这个操作：

\[A(x)=B(x)D(x)+R(x) \\ A(\frac{1}{x})=B(\frac{1}{x})D(\frac{1}{x})+R(\frac{1}{x}) \\ x^nA(\frac{1}{x})=x^nB(\frac{1}{x})D(\frac{1}{x})+x^nR(\frac{1}{x}) \\ x^nA(\frac{1}{x})=x^mB(\frac{1}{x})x^{n-m}D(\frac{1}{x})+x^{n-m+1}x^{m-1}R(\frac{1}{x}) \\ A^R(x)=B^R(x)D^R(x)+x^{n-m+1}R^R(x) \]

注意到\(x^{n-m+1}R^R(x)\)的最低次项至少为\(n-m+1\)次，我们对整个式子取模：

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

可以这样做是因为\(A\)和\(B\)是直接给出的，不会出问题，而\(deg\ D\le n-m\)，\(deg\ D^R\le n-m\)，所以也不会出问题。

这样就可以轻松地算出$$D^R=A_R(BR)^{-1}$$，运用多项式求逆元的方法。之后把\(D^R\)转换回\(D\)，带入减一减就能得到\(R\)的值。

整个做法的关键就是消去余数。

代码

#include<cstdio>
#include<cctype>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long giant;
const int q=998244353;
const int basic=3;
const int maxn=1<<17|1;
const int maxj=18;
int read() {
	int x=0,f=1;
	char c=getchar();
	for (;!isdigit(c);c=getchar()) if (c=='-') f=-1;
	for (;isdigit(c);c=getchar()) x=x*10+c-'0';
	return x*f;
}
int A[maxn],B[maxn],C[maxn],D[maxn],R[maxn];
int AR[maxn],BR[maxn],IBR[maxn],DR[maxn],T[maxn];
int rev[maxn],Wn[maxj][2];
void print(int a[],int len) {
	len=max(len,0);
	for (int i=0;i<=len;++i) printf("%d ",a[i]);
	puts("");
}
int mi(int x,int y) {
	int ret=1;
	for (;y;y>>=1,x=(giant)x*x%q) if (y&1) ret=(giant)ret*x%q;
	return ret;
}
int inv(int x) {
	return mi(x,q-2);
}
void ntt(int a[],int len,bool op=0) {
	for (int i=0;i<len;++i) if (i<rev[i]) swap(a[i],a[rev[i]]);
	for (int i=2,the=1;i<=len;i<<=1,++the) {
		int wn=Wn[the][op];
		for (int j=0;j<len;j+=i) {
			int w=1;
			for (int k=j;k<j+i/2;++k,w=(giant)w*wn%q) {
				int u=a[k]%q,v=(giant)a[k+i/2]*w%q;
				a[k]=(u+v)%q,a[k+i/2]=((u-v)%q+q)%q;
			}
		}
	}
	if (op) {
		int iv=inv(len);
		for (int i=0;i<len;++i) a[i]=(giant)a[i]*iv%q;
	}
}
void TIME(int a[],int la,int b[],int lb,int c[]) {
	int len=max(la,lb),M,xj;
	for (M=1,xj=0;M<=(len<<1);M<<=1,++xj);
	for (int i=0;i<M;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(xj-1));
	ntt(a,M);
	ntt(b,M);
	for (int i=0;i<M;++i) c[i]=(giant)a[i]*b[i]%q;
	ntt(c,M,1);
	ntt(a,M,1);
	ntt(b,M,1);
}
void INV(int a[],int c[],int len,int mod) {
	if (mod==1) {
		c[0]=inv(a[0]);
		return;
	}
	INV(a,c,len,(mod+1)>>1);
	int M,xj,chang=max(len,mod);
	for (M=1,xj=0;M<=(chang<<1);M<<=1,++xj);
	for (int i=0;i<M;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(xj-1));
	memset(T,0,sizeof T);
	copy(c,c+chang,T);
	ntt(a,M);
	ntt(T,M);
	for (int i=0;i<M;++i) {
		c[i]=(2ll-(giant)a[i]*T[i]%q)%q;
		if (c[i]<0) c[i]+=q;
		c[i]=(giant)c[i]*T[i]%q;
		if (c[i]<0) c[i]+=q;
	}
	ntt(c,M,1);
	ntt(a,M,1);
	for (int i=mod;i<M;++i) c[i]=0;
}
int main() {
#ifndef ONLINE_JUDGE
	freopen("7.in","r",stdin);
	freopen("my.out","w",stdout);
#endif
	Wn[0][0]=Wn[0][1]=1;
	for (int i=1;i<maxj;++i) {
		Wn[i][0]=mi(basic,(q-1)>>i);
		Wn[i][1]=inv(Wn[i][0]);
	}
	int n=read();
	for (int i=0;i<=n;++i) A[i]=read();
	int m=read();
	for (int i=0;i<=m;++i) B[i]=read();
	if (n<m) {
		puts("0");
		print(A,n);
		return 0;
	}
	for (int i=0;i<=n-m;++i) AR[i]=A[n-i],BR[i]=B[m-i];
	INV(BR,IBR,m,n-m+1);
	TIME(AR,n-m,IBR,n-m,DR);
	for (int i=0;i<=n-m;++i) D[i]=DR[n-m-i];
	TIME(B,m,D,n-m,C);
	for (int i=0;i<m;++i) R[i]=((giant)(A[i]-C[i])%q+q)%q;
	print(D,n-m);
	print(R,m-1);
	return 0;
}