多项式除法
给出一个 n n n 次多项式 F ( x ) F(x) F(x) 和一个 m m m 次多项式 G ( x ) G(x) G(x),求出一个 n − m n-m n−m 次多项式 Q ( x ) Q(x) Q(x) 和一个小于 m m m 次的多项式 R ( x ) R(x) R(x),满足 F = G Q + R F= GQ+R F=GQ+R。
对于一个 n n n 次多项式 A ( x ) A(x) A(x), x n A ( 1 x ) x^nA(\frac 1 x) xnA(x1) 相当于将 A A A 的系数翻转后得到的多项式,即第 0 0 0 项系数变第 n n n 项系数,第 1 1 1 项系数变第 n − 1 n-1 n−1 项系数……不妨将 A ( x ) A(x) A(x) 翻转后的多项式记为 A R ( x ) A_R(x) AR(x)
将上面那个柿子左右同时翻转前
n
n
n 项系数,得到:
x
n
F
(
1
x
)
=
x
n
(
G
(
1
x
)
Q
(
1
x
)
+
R
(
1
x
)
)
x
n
F
(
1
x
)
=
x
m
G
(
1
x
)
x
n
−
m
Q
(
1
x
)
+
x
n
−
m
+
1
x
m
−
1
R
(
1
x
)
F
R
(
x
)
=
G
R
(
x
)
Q
R
(
x
)
+
x
n
−
m
+
1
R
R
(
x
)
\begin{aligned} x^nF(\frac1 x)&=x^n(G(\frac 1 x)Q(\frac 1 x)+R(\frac 1 x))\\ x^nF(\frac1 x)&=x^mG(\frac 1 x)x^{n-m}Q(\frac 1 x)+x^{n-m+1}x^{m-1}R(\frac 1 x)\\ F_R(x)&=G_R(x)Q_R(x)+x^{n-m+1}R_R(x)\\ \end{aligned}
xnF(x1)xnF(x1)FR(x)=xn(G(x1)Q(x1)+R(x1))=xmG(x1)xn−mQ(x1)+xn−m+1xm−1R(x1)=GR(x)QR(x)+xn−m+1RR(x)
假如让这个柿子对
x
n
−
m
+
1
x^{n-m+1}
xn−m+1 取模,那么后面的
x
n
−
m
+
1
R
R
(
x
)
x^{n-m+1}R_R(x)
xn−m+1RR(x) 就没了,然后
F
R
,
G
R
,
Q
R
F_R,G_R,Q_R
FR,GR,QR 保留小于
x
n
−
m
+
1
x^{n-m+1}
xn−m+1 的项,变成:
F
R
(
x
)
≡
G
R
(
x
)
Q
R
(
x
)
(
m
o
d
x
n
−
m
+
1
)
F_R(x)\equiv G_R(x)Q_R(x)\pmod {x^{n-m+1}}
FR(x)≡GR(x)QR(x)(modxn−m+1)
诶,我们要求的 Q ( x ) Q(x) Q(x) 就是一个 n − m n-m n−m 次的多项式,这不就刚好了嘛。
转化一下,就是:
Q
R
(
x
)
≡
F
R
(
x
)
G
R
(
x
)
(
m
o
d
x
n
−
m
+
1
)
Q_R(x)\equiv \frac {F_R(x)} {G_R(x)}\pmod {x^{n-m+1}}
QR(x)≡GR(x)FR(x)(modxn−m+1)
于是对 G R ( x ) G_R(x) GR(x) 求个逆,再乘上 F R ( x ) F_R(x) FR(x) 得到 Q R ( x ) Q_R(x) QR(x),翻转回去就得到了 Q ( x ) Q(x) Q(x),然后又可以把 R ( x ) R(x) R(x) 求出来了。
代码如下:
#include <cstdio>
#include <cstring>
#include <vector>
#include <cmath>
#include <algorithm>
using namespace std;
#define maxn 400010
#define mod 998244353
#define bin(x) (1<<(x))
int n,m;
int inv[maxn];
int ksm(int x,int y){int re=1;for(;(y&1?re=1ll*re*x%mod:0),y;y>>=1,x=1ll*x*x%mod);return re;}
#define INV(x) ksm(x,mod-2)
struct NTT{
vector<int> w[30];NTT(){
inv[1]=1;for(int i=2;i<=maxn-10;i++)inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
for(int i=1,wn;i<=19;i++){
w[i].resize(bin(i));w[i][0]=1;wn=ksm(3,(mod-1)/bin(i));
for(int j=1;j<bin(i-1);j++)w[i][j]=1ll*w[i][j-1]*wn%mod;
}
}
int limit,r[maxn];void dft(int *f,int lg,int type=0){
limit=bin(lg);if(type)reverse(f+1,f+limit);
for(int i=1;i<limit;i++){r[i]=(r[i>>1]>>1)|((i&1)<<(lg-1));if(i<r[i])swap(f[i],f[r[i]]);}
for(int mid=1,Lg=1;mid<limit;mid<<=1,Lg++)for(int j=0;j<limit;j+=(mid<<1))for(int i=0;i<mid;i++)
{int t=1ll*f[j+i+mid]*w[Lg][i]%mod;f[j+i+mid]=(f[j+i]-t+mod)%mod;f[j+i]=(f[j+i]+t)%mod;}
}
}ntt;
int A[maxn],B[maxn],M;
struct POLY{
vector<int> a;int len;void rs(int N){a.resize(len=N);}POLY(){rs(M);}
int &operator [](int x){return a[x];}
void dft(int *A_,int lg,int ln){for(int i=0;i<bin(lg);i++)A_[i]=(i<min(ln,len)?a[i]:0);ntt.dft(A_,lg);}
void idft(int *A_,int lg,int ln){ntt.dft(A_,lg,1);rs(ln);for(int i=0;i<ln;i++)a[i]=1ll*A_[i]*inv[bin(lg)]%mod;}
const POLY Mul(POLY b,int ln){
int lg=ceil(log2(ln*2-1));dft(A,lg,ln);b.dft(B,lg,ln);
for(int i=0;i<bin(lg);i++)B[i]=1ll*A[i]*B[i]%mod;b.idft(B,lg,ln);return b;
}
}F,G,Q,R;
void getinv(POLY &f,POLY &g,int ln)
{
if(ln==1){g.rs(1);g[0]=INV(f[0]);return;}getinv(f,g,(ln+1)>>1);
int lg=ceil(log2(ln*2-1));f.dft(A,lg,ln);g.dft(B,lg,ln);
for(int i=0;i<bin(lg);i++)B[i]=1ll*(2-1ll*A[i]*B[i]%mod)%mod*B[i]%mod;g.idft(B,lg,ln);
}
POLY rev(POLY f){reverse(f.a.begin(),f.a.end());return f;}
POLY getinv(POLY f,int ln=M){POLY g;getinv(f,g,ln);return g;}
POLY getdiv(POLY &f,POLY &g){return rev(rev(f).Mul(getinv(rev(g),f.len-g.len+1),f.len-g.len+1));}
POLY getmod(POLY &f,POLY &g,POLY &q){POLY r=g.Mul(q,f.len);for(int i=0;i<f.len;i++)r[i]=(f[i]-r[i]+mod)%mod;return r;}
int main()
{
scanf("%d %d",&n,&m);F.rs(n+1);G.rs(m+1);M=n+1;
for(int i=0;i<=n;i++)scanf("%d",&F[i]);
for(int i=0;i<=m;i++)scanf("%d",&G[i]);
Q=getdiv(F,G);R=getmod(F,G,Q);
for(int i=0;i<=n-m;i++)printf("%d ",Q[i]);printf("\n");
for(int i=0;i<m;i++)printf("%d ",R[i]);
}