[BZOJ3992][SDOI2015]序列统计(DP+原根+NTT)

3992: [SDOI2015]序列统计

Time Limit: 30 Sec  Memory Limit: 128 MB
Submit: 1888  Solved: 898
[Submit][Status][Discuss]

Description

小C有一个集合S，里面的元素都是小于M的非负整数。他用程序编写了一个数列生成器，可以生成一个长度为N的数

列，数列中的每个数都属于集合S。小C用这个生成器生成了许多这样的数列。但是小C有一个问题需要你的帮助：

给定整数x，求所有可以生成出的，且满足数列中所有数的乘积mod M的值等于x的不同的数列的有多少个。小C认为

，两个数列{Ai}和{Bi}不同，当且仅当至少存在一个整数i，满足Ai≠Bi。另外，小C认为这个问题的答案可能很大

，因此他只需要你帮助他求出答案mod 1004535809的值就可以了。

Input

一行，四个整数，N、M、x、|S|，其中|S|为集合S中元素个数。

第二行，|S|个整数，表示集合S中的所有元素。

1<=N<=10^9，3<=M<=8000，M为质数

0<=x<=M-1，输入数据保证集合S中元素不重复x∈[1,m-1]

集合中的数∈[0,m-1]

Output

一行，一个整数，表示你求出的种类数mod 1004535809的值。

Sample Input

4 3 1 2
1 2

Sample Output

8
【样例说明】
可以生成的满足要求的不同的数列有(1,1,1,1)、(1,1,2,2)、(1,2,1,2)、(1,2,2,1)、
(2,1,1,2)、(2,1,2,1)、(2,2,1,1)、(2,2,2,2)

HINT

Source

Round 1 感谢yts1999上传

[Submit][Status][Discuss]

好久没有调过这么痛苦了。。

首先有一个简单的DP方程:$dp[i+j][(x*y)\%p]+=dp[i][x]*dp[j][y]$，dp[i][j]表示前i个数凑成余数为j的方案数。

$n^2$转移很简单，然后可以矩阵优化，但M的范围仍然过不了。

这时候就要敢往FFT方面去想。

现在这个方程之所以不能用FFT的原因在于x,y转移到的不是x+y而是x*y%mod，我们可以想到，乘法运算在作为指数的时候就是加法，又发现M是质数，于是我们考虑用原根的次幂替代S中的每个数。

设$g^{x'} equiv x (mod p)$，$g^{y'} equiv y (mod p)$，这样方程就变为$dp[i+j][(x'+y')\% phi(p)]+=dp[i][x']*dp[j][y']$，这个就是标准的循环卷积了。

注意循环卷积次数界要再放大一倍！！还有S中要忽略0！！两个问题调到崩溃。

#include<cstdio>
#include<algorithm>
#define rep(i,l,r) for (int i=l; i<=r; i++)
typedef long long ll;
using namespace std;

const int N=20100,mod=1004535809,G=3;
int k,p,x,n,s,g,inv,cnt[N],f[N],pw[N],ind[N],rev[N],c[N];

int ksm(int a,int b,int p){
    int res;
    for (res=1; b; a=(1ll*a*a)%p,b>>=1)
        if (b & 1) res=(1ll*res*a)%p;
    return res;
}

bool chk(){
    for (int i=2; i*i<=p; i++) if ((p-1)%i==0 && ksm(g,(p-1)/i,p)==1) return 0;
    return 1;
}

void getroot(){
    if (p==2) g=1; else for (g=2; !chk(); g++);
    ind[1]=0; pw[0]=1;
    for (int i=1; i<p-1; i++) pw[i]=pw[i-1]*g%p,ind[pw[i]]=i;
}

namespace NTT{
    int n,L,rev[N];
    void init(int m){
        n=1; L=0;
        for (; n<=m; n<<=1) L++;
        n<<=1; L++; inv=ksm(n,mod-2,mod);
        for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    }
    void DFT(int a[],int n,int f){
        for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
        for (int i=1; i<n; i<<=1){
            int wn=ksm(G,(f==1) ? (mod-1)/(i<<1) : (mod-1)-(mod-1)/(i<<1),mod);
            for (int p=i<<1,j=0; j<n; j+=p){
                int w=1;
                for (int k=0; k<i; k++,w=1ll*w*wn%mod){
                    int x=a[j+k],y=1ll*w*a[i+j+k]%mod;
                    a[j+k]=(x+y)%mod; a[i+j+k]=(x-y+mod)%mod;
                }
            }
        }
        if (f==1) return;
        for (int i=0; i<n; i++) a[i]=1ll*a[i]*inv%mod;
    }
    void mul(int a[],int b[]){
        for (int i=0; i<n; i++) c[i]=b[i];
        DFT(a,n,1); DFT(c,n,1);
        for (int i=0; i<n; i++) a[i]=1ll*a[i]*c[i]%mod;
        DFT(a,n,-1);
        for (int i=n-1; i>=p-1; i--) a[i-p+1]=(a[i-p+1]+a[i])%mod,a[i]=0;
    }
};

int main(){
    freopen("bzoj3992.in","r",stdin);
    freopen("bzoj3992.out","w",stdout);
    scanf("%d%d%d%d",&k,&p,&x,&n); getroot(); NTT::init(p);
    rep(i,1,n) { scanf("%d",&s); if (s) cnt[ind[s]]++; }
    for (f[0]=1; k; k>>=1,NTT::mul(cnt,cnt)) if (k & 1) NTT::mul(f,cnt);
    printf("%d
",f[ind[x]]);
    return 0;
}