「SDOI 2015」序列统计

Problem

Description

给定一个集合 S ,其元素的值域为 \left[0,m\right) ,求满足 \displaystyle\prod_{i=1}^{n}a_{i}\equiv x\pmod{m}, \forall i : a_{i} \in S 的长为 n 的数组 a 的个数。

两个数组 a b 不同当且仅当 \exists i \in [1, n] : a_{i}\ne b_{i}

答案对 1004535809 取模。

Constraints

1\le n\le 10^{9}, 3\le m\le 8\times 10^{3}, 1\le x\le m-1

数据保证 S 中元素不重复且 m 为质数。

Solution

Analysis

由于有

\prod_{i=1}^{n}a_{i}\equiv x\pmod{m}

\sum_{i=1}^{n}\log_{g}a_{i}\equiv \log_{g}x\pmod{m-1}

其中 g m 的原根。

故可以通过对集合 S 和数 x 取离散对数将 a_{i} 的乘积转化为 \log_{g}{a_{i}} 的和,此时有 OGF:

f(x)=\sum_{s\in S}x^{\log_{g}s}

又因集合 S 中的数互不相同,有:

f(x)=\sum_{i=0}^{m-2}\left[g^{i}\bmod{m}\in S\right]x^{i}

则:

ans_{n}(x)=f^{n}(x)

使用 NTT + 快速幂即可。

Notice: 由于转化后条件是在 \mathbb{Z}_{m-1} 下成立,故 NTT 时应将 [x^{k}]f(x) 叠加到 \left[x^{k\bmod{\left(m-1\right)}}\right]f(x) 上。

时间复杂度 O\left(m\log{m}\log{n}\right)

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
#include<cstdio>
#include<algorithm>

using i64=long long;
using uchar=unsigned char;

constexpr int maxm(8000);
constexpr int mxdg(16384);
constexpr int p(1004535809);
constexpr int proot(3);

namespace IOManager{
constexpr int FILESZ(131072);
char buf[FILESZ];
const char*ibuf=buf,*tbuf=buf;

struct IOManager{
inline char gc()
{return(ibuf==tbuf)&&(tbuf=(ibuf=buf)+fread(buf,1,FILESZ,stdin),ibuf==tbuf)?EOF:*ibuf++;}

template<typename _Tp>
inline operator _Tp(){
_Tp s=0u;char c=gc();
for(;c<48;c=gc());
for(;c>47;c=gc())
s=(_Tp)(s*10u+c-48u);
return s;
}
};
}IOManager::IOManager io;

namespace math{
using Z=int;
using mZ=i64;

template<typename _Tp>
inline Z fpow(_Tp v,int n,const Z&fp=p){
Z pw=1;
for(;n;n>>=1,v=(mZ)v*v%fp)
if(n&1)pw=(mZ)pw*v%fp;
return pw;
}

inline Z calcg(const Z&fp){
static Z divs[50],*ps;ps=divs;

for(Z i=2,mp=fp-1;i*i<=mp;++i)
if(mp%i==0)
for(*++ps=(fp-1)/i,mp/=i;mp%i==0;mp/=i);
for(Z g=2;g!=fp;++g){
uchar t=1;
for(Z*s=divs+1;s<=ps;++s)
if(fpow(g,*s,fp)==1)
{t=0;break;}
if(t)
return g;
}return-1;
}

template<typename _Tp>
inline void swap(_Tp&x,_Tp&y)
{_Tp z=x;x=y;y=z;}

template<typename x_tp,typename y_tp>
inline void inc(x_tp&x,const y_tp&y)
{x+=y;(p<=x)&&(x-=p);}
template<typename x_tp,typename y_tp>
inline void dec(x_tp&x,const y_tp&y)
{x-=y;(x<0)&&(x+=p);}

template<typename x_tp,typename y_tp>
inline Z pls(x_tp&x,const y_tp&y)
{return(x+y<p?x+y:(x+y-p));}
template<typename x_tp,typename y_tp>
inline Z sub(x_tp&x,const y_tp&y)
{return(x<y?x-y+p:(x-y));}

template<typename _Tp>
inline Z sqr(const _Tp&v)
{return(mZ)v*v%p;}
}

namespace poly{
using namespace math;

using poly_t=Z[mxdg];
using poly=Z*const;

inline int calcpw2(const int&n){
int t=1;
for(;t<n;t<<=1);
return t;
}

poly_t wn,iwn;

inline void polyinit(){
const Z wnv=fpow(proot,(p-1)/mxdg),
iwnv=fpow(wnv,p-2);
wn[0]=iwn[0]=1;
for(int i=1,w=1,iw=1;i!=mxdg;++i){
w=(mZ)w*wnv%p;wn[i]=w;
iw=(mZ)iw*iwnv%p;iwn[i]=iw;
}
}

void DFT(poly&A,const int&n){
for(int i=0,j=0;i!=n;++i){
if(i<j)swap(A[i],A[j]);
for(int k=n>>1;(j^=k)<k;k>>=1);
}
Z*const&td=A+n;
for(int l=1,tp=mxdg>>1;l!=n;l+=l,tp>>=1)
for(Z*i=A,*w,z;i!=td;i+=l+l)
for(int j=(w=wn,0);j!=l;++j,w+=tp)
z=(mZ)i[j+l]**w%p,
i[j+l]=sub(i[j],z),
inc(i[j],z);
}
void IDFT(poly&A,const int&n){
for(int i=0,j=0;i!=n;++i){
if(i<j)swap(A[i],A[j]);
for(int k=n>>1;(j^=k)<k;k>>=1);
}
Z*const&td=A+n;
for(int l=1,tp=mxdg>>1;l!=n;l+=l,tp>>=1)
for(Z*i=A,*w,z;i!=td;i+=l+l)
for(int j=(w=iwn,0);j!=l;++j,w+=tp)
z=(mZ)i[j+l]**w%p,
i[j+l]=sub(i[j],z),
inc(i[j],z);

const Z invn=fpow(n,p-2);
for(Z*i=A;i!=td;++i)
*i=(mZ)*i*invn%p;
}
}

const int n=io,m=io,x=io;

inline void mul(poly::poly&f,const poly::poly&g){
static int deg=poly::calcpw2(m+m);
static poly::poly_t mul_t;

int*const&td=f+deg,*const&md=f+m-1,*const&td1=f+m+m-3;

poly::DFT(f,deg);
if(f==g)
for(int*i=f;i!=td;++i)
*i=(i64)*i**i%p;
else{
std::copy(g,g+m-1,mul_t);
std::fill(mul_t+m-1,mul_t+deg,0);

poly::DFT(mul_t,deg);
for(int i=0;i!=deg;++i)
f[i]=(i64)f[i]*mul_t[i]%p;
}

poly::IDFT(f,deg);
for(int*i=f+m-1,*j=f;i!=td1;++i,(++j==md)&&(j=f))
*i&&(math::inc(*j,*i),*i=0);
}

int lg[maxm];
poly::poly_t s,pw;

int main(){
poly::polyinit();

const int g=math::calcg(m);

for(int i=1,pwg=g;i!=m;++i,pwg=(i64)pwg*g%m)
lg[pwg]=i;
lg[1]=0;

for(int i=io,v;i!=0;--i)
if(v=io)
s[lg[v]]=1;

pw[0]=1;
for(int pn=n;pn;pn>>=1,mul(s,s))
if(pn&1)mul(pw,s);
printf("%d",pw[lg[x]]);

return 0;
}