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\times 2^{21}+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内元素为a_{1},a_{2},...,a_{|S|}，用M的原根rt分别对应表示为rt^{b1},rt^{b2},...,rt^{b|S|}，x对应表示为rt^{B}。现在中放回地取出N个数，求取出的所有数的指数和与B在模(M-1)意义下同余的方案总数。（大概意会就好，我也不知道说错没）。$。

简单概括下就是把原数的乘积变成了用原根表示的指数相加。

注意！这里有本题的一个坑点：S中的元素是非负的，也就是说可能会有0。然而x是正数，所以遇到要稍微处理一下，否则会WA

这样就可以转换为多项式乘积的形式了。

设$F(y)=a_{0}y_{0}+a_{1}y_{1}+...+a_{M-2}y_{M-2}$，其中$a_{i}$的含义是：将所有数用原根表示后，指数i出现的次数。那么答案就是$(F(y))^{N}$中，满足次数在模$(M-1)$意义下与$B$同余的系数之和。（还是意会即可）

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

时间复杂度：$O(MlogMlogN)$

代码

首先是求原根的算法：

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