欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页  >  IT编程

快速数论变换(NTT)小结

程序员文章站 2022-05-24 21:24:28
NTT 在FFT中,我们需要用到复数,复数虽然很神奇,但是它也有自己的局限性——需要用double类型计算,精度太低 那有没有什么东西能够代替复数且解决精度问题呢? 这个东西,叫 原根 原根 原根的定义 设$m$是正整数,$a$是整数,若$a$模$m$的阶等于$\phi(m)$,则称$a$为模$m$ ......

NTT

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

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

这个东西,叫原根

原根

原根的定义

\(m\)是正整数,\(a\)是整数,若\(a\)\(m\)的阶等于\(\phi(m)\),则称\(a\)为模\(m\)的一个原根

定义中用到了群论的一些知识,不过不会也没关系,不影响接下来的学习

我们定义\(P\)为素数,\(g\)\(P\)的原根

接下来不加证明的扔出一个很重要定理

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

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




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

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

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

1 . 对于所有\(\omega_n ^ t (0 \leq t \leq n - 1)\)均不相同

这一条可以由上面的定理得到

2 . \(\omega_{2n} ^ {2k} = \omega_n ^ k\)

通过代换可以得到

3 . \(\omega_n ^ { k + \frac{n}{2} } = -\omega_n ^ k\)

根据费马小定理和性质1可以得到

4 . $1 + \omega_n ^ k + (\omega_n ^ k) ^ 2 + \dots + (\omega_n ^ k) ^ {n - 1} = 0 $

由性质3和FFT中傅里叶逆变换的定理可以得到




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

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

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

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

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

可以证明满足\(g^r \equiv 1(\mod p)\)的最小的\(r\)一定是\(p-1\)的约数

对于质数\(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;
}