快速数论变换 NTT

前言：
在学习NTT之前，应当先掌握FFT（快速傅立叶变换）的基本知识，并能动手完成代码实现。如果有时间（心情）我会写一篇FFT的算法介绍。

在FFT中起到相当重要的作用的就是那个主n次单位根\(w_n=e^{\frac{2i\pi}{n}}\)，一切的一切都围绕这个神奇的复数展开。但是复数的运算设计到两个double变量，浮点数的计算不仅很慢，而且经常出现精度误差。

为此我们找到了NTT作为代替，其中NTT本身是在模一个费马素数\(P\)的意义下展开的，我们设\(G\)为\(P\)的原根，会发现\(G^{\frac{P-1}{n}}\)具有和\(w_n=e^{\frac{2i\pi}{n}}\)相似的性质，可以代替\(w_n\)作为单位根。这个就不写详细证明了，反正也没什么要紧，知道这个性质的存在就好。

所以，我们现在用\(G^{\frac{P-1}{n}}\)代替原来的\(w_n\)，继续做原来的分治过程。原本的除法全部改成乘以逆元即可。

在51Nod上1027和1028都可以拿来练习模板，好想Uoj的#37也是高精度乘法的样子。

下面给出我的NTT代码实现。

#include <cstdio>
#include <cstring>

template <class T>
inline void swap(T &a, T &b)
{
    T c;
    c = a;
    a = b;
    b = c;
}

const int siz = 500005;

const int P = 998244353, G = 3;

inline int pow(int a, int b)
{
    int r = 1;
    
    while (b)
    {
        if (b & 1)
            r = 1LL * r * a % P;
        
        b >>= 1, a = 1LL * a * a % P;
    }
    
    return r;
}

inline void calculateNTT(int *s, int n, int f)
{
    {
        int cnt = 0;
        
        static int rev[siz];
        
        while (n >> cnt)++cnt; --cnt;
        
        memset(rev, 0, sizeof rev);
        
        for (int i = 0; i < n; ++i)
        {
            rev[i] |= rev[i >> 1] >> 1;
            rev[i] |= (i & 1) << (cnt - 1);
        }
        
        for (int i = 0; i < n; ++i)if (i < rev[i])swap(s[i], s[rev[i]]);
    }
    
    {
        for (int i = 1; i < n; i <<= 1)
        {
            int wn = pow(G, (P - 1) / (i * 2));
            
            if (f == -1)wn = pow(wn, P - 2);
            
            for (int j = 0; j < n; j += (i << 1))
            {
                int wk = 1;
                
                for (int k = 0; k < i; ++k, wk = 1LL * wk * wn % P)
                {
                    int x = s[j + k];
                    int y = 1LL * s[i + j + k] * wk % P;
                    
                    s[j + k] = x + y;
                    s[i + j + k] = x - y;
                    
                    s[j + k] = (s[j + k] % P + P) % P;
                    s[i + j + k] = (s[i + j + k] % P + P) % P;
                }
            }
        }
    }
    
    {
        if (f == -1)
        {
            int inv = pow(n, P - 2);
            
            for (int i = 0; i < n; ++i)
                s[i] = 1LL * s[i] * inv % P;
        }
    }
}

signed main(void)
{
    static char sa[siz]; 
    static char sb[siz];
    
    scanf("%s", sa);
    scanf("%s", sb);
    
    static int la, a[siz];
    static int lb, b[siz];
    
    la = strlen(sa);
    lb = strlen(sb);
    
    for (int i = 0; i < la; ++i)a[i] = sa[la - i - 1] - '0';
    for (int i = 0; i < lb; ++i)b[i] = sb[lb - i - 1] - '0';
    
    int len; for (len = 1; len < la || len < lb; len <<= 1);
    
    calculateNTT(a, len << 1, +1);
    calculateNTT(b, len << 1, +1);
    
    
    for (int i = 0; i < len << 1; ++i)a[i] = 1LL * a[i] * b[i] % P;
    
    calculateNTT(a, len << 1, -1); 
    
    for (int i = 0; i < len << 1; ++i)a[i + 1] += a[i] / 10, a[i] = a[i] % 10;
    
    len <<= 1; while (!a[len])--len;
    
    for (int i = len; ~i; --i)printf("%d", a[i]); puts("");
}

@Author: YouSiki

转载于:https://www.cnblogs.com/yousiki/p/6426556.html