洛谷P3321 [SDOI2015]序列统计(NTT)

传送门

题意:$a_iin S$,求$prod_{i=1}^na_iequiv xpmod{m}$的方案数

这题目太珂怕了……数学渣渣有点害怕……kelin大佬TQL

设$f[i][j]$表示$prod_{k=1}^ia_kequiv jpmod{m}$的方案数

那么$$f[2*i][j]=sum_{abequiv jpmod{m}}f[i][a]f[i][b]$$

然后因为$m$是质数。质数有一个叫做原根的东西,质数$p$的原根$g$满足$g^i mod p$在$i$为$[1,p-1]$时互不相同,然后根据费马小定理$g^{p-1}equiv 1pmod{p}$,所以$i$在$[0,p-2]$的范围内$g^{i}$能取遍[1,p-1]的所有数(因为原根一般都特别小,所以可以直接暴力求,这个可以看代码)

因为$g^Aequiv apmod{m}$,每一个$A$和$a$互相对应,那么我们可以用$A$来代替$a$

则上式变为$$f[2*i][K]equiv sum_{g^Ag^Bequiv g^Kpmod{m}}f[i][A]f[i][B]$$

$$f[2*i][K]equiv sum_{g^{A+B}equiv g^Kpmod{m}}f[i][A]f[i][B]$$

根据费马小定理$a^bequiv a^{b mod p-1}pmod{p}$,可得$$f[2*i][K]equiv sum_{g^{(A+B) mod m-1}=g^K}f[i][A]f[i][B]$$

$$f[2*i][K]equiv sum_{A+Bequiv Kpmod{m-1}}f[i][A]f[i][B]$$

设$g[K]=sum_{A+B=K}f[i][A]f[i][B]$,则$$f[i*2][K]=g[K]+g[K+m-1]$$

然后这个每一层的转移都是一样的,所以用多项式快速幂来计算

 1 //minamoto
 2 #include<iostream>
 3 #include<cstdio>
 4 #include<algorithm>
 5 #define ll long long
 6 #define swap(x,y) (x^=y,y^=x,x^=y)
 7 #define mul(x,y) (1ll*x*y%P)
 8 #define add(x,y) (x+y>=P?x+y-P:x+y)
 9 #define dec(x,y) (x-y<0?x-y+P:x-y)
10 using namespace std;
11 #define getc() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
12 char buf[1<<21],*p1=buf,*p2=buf;
13 inline int read(){
14     #define num ch-'0'
15     char ch;bool flag=0;int res;
16     while(!isdigit(ch=getc()))
17     (ch=='-')&&(flag=true);
18     for(res=num;isdigit(ch=getc());res=res*10+num);
19     (flag)&&(res=-res);
20     #undef num
21     return res;
22 }
23 const int N=20005,P=1004535809;
24 int n,m,vis[N],F[N],G[N],A[N],B[N],O[N];
25 int l,limit=1,r[N],x,S,pr;
26 inline int ksm(int a,int b,int P){
27     int res=1;a%=P;
28     while(b){
29         if(b&1) res=mul(res,a);
30         a=mul(a,a),b>>=1;
31     }
32     return res;
33 }
34 inline int calc(int x){
35     //求原根,若只有j=x-1时i^j=1(mod x),i是x的原根 
36     if(x==2) return 1;
37     for(int i=2;;++i){
38         bool flag=1;
39         for(int j=2;j*j<x;++j)
40         if(ksm(i,(x-1)/j,x)==1){flag=false;break;}
41         if(flag) return i;
42     }
43 }
44 void NTT(int *A,int type){
45     for(int i=0;i<limit;++i)
46     if(i<r[i]) swap(A[i],A[r[i]]);
47     for(int mid=1;mid<limit;mid<<=1){
48         int R=mid<<1,Wn=ksm(3,(P-1)/R,P);O[0]=1;
49         for(int j=1;j<mid;++j) O[j]=mul(O[j-1],Wn);
50         for(int j=0;j<limit;j+=R){
51             for(int k=0;k<mid;++k){
52                 int x=A[j+k],y=mul(O[k],A[j+k+mid]);
53                 A[j+k]=add(x,y),A[j+k+mid]=dec(x,y);
54             }
55         }
56     }
57     if(type==-1){
58         reverse(A+1,A+limit);
59         for(int i=0,inv=ksm(limit,P-2,P);i<limit;++i)
60         A[i]=mul(A[i],inv);
61     }
62 }
63 void Mul(int *G,int *F){
64     for(int i=0;i<limit;++i) A[i]=G[i],B[i]=F[i];
65     NTT(A,1),NTT(B,1);
66     for(int i=0;i<limit;++i) A[i]=mul(A[i],B[i]);
67     NTT(A,-1);
68     for(int i=0;i<m-1;++i) G[i]=add(A[i],A[i+m-1]);
69     
70 }
71 void solve(int b){
72     G[0]=1;
73     while(b){
74         if(b&1) Mul(G,F);
75         b>>=1,Mul(F,F);
76     }
77 }
78 int main(){
79 //    freopen("testdata.in","r",stdin);
80     n=read(),m=read(),x=read(),S=read();
81     pr=calc(m);
82     for(int i=1,v;i<=S;++i)
83     vis[v=read()]=1;
84     int pos=-1;
85     for(int i=0,j=1;i<m-1;++i,j=j*pr%m){
86         if(vis[j]) F[i]=1;
87         if(j==x) pos=i;
88     }
89     if(pos==-1) return puts("0"),0;
90     while(limit<=((m-1)<<1)) limit<<=1,++l;
91     for(int i=0;i<limit;++i)
92     r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
93     solve(n);
94     printf("%d
",G[pos]%P);
95     return 0;
96 }
原文地址:https://www.cnblogs.com/bztMinamoto/p/9745412.html