[SDOI2015]序列统计 （NTT：快速数论变换 + 快速幂）

调试了一个多小时可算是调试出来了，我的NTT真是漏洞百出，还好最后AC了…

题目描述

小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<=109，3<=M<=8000，M为质数，1<=x<=M-1，输入数据保证集合S中元素不重复

思路

在想出标算之前，先想一些显然会超时但是具有正确性的算法，比如DP。

令fi,j表示一个具有i个数的序列，其所有元素的乘积除以M的余数为j的方案数。

那么就有状态转移方程（其中k−1表示k的数论倒数）：

fi,j=∑k∈Sfi−1,(j⋅k−1modM)

乘法运算不利于状态转移，考虑如何把乘法降为加法，这不由得让我想起了计算尺的使用原理（好像现在的孩子都没用过计算尺…）——就是把两个数都表示成同一个数的若干次幂的形式，求着两个数的乘积，就只需要把两个指数加和即可。即：

ga×gb=ga+b

然而，在模意义下，恰好有这样的一个数，能使得剩余系中的每一个元素恰都可以表示为这个数的整数次幂的形式，这个数被称为模意义下的“原根”。

这样的话我们，可以先求出模M意义下的原根g。再把剩余系中的没一个数，都表示成ga的形式，再看看DP转移会变成什么样。

fi,j=∑k∈Sfi−1,(j⋅k−1modM)

fi,gj=∑gk∈Sfi−1,(gj⋅(gk)−1modM)=∑gk∈Sfi−1,(gj−kmodM)

不妨用j替换gj，用k替换gk：

p.s.本文不规范的表示：定义：logg(S)={x|gx∈SmodM}。

fi,j=∑k∈logg(S)fi−1,(j−kmodM−1)

/*
注：因为M是质数，所以根据费马小定理有：

gcd(a,M)=1→aM−1≡1modM

也就是说：

ab≡a(bmodM−1)modM

*/

但是上面的式子是有问题的，因为不同的数在模意义下求其以原根为底的离散对数可能相同。所以，上式实际上应该是：

fi,j=∑k=0M−2fi−1,(j−kmodM−1)⋅Cntk

其中Cntk表示k在logg(S)中“理应”出现的次数。

考虑到N≤109，空间一定是存不下f的，考虑使用滚动数组。令A=fi，B=fi−1，替换到原式中：

Aj=∑k=0M−2Bj−kmodM−1⋅Cntk=∑k=0M−2Bj−k⋅Cntk    (1)

貌似发现了一些神奇的事情：A恰好是B和Cnt的卷积，即：A=B×Cnt。

/*

注：（1）处之所以直接删去了“modM−1”，是因为在进行卷积的过程中，我们可以把所有的大于等于M-1一 一对应到，0~M-2之间，并把结果累加到相应的位置。

*/

最后，考虑递推的初始状态：假定在最开始的时候有一个数字“1”，这对结果并没有影响。令：Base0=1，Basei=0(i≠0)(Base是一个向量)。答案即为：

(Base×CntN)logg(x)

CntN可以用“NTT + 快速幂”求得。

附上AC代码：

#include<stdio.h>#include<stdlib.h>#include<string.h>typedef long long LL;void exgcd(int a,int b,int& x,int& y){    if(b==0){        x=1;        y=0;    }else{        exgcd(b,a%b,y,x);        y-=int(a/b)*x;    }}int Inv(int a,int p){//求逆元    int x,y;    exgcd(a,p,x,y);    return (x%p+p)%p;}int Qpow(int a,int b,int p){//快速幂    if(b<0){        b=-b;        a=Inv(a,p);    }    LL ans=1,mul=a%p;    while(b){        if(b&1)            ans=ans*mul%p;        b>>=1;        mul=mul*mul%p;    }    return ans;}const int maxn=8193*2;const int MOD=479*(1<<21)+1,G=3;//NTT的必要常量LL a[maxn],b[maxn];int rev[maxn];void get_rev(int bit){    for(int i=0;i<(1<<bit);i++)        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));}void swap(LL& a,LL& b){    LL t=a;a=b;b=t;}void ntt(LL* a,int n,int dft){//NTT主过程    for(int i=0;i<n;i++){        if(i<rev[i])            swap(a[i],a[rev[i]]);    }    for(int step=1;step<n;step<<=1){        LL wn=Qpow(G,dft*(MOD-1)/(2*step),MOD);        for(int j=0;j<n;j+=step<<1){            LL wnk=1;            for(int k=j;k<j+step;k++){                LL x=a[k]%MOD;                LL y=(a[k+step]*wnk)%MOD;                a[k]=(x+y)%MOD;                a[k+step]=((x-y)%MOD+MOD)%MOD;                wnk=(wnk*wn)%MOD;            }        }    }    if(dft==-1){        int nI=Inv(n,MOD);        for(int i=0;i<n;i++)            a[i]=a[i]*nI%MOD;    }}int N,M,bit,s;struct Item{//定义"向量"    LL cnt[maxn];    void init(){        memset(cnt,0,sizeof(cnt));    }};Item operator*(Item A,Item B){//定义向量卷积    ntt(A.cnt,s,1);    ntt(B.cnt,s,1);    for(int i=0;i<s;i++)        A.cnt[i]=A.cnt[i]*B.cnt[i]%MOD;    ntt(A.cnt,s,-1);    for(int i=M-1;i<s;i++){        (A.cnt[i%(M-1)]+=A.cnt[i])%=MOD;        A.cnt[i]=0;    }    return A;}Item Qpow(Item a,int b){//卷积快速幂    Item ans,mul=a;    ans.init();ans.cnt[0]=1;    while(b){        if(b&1)            ans=ans*mul;        b>>=1;        mul=mul*mul;    }    return ans;}int pcnt,ptab[maxn];bool checkG(int a,int p){//判断一个数是否是原根    for(int i=1;i<=pcnt;i++){        if(Qpow(a,(p-1)/ptab[i],p)==1)            return false;    }    return true;}int solveG(int p){    int n=p-1;pcnt=0;//对p-1进行分解质因数    for(int i=2;i*i <=p-1 && n!=1;i++){        if(n%i==0){            ptab[++pcnt]=i;            while(n%i==0){                n/=i;            }        }    }    if(n!=1)ptab[++pcnt]=n;    for(int i=2;i<p;i++){//一个一个的去试        if(checkG(i,p)){            return i;        }    }    return -1;}int zcnt=0;//统计0出现的次数(在本题数据中没有意义)int trans[maxn];Item C;//即上文中所说的Cnt向量int main(){    int x,cnt;    scanf("%d%d%d%d",&N,&M,&x,&cnt);    for(bit=1,s=2;(1<<bit)<(M-1)*2-1;bit++)s<<=1;    get_rev(bit);    int g=solveG(M),nw=1;//求g<-原根    for(int i=0;i<(M-1);i++){//init trans        trans[nw]=i;//即g^trans[nw]=nw        nw=(nw*g)%M;    }    for(int i=1;i<=cnt;i++){//input set        int nx;scanf("%d",&nx);//预处理Cnt向量        nx%=M;        if(nx==0){            zcnt++;        }else{            C.cnt[trans[nx]]++;        }    }    C=Qpow(C,N);    if(x!=0)        printf("%lld",C.cnt[trans[x]]);//输出结果    else ;//其实当x=0的情况需要特判,但好像并没有这样的数据...    return 0;}

NTT可以用来计算模意义下的卷积，效率和空间相对于FFT都有所提升，但只适用于整数，且被模的质数一般只能为“费马素数”（本题中的479⋅221+1=1004535809就是一个费马素数）。