BZOJ 3992: [SDOI2015]序列统计 快速幂+NTT（离散对数下）

3992: [SDOI2015]序列统计

Description

小C有一个集合S，里面的元素都是小于M的非负整数。他用程序编写了一个数列生成器，可以生成一个长度为N的数列，数列中的每个数都属于集合S。

小C用这个生成器生成了许多这样的数列。但是小C有一个问题需要你的帮助：给定整数x，求所有可以生成出的，且满足数列中所有数的乘积mod M的值等于x的不同的数列的有多少个。小C认为，两个数列{Ai}和{Bi}不同，当且仅当至少存在一个整数i，满足Ai≠Bi。另外，小C认为这个问题的答案可能很大，因此他只需要你帮助他求出答案mod 1004535809的值就可以了。

Input

一行，四个整数，N、M、x、|S|，其中|S|为集合S中元素个数。第二行，|S|个整数，表示集合S中的所有元素。

Output

一行，一个整数，表示你求出的种类数mod 1004535809的值。

Sample Input

4 3 1 2
1 2

Sample Output

8

HINT

【样例说明】

可以生成的满足要求的不同的数列有(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中元素不重复

Source

题解：

　　利用原根G

　　将序列的乘积替换成加法

　　想快速幂那般的计算，得到长度为n的序列

　　每次NTT加速

#include<bits/stdc++.h>
using namespace std;
#pragma comment(linker, "/STACK:102400000,102400000")
#define ls i<<1
#define rs ls | 1
#define mid ((ll+rr)>>1)
#define pii pair<int,int>
#define MP make_pair
typedef long long LL;
typedef unsigned long long ULL;
const long long INF = 1e18+1LL;
const double pi = acos(-1.0);
const int N = 6e4+10, M = 2e2+20,inf = 2e9;
int MOD;
inline int muls(int a, int b){
    return (long long)a * b % MOD;
}
int power(int a, int b){
    int ret = 1;
    for (int t = a; b; b >>= 1){
        if (b & 1)ret = muls(ret, t);
        t = muls(t, t);
    }
    return ret;
}
int cal_root(int mod)
{
    int factor[20], num = 0, s = mod - 1;
    MOD = mod--;
    for (int i = 2; i * i <= s; i++){
        if (s % i == 0){
            factor[num++] = i;
            while (s % i == 0)s /= i;
        }
    }
    if (s != 1)factor[num++] = s;
    for (int i = 2;; i++){
        int j = 0;
        for (; j < num && power(i, mod / factor[j]) != 1; j++);
        if (j == num)return i;
    }
}

const LL G = 3LL, P = 1004535809LL;
LL mul(LL x,LL y){
    return (x*y-(LL)(x/(long double)P*y+1e-3)*P+P)%P;
}
LL qpow(LL x,LL k,LL p){
    LL ret=1;
    while(k){
        if(k&1) ret=mul(ret,x);
        k>>=1;
        x=mul(x,x);
    }
    return ret;
}
LL wn[50];
void getwn(){
    for(int i=1; i<=35; ++i){
        int t=1<<i;
        wn[i]=qpow(G,(P-1)/t,P);
    }
}int len;
void NTT_init() {
    getwn();
}

void NTT(LL y[],int op){
    for(int i=1,j=len>>1,k; i<len-1; ++i){
        if(i<j) swap(y[i],y[j]);
        k=len>>1;
        while(j>=k){
            j-=k;
            k>>=1;
        }
        if(j<k) j+=k;
    }
    int id=0;
    for(int h=2; h<=len; h<<=1) {
        ++id;
        for(int i=0; i<len; i+=h){
            LL w=1;
            for(int j=i; j<i+(h>>1); ++j){
                LL u=y[j],t=mul(y[j+h/2],w);
                y[j]=u+t;
                if(y[j]>=P) y[j]-=P;
                y[j+h/2]=u-t+P;
                if(y[j+h/2]>=P) y[j+h/2]-=P;
                w=mul(w,wn[id]);
            }
        }
    }
    if(op==-1){
        for(int i=1; i<len/2; ++i) swap(y[i],y[len-i]);
        LL inv=qpow(len,P-2,P);
        for(int i=0; i<len; ++i) y[i]=mul(y[i],inv);
    }
}
LL mo[N],fmo[N],now[N],xx[N],yy[N];
int m;
void cal(LL *x,LL *y,LL *z) {
    for(int i = 0; i < len; ++i) xx[i] = x[i],yy[i] = y[i];
    NTT(xx,1);NTT(yy,1);
    for(int i = 0; i < len; ++i) z[i] = xx[i] * yy[i] % P,now[i] = 0;
    NTT(z,-1);

    for(int i = 0; i <= 2*m-2; ++i)
        now[mo[i%(m-1)]] += z[i],z[i] = 0;

    for(int i = 0; i < len; ++i) z[i] = 0;
    for(int i = 0; i <= m-1; ++i) {
        z[fmo[i]] += now[i];
    }
}
LL y[N],ans[N],ny[N],num[N],an[N];
int n,x,s,root,cnt0;

int main() {
    NTT_init();
    scanf("%d%d%d%d",&n,&m,&x,&s);
    root = cal_root(m);
    for(int i = 0, t = 1; i < m-1; ++i,t = t*root%m)
        mo[i] = t,fmo[t] = i;
    for(int i = 1; i <= s; ++i) {
        int xx;
        scanf("%d",&xx);
        if(xx == 0)
            cnt0+=1;
        else
            num[fmo[xx]]+=1;
    }
    if(x == 0) {
        LL L = (qpow(2LL,cnt0,P) - 1 + P) % P + qpow(2LL,s-cnt0,P)%P;
        printf("%lld\n",L%P);
        return 0;
    }
    len = 1;
    int f = 1;
    while(len <= 2*m) len<<=1;
    for(int i = 0; i < len; ++i) y[i] = num[i],ans[i] = 0;
    int j = 1;
    int chan = 0;
    while(n) {
        if(n&1) {
            if(f) {
                f = 0;
                for(int i = 0; i < len; ++i) ans[i] = y[i];
            }
            else {
                cal(ans,y,ans);
            }
        }
        cal(y,y,y);
        n>>=1;
    }
    int l = 0;
    for(int i = 0; i < len; ++i) {
        an[mo[i%(m-1)]] += ans[i],an[mo[i%(m-1)]]%=P;
    }
    cout<<an[x]<<endl;
    return 0;
}