快速数论变换(NTT)小结

 


NTT

在FFT中,我们需要用到复数,复数虽然很神奇,但是它也有自己的局限性——需要用double类型计算,精度太低

那有没有什么东西能够代替复数且解决精度问题呢?

这个东西,叫原根

原根

a,p互素,且p>1

对于an1(modp)最小n,我们称之为ap的阶,记做δp(a)

例如:

δ7(2)=3

212(mod7)

224(mod7)

231(mod7)

原根

原根的定义

p是正整数,a是整数,若δp(a)等于ϕ(p),则称a为模p的一个原根

δ7(3)=6=ϕ(7),因此3是模7的一个原根

注意原根的个数是不唯一的

如果模数p有原根,那么它一定有ϕ(ϕ(p))个原根

原根存在的重要条件为m=2,4,pa,2pa,其中p为奇素数a1

例如

原根有一个非常重要的定理:

  • P为素数,假设一个数gP的原根,那么gimodP(1<g<P,0<i<P)的结果两两不同

不要问我为什么,因为我也不知道。。



考虑原根为什么能代替单位根进行运算,(这部分可以跳过)

原因很简单,因为它具有和单位根相同的性质

在FFT中,我们用到了单位根的四条性质,而原根也满足这四条性质

这样我们最终可以得到一个结论

ωngp1nmodp

然后把FFT中的ωn都替换掉就好了

p建议取998244353,它的原根为3

如何求任意一个质数的原根呢?

对于质数p,质因子分解p1,若gp1pi1(modp)恒成立,gp的原根

实现

NTT求卷积代码:

确实比FFT快了不少

#include<cstdio>
#define getchar() (p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 1<<21, stdin), p1 == p2) ? EOF : *p1++)
#define swap(x,y) x ^= y, y ^= x, x ^= y
#define LL long long 
const int MAXN = 3 * 1e6 + 10, P = 998244353, G = 3, Gi = 332748118; 
char buf[1<<21], *p1 = buf, *p2 = buf;
inline int read() { 
    char c = getchar(); int x = 0, f = 1;
    while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();}
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x * f;
}
int N, M, limit = 1, L, r[MAXN];
LL a[MAXN], b[MAXN];
inline LL fastpow(LL a, LL k) {
	LL base = 1;
	while(k) {
		if(k & 1) base = (base * a ) % P;
		a = (a * a) % P;
		k >>= 1;
	}
	return base % P;
}
inline void NTT(LL *A, int type) {
	for(int i = 0; i < limit; i++) 
		if(i < r[i]) swap(A[i], A[r[i]]);
	for(int mid = 1; mid < limit; mid <<= 1) {	
		LL Wn = fastpow( type == 1 ? G : Gi , (P - 1) / (mid << 1));
		for(int j = 0; j < limit; j += (mid << 1)) {
			LL w = 1;
			for(int k = 0; k < mid; k++, w = (w * Wn) % P) {
				 int x = A[j + k], y = w * A[j + k + mid] % P;
				 A[j + k] = (x + y) % P,
				 A[j + k + mid] = (x - y + P) % P;
			}
		}
	}
}
int main() {
	N = read(); M = read();
	for(int i = 0; i <= N; i++) a[i] = (read() + P) % P;
	for(int i = 0; i <= M; i++) b[i] = (read() + P) % P;
	while(limit <= N + M) limit <<= 1, L++;
	for(int i = 0; i < limit; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));	
	NTT(a, 1);NTT(b, 1);	
	for(int i = 0; i < limit; i++) a[i] = (a[i] * b[i]) % P;
	NTT(a, -1);	
	LL inv = fastpow(limit, P - 2);
	for(int i = 0; i <= N + M; i++)
		printf("%d ", (a[i] * inv) % P);
	return 0;
}
posted @   自为风月马前卒  阅读(9120)  评论(14编辑  收藏  举报
编辑推荐:
· 从 HTTP 原因短语缺失研究 HTTP/2 和 HTTP/3 的设计差异
· AI与.NET技术实操系列:向量存储与相似性搜索在 .NET 中的实现
· 基于Microsoft.Extensions.AI核心库实现RAG应用
· Linux系列:如何用heaptrack跟踪.NET程序的非托管内存泄露
· 开发者必知的日志记录最佳实践
阅读排行:
· TypeScript + Deepseek 打造卜卦网站:技术与玄学的结合
· Manus的开源复刻OpenManus初探
· AI 智能体引爆开源社区「GitHub 热点速览」
· 从HTTP原因短语缺失研究HTTP/2和HTTP/3的设计差异
· 三行代码完成国际化适配,妙~啊~

Contact with me

点击右上角即可分享
微信分享提示