lp3803 【模板】多项式乘法(FFT)(NTT实现)

程序运行效率。(评测系统为duck.ac)

我们现在已经有一种快速完成多项式卷积的方案了。然而这种方案还存在两个缺陷。
首先,是浮点型运算带来的精度误差。这在多次乘积之后会不断扩大,最终到一个无法解决的地步。
第二,就是结构体运算带来的大常数。这也使得这道题的时间限制变得比较紧张。
于是,我们考虑另辟蹊径,使用某一种大模数来完成整数多项式卷积。
快速数论变换(Fast Number Transform)就是一种这样的算法。

回到FFT的本质,FFT是一种精确选取取值以后代值优化系数表达和点值表达相互转化的方案。我们尝试在模意义下找到一些满足与单位复根性质类似的数。
一种在数论上被称为原根的数可以满足这一族性质。
首先我们有原根的定义。
形式化的说,若\(g\)模\(p\)的阶等于\(\phi(p)\),则称\(g\)为\(p\)的一个原根。
具体来说,我们认为,若\(g^x \equiv 1\ (mod\ p) \)总是有解,那么我们称\(g\)是质数\(p\)的一个原根。
对于\(\forall g\in Z,p\in P\),总是有\(g^{p-1}\equiv 1\ (mod\ p) \)
那么,对于一个\(p=t2^{k}+1\)我们有,\(g^{\frac{p-1}{2^i}}\ (0\le i\le k)\)在模$latexp$意义下始终有正整数解。这一点和单位复根类似。
通过选取满足上述条件的质数,即\(p=t2^{k}+1\),我们还有这样的性质:
消去引理:
$$g^{\frac{dk(p-1)}{dn}}\equiv g^{\frac{k(p-1)}{n}}$$
折半引理:
$$(g^{\frac{(k+\frac{n}{2})(p-1)}{n}})^2\equiv (g^{\frac{k(p-1)}{n}}g^{\frac{p-1}{2}})^2\equiv (g^{\frac{k(p-1)}{n}}(p-1))^2\equiv (g^{\frac{k(p-1)}{n}})^2$$
求和引理:
$$\sum_{j=0}^{n-1}(g^{\frac{k(p-1)}{n}})^j\equiv \frac{(g^{\frac{k(p-1)}{n}})^n-1}{g^{\frac{k(p-1)}{n}}-1}\equiv \frac{g^{k(p-1)}-1}{g^{\frac{k(p-1)}{n}}-1}\equiv \frac{1^k-1}{g^{\frac{k(p-1)}{n}}}\equiv 0$$
折半引理的一个推论:
$$\sum_{i=0}^{n-1}g^{\frac{i(p-1)}{n}}\equiv 0$$
我们惊讶地发现原根具有单位复根具有的三个定理和一个推论。这说明,我们可以将快速傅里叶变换的操作以原根的形式拓展到模意义下。
我们令:
$$\omega_{n}^{k}=g^{\frac{k(p-1)}{n}}$$
那么我们就可以简单地将快速傅里叶变换中的单位复根全部替换为原根,这样就实现了快速数论变换。
对于这里的模数\(p\),我们可以选取998244353,它是\(2^{23}*119+1\)。

#include<iostream>
#include<cstdio>
#define Swap(A,B) (A^=B^=A^=B)
const long long P=998244353;
const long long g0=3,gi=332748118;
//*****非常重要:gi是332748118!!!!!!! 
int L=1,R[1<<21|1];
long long invL;
inline int pw(int A,int X){
	long long BS=A,RT=1;
	while(X){
		if(X&1){
			RT=RT*BS%P;
		}
		BS=BS*BS%P;
		X>>=1;
	}
	return RT;
}
inline void prpr(int LEN){
	int B=0;
	while(L<=LEN){
		L<<=1;
		++B;
	}
	invL=P-(P-1)/L;
	for(int i=0;i<L;++i){
		R[i]=R[i>>1]>>1|(i&1)<<(B-1);
	}
}
inline void FNTT(int *A,int typ){
	for(int i=0;i<L;++i){
		if(R[i]<i){
			Swap(A[R[i]],A[i]);
		}
	}
	int gn,g,X,Y,M;
	for(int i=2;i<=L;i<<=1){
		M=i>>1;
		gn=pw(~typ?g0:gi,(P-1)/i);
		for(int j=0;j<L;j+=i){
			g=1;
			for(int k=0;k<M;++k,g=1ll*g*gn%P){
				X=A[j+k],Y=1ll*g*A[M+j+k]%P;
				A[j+k]=(X+Y)%P,A[M+j+k]=(X-Y)%P;
			}
		}
	}
}
int n,m,a[1<<21|1],b[1<<21|1];
void init(){
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;++i){
		scanf("%d",a+i);
	}
	for(int i=0;i<=m;++i){
		scanf("%d",b+i);
	}
	prpr(n+m);
	FNTT(a,1);
	FNTT(b,1);
	for(int i=0;i<L;++i){
		a[i]=1ll*a[i]*b[i]%P;
	}
	FNTT(a,-1);
	for(int i=0;i<=n+m;++i){
		printf("%d ",1ll*(a[i]+P)*invL%P);
	}
}

int main(){
	init();
	return 0;
}

发表评论

您的电子邮箱地址不会被公开。