快速数论变换(NTT)小结

时间：2021-06-09 00:22:10

NTT

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

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

这个东西，叫原根

原根

阶

若\(a,p\)互素，且\(p>1\)，

对于\(a^n \equiv 1 \pmod{p}\)最小的\(n\)，我们称之为\(a\)模\(p\)的阶,记做\(\delta_p(a)\)

例如：

\(\delta_7(2)=3\)，

\(2^1 \equiv 2 \pmod{7}\)

\(2^2 \equiv 4 \pmod{7}\)

\(2^3 \equiv 1 \pmod{7}\)

原根

原根的定义

设\(p\)是正整数，\(a\)是整数，若\(\delta_p(a)\)等于\(\phi(p)\)，则称\(a\)为模\(p\)的一个原根

\(\delta_7(3)=6=\phi(7)\)，因此\(3\)是模\(7\)的一个原根

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

如果模数\(p\)有原根，那么它一定有\(\phi(\phi(p))\)个原根

原根存在的重要条件为\(m = 2,4,p^a,2p^{a}\)，其中\(p\)为奇素数\(a \ge 1\)

例如

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

  • 若\(P\)为素数，假设一个数\(g\)是\(P\)的原根，那么\(g^i \mod P (1<g<P,0<i<P)\)的结果两两不同

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

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

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

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

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

\[\omega_n \equiv g^\frac{p-1}{n} \mod p
\]

然后把FFT中的\(\omega_n\)都替换掉就好了

\(p\)建议取\(998244353\)，它的原根为\(3\)。

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

对于质数\(p\)，质因子分解\(p−1\)，若\(g^{\frac{p-1}{p_i}} \neq 1 \pmod p\)恒成立，\(g\)为\(p\)的原根

实现

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;

}

相关文章