【洛谷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;
}
败得义无反顾,弱得一无是处
原文地址:https://www.cnblogs.com/chenxiaoran666/p/MTT.html