BZOJ 3992: [SDOI2015]序列统计（NTT+快速幂+DP）

写在前面

最近入坑LL，带给我很多感动。如果有那么一群人陪着你共同努力，共同追求一个目标，完成一个梦想真是太好了。

一个人可以走得很快，但一群人可以走得更远。就像Google一样的大企业，往往聚集着众多有共同目标的人才，单靠一个人是不行的，是无法完成改变世界的理想的。（好像真的走得有点远了）

学过FFT单位复数根记得不太熟，下面我想详细讲一下原根和NTT。为了以后能巩固知识和温故知新。（然而我几乎不看自己的博客）

---

快速数论变换(NTT)

与快速傅里叶变换的区别：FFT采用单位复数根进行点值和插值，而NTT用原根代替。

原根的性质

设$n$是$2$的正整数次幂，令

$$g_n=g^{(p-1)/n}$$
其中 $p$为素数且$n|(p-1)$， $g$是模$p$意义下的原根。

对比单位复数根的$\omega^n_n=1,\omega_n^{n/2}=-1$，原根也有$g_n^n=g^{p-1}=g^{\phi(p)}=1$，$g_n^{n/2}=g^{(p-1)/2}=\pm\sqrt1$（当然都是在模$p$的环境下）。

而原根$g$的$0..p-2$次幂对应表示了$1..p-1$的所有数，且$g_n^n=g^{p-1}=1$，所以$g_n^{n/2}=-1$。

$\omega_{dn}^{dk}=\omega_n^k:g_{dn}^{dk}=g^{k(p-1)/n}=g_n^k$

$(\omega_n^k)^2=\omega^k_{n/2}:(g_n^k)^2=g_n^{2k}=g^{2k(p-1)/n}=g^{k(p-1)/(n/2)}=g_{n/2}^{k}$

$(\omega_n^{k+n/2})^2=\omega^k_{n/2}:(g_n^{k+n/2})^2=(-g_n^k)^2=g_{n/2}^k$

所以总之，单位复数根拥有的性质，什么折半、消去、求和，原根都有！原根可以很好地代替单位根参与DFT！

可以认为：$\omega_n= g_n$。

$\omega_n^0=g_n^0=1,\omega_n^{n/2}=g_n^{n/2}=-1,\omega_n^n=g_n^n=1$，单位复数根的旋转就等同原根的一个循环，它们都是一个循环群。而我们做DFT需要的就是这个。

NTT的实现

直接将单位复数根换成原根，整个运算都要在模$p$的环境下进行，所以一堆模。IDFT时除以$n$就改成乘逆元。

一般FFT能解决的NTT都能解决，但NTT没有精度问题，而且题目要求取模就只能NTT了。

至于$g$和$p$的选取，由于$n$都是$2$的幂，所以$p-1$要含有一个$2$的幂，而且要大于$n$。通常选形如$a·2^k+1$的素数。一般题目会给模数，有一些模数一看到就知道很可能是NTT。像$1004535809=479×2^{21}+1$，它的原根是$3$。另外还有$998244353=119×2^{23}+1$，原根是$3$，但它经常出现在非NTT的题目中-_-||

NTT的时间复杂度：$O(nlogn)$，但常数比FFT大一些。

---

题解

题面

NTT第一道题，当然是板子题，不过这题涉及到求原根和快速幂优化等套路。

首先DP模型很明显，以$n$为阶段，记录当前取模结果为$j$的方案数，枚举取哪个数$k$转移。一位一位的转移肯定不行，考虑两个DP数组代表前、后半段，合并然后长度相加。只能将$n$分解做一个快速幂，但这样依然会T。

因为将数字当下标计算答案是乘不是加，但化乘为加是FFT的套路，将数字放在指数，如果能做到就显然可以卷积了。

