【算法】快速数论变换(NTT)初探

【简介】

　　快速傅里叶变换(FFT)运用了单位复根的性质减少了运算，但是每个复数系数的实部和虚部是一个余弦和正弦函数，因此系数都是浮点数，而浮点数的运算速度较慢且可能产生误差等精度问题，因此提出了以数论为基础的具有循环卷积性质的快速数论变换(NTT)。

　　在FFT中，通过$n$次单位复根即$\omega^n=1$的$\omega$来运算，而对于NTT来说，则是运用了素数的原根来运算。

【原根】

【定义】

　　对于两个正整数$a,m$满足$gcd(a, m)=1$，由欧拉定理可知，存在正整数$d\leq m-1$，如$d=\varphi(m)$，使得$a^d\equiv 1(mod\ m)$。

　　因此，在$gcd(a, m)=1$时，定义$a$对模$m$的指数$\delta m(a)$为使$a^d\equiv 1(mod\ m)$成立的最小正整数$d$。若$\delta m(a)=\varphi (m)$，则称$a$是模$m$的原根。

【性质/定义2】

　　若一个数$g$是对于$P$的原根，那么$g^i\ mod\ P,1\leq i<P$的结果互不相同。

【求原根方法】

　　对质数$P-1$分解质因数得到不同的质因子$p_1,p_2,p_3,...,p_n$，对于任何$2\leq a\leq P-1$，判定$a$是否为$P$的原根，只需要检验$a^{\frac{P-1}{p_1}},a^{\frac{P-1}{p_2}},...,a^{\frac{P-1}{p_n}}$这n个数中，是否存在一个数$mod\ P$为$1$，若存在，则$a$不是$P$的原根，否则$a$是$P$的原根。

【正确性证明】

　　假设存在一个$t<\varphi(P)=P-1$使得$a^t\equiv 1(mod\ P)$且$\forall i\subseteq [1,P),a^{\frac{P-1}{p_i}}\ mod\ P \neq 1$。

　　由裴蜀定理得，一定存在一组$k,x$使得$kt=(P-1)x+gcd(t, P-1)$

　　由欧拉定理/费马小定理得，$a^{P-1}\equiv 1(mod\ P)$

　　于是$1\equiv a^{kt}\equiv a^{(P-1)x+gcd(t, P-1)}\equiv a^{gcd(t, P-1)}(mod\ P)$

　　$t<P-1$故$gcd(t, P-1)<P-1$

　　又$gcd(t, P-1)|P-1$，于是$gcd(t, P-1)$必整除$a^{\frac{P-1}{p_1}},a^{\frac{P-1}{p_2}},...,a^{\frac{P-1}{p_n}}$中至少一个，设$gcd(t, P-1)|a^{\frac{P-1}{p_i}}$，则$a^{\frac{P-1}{p_i}}\equiv a^{gcd(t, P-1)}\equiv 1(mod\ P)$。

　　故假设不成立。

【用途】

　　我们可以发现原根$g$拥有所有FFT所需的$\omega$的性质，于是如果我们用$g^{\frac{P-1}{N}}(mod\ P)$来代替$\omega_n=e^{\frac{2\pi i}{N}}$，就能把复数对应成一个整数，在$(mod\ P)$意义下做快速变换了。

【NTT模数】

　　显然在上述的用途中，$P$必须是素数且$N$必须是$P-1$的因数，因为$N$是$2$的幂，所以可以构造形如$P=c\cdot 2^k+1$的素数。

　　常见的形如$P=c\cdot 2^k+1$的素数有$998244353=119\cdot 2^{23}+1,1004535809=479\cdot 2^{21}+1$，它们的原根都为$3$

　　如果题目的模数$P$任意怎么办？我们取的模数必须超过$n(P-1)^2$

　　那么我们可以取多个模数（乘积$> n(P-1)^2$）做完NTT之后用CRT合并...

$1004535809 = 479 \cdot 2^{21} + 1$

#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<algorithm>
#define ll long long
#define MOD(x) ((x)>=mod?(x)-mod:(x))
using namespace std;
const int maxn=4000010, inf=1e9, mod=998244353;
int n, N;
int a[maxn], b[maxn], c[maxn];
char s[maxn];
inline int power(int a, int b)
{
    int ans=1;
    for(;b;b>>=1, a=1ll*a*a%mod)
    if(b&1) ans=1ll*ans*a%mod;
    return ans;
}
inline void ntt(int *a, int f)
{
    for(int i=0, j=0;i<N;i++)
    {
        if(i<j) swap(a[i], a[j]);
        for(int k=N>>1;(j^=k)<k;k>>=1);
    }
    for(int i=2;i<=N;i<<=1)
    {
        int nw=power(3, (mod-1)/i);
        if(f==-1) nw=power(nw, mod-2);
        for(int j=0, m=i>>1;j<N;j+=i)
        for(int k=0, w=1;k<m;k++)
        {
            int t=1ll*a[j+k+m]*w%mod;
            a[j+k+m]=MOD(a[j+k]-t+mod);
            a[j+k]=MOD(a[j+k]+t);
            w=1ll*w*nw%mod;
        }
    }
    if(f==-1) for(int i=0, inv=power(N, mod-2);i<N;i++) a[i]=1ll*a[i]*inv%mod;
}
int main()
{
    scanf("%d", &n); for(N=1;N<(n<<1);N<<=1); 
    scanf("%s", s); for(int i=0;i<n;i++) a[n-i-1]=s[i]-'0';
    scanf("%s", s); for(int i=0;i<n;i++) b[n-i-1]=s[i]-'0';
    ntt(a, 1); ntt(b, 1);
    for(int i=0;i<N;i++) c[i]=1ll*a[i]*b[i]%mod;
    ntt(c, -1); 
    for(int i=0;i<N;i++)
    if(c[i]>=10)
    {
        c[i+1]+=c[i]/10; c[i]%=10;
        if(i==N-1) N++;
    }
    N--; while(!c[N] && N>0) N--;
    for(int i=N;~i;i--) printf("%d", c[i]);
}

View Code