[BZOJ3992] [SDOI2015] [NTT] 序列统计

最新推荐文章于 2019-03-15 15:44:51 发布

原创最新推荐文章于 2019-03-15 15:44:51 发布 · 518 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#NTT

FFT & FWT 专栏收录该内容

22 篇文章

订阅专栏

本文探讨了一个特定模数下组合计数的问题，通过将乘法转化为加法，利用NTT（Number Theoretic Transform）和快速幂算法进行优化解决。文章提供了完整的C++实现代码。

先写出DP方程
$f(i,j)$ 表示放了 $i$ 个数，乘积模 $m$ 为 $j$ 的方案数
那么 $f(i,j)=\sum _{ab=j,a\in S,b\in S} f(i-1,a)*f(f-1,b)$
这个可以用矩阵+快速幂优化但是不够优
发现模数比较特别，如果可以把乘号变成加号，就可以用NTT+快速幂来做
用 $m$ 的原根 $g$ 把每个数表示成 $g^i$ 就可以啦

#include <cstdio>
#include <iostream>
#include <algorithm>

using namespace std;

const int N=8010<<3,P=1004535809,G=3;

inline char nc(){
    static char buf[100000],*p1=buf,*p2=buf;
    return p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++;
}

inline void rea(int &x){
    char c=nc(); x=0;
    for(;c>'9'||c<'0';c=nc());for(;c>='0'&&c<='9';x=x*10+c-'0',c=nc());
}

int n,m,x,s,g;
int a[N],numb[N];

inline bool prim(int x){
    static int vis[N];
    for(int i=0;i<m;i++) vis[i]=0;
    for(int i=1,prod=x;i<m;i++,prod=1LL*prod*x%m)
        if(vis[prod]) return false;
        else vis[prod]=1;
    return g=x;
}

int _x[N],rev[N],w[2][N];
int num,L,M;

inline int Pow(int x,int y){
    int ret=1;
    for(;y;y>>=1,x=1LL*x*x%P)
        if(y&1) ret=1LL*ret*x%P;
    return ret;
}

inline void Pre(int n){
    num=n;
    int g=Pow(G,(P-1)/num),invg=Pow(g,P-2);
    w[0][0]=w[1][0]=1;
    for(int i=1;i<num;i++)
        w[0][i]=1LL*w[0][i-1]*invg%P,w[1][i]=1LL*w[1][i-1]*g%P;
}

inline void NTT(int *a,int n,int r){
    for(int i=1;i<n;i++) if(rev[i]>i) swap(a[rev[i]],a[i]);
    for(int i=1;i<n;i<<=1)
        for(int j=0;j<n;j+=(i<<1))
            for(int k=0;k<i;k++){
                int x=a[j+k],y=1LL*a[j+k+i]*w[r][num/(i<<1)*k]%P;
                a[j+k]=(x+y)%P; a[j+k+i]=(x+P-y)%P;
            }
    if(!r) for(int i=0,inv=Pow(n,P-2);i<n;i++) a[i]=1LL*a[i]*inv%P;
}

inline void mul(int *a,int *b){
    static int tmp[N];
    for(int i=0;i<m-1;i++) tmp[i]=b[i];
    NTT(a,M<<1,1); NTT(tmp,M<<1,1); 
    for(int i=0;i<(M<<1);i++) a[i]=1LL*a[i]*tmp[i]%P;
    NTT(a,M<<1,0);
    for(int i=0;i<m-1;i++) (a[i]+=a[i+m-1])%=P;
    for(int i=m-1;i<(M<<1);i++) a[i]=tmp[i]=0;
}

inline int *Pow(int *x,int y){
    int *ret=_x; ret[0]=1;
    for(;y;y>>=1,mul(x,x))
        if(y&1) mul(ret,x);
    return ret;
}

int main(){
    rea(n); rea(m); rea(x); rea(s);
    M=1; for(;M<=(m<<1);M<<=1); Pre(M<<1);
    L=0; while(!(M>>L&1)) L++;
    for(int i=1;i<(M<<1);i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L);
    for(int i=2;i<m;i++)
        if(prim(i)) break;
    for(int i=1,prod=g;i<m-1;i++,prod=1LL*prod*g%m) numb[prod]=i;
    for(int i=1;i<=s;i++){
        int x; rea(x); 
        if(!x) continue;
        a[numb[x%m]]++;
    }
    printf("%d\n",Pow(a,n)[numb[x]]);
    return 0;
}