而根据上面的性质，原根刚好可以做到这一点，$g^{0..M-2}$刚好代替原来的$1..M-1$（一一对应），而质数必有原根，原根很小直接枚举，检验的话就是将$M-1$分解质因数，看看是否存在$g^{(M-1)/n}=1$就行了，因出现$1$代表出现循环，就有数字重复，不符合原根的性质。质因数都可以则$M-1$的因数都可以。

每次做NTT转移就使总时间优化到$O(MlogMlogN)$，注意每次NTT后超过$M-2$的那部分要手动累加到次数界内，然后清掉。

在膜了某大神的模板后终于写对NTT了。1A留念。

---

代码

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <iostream>
#define MAXN 8010

using namespace std;

const int MOD = 1004535809;
int n = 1, N, M, X, S, g, fX;
int f[MAXN<<2], muse[MAXN<<2], ans[MAXN<<2], W[MAXN];

int Pow(int x, int y, int Mod){
    int res = 1;
    while(y){
        if(y & 1)  res = 1LL * res * x % Mod;
        x = 1LL * x * x % Mod;
        y >>= 1;
    }
    return res;
}

void Expd(int x){
    for(int i = 1; i <= x; i++)  if(n < i)  n <<= 1;
    n <<= 1;
}

void Reverse(int *A){
    for(int i = 0; i < n-1; i++){
        int j = 0;
        for(int k = 1, tmp = i; k < n; k <<= 1, tmp >>= 1)
            j = ((j << 1) | (tmp & 1));
        if(j > i)  swap(A[i], A[j]);
    }
}

void NTT(int *A, int DFT){
    Reverse(A);
    for(int s = 1; (1<<s) <= n; s++){
        int m = 1 << s, wn = Pow(3, (DFT == 1) ? (MOD-1)/m : MOD-1-(MOD-1)/m, MOD);
        for(int k = 0; k < n; k += m)
            for(int j = 0, w = 1; j < (m>>1); j++){
                int u = A[k+j], t = 1LL * w * A[k+j+(m>>1)] % MOD;
                A[k+j] = (u+t) % MOD;
                A[k+j+(m>>1)] = (u-t+MOD) % MOD;
                w = 1LL * w * wn % MOD;
            }
    }
    if(DFT == -1){
        for(int i = 0, inv = Pow(n, MOD-2, MOD); i < n; i++)
            A[i] = 1LL * A[i] * inv % MOD;
    }
}

int Find(){
    int res = 2;
    while(true){
        int temp = M - 1, ok = 1;
        for(int i = 2; i * i <= temp && ok; i++){
            if(temp % i == 0){
                if(Pow(res, (M-1)/i, M) == 1)  ok = 0;
                while(temp % i == 0)  temp /= i;
            }
        }
        if(ok)  return res;
        res ++;
    }
}

void Work(int *A, int *B){

    for(int i = 0; i < n; i++)  muse[i] = B[i];

    NTT(A, 1);  NTT(muse, 1);

    for(int i = 0; i < n; i++)  A[i] = 1LL * A[i] * muse[i] % MOD;

    NTT(A, -1);

    for(int i = M-1; i < n; i++){
        A[i-M+1] = (A[i-M+1] + A[i]) % MOD;
        A[i] = 0;
    }
}

int main(){

    scanf("%d%d%d%d", &N, &M, &X, &S);

    g = Find();
    Expd(M-1);

    for(int i = 1, x; i <= S; i++){
        scanf("%d", &x);
        if(x)  W[x] = 1;
    }

    for(int i = 0; i < M-1; i++){
        int temp = Pow(g, i, M);
        if(W[temp])  f[i] = 1;
        if(temp == X)  fX = i;
    }

    ans[0] = 1;
    while(N){
        if(N & 1)  Work(ans, f);
        Work(f, f);
        N >>= 1;
    }

    printf("%d\n", ans[fX]);

    return 0;
} 

---

写在后面

蒟蒻我现在才学NTT，大佬们都投来轻蔑的目光。不说了，快去补生物作业吧……

---