【BZOJ3992/SDOI2015】序列统计 NTT

原题走这里

一道非常奇妙的题，操作十分神奇。
首先原题中的式子是 $a_1a_2a_3...a_{|S|}\equiv p (mod m) ， a_i\in S$
然而乘法是没有办法直接做FFT的，所以我们要把式子中的连乘改为连加
也就是说，给上面式子“取对数”

求出m的原根g，将上述式子中的$a_i$和p还有原题描述中的x全部表示成g的指数
离散对数一波（其实数据很小暴力就可以了）。
上述式子就可以转化成：$g^{A_1}...g^{A_n}\equiv g^P (mod m)$ 此外 $x=g^X$
根据费马小定理：$A_1+...+A_n\equiv P (mod {m-1})$
那么就相当于多项式N次幂，同时对$x^{m-1}$取模，NTT一波即可
（这也许可以算上是多项式快速幂？）

简要介绍一下NTT，说白了就是带模FFT
但是难点在于：我们要把FFT中的一切运算一切数值全部投射到模P域里面
其中最大的问题就是单位根，我们要在模P域下找到对应物，也就是原根
原根就是这样的g使得g的0到P-2次幂各不相同
此外费马小定理：$g^{P-1}\equiv 1(mod P)$
这两点使得原根和单位根非常相似，而在NTT的题目中，模大多是$m·2^k+1$的形式
因此，我们可以用$g^{\frac{x(P-1)}{n}}$来代替FFT中的$e^{\frac{2\pi ix}{n}}$
其他的基本和FFT一样，保持运算取模就是了。
值得一提的是：本题中的模数1004535809的原根是3

此外还有一个小问题
本题的输入文件结尾似乎是没有换行符的，至少在洛谷上是这样的，这种情况下scanf和cin还有getchar都读不到，但可以用fread

代码如下：

#include <bits/stdc++.h>
#define MOD 1004535809
#define LL long long
#define getchar() *(pp++);
using namespace std;
char buf[200010],*pp=buf; 
inline int read()
{
    int x=0;char ch=getchar();
    while(ch>'9'||ch<'0')ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-'0',ch=getchar();
    return x;
}
LL n,m,im,x,X,p,M,g,s[10010],a[20010],b[20010],A[20010],B[20010],w[20010];
LL qpow(int xx,int y)
{
    LL ans=1,temp=xx;
    while(y)
    {
        if(y&1)(ans*=temp)%=MOD;
        (temp*=temp)%=MOD;
        y>>=1;
    }
    return ans;
}
void init()
{
    M=m-1;
    for(LL i=1;i!=M;i<<=1,M=max(M,i));
    M<<=1;
    im=qpow(M,MOD-2);
    for(g=2;g<=m;g++)
    {
        bool flag=1;
        for(int i=1,t=g;i<m-1;i++,t=t*g%m)
        {
            if(t==1)
            {
                flag=0;
                break;
            }
        }
        if(flag)break;
    }
    for(int i=0;i<M;i++)
    {
        w[i]=qpow(3,(MOD-1)/M*i);
    }
}
void transform()
{
    for(int i=1,j=0,k=1;i<=p;i++)
    {
//      cout<<i<<endl;
        if(!s[i])continue;
        for(j=0,k=1;k!=s[i];k=k*g%m,j++);
        b[j]++;
    }
    for(int i=1;i!=x;(i*=g)%=m,X++);
}
void FFT(LL *a)
{
    for(int i=0,j=0;i!=M;i++)
    {
        if(i>j)swap(a[i],a[j]);
        for(int k=M>>1;(j^=k)<k;k>>=1);
    }
    for(int i=2;i<=M;i<<=1)
    {
        LL r=i>>1,l=M/i;
        for(int j=0;j<M;j+=i)
        {
            for(int k=0;k<r;k++)
            {
                LL t=1LL*a[j+r+k]*w[l*k]%MOD;
                a[j+r+k]=(a[j+k]-t+MOD)%MOD;
                (a[j+k]+=t)%=MOD;
            }
        }
    }
}
void MUL(LL *a,LL *b)
{
    memcpy(A,a,M*sizeof(LL));
    memcpy(B,b,M*sizeof(LL));
    FFT(A);
    FFT(B);
    for(int i=0;i<M;i++)(A[i]*=B[i])%=MOD;
    reverse(w+1,w+M);
    FFT(A);
    for(int i=0;i<M;i++)a[i]=(A[i]*im%MOD);
    reverse(w+1,w+M);
    for(int i=0;i<m-1;i++)
    {
        (a[i]+=a[i+m-1])%=MOD;
        a[i+m-1]=0;
    }
}
void qpow_for_polynomial()
{
    a[0]=1;
    while(n)
    {
        if(n&1)MUL(a,b);
        MUL(b,b);
        n>>=1;
    }
}
int main()
{
//  freopen("sdoi2015_sequence.in","r",stdin);
//  freopen("sdoi2015_sequence.out","w",stdout);
    fread(buf,sizeof(char),sizeof(buf),stdin);
    n=read(),m=read(),x=read(),p=read();
    init();
    for(int i=1;i<=p;i++)
    {
        s[i]=read();
    }
    transform();
    qpow_for_polynomial();
    cout<<a[X]<<endl;
    return 0;
}