[bzoj3992][SDOI2015]序列统计

题目描述

小C有一个集合S，里面的元素都是小于M的非负整数。他用程序编写了一个数列生成器，可以生成一个长度为N的数列，数列中的每个数都属于集合S。
小C用这个生成器生成了许多这样的数列。但是小C有一个问题需要你的帮助：给定整数x，求所有可以生成出的，且满足数列中所有数的乘积mod M的值等于x的不同的数列的有多少个。小C认为，两个数列{Ai}和{Bi}不同，当且仅当至少存在一个整数i，满足Ai≠Bi。另外，小C认为这个问题的答案可能很大，因此他只需要你帮助他求出答案mod 1004535809的值就可以了。
1<=N<=109，3<=M<=8000，M为质数，1<=x<=M-1，输入数据保证集合S中元素不重复。

DP

我们设f[i,j]表示放到第i个数，此时乘积为j。
$F[i,j]=\sum F[i-1,k]$
满足k可以与集合内某数相乘模mo的结果为j。
N太大，但是我们显然可以用矩阵乘法优化。
我们把它化为另一种形式：
$F[i+1,j*k\% mo]=\sum F[i,j]*C[k]$
其中C[k]=0或1表示k这个数在集合里有没有出现。
这样做有什么用呢？

往FFT方面进军

M是质数，我们就可以找到它的原根r。
根据原根定义，我们可以把每一个m以内的数进行映射对应，我们设
ind[i]=j表示$r^j\equiv i(Mod \ m)$
那对于任何多项式A我们令A[i]=A[ind[i]]
DP式变为
$F[i+1,(j+k)\% (m-1)]=\sum F[i,j]*C[k]$
“标准”卷积形式。
因此我们可以对多项式C进行快速幂，多项式乘法用FFT加速。
注意FFT最后大于m-2的项要加到其下标对于m-1取模的下标里，详见代码。

参考代码

#include<cstdio>
#include<algorithm>
#include<cmath>
#define fo(i,a,b) for(i=a;i<=b;i++)
#define fd(i,a,b) for(i=a;i>=b;i--)
using namespace std;
typedef long long ll;
const ll maxm=8000+10,mo=1004535809,GG=3;
ll ind[maxm],a[maxm*3],b[maxm*3],c[maxm*3],e[maxm*3],f[maxm*3],s[maxm*3],w[maxm*3],tt[maxm*3];
double ce;
ll i,j,k,l,r,t,n,m,q,top,len,ni;
ll quicksortmi(ll x,ll y,ll z){
    if (!y) return 1;
    ll t=quicksortmi(x,y/2,z);
    t=t*t%z;
    if (y%2) t=t*(x%z)%z;
    return t;
}
void DFT(ll *a,ll sig){
    fo(i,0,len-1){
        ll p=0;
        for(ll j=0,tp=i;j<ce;j++,tp/=2) p=(p<<1)+(tp%2);
        tt[p]=a[i];
    }
    for(ll m=2;m<=len;m*=2){
        ll half=m/2,bei=len/m;
        fo(i,0,half-1){
            ll wi=(sig>0)?w[i*bei]:w[len-i*bei];
            for(ll j=i;j<len;j+=m){
                ll u=tt[j],v=tt[j+half]*wi%mo;
                tt[j]=(u+v)%mo;
                tt[j+half]=((u-v)%mo+mo)%mo;
            }
        }
    }
    if (sig==-1)
        fo(i,0,len-1) tt[i]=tt[i]*ni%mo;
    fo(i,0,len-1) a[i]=tt[i];
}
void FFT(ll *a,ll *b,ll *c){
    ll i;
    fo(i,0,len-1) e[i]=a[i],f[i]=b[i];
    DFT(e,1);DFT(f,1);
    fo(i,0,len-1) e[i]=e[i]*f[i]%mo,c[i]=0;
    DFT(e,-1);
    fo(i,0,len-1) c[i%(m-1)]=(c[i%(m-1)]+e[i])%mo;
    //fo(i,m-1,len-1) c[i]=0;
}
void solve(ll n){
    ll i;
    if (n==1){
        fo(i,0,len-1) b[i]=c[i];
        return;
    }
    solve(n/2);
    FFT(b,b,b);
    if (n%2) FFT(b,c,b);
}
int main(){
    scanf("%lld%lld%lld%lld",&n,&m,&q,&top);
    fo(i,1,top){
        scanf("%lld",&j);
        s[i]=j;
    }
    len=1;
    while (len<m*2) len*=2;
    ni=quicksortmi(len,mo-2,mo);
    ce=double(log(len)/log(2));
    w[0]=1;w[1]=quicksortmi(GG,(mo-1)/len,mo);
    fo(i,2,len) w[i]=w[i-1]*w[1]%mo;
    fo(i,1,m-1){
        fo(j,1,m-1){
            k=quicksortmi(i,j,m);
            if (k==1) break;
        }
        if (j==m-1){
            r=i;
            break;
        }
    }
    fo(j,0,m-2) ind[quicksortmi(i,j,m)]=j;
    a[ind[1]]=1;
    fo(i,1,top) 
        if (s[i]) c[ind[s[i]]]=1;
    solve(n);
    printf("%lld\n",b[ind[q]]);
}