【瞎口胡】快速数论变换 NTT

在 FFT 中，因为是浮点数计算因此会掉精度。如果你不知道 FFT 是什么，请阅读这里。

如果在模意义下，我们可以选择不使用复平面的单位根，而是模意义下的单位根。

考虑单位根的性质：

• \(\omega_{n}^{n}=\omega_{n}^{0}=1\)
• 对于 \(n=2m\)，\(\omega_{n}^{k}=-\omega_{n}^{k+m}\)
• \(\omega_{n}^{k} = \omega_{2n}^{2k}\)

一个常见的 NTT 模数是 \(998244353=7 \times 17 \times 2^{23}+1\)（它是一个质数）。它有一个原根 \(3\)。别忘了，模 \(p\) 意义下的原根 \(a\) 满足 \(a^{\varphi(p)}\equiv 1 \pmod p\)，并且不存在更小的指数 \(x\) 使得 \(a^{x} \equiv 1 \pmod p\)。

因为 \(\varphi(p = 998244353) = 998244352 = 7 \times 17 \times 2^{23}\)，容易发现模意义下的 \(\omega_2 = a^{\frac 12 \varphi(p)}\)，\(\omega_{4} = a^{\frac 14 \varphi(p)}\)，以此类推。由 \(\varphi(p)\) 的唯一分解式可以发现，它最多有 \(2^{23}\) 次单位根。于是用 \(998244353\) 作为模数时，最多可以做 长度为 \(2^{23}\) 的多项式乘法。

模板 Luogu P1919 A*B Problem 升级版

题意

给定两个正整数 \(a,b\)，输出它们的乘积。

\(1 \leq a,b \leq 10^{10^6}\)

题解

把数看成多项式，做 NTT 即可。

# include <bits/stdc++.h>

const int N=4000010,INF=0x3f3f3f3f,MOD=998244353;

typedef long long ll;

int w[30],winv[30]; // 998244353 = 7*17*2^23 w[i] = 2^i 等分单位根 
typedef std::vector <int> Poly;
int rev[N];
char A[N],B[N];
Poly F,G;
int ans[N<<1],anslen;

inline int read(void){
	int res,f=1;
	char c;
	while((c=getchar())<'0'||c>'9')
		if(c=='-')f=-1;
	res=c-48;
	while((c=getchar())>='0'&&c<='9')
		res=res*10+c-48;
	return res*f;
}
inline int qpow(ll d,ll p){
	ll ans=1;
	while(p){
		if(p&1)
			ans=ans*d%MOD;
		p>>=1,d=d*d%MOD;
	}
	return ans;
}
inline void init(void){
	for(int i=0;i<=23;++i){
		w[i]=qpow(3,(MOD-1)/(1<<i)),winv[i]=qpow(w[i],MOD-2);
        // 2^i 次单位根及其逆元
	}
	return;
}
inline void Revchange(Poly &f,int len){ // 蝴蝶变换
	for(int i=1;i<len;++i){
		rev[i]=(rev[i>>1]>>1)|((i&1)*(len>>1));
	}
	for(int i=0;i<len;++i){
		if(i<rev[i])
			std::swap(f[i],f[rev[i]]);
	}
	return;
}
void NTT(Poly &f,int len,int op){
	Revchange(f,len);
	for(int h=2,hpow=1;h<=len;h<<=1,++hpow){
		int wh=((op==1)?w[hpow]:winv[hpow]); // idft 的时候记得取逆元
		for(int j=0;j<len;j+=h){
			int noww=1;
			for(int k=j;k<j+h/2;++k){
				int u=f[k],v=1ll*noww*f[k+h/2]%MOD;
				f[k]=(u+v)%MOD,f[k+h/2]=((u-v)%MOD+MOD)%MOD;
				noww=1ll*noww*wh%MOD;
			}
		}
	}
	return;
}
int main(void){
	init();
	scanf("%s",A),scanf("%s",B);
	int n=strlen(A),m=strlen(B),maxlen=1;
	while(maxlen<=n+m-2){
		maxlen<<=1;
	}
	F.resize(maxlen+5),G.resize(maxlen+5);
	for(int i=0;i<n;++i){
		F[i]=A[n-i-1]-'0';
	}
	for(int i=0;i<m;++i){
		G[i]=B[m-i-1]-'0';
	}
	NTT(F,maxlen,1),NTT(G,maxlen,1);
	for(int i=0;i<maxlen;++i){
		F[i]=1ll*F[i]*G[i]%MOD;
	}
	NTT(F,maxlen,-1);
	for(int i=0;i<maxlen;++i)
		F[i]=1ll*F[i]*qpow(maxlen,MOD-2)%MOD,ans[i]=F[i];
	for(int i=0;i<maxlen+200;++i){
		ans[i+1]+=ans[i]/10,ans[i]%=10;
	}
	for(int i=maxlen+200;i>=0;--i){
		if(ans[i]){
			for(int j=i;j>=0;--j)
				printf("%d",ans[j]);
			return 0;
		}
	}
	return 0;
}