zoukankan      html  css  js  c++  java
  • 【洛谷4245】【模板】任意模数多项式乘法

    点此看题面

    • 给定两个多项式(F(x),G(x)),求(F(x)*G(x))每一项系数在模(p)意义下的值。(​(p)可能是任意模数)
    • (n,mle10^5,ple10^9+9)

    中国剩余定理

    首先发现答案的上界为(O(np^2)),即(10^{23}),只要用三个质数就可以保证它们的乘积大于(10^{23})

    考虑我们分别在三个不同的(NTT)模数下做三次卷积,对于每一个位置上的值,都可以得出三个式子:

    [egin{cases}xequiv r_1(mod p_1),\xequiv r_2(mod p_2),\xequiv r_3(mod p_3)end{cases} ]

    然后我们可以利用这三个式子推出模(p_1p_2p_3)意义下的答案。

    (p_1,p_2)的合并为例:

    [k imes p_1+r_1equiv r_2(mod p_2)\ kequivfrac{r_2-r_1}{p_1}(mod p_2) ]

    既然得出了(k),那么就有(r_{1+2}=k imes p_1+r_1)

    同理我们可以合并(xequiv r_{1+2}(mod p_1p_2))(xequiv r_{3}(mod p_3))

    我们计算出一个新的(k),然后就可以得出(r_{1+2+3}=k imes p_1p_2+r_{1+2}),即(x)在模(p_1p_2p_3)意义下的值。

    然而实际的(x)小于(p_1p_2p_3),所以(x)就等于(k imes p_1p_2+r_{1+2})

    直接计算(k imes p_1p_2+r_{1+2})显然会爆(long long)(就和直接不取模(NTT)没区别了),我们可以在计算之前先给(p_1p_2)模上(p),就能解决这个问题了。

    代码:(O(nlogn))

    #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 LL long long
    using namespace std;
    int n,m,XX,a[2*N+5],b[N+5];
    I int QP(RI x,RI y,CI p=XX) {RI t=1;W(y) y&1&&(t=1LL*t*x%p),x=1LL*x*x%p,y>>=1;return t;}
    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--]);pc(' ');}
    		I void clear() {fwrite(FO,1,C-FO,stdout),C=FO;}
    }F;
    namespace Poly
    {
    	int P,L,R[N<<2];struct Poly_NTT
    	{
    		private:
    			int A[N<<2],B[N<<2];
    			I void T(int* s,CI op)
    			{
    				RI i,j,k,x,y,U,S;for(i=0;i^P;++i) i<R[i]&&(swap(s[i],s[R[i]]),0);
    				for(i=1;i^P;i<<=1) for(U=QP(QP(3,op,X),(X-1)/(i<<1),X),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;
    			}
    		public:
    			int X;I int operator [] (CI x) Con {return A[x];}
    			I void Mul(CI n,int* a,CI m,int* b)
    			{
    				RI i;for(i=0;i^P;++i) A[i]=B[i]=0;for(i=0;i<=n;++i) A[i]=a[i];for(i=0;i<=m;++i) B[i]=b[i];
    				for(T(A,1),T(B,1),i=0;i^P;++i) A[i]=1LL*A[i]*B[i]%X;
    				RI t=QP(P,X-2,X);for(T(A,X-2),i=0;i<=n+m;++i) A[i]=1LL*A[i]*t%X;
    			}
    	}Q[3];
    	I LL CRT(Con LL& r1,Con LL& p1,CI r2,CI p2,CI fg)//合并
    	{
    		RI k=1LL*((r2-r1)%p2+p2)*QP(p1%p2,p2-2,p2)%p2;//计算系数
    		if(!fg) return (p1*k+r1)%(p1*p2);return ((p1%XX)*k+r1)%XX;//fg=0直接计算,fg=1取模
    	}
    	I void Mul(CI n,int* a,CI m,int* b)//MTT
    	{
    		RI i;P=1,L=0;W(P<=n+m) P<<=1,++L;for(i=0;i^P;++i) R[i]=(R[i>>1]>>1)|((i&1)<<L-1);
    		Q[0].Mul(n,a,m,b),Q[1].Mul(n,a,m,b),Q[2].Mul(n,a,m,b);//三个NTT
    		for(i=0;i<=n+m;++i) a[i]=CRT(CRT(Q[0][i],Q[0].X,Q[1][i],Q[1].X,0),1LL*Q[0].X*Q[1].X,Q[2][i],Q[2].X,1);//每个位置合并答案
    	}
    }
    int main()
    {
    	Poly::Q[0].X=998244353,Poly::Q[1].X=469762049,Poly::Q[2].X=1004535809;//三个NTT模数
    	RI i;for(F.read(n),F.read(m),F.read(XX),i=0;i<=n;++i) F.read(a[i]);
    	for(i=0;i<=m;++i) F.read(b[i]);for(Poly::Mul(n,a,m,b),i=0;i<=n+m;++i) F.write(a[i]);return F.clear(),0;
    }
    
    败得义无反顾,弱得一无是处
  • 相关阅读:
    Java面试题:栈和队列的实现
    Java面试题:如何对HashMap按键值排序
    经典的Java基础面试题集锦
    9个Java初始化和回收的面试题
    20个高级Java面试题汇总
    Spring、Spring MVC、MyBatis整合文件配置详解2
    Spring、Spring MVC、MyBatis整合文件配置详解
    Spring:基于注解的Spring MVC
    margin百分比的相对值--宽度!
    jquery.cxSelect插件,城市没单位
  • 原文地址:https://www.cnblogs.com/chenxiaoran666/p/MTT.html
Copyright © 2011-2022 走看看