Loading

https://blog.csdn.net/jwg2732/article/details/123179632?spm=1001.2014.3001.5501

任意模数NTT

已知多项式\(F(x),G(x)\),求\(H(x)=F(x)G(x)(mod\;P)\)
其中\(P\)是给定模数,不保证是NTT模数
我们先取三个NTT模数
常用的是(998244353,469762049,1004535809)
因为他们的原根都是3
我们分别用这三个模数做NTT,然后再用中国剩余定理合并
不过这三个模数乘起来爆longlong
所以我们必须一步一合并
我们可以对于每一位分别考虑
\(A=998244353,B=469762049,C=1004535809\)
设这一位的真实值是\(x\)
我们已知

\[\left\{\begin{matrix}x\equiv x_1(mod\;A)\\x\equiv x_2(mod\;B)\\x\equiv x_3(mod\;C)\end{matrix}\right. \]

也就是
第一个方程的通解是\(x_1+k_1A\)
也就是

\[x_1+k_1A\equiv x_2(mod \;B) \]

移项可得

\[k_1\equiv\frac{x_2-x_1}{A} \]

我们代入可以的到\(X=x_1+k_1A\)
那么前两个方程的通解是

\[X+K(AB) \]

类似的代入第三个方程

\[X+K(AB)\equiv x_3(mod\;C) \]

移项

\[K\equiv \frac{x_3-X}{AB} \]

那么我们可以算出一个特解\(x_0\)
总的通解式就是

\[x_0+k_0ABC(mod\;ABC) \]

因为\(ABC\)一定大于\(x_0\)所以答案就是\(x_0\)
因为要做\(9\)次NTT,常数很大

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL A = 998244353;
const LL B = 1004535809;
const LL C = 469762049; 
#define Poly vector<LL>
#define len(x) (int)x.size()
#define turn(x,v) x.resize(v)
const int N = 4e5+7;
LL Pow(LL a,LL b,LL mod)
{
	LL res=1;
	while(b)
	{
		if(b&1) res=res*a%mod;
		a=1ll*a*a%mod;
		b>>=1;
	}
	return res;
}
int tr[N];
void Retry(int n)
{
	for(int i=0;i<n;i++)
	tr[i]=((tr[i>>1]>>1)|((i&1)?(n>>1):0));
}
void NTT(Poly &f,int opt,LL mod)
{
	int n=len(f);
	LL G=3,IG=Pow(G,mod-2,mod);
	for(int i=0;i<n;i++)
	if(i<tr[i]) swap(f[i],f[tr[i]]);
	for(int k=2;k<=n;k<<=1)
	{
		int l=(k>>1);
		LL w=Pow(opt==1?G:IG,(mod-1)/k,mod);
		for(int i=0;i<n;i+=k)
		{
			LL g=1;
			for(int j=i;j<i+l;j++)
			{
				LL v=f[j+l]*g%mod;
				f[j+l]=(f[j]-v+mod)%mod;
				f[j]=(f[j]+v)%mod;
				g=g*w%mod;
			}
		}
	}
}
Poly Times(Poly f,Poly g,int len,LL mod)
{
	NTT(f,1,mod);NTT(g,1,mod);
	for(int i=0;i<len;i++)
	f[i]=f[i]*g[i]%mod;
	NTT(f,-1,mod);
	LL invn=Pow(len,mod-2,mod);
	for(int i=0;i<len;i++)
	f[i]=f[i]*invn%mod;
	return f;
}
LL mul(LL a,LL b,LL mod)
{
	LL res=0;
	while(b)
	{
		if(b&1) res=(res+a)%mod;
		a=(a+a)%mod;
		b>>=1; 
 	}
	return res; 
}
LL Inv(LL x,LL mod)
{
	LL a=x,b=mod-2,res=1;
	while(b)
	{
		if(b&1) res=res*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return res;
}
LL Merge(LL x1,LL x2,LL x3,LL mod)
{
	LL k1=((x2-x1)%B+B)%B*Inv(A,B)%B;
	LL x0=(x1+k1*A%(A*B))%(A*B);
	LL k0=((x3-x0)%C+C)%C*Inv(A*B%C,C)%C; 
	return (x0+k0*A%mod*B%mod)%mod;
}
inline LL read() {
  LL x = 0, w = 1;
  char ch = 0;
  while (ch < '0' || ch > '9') {  // ch 不是数字时
    if (ch == '-') w = -1;        // 判断是否为负
    ch = getchar();               // 继续读入
  }
  while (ch >= '0' && ch <= '9') {  // ch 是数字时
    x = x * 10 + (ch - '0');  // 将新读入的数字’加’在 x 的后面
    // x 是 int 类型,char 类型的 ch 和 ’0’ 会被自动转为其对应的
    // ASCII 码,相当于将 ch 转化为对应数字
    // 此处也可以使用 (x<<3)+(x<<1) 的写法来代替 x*10
    ch = getchar();  // 继续读入
  }
  return x * w;  // 数字 * 正负号 = 实际数值
} 
void Put(Poly f,int len)
{
	for(int i=0;i<len;i++)
	printf("%lld ",f[i]);
	printf("\n");
}
Poly FreeMul(Poly f,Poly g,int len,LL mod)
{
	Retry(len);
	Poly X1=Times(f,g,len,A);
	Poly X2=Times(f,g,len,B);
	Poly X3=Times(f,g,len,C);
//	Put(X1,len);
	Poly res;
	turn(res,len);
	for(int i=0;i<len;i++)
	res[i]=Merge(X1[i],X2[i],X3[i],mod);
	return res;
}
int n,m,p;
int main()
{
	cin>>n>>m>>p;
	n++;m++;
	Poly f,g;
	int k=n+m;
	int T;
	for(T=1;T<k;T<<=1);
	turn(f,T);turn(g,T);
	for(int i=0;i<n;i++)
	f[i]=read();
	for(int i=0;i<m;i++)
	g[i]=read();
	Poly res=FreeMul(f,g,T,p);
	for(int i=0;i<n+m-1;i++)
	printf("%lld ",res[i]);
	return 0;
}
posted @ 2022-03-01 13:19  Larunatrecy  阅读(244)  评论(0)    收藏  举报