多项式除法

问题

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

即求多项式(A(x))(B(x))带余除法

做法

首先,对于任意多项式$$P(x)=sum _{i=0}^na_ix^i$$,显然有$$P^R(x)=x^nP(frac{1}{x})=sum {i=0}^na{n-i}x^i$$。

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

[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 Dle n-m)(deg D^Rle n-m),所以也不会出问题。

这样就可以轻松地算出$$D^R=A_R(B^R)^{-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;
}
原文地址:https://www.cnblogs.com/owenyu/p/6724611.html