【洛谷4245】【模板】任意模数多项式乘法

点此看题面

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

中国剩余定理

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

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

\[\begin{cases}x\equiv r_1(mod\ p_1),\\x\equiv r_2(mod\ p_2),\\x\equiv r_3(mod\ p_3)\end{cases} \]

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

以\(p_1,p_2\)的合并为例：

\[k\times p_1+r_1\equiv r_2(mod\ p_2)\\ k\equiv\frac{r_2-r_1}{p_1}(mod\ p_2) \]

既然得出了\(k\)，那么就有\(r_{1+2}=k\times p_1+r_1\)。

同理我们可以合并\(x\equiv r_{1+2}(mod\ p_1p_2)\)和\(x\equiv r_{3}(mod\ p_3)\)。

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

然而实际的\(x\)小于\(p_1p_2p_3\)，所以\(x\)就等于\(k\times p_1p_2+r_{1+2}\)。

直接计算\(k\times 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;
}