[多项式] 任意模数NTT

前言

两种解决方法
MTT & 拆系数FFT

FFT很讨厌了,卡精度,弃
于是学习易卡常的MTT
钻研一下用vector写不会炸裂的

总览

利用中国剩余定理得到一个超大的结果( > 1 0 23 >10^{23} >1023)
然后一模,一切都解决了

  • step1
    在三个模数(两个不够 > 1 0 23 >10^{23} >1023)下分别得到三个多项式
  • step2
    合并三个多项式

合并

考虑合并多项式
对于每一位,有
x ≡ a 1 ( m o d p 1 ) x ≡ a 2 ( m o d p 2 ) x ≡ a 3 ( m o d p 3 ) x\equiv a_1 \pmod {p_1}\\ x\equiv a_2 \pmod {p_2}\\ x\equiv a_3 \pmod {p_3} xa1(modp1)xa2(modp2)xa3(modp3)
合并三项发现炸long long
考虑合并前两项

x ≡ a 2 p 1 p 1 ′ + a 1 p 2 p 2 ′ ( m o d p 1 p 2 ) x\equiv a_2p_1p_1'+a_1p_2p_2'\pmod {p_1p_2} xa2p1p1+a1p2p2(modp1p2)
其中 p 1 p 1 ′ ≡ 1 ( m o d p 2 ) p_1p_1'\equiv1\pmod{p_2} p1p11(modp2) p 2 p 2 ′ ≡ 1 ( m o d p 1 ) p_2p_2'\equiv1\pmod{p_1} p2p21(modp1)
a 2 p 1 p 1 ′ + a 1 p 2 p 2 ′ = a 4 a_2p_1p_1'+a_1p_2p_2'=a_4 a2p1p1+a1p2p2=a4 x = a 5 p 1 p 2 + a 4 x=a_5p_1p_2+a_4 x=a5p1p2+a4

考虑第三个同余式
a 5 p 1 p 2 + a 4 ≡ a 3 ( m o d p 3 ) a 5 = a 3 − a 4 p 1 p 2 ( m o d p 3 ) a_5p_1p_2+a_4\equiv a_3\pmod{p_3}\\ a_5=\frac{a_3-a_4}{p_1p_2}\pmod{p_3} a5p1p2+a4a3(modp3)a5=p1p2a3a4(modp3)

考虑是否被模掉一些东西
发现 a 4 , a 5 a_4,a_5 a4,a5都是真实值,可以代入 x = a 5 p 1 p 2 + a 4 x=a_5p_1p_2+a_4 x=a5p1p2+a4得到答案

太菜了
肝不动了

#include<bits/stdc++.h>
using namespace std;
#define in Read()
#define int long long
int in{
	int i=0,f=1;char ch=0;
	while(!isdigit(ch)&&ch!='-') ch=getchar();
	if(ch=='-') ch=getchar(),f=-1;
	while(isdigit(ch)) i=(i<<1)+(i<<3)+ch-48,ch=getchar();
	return i*f;
}

const int N=1e6+5;
typedef long long ll;
typedef vector<int> poly;
int n,m,pmod,lim,rev[N];
int mod[5]={496762049,998244353,1004535809};
poly A[3],B[3];
ll a1,a2;

void print(poly f){
	for(int i=0;i<f.size();++i) printf("%lld ",f[i]);
	puts("");
	return ;
}

int add(int a,int b,int id){return a+b>=mod[id]?a+b-mod[id]:a+b;}
int dec(int a,int b,int id){return a<b?a-b+mod[id]:a-b;}
int mul(int a,int b,int id){return 1ll*a*b%mod[id];}
int qpw(int a,int b,int id){
	int res=1;
	while(b){
		if(b&1) res=mul(res,a,id);
		a=mul(a,a,id);
		b>>=1;
	}
	return res;
}

void init(int deg){
	lim=1;
	while(lim<deg) lim<<=1;
	for(int i=0;i<lim;++i)
		rev[i]=(rev[i>>1]>>1)|((i&1)?lim>>1:0);
	return;
}

void NTT(poly &f,bool sgn,int id){
	for(int i=0;i<lim;++i)
		if(i<rev[i]) swap(f[i],f[rev[i]]);
	for(int len=1;len<lim;len<<=1){
		int siz=len<<1;
		int wn=qpw(3,(mod[id]-1)/siz,id);
		for(int l=0;l<lim;l+=siz){
			int w=1;
			for(int i=l;i<l+len;++i){
				int a=f[i],b=mul(f[i+len],w,id);
				f[i]=add(a,b,id);
				f[i+len]=dec(a,b,id);
				w = mul(w, wn, id);
			}
		}
	}
	if(sgn) return;
	int INV=qpw(lim,mod[id]-2,id);
	reverse(f.begin()+1,f.end());
	for(int i=0;i<lim;++i)
		f[i] = mul(f[i], INV, id);
	return;
}

poly mul(poly a,poly b,int id){
	int deg=a.size()+b.size()-1;
	if(a.size()<=32||b.size()<=32){
		poly c(deg,0);
		for(int i=0;i<=a.size();++i)
			for(int j=0;j<b.size();++j)
				c[i + j] = add(c[i + j], mul(a[i], b[j],id),id);
		return c;
	}
	init(deg);
	a.resize(lim),NTT(a,1,id);
	b.resize(lim),NTT(b,1,id);
	for(int i=0;i<lim;++i)
		a[i] = mul(a[i], b[i],id);
	NTT(a,0,id);
	a.resize(deg);
	return a;
}

ll sqpw(ll a,ll b,ll p){
	ll res=1;
	while(b){
		if(b&1) res=res*a%p;
		a=a*a%p;
		b>>=1;
	}
	return res;
}

ll inv(ll x,ll p){
	return sqpw(x%p,p-2,p);
}

signed main(){
	n=in,m=in,pmod=in;
	init(n+m);
	for(int i=0;i<=n;++i){
		int x=in;
		A[0].push_back(x);
		A[1].push_back(x);
		A[2].push_back(x);
	}
	for(int i=0;i<=m;++i) {
		int x=in;
		B[0].push_back(x);
		B[1].push_back(x);
		B[2].push_back(x);
	}
	for(int id=0;id<3;++id) A[id]=mul(A[id],B[id],id);
	const ll p1=mod[0],p2=mod[1],p3=mod[2],p1p2=1ll*mod[0]*mod[1];
	const ll inv1=inv(p1,p2),inv2=inv(p2,p1),inv12=inv(p1p2,p3);
	for(int i=0;i<=n+m;++i){
		ll a=A[0][i],b=A[1][i],c=A[2][i];
		ll d=(inv1*p1%p1p2*b%p1p2+inv2*p2%p1p2*a%p1p2+p1p2)%p1p2;
		ll e=(c-d+p3)%p3*inv12%p3;
		printf("%lld ",(p1p2%pmod*e%pmod+d+pmod)%pmod);
	}
	return 0;
}
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值