进军多项式(二):任意模数NTT

任意模数NTT

已知多项式 F ( x ) , G ( x ) F(x),G(x) F(x),G(x),求 H ( x ) = F ( x ) G ( x ) ( m o d    P ) H(x)=F(x)G(x)(mod\;P) H(x)=F(x)G(x)(modP)
其中 P P P是给定模数,不保证是NTT模数
我们先取三个NTT模数
常用的是(998244353,469762049,1004535809)
因为他们的原根都是3
我们分别用这三个模数做NTT,然后再用中国剩余定理合并
不过这三个模数乘起来爆longlong
所以我们必须一步一合并
我们可以对于每一位分别考虑
A = 998244353 , B = 469762049 , C = 1004535809 A=998244353,B=469762049,C=1004535809 A=998244353,B=469762049,C=1004535809
设这一位的真实值是 x x x
我们已知
{ x ≡ x 1 ( m o d    A ) x ≡ x 2 ( m o d    B ) x ≡ x 3 ( m o d    C ) \left\{\begin{matrix}x\equiv x_1(mod\;A)\\x\equiv x_2(mod\;B)\\x\equiv x_3(mod\;C)\end{matrix}\right. xx1(modA)xx2(modB)xx3(modC)
也就是
第一个方程的通解是 x 1 + k 1 A x_1+k_1A x1+k1A
也就是
x 1 + k 1 A ≡ x 2 ( m o d    B ) x_1+k_1A\equiv x_2(mod \;B) x1+k1Ax2(modB)
移项可得
k 1 ≡ x 2 − x 1 A k_1\equiv\frac{x_2-x_1}{A} k1Ax2x1
我们代入可以的到 X = x 1 + k 1 A X=x_1+k_1A X=x1+k1A
那么前两个方程的通解是
X + K ( A B ) X+K(AB) X+K(AB)
类似的代入第三个方程
X + K ( A B ) ≡ x 3 ( m o d    C ) X+K(AB)\equiv x_3(mod\;C) X+K(AB)x3(modC)
移项
K ≡ x 3 − X A B K\equiv \frac{x_3-X}{AB} KABx3X
那么我们可以算出一个特解 x 0 x_0 x0
总的通解式就是
x 0 + k 0 A B C ( m o d    A B C ) x_0+k_0ABC(mod\;ABC) x0+k0ABC(modABC)
因为 A B C ABC ABC一定大于 x 0 x_0 x0所以答案就是 x 0 x_0 x0
因为要做 9 9 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;
}
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值