SDOI2015 序列统计 NTT+原根+快速幂

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=rt0，2=rt1

建立好对应关系之后，问题转变为了这样：

设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加速，多项式的乘方可以用快速幂的思想解决。注意：(F(y))N的项数可以非常多，为了不超过空间限制，在每一次NTT后我们都要把次数高于(M−2)的项的系数加到前面对应的项的系数里面去。这样就只需要2M左右的空间了。

时间复杂度：O(MlogMlogN)

代码

首先是求原根的算法：

设要进行操作的数为X
1.判断该数是否有原根（在本题中不用判断，下面的代码也不对此进行讨论）
2.找出(X−1)的所有质因数，用数列Fi表示
3.从2开始枚举，依次检验该数是否为原根。设检验的数为y,检验方法如下： a.依次算出y(X−1)/Fi对X取模后的值 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 longusing 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]]);}