[原根 + NTT] LOJ#2183 BZOJ3992:「SDOI2015」序列统计

做法比较显然,应该就是这样的 DP D P :

f(i)=jki(modm)f(j)g(k) f ( i ) = ∑ j ∗ k ≡ i ( mod m ) f ′ ( j ) g ( k )

可以用原根转化为加法,就变成
f(i)=j+ki(modm)f(j)g(k) f ( i ) = ∑ j + k ≡ i ( mod m ) f ′ ( j ) g ( k )

是个循环卷积,循环卷积的求法就是先正常乘一次,然后把后面超出 m m <script type="math/tex" id="MathJax-Element-567">m</script> 的部分加到前面去,并清0。
乘一次会算了,然后直接快速幂就行了。

#include<cstdio>
#include<cctype>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long LL;
const LL P=1004535809ll,maxn=1000005,P_G=3,maxm=8005;
inline char gc(){
  static char buf[100000],*p1=buf,*p2=buf;
  return p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++;
}
inline LL getint(){
  char ch=gc(); LL res=0,ff=1;
  while(!isdigit(ch)) ch=='-'?ff=-1:0, ch=gc();
  while(isdigit(ch)) res=(res<<3)+(res<<1)+ch-'0', ch=gc();
  return res*ff;
}
LL Pow(LL a,LL b,LL MOD){
  LL res=1; a%=MOD;
  for(;b;b>>=1,a=a*a%MOD) if(b&1) res=res*a%MOD;
  return res;
}
int rev[maxn];
void get_rev(int n){
  int t=0; while((1<<t)<n) t++;
  rev[0]=0; for(int i=1;i<=n-1;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<t-1);
}
void NTT(LL a[],int n,int k){
  for(int i=0;i<=n-1;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
  for(int m=2;m<=n;m<<=1){
    int wm=Pow(P_G,(P-1)/m,P); if(k==-1) wm=Pow(wm,P-2,P);
    for(int i=0;i<=n-1;i+=m){
      LL w=1,t0,t1;
      for(int j=0;j<=(m>>1)-1;j++,w=w*wm%P)
    t0=a[i+j], t1=w*a[i+j+(m>>1)]%P, a[i+j]=(t0+t1)%P, a[i+j+(m>>1)]=((t0-t1)%P+P)%P;
    }
  }
  if(k==-1){
    int t=Pow(n,P-2,P);
    for(int i=0;i<=n-1;i++) a[i]=(LL)a[i]*t%P;
  }
}

int num,m,X,S,g,N;
int gp[maxm],ind[maxm];
void get_mG(int m){
  for(g=2;;g++){
    bool pd=true;
    for(int j=1;j<=m-2&&pd;j++) if(Pow(g,j,m)==1) pd=false;
    if(pd) return;
  }           
}
LL A[maxn],Ans[maxn];
void Calc(int b){
  Ans[0]=1;
  for(;b;b>>=1){
    NTT(A,N,1);
    if(b&1){
      NTT(Ans,N,1); 
      for(int i=0;i<=N-1;i++) Ans[i]=Ans[i]*A[i]%P;
      NTT(Ans,N,-1);
      for(int i=N-1;i>=m-1;i--) Ans[i-m+1]=(Ans[i-m+1]+Ans[i])%P, Ans[i]=0;
    }
    for(int i=0;i<=N-1;i++) A[i]=A[i]*A[i]%P;
    NTT(A,N,-1);
    for(int i=N-1;i>=m-1;i--) A[i-m+1]=(A[i-m+1]+A[i])%P, A[i]=0;
  }
}
int main(){ 
  freopen("loj2183.in","r",stdin);
  freopen("loj2183.out","w",stdout);
  scanf("%d%d%d%d",&num,&m,&X,&S); 
  get_mG(m);
  gp[0]=1; for(int i=1;i<=m-2;i++) gp[i]=gp[i-1]*g%m;
  for(int i=0;i<=m-2;i++) ind[gp[i]]=i;
  while(S--){
    int x; scanf("%d",&x);
    if(x) A[ind[x]]=1;
  }
  N=1; while(N<m) N<<=1;  N<<=1; get_rev(N);
  Calc(num);
  printf("%d\n",Ans[ind[X]]);
  return 0;
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值