zoukankan      html  css  js  c++  java
  • 【洛谷4512】【模板】多项式除法

    点此看题面

    大致题意: 给定多项式(F(x),G(x)),求(Q(x),R(x))满足(F(x)=Q(x)*G(x)+R(x))

    前言

    公式的推导看起来十分自然,然而远不是我这种蒟蒻所能想到的。

    推式子

    定义(F_R(x))表示将(F(x))系数翻转得到的式子(即(F_R(x)[i]=F(x)[n-i])),显然有:

    [F_R(x)=x^n*F(frac1x) ]

    然后对于原式我们进行转化:

    [F(x)=Q(x)*G(x)+R(x) ]

    [F(frac1x)=Q(frac1x)*G(frac1x)+R(frac1x) ]

    [x^n*F(frac1x)=(x^{n-m}*Q(frac1x))*(x^m*G(frac1x))+x^{n-m+1}*(x^{m-1}*R(frac1x)) ]

    [F_R(x)=Q_R(x)*G_R(x)+x^{n-m+1}*R_R(x) ]

    我们把这个式子两边同时向(x^{n-m+1})取模,这样就可以消去(R_R(x))这一项,同时依然能保留(Q_R(x))的有效系数,得到:

    [F_R(x)equiv Q_R(x)*G_R(x)(mod x^{n-m+1}) ]

    [Q_R(x)equiv F_R(x)*G_R^{-1}(x)(mod x^{n-m+1}) ]

    于是只要通过多项式乘法逆就能求出(Q(x))了,而(R(x))的计算是很简单的:

    [R(x)=F(x)-Q(x)*G(x) ]

    代码

    #include<bits/stdc++.h>
    #define Tp template<typename Ty>
    #define Ts template<typename Ty,typename... Ar>
    #define Reg register
    #define RI Reg int
    #define Con const
    #define CI Con int&
    #define I inline
    #define W while
    #define N 100000
    #define X 998244353
    using namespace std;
    int n,m,f[N+5],g[N+5],q[N+5],r[N+5];
    class FastIO
    {
    	private:
    		#define FS 100000
    		#define tc() (A==B&&(B=(A=FI)+fread(FI,1,FS,stdin),A==B)?EOF:*A++)
    		#define pc(c) (C==E&&(clear(),0),*C++=c)
    		#define D isdigit(c=tc())
    		int T;char c,*A,*B,*C,*E,FI[FS],FO[FS],S[FS];
    	public:
    		I FastIO() {A=B=FI,C=FO,E=FO+FS;}
    		Tp I void read(Ty& x) {x=0;W(!D);W(x=(x<<3)+(x<<1)+(c&15),D);}
    		Tp I void write(Ty x) {W(S[++T]=x%10+48,x/=10);W(T) pc(S[T--]);}
    		Tp I void write(Con Ty& x,Con char& y) {write(x),pc(y);}
    		I void clear() {fwrite(FO,1,C-FO,stdout),C=FO;}
    }F;
    I int QP(RI x,RI y) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
    #define QI(x) QP(x,X-2)
    namespace Poly
    {
    	#define Init(n) P=1,L=0;W(P<=2*(n)) P<<=1,++L;
    		for(i=0;i^P;++i) A[i]=B[i]=0,R[i]=(R[i>>1]>>1)|((i&1)<<L-1);
    	int PR=3,IPR=QI(3),P,L,R[4*N+5],A[4*N+5],B[4*N+5],p[N+5];
    	I void NTT(int *s,CI t)//NTT
    	{
    		RI i,j,k,x,y,U,S;for(i=0;i^P;++i) i<R[i]&&(x=s[i],s[i]=s[R[i]],s[R[i]]=x);
    		for(i=1;i^P;i<<=1) for(U=QP(t,(X-1)/(i<<1)),j=0;j^P;j+=i<<1) for(S=1,k=0;k^i;
    			++k,S=1LL*S*U%X) s[j+k]=((x=s[j+k])+(y=1LL*S*s[i+j+k]%X))%X,s[i+j+k]=(x-y+X)%X;
    	}
    	I void Mul(CI n,int *a,int *b,int *c)//多项式乘法
    	{
    		RI i;Init(n);for(i=0;i<=n;++i) A[i]=a[i],B[i]=b[i];
    		for(NTT(A,PR),NTT(B,PR),i=0;i^P;++i) A[i]=1LL*A[i]*B[i]%X;
    		RI t=QI(P);for(NTT(A,IPR),i=0;i<=n;++i) c[i]=1LL*A[i]*t%X;
    	}
    	I void Inv(CI n,int *a,int *b)//多项式求逆
    	{
    		if(!n) return (void)(b[0]=QI(a[0]));RI i;Inv(n>>1,a,b);
    		Init(n);for(i=0;i<=n;++i) A[i]=a[i],B[i]=b[i];
    		for(NTT(A,PR),NTT(B,PR),i=0;i^P;++i) B[i]=(2LL*B[i]-1LL*A[i]*B[i]%X*B[i]%X+X)%X;
    		RI t=QI(P);for(NTT(B,IPR),i=0;i<=n;++i) b[i]=1LL*B[i]*t%X;
    	}
    	I void Rev(CI n,int *a)//翻转系数
    	{
    		RI i;for(i=0;i<=n;++i) A[i]=a[n-i];for(i=0;i<=n;++i) a[i]=A[i];
    	}
    	I void Div(int *f,int *g,int *q,int *r)//多项式除法
    	{
    		RI i;Rev(n,f),Rev(m,g),Inv(n-m,g,p),Mul(n-m,f,p,q);Rev(n-m,q);//通过多项式求逆得到q
    		Rev(n,f),Rev(m,g),Mul(n,g,q,p);for(i=0;i^m;++i) r[i]=(f[i]-p[i]+X)%X;//余数(r)=被除数(f)-除数(g)×商(q)
    	}
    }
    int main()
    {
    	RI i;for(F.read(n),F.read(m),i=0;i<=n;++i) F.read(f[i]);for(i=0;i<=m;++i) F.read(g[i]);//读入
    	for(Poly::Div(f,g,q,r),i=0;i<=n-m;++i) F.write(q[i]," 
    "[i==n-m]);//输出商
    	for(i=0;i^m;++i) F.write(r[i]," 
    "[i==m-1]);return F.clear(),0;//输出余数
    }
    
  • 相关阅读:
    POJ 3253 Fence Repair
    POJ 2431 Expedition
    NYOJ 269 VF
    NYOJ 456 邮票分你一半
    划分数问题 DP
    HDU 1253 胜利大逃亡
    NYOJ 294 Bot Trust
    NYOJ 36 最长公共子序列
    HDU 1555 How many days?
    01背包 (大数据)
  • 原文地址:https://www.cnblogs.com/chenxiaoran666/p/PolyDiv.html
Copyright © 2011-2022 走看看