SDOI2015 序列统计
问题描述
小C有一个集合S,里面的元素都是小于M的非负整数。他用程序编写了一个数列生成器,可以生成一个长度为N的数列,数列中的每个数都属于集合S。
小C用这个生成器生成了许多这样的数列。但是小C有一个问题需要你的帮助:给定整数x,求所有可以生成出的,且满足数列中所有数的乘积mod
M的值等于x的不同的数列的有多少个。小C认为,两个数列{Ai}和{Bi}不同,当且仅当至少存在一个整数i,满足Ai≠Bi。另外,小C认为这个问题的答案可能很大,因此他只需要你帮助他求出答案mod
1004535809的值就可以了。
输入格式
一行,四个整数,N、M、x、|S|,其中|S|为集合S中元素个数。第二行,|S|个整数,表示集合S中的所有元素。
输出格式
一行,一个整数,表示你求出的种类数mod 1004535809的值。
样例输入
4 3 1 2
1 2
样例输出
8
提示
【样例说明】
可以生成的满足要求的不同的数列有(1,1,1,1)、(1,1,2,2)、(1,2,1,2)、(1,2,2,1)、(2,1,1,2)、(2,1,2,1)、(2,2,1,1)、(2,2,2,2)。
【数据规模和约定】
对于10%的数据,1<=N<=1000;
对于30%的数据,3<=M<=100;
对于60%的数据,3<=M<=800;
对于全部的数据,1<=N<=10^9,3<=M<=8000,M为质数,1<=x<=M-1,输入数据保证集合S中元素不重复
目前为止做过的最难的题之一(我还是太菜了), 做了一天半的样子。最后与六号@Mogician_Evian 讨论了一下,终于搞出来了。
标题中差不多把做法说出来了。其中找M的原根是关键,NTT和快速幂是用来加速的。
讲得有点啰嗦,可能还有错误,直接看代码问题不大。
读题
1.注意到模数1004535809是一个良好的质数,即 1004535809=479×221+1 ,那么我们差不多可以料到这道题多半要用NTT,当然这是后话。
2.M也是质数。用一个质数作为模数,虽然有的题目中确实与算法没多大关系,但在本题中确实一个关键。我们可以利用质数的很多有用的性质。
3.“数列中所有数的乘积mod M的值等于x“。这是本题最难操作的地方,对我们来说,乘法确实难以操作,要想办法把它转化一下。
分析
接下来讨论如何利用原根,将乘法转化为加法。
首先简要介绍一下原根的几点基本知识(摘自nodgd董神的ppt):
注意性质2和性质4。这对本题非常重要。
1.性质4表明,任何质数都有原根。因此,对M是能够把原根求出来的。
2.性质2表明,若求出了M的原根rt,那么一定可以用rt的0~M-2次幂中的某一个对M取模后把集合S中的任意一个元素表示出来。
以样例为例:
M=3,原根rt=2。1=rt
0
,2=rt
1
建立好对应关系之后,问题转变为了这样:
设S内元素为a1,a2,...,a|S|,用M的原根rt分别对应表示为rtb1,rtb2,...,rtb|S|,x对应表示为rtB。现在中放回地取出N个数,求取出的所有数的指数和与B在模(M−1)意义下同余的方案总数。(大概意会就好,我也不知道说错没)。 。
简单概括下就是把原数的乘积变成了用原根表示的指数相加。
注意!这里有本题的一个坑点:S中的元素是非负的,也就是说可能会有0。然而x是正数,所以遇到要稍微处理一下,否则会WA
这样就可以转换为多项式乘积的形式了。
设 F(y)=a0y0+a1y1+...+aM−2yM−2 ,其中 ai 的含义是:将所有数用原根表示后,指数i出现的次数。那么答案就是 (F(y))N 中,满足次数在模 (M−1) 意义下与 B 同余的系数之和。(还是意会即可)
其中,多项式乘法可以用NTT加速,多项式的乘方可以用快速幂的思想解决。注意:
时间复杂度: O(MlogMlogN)
代码
首先是求原根的算法:
设要进行操作的数为 X
1.判断该数是否有原根(在本题中不用判断,下面的代码也不对此进行讨论)
2.找出(X−1) 的所有质因数,用数列 Fi 表示
3.从2开始枚举,依次检验该数是否为原根。设检验的数为 y ,检验方法如下: a.依次算出y(X−1)/Fi 对 X <script type="math/tex" id="MathJax-Element-18">X</script>取模后的值 b.如果这些数里面有一个是1,那么这个数不是原根;反之则是原根。
质数的原根是较多的。我们一般不需要求出所有原根,而是拿最小的一个进行操作。很好的一点是,很多质数最小的原根是很小的(如本题中1004535809的最小原根是3),所以这个算法一般是很快的。
int ksm(int a,int b,int c)
{
int ans=1;a%=c;
while(b)
{
if(b&1)ans=1ll*ans*a%c;
b>>=1;
a=1ll*a*a%c;
}
return ans;
}
int Prime[MAXN];
bool mark[MAXN];
void Euler(int x)
{
int i,j;
mark[1]=true;
for(i=2;i<=x;i++)
{
if(!mark[i])Prime[++Prime[0]]=i;
for(j=1;j<=Prime[0]&&i*Prime[j]<=x;j++)
{
mark[i*Prime[j]]=true;
if(i%Prime[j]==0)break;
}
}
}
int F[MAXN];
void Factor(int x)
{
x--;
int i;
for(i=1;i<=Prime[0];i++)if(x%Prime[i]==0)F[++F[0]]=Prime[i];
}
bool Judge(int x,int p)
{
int i;
for(i=1;i<=F[0];i++)if(ksm(x,(p-1)/F[i],p)==1)return false;
return true;
}
int Primitive_Root(int x)
{
int i;
Euler(x);
Factor(x);
for(i=2;i<x;i++)if(Judge(i,x))return i;
}
本题代码,写得很丑:
#include<stdio.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#define ll long long
using namespace std;
const int MAXN=23456;
ll N,M,X,S,rt,data[MAXN],f[MAXN],AA[MAXN],Ans[MAXN],Len=1;
ll ksm(ll a,ll b,ll c)
{
ll ans=1;a%=c;
while(b)
{
if(b&1)ans=1ll*ans*a%c;
b>>=1;
a=1ll*a*a%c;
}
return ans;
}
ll Prime[MAXN];
bool mark[MAXN];
void Euler(ll x)
{
ll i,j;
mark[1]=true;
for(i=2;i<=x;i++)
{
if(!mark[i])Prime[++Prime[0]]=i;
for(j=1;j<=Prime[0]&&i*Prime[j]<=x;j++)
{
mark[i*Prime[j]]=true;
if(i%Prime[j]==0)break;
}
}
}
ll F[MAXN];
void Factor(ll x)
{
x--;
ll i;
for(i=1;i<=Prime[0];i++)if(x%Prime[i]==0)F[++F[0]]=Prime[i];
}
bool Judge(ll x,ll p)
{
ll i;
for(i=1;i<=F[0];i++)if(ksm(x,(p-1)/F[i],p)==1)return false;
return true;
}
ll Primitive_Root(ll x)
{
ll i;
Euler(x);
Factor(x);
for(i=2;i<x;i++)if(Judge(i,x))return i;
}
ll ntt[MAXN],G=3,P=1004535809;
void NTT(ll A[MAXN],ll n,ll ty)
{
ll i,j,k,m,t0,t1;
for(i=0;i<n;i++)
{
for(j=0,m=1,k=i;m<n;m<<=1,j=(j<<1)|(k&1),k>>=1);
if(i<j)t0=A[i],A[i]=A[j],A[j]=t0;
}
ntt[0]=1;
for(m=1;m<n;m<<=1)
{
t0=ksm(G,(P-1)+ty*(P-1)/(m<<1),P);
for(i=1;i<m;i++)ntt[i]=1ll*ntt[i-1]*t0%P;
for(k=0;k<n;k+=(m<<1))
for(i=k;i<k+m;i++)
{
t0=A[i];
t1=1ll*A[i+m]*ntt[i-k]%P;
A[i]=t0+t1;
A[i+m]=t0-t1;
A[i]-=A[i]>=P?P:0;
A[i+m]+=A[i+m]<0?P:0;
}
}
if(ty==1)return;
t0=ksm(n,P-2,P);
for(i=0;i<n;i++)A[i]=1ll*A[i]*t0%P;
}
ll a[MAXN],b[MAXN];
void Multiple(ll x[MAXN],ll y[MAXN])
{
ll i;
for(i=0;i<Len;i++)a[i]=x[i];
for(i=0;i<Len;i++)b[i]=y[i];
NTT(a,Len,1);NTT(b,Len,1);
for(i=0;i<Len;i++)a[i]=1ll*a[i]*b[i]%P;
NTT(a,Len,-1);
for(i=Len-1;i>=M-1;i--)a[i-M+1]=(a[i-M+1]+a[i])%P,a[i]=0;//保证项数不会太大
memcpy(x,a,sizeof(a));
}
void KSM(ll x)
{
Ans[0]=1;
//为什么这样赋初值呢?可以这样理解:X^0的系数为1 ,其余系数为0。X^0就是1
while(x)
{
if(x&1)Multiple(Ans,AA);
x>>=1;
Multiple(AA,AA);
}
}
main()
{
ll i,tmp;
scanf("%lld%lld%lld%lld",&N,&M,&X,&S);
rt=Primitive_Root(M);//求M原根
for(i=1;i<=S;i++)scanf("%lld",&data[i]);
tmp=1;
for(i=0;i<M-1;i++)f[tmp]=i,tmp=tmp*rt%M; //M-1=phi[M],找对应关系
for(i=1;i<=S;i++)if(data[i])AA[f[data[i]]]++;
while(Len<=(M-1)*2)Len<<=1;//满足NTT要求
KSM(N);
printf("%lld",Ans[f[X]]);
}