在FFT中,咱们须要用到复数,复数虽然很神奇,可是它也有本身的局限性——须要用double类型计算,精度过低ui
那有没有什么东西可以代替复数且解决精度问题呢?spa
这个东西,叫原根code
若\(a,p\)互素,且\(p>1\),get
对于\(a^n \equiv 1 \pmod{p}\)最小的\(n\),咱们称之为\(a\)模\(p\)的阶,记作\(\delta_p(a)\)it
例如:io
\(\delta_7(2)=3\),ast
\(2^1 \equiv 2 \pmod{7}\)class
\(2^2 \equiv 4 \pmod{7}\)im
\(2^3 \equiv 1 \pmod{7}\)di
原根的定义
设\(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\)
例如
原根有一个很是重要的定理:
不要问我为何,由于我也不知道。。
考虑原根为何能代替单位根进行运算,(这部分能够跳过)
缘由很简单,由于它具备和单位根相同的性质
在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; }