[SDOI2015] 序列统计

Description

给你一个模 \(m\) 意义下的数集，需要用这个数集生成一个数列，使得这个数列在的乘积为 \(x\)。
问方案数模 \(1004535809\)。 保证 \(m\) 是质数。

Solution

首先要知道 \(1004535809\) 是个 \(NTT\) 常用模数，原根为 \(3\)。

然后推一推 \(DP\) 式子，设 \(f[i][j]\) 表示填好了前 \(i\) 个数，乘积为 \(j\) 的方案数。

转移显然 \(f[i][j]=\sum\limits_{p\times k=j}f[i-1][p]\times sum[k]\) 。其中 \(sum[i]\) 表示数集中 \(i\) 是否出现。

然鹅这样做是 \(O(m^2n)\) 的

观察到如果把转移写成矩阵，那每一次的转移矩阵都是相等的，这启发我们进行快速幂，复杂度变为 \(O(m^2\log n)\)。

发现这个 \(p\times k=j\) 不好处理，如果把它变成 \(p+k=j\) 的形式就可以 \(NTT\) 优化了。

\(m\) 是质数，它的原根数量为 \(\phi(m-1)\) 个。

\(m\ge 3\)。

相当于保证 \(m\) 一定有原根。

根据原根的性质，\(g^0,g^1,g^2\dots g^{m-2}\) 两两互不相同。

那 \(DP\) 式子可以改写为 \[f[i][g^a]=\sum_{g^b\times g^c=g^a} f[i-1][g^b]\times sum[g^c]\]

\[f[i][g^a]=\sum\limits_{g^{b+c}=g^a} f[i-1][g^b]\times sum[g^c]\]

\[f[i][a]=\sum\limits_{b+c=a} f[i-1][b]\times sum[c]\]

快速幂的时候 \(NTT\) 优化即可。

哦对了，因为这个式子只有 \(g^0\sim g^{m-2}\) 总共 \(m-1\) 项，所以下标实际上都应该对 \(m-1\) 取模，因为 \(0\) 是没有对应的项的。

时间复杂度 \(O(m\log n\log m)\)

求原根的时候对 \(m-1\) 因数分解，然后暴力判是否可行即可。

Code

#include<bits/stdc++.h>
using std::min;
using std::max;
using std::swap;
using std::vector;
typedef double db;
typedef long long ll;
#define pb(A) push_back(A)
#define pii std::pair<int,int>
#define all(A) A.begin(),A.end()
#define mp(A,B) std::make_pair(A,B)
#define int long long
const int G=3;
const int N=8005*8;
const int mod=1004535809;

int lim,rev[N];
int n,m,X,g;
int p[N],pc;
int pos[N],in;
int a[N],b[N],c[N];

int ksm(int a,int b,int mod,int ans=1){
    while(b){
        if(b&1) ans=ans*a%mod;
        a=a*a%mod;b>>=1;
    } return ans;
}

int getg(){
    int Phi=m-1,sq=sqrt(Phi);
    for(int i=2;i<=sq;i++){
        if(Phi%i==0) p[++pc]=i,i*i!=Phi?p[++pc]=Phi/i:0;
    } for(int i=2;i<=Phi;i++){
        int flag=0;
        for(int j=1;j<=pc;j++){
            if(ksm(i,p[j],m)==1) {
                flag=1;
                break;
            }
        } if(!flag) return i;
    }
}

int inv(int x){
    return ksm(x,mod-2,mod);
}

int getint(){
    int X=0,w=0;char ch=getchar();
    while(!isdigit(ch))w|=ch=='-',ch=getchar();
    while( isdigit(ch))X=X*10+ch-48,ch=getchar();
    if(w) return -X;return X;
}

void ntt(int *f,int opt){
    for(int i=0;i<lim;i++) if(i<rev[i]) swap(f[i],f[rev[i]]);
    for(int mid=1;mid<lim;mid<<=1){
        int tmp=ksm(G,(mod-1)/(mid<<1),mod);
        if(opt==-1) tmp=inv(tmp);
        for(int R=mid<<1,j=0;j<lim;j+=R){
            int w=1;
            for(int k=0;k<mid;k++,w=w*tmp%mod){
                int x=f[j+k],y=f[j+k+mid]*w%mod;
                f[j+k]=(x+y)%mod;f[j+k+mid]=(mod+x-y)%mod;
            }
        }
    }
}

void mul(int *a,int *b){
    for(int i=0;i<m;i++) c[i]=b[i];
    ntt(a,1),ntt(c,1);
    for(int i=0;i<lim;i++) a[i]=a[i]*c[i]%mod;
    ntt(a,-1);
    for(int i=0;i<lim;i++) a[i]=a[i]*in%mod;
    for(int i=0;i<m;i++) a[i]=(a[i]+a[i+m])%mod;
    for(int i=m;i<lim;i++) a[i]=c[i]=0;
}

void sqr(int *a){
    ntt(a,1);
    for(int i=0;i<lim;i++) a[i]=a[i]*a[i]%mod;
    ntt(a,-1);
    for(int i=0;i<lim;i++) a[i]=a[i]*in%mod;
    for(int i=0;i<m;i++) a[i]=(a[i]+a[i+m])%mod;
    for(int i=m;i<lim;i++) a[i]=0;
}

void ksm(int *ans,int *a,int b){
    while(b){
        if(b&1) mul(ans,a);
        sqr(a);b>>=1;
    }
}

signed main(){
    n=getint(),m=getint(),X=getint();
    g=getg();int now=1;
    pos[1]=0;for(int i=1;i<m-1;i++) pos[(now*=g)%=m]=i;
    lim=1;while(lim<=m) lim<<=1;lim<<=1;in=inv(lim);
    for(int i=0;i<lim;i++) rev[i]=(rev[i>>1]>>1)|(i&1?lim>>1:0);
    int len=getint();
    for(int i=1;i<=len;i++){
        int x=getint();
        if(x) b[pos[x]]++;
    } a[0]=1;m--;
    ksm(a,b,n);
    printf("%lld\n",a[pos[X]]);
    return 0;
}

转载于:https://www.cnblogs.com/YoungNeal/p/10162307.html