BZOJ3992：[SDOI2015]序列统计 （NTT+倍增+DP）

题目传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=3992

---

题目分析：“如果一道题的模数是998开头的那个费马素数，那它多半不是NTT。如果模数是1004开头的，那它基本上就是NTT了。”——tututu

我们可以将长度为n的这个序列分成两段，分别算出前半段和后半段乘积模m等于0~m-1的方案数，记进f和g数组。然后令$f[i]*g[j]$向$ans[i*j\mod m]$贡献答案。

然而这样是$O(m^2\log(n))$的，明显超时。这个式子很像卷积，能不能把乘法换成加法呢？由于m为质数，它必定存在一个原根g。这样1~m-1就可以换成原根的0~m-2次方，我们就可以将模M意义下的乘法换成原根指数模m-1意义下的加法。至于原根的求法，只要找到一个$x(0<x<M)$，使得对于所有整除M-1的$y(0<y<M-1)$，$x^y$都不为1，x即为M的原根。一般质数的原根都不会很大，暴力检验即可。最后用NTT优化卷积，时间复杂度为$O(m\log(m)\log(n))$。

其实这题太水了，就是我用来练习一下求原根和NTT的。

---

CODE：

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
using namespace std;

const int maxn=300000;
const long long M=1004535809;
const long long g=3;
typedef long long LL;

LL A[maxn];
LL B[maxn];

int Rev[maxn];
int N,Lg;

int p[maxn];
int a[maxn];
int n,m,val,num,gm;

LL Pow(LL x,LL y,LL Mod)
{
    if (!y) return 1;
    LL temp=Pow(x,y>>1,Mod);
    temp=temp*temp%Mod;
    if (y&1) temp=temp*x%Mod;
    return temp;
}

int Find()
{
    int x=m-1,y=(int)floor( sqrt( (double)x )+0.5 );
    for (int i=2; i<m; i++)
    {
        bool sol=true;
        for (int j=2; j<=y; j++)
            if (x%j==0)
            {
                if (Pow(i,j,m)==1) sol=false;
                if (Pow(i,x/j,m)==1) sol=false;
                if (!sol) break;
            }
        if (sol) return i;
    }
}

void DFT(LL *a,int f)
{
    for (int i=0; i<N; i++)
        if (i<Rev[i]) swap(a[i],a[ Rev[i] ]);

    for (int len=2; len<=N; len<<=1)
    {
        int mid=(len>>1);
        LL e=Pow(g,(M-1)/len,M);
        if (f==-1) e=Pow(e,M-2,M);

        for (LL *p=a; p!=a+N; p+=len)
        {
            LL wn=1;
            for (int i=0; i<mid; i++)
            {
                LL temp=wn*p[mid+i]%M;
                p[mid+i]=(p[i]-temp+M)%M;
                p[i]=(p[i]+temp)%M;
                wn=wn*e%M;
            }
        }
    }
}

void NTT()
{
    DFT(A,1);
    DFT(B,1);
    for (int i=0; i<N; i++) A[i]=A[i]*B[i]%M;
    DFT(A,-1);
    LL inv=Pow(N,M-2,M);
    for (int i=0; i<N; i++) A[i]=A[i]*inv%M;
    for (int i=m-1; i<N; i++) A[i%(m-1)]+=A[i],A[i%(m-1)]%=M,A[i]=0;
}

void Work(LL x)
{
    if (x==1)
    {
        for (int i=0; i<N; i++) A[i]=a[i];
        return;
    }

    Work(x>>1);
    for (int i=0; i<N; i++) B[i]=A[i];
    NTT();

    if (x&1)
    {
        for (int i=0; i<N; i++) B[i]=a[i];
        NTT();
    }
}

int main()
{
    freopen("3992.in","r",stdin);
    freopen("3992.out","w",stdout);

    scanf("%d%d%d%d",&n,&m,&val,&num);

    gm=Find();
    p[1]=0;
    int v=1;
    for (int i=1; i<m-1; i++) v=v*gm%m,p[v]=i;

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

    N=1,Lg=0;
    while (N<2*m) N<<=1,Lg++;
    for (int i=0; i<N; i++)
        for (int j=0; j<Lg; j++)
            if (i&(1<<j)) Rev[i]|=(1<<(Lg-j-1));

    Work(n);
    printf("%lld\n",A[ p[val] ]);

    return 0;
}