【算法】快速數論變換(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合並...

　　模數表：http://blog.miskcoo.com/2014/07/fft-prime-table 

【例題】

例1：bzoj2179: FFT快速傅立葉 

　　原來的配方，熟悉的味道...（怎么比FFT慢了

#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]);
}

例2：bzoj3992: [SDOI2015]序列統計 

　　DP。

　　$f[i][j]$為前$i$個數總和$mod\ m$后為$j$的方案數，$B[i]$為集合內數$mod\ m$后為$i$的數的個數。$$f[i][j*k\%m]=\sum f[i-1][j]*B[k]$$

　　這個式子顯然是可以矩陣快速冪優化的，效率$O(m^2logn)$，還是會TLE。

　　對所有數$x$，設$g^i(mod\ m)=x$，則將$x$映射為$i$，即取$x$在模$m$意義下的離散對數，$c$為$f[L]$，$a$為$f[L-1]$。

　　原式改為：$$c[(i+j)\%(\varphi(m)=m-1)]=\sum_{k=0}^{i+j} a[i]\cdot b[i+j-k]$$

　　卷積形式，上NTT

　　注意乘出來超過$m-1$項，而次數超過$m-1$項的貢獻應加回次數減去$m-1$的項中，即上方式子的取模操作。

　　NTT次數較多預處理出$g$的所有次冪會快一些，而且因為不明原因直接取模會比快速取模快，CPU你怎么這么難伺候...

#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#define MOD(x) ((x)>=mod?(x)-mod:(x))
#include<algorithm>
#define ll long long
using namespace std;
const int maxn=20010, inf=1e9, mod=1004535809;
int n, m, x, s, N, g, X;
int a[maxn], b[maxn], A[maxn], B[maxn], ind[maxn], p[maxn], g1[maxn], g2[maxn];
inline void read(int &k)
{
    int f=1; k=0; char c=getchar();
    while(c<'0' || c>'9') c=='-'&&(f=-1), c=getchar();
    while(c<='9' && c>='0') k=k*10+c-'0', c=getchar();
    k*=f;   
} 
inline int power(int a, int b, int p)
{
    int ans=1;
    for(;b;b>>=1, a=1ll*a*a%p)
    if(b&1) ans=1ll*ans*a%p;
    return ans;
}
inline bool judge(int x)
{
    for(int i=1;i<=p[0];i++)
    if(power(x, (m-1)/p[i], m)==1) return 0;
    return 1;
}
inline int yg(int x)
{
    x--;
    for(int i=2;i*i<=x;i++)
    if(x%i==0)
    {
        p[++p[0]]=i;
        while(x%i==0) x/=i;
    }
    if(x>1) p[++p[0]]=x;
    for(int i=2;;i++) if(judge(i)) return i;
}
inline void init()
{
    int base1=power(3, (mod-1)/N, mod), base2=power(base1, mod-2, mod);
    g1[0]=g2[0]=1;
    for(int i=1;i<=N;i++) g1[i]=1ll*g1[i-1]*base1%mod, g2[i]=1ll*g2[i-1]*base2%mod;
}
inline void ntt(int *a, int *w, 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)
    for(int j=0, m=i>>1;j<N;j+=i)
    for(int k=0;k<m;k++)
    {
        int t=1ll*a[j+k+m]*w[N/i*k]%mod;
        a[j+k+m]=(a[j+k]-t+mod)%mod;
        a[j+k]=(a[j+k]+t)%mod;
    }
    if(f==-1) for(int i=0, inv=power(N, mod-2, mod);i<N;i++) a[i]=1ll*a[i]*inv%mod;
}
inline void mul(int *a, int *b, int *c)
{
    memcpy(A, a, N<<2); memcpy(B, b, N<<2);
    ntt(A, g1, 1); ntt(B, g1, 1);
    for(int i=0;i<N;i++) c[i]=1ll*A[i]*B[i]%mod;
    ntt(c, g2, -1);
    for(int i=0;i<m-1;i++) c[i]=MOD(c[i]+c[i+m-1]), c[i+m-1]=0;
}
int main()
{
    read(n); read(m); read(X); read(s); g=yg(m);
    for(int i=0, j=1;i<m-1;i++, j=1ll*j*g%m) ind[j]=i; X=ind[X];
    for(N=1;N<m;N<<=1); N<<=1; init();
    for(int i=1;i<=s;i++) read(x), x && (a[ind[x]]++); b[0]=1;
    for(;n;n>>=1, mul(a, a, a)) if(n&1) mul(a, b, b);
    printf("%d\n", b[X]);
}

例3：