[BZOJ 3992][SDOI2015]序列统计

3992: [SDOI2015]序列统计

Time Limit: 30 Sec  Memory Limit: 128 MB
Submit: 2275  Solved: 1090
[Submit][Status][Discuss]

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中的所有元素。

1<=N<=10^9，3<=M<=8000，M为质数

0<=x<=M-1，输入数据保证集合S中元素不重复x∈[1,m-1]

集合中的数∈[0,m-1]

Output

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

Sample Input

4 3 1 2
1 2

Sample Output

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)

 

题解

这题的要求非常类似背包, 但是值之间做的贡献是乘积的形式. 我们考虑把它转成加法.

怎么转加法呢? 当然是取对数了!

离散对数哪家强? SDOI找...(划掉)

注意到如果数列里有 $0$ 那么贡献就是 $0$, 但是题目规定要求的目标值非 $0$, 所以我们首先先把 $0$ 扔掉

然后暴力计算出 $m$ 的一个原根, 同时处理出每个数的对数

那么这就是一个长度为 $m-1$ 的循环卷积辣!

然而循环卷积需要膜数有 $m-1$ 次单位根(吧), 于是只能NTT倍增爆算了...(或者博主naive?)

参考代码

多项式左右移操作现已加入豪华午餐

#include <bits/stdc++.h>

const int G=3;
const int DFT=1;
const int IDFT=-1;
const int MAXN=550000;
const int MOD=1004535809;
const int INV2=(MOD+1)>>1;
const int PHI=MOD-1;

typedef std::vector<int> Poly;

Poly Sqrt(Poly);
void Read(Poly&);
Poly Inverse(Poly);
Poly Ln(const Poly&);
Poly Exp(const Poly&);
void Print(const Poly&);
void NTT(Poly&,int,int);
Poly Pow(const Poly&,int);
Poly Integral(const Poly&);
Poly Derivative(const Poly&);
Poly operator*(Poly,Poly);
Poly operator/(Poly,Poly);
Poly operator%(Poly,Poly);
Poly operator<<(const Poly&,int);
Poly operator>>(const Poly&,int);
Poly operator+(const Poly&,const Poly&);
Poly operator-(const Poly&,const Poly&);

int rev[MAXN];

int FindG(int);
int NTTPre(int);
int Sqrt(int,int);
int Pow(int,int,int);
int Log(int,int,int);
int ExGCD(int,int,int&,int&);

int k,m,x,n;
int lg[MAXN];

int main(){
    scanf("%d%d%d%d",&k,&m,&x,&n);
    Poly a(m-1);
    int g=FindG(m);
    for(int i=0,pw=1;i<m-1;i++,pw=pw*g%m)
        lg[pw]=i;
    for(int i=0;i<n;i++){
        int x;
        scanf("%d",&x);
        if(x)
            a[lg[x]]=1;
    }
    --m;
    Poly ans(a);
    --k;
    while(k>0){
        if(k&1){
            ans=a*ans;
            ans=ans+(ans>>m);
            ans.resize(m);
        }
        a=a*a;
        a=a+(a>>m);
        a.resize(m);
        k>>=1;
    }
    printf("%d
",ans[lg[x]]);
    return 0;
}

Poly operator<<(const Poly& a,int x){
    Poly ans(a.size()+x);
    for(size_t i=0;i<a.size();i++)
        ans[i+x]=a[i];
    return ans;
}

Poly operator>>(const Poly& a,int x){
    Poly ans(a.size()-x);
    for(size_t i=x;i<a.size();i++)
        ans[i-x]=a[i];
    return ans;
}

int FindG(int mod){
    for(int g=2;g<mod;g++){
        bool flag=true;
        for(int i=1,pw=g;i<mod-1;i++,pw=pw*g%mod){
            if(pw==1){
                flag=false;
                break;
            }
        }
        if(flag)
            return g;
    }
    return -1;
}

void Read(Poly& a){
    for(auto& i:a)
        scanf("%d",&i);
}

void Print(const Poly& a){
    for(auto i:a)
        printf("%d ",i);
    puts("");
}

Poly Pow(const Poly& a,int k){
    Poly log=Ln(a);
    for(auto& i:log)
        i=1ll*i*k%MOD;
    return Exp(log);
}

Poly Sqrt(Poly a){
    int len=a.size();
    if(len==1)
        return Poly(1,Sqrt(a[0],MOD));
    else{
        Poly b=a;
        b.resize((len+1)>>1);
        b=Sqrt(b);
        b.resize(len);
        Poly inv=Inverse(b);
        int bln=NTTPre(inv.size()+a.size());
        NTT(a,bln,DFT);
        NTT(inv,bln,DFT);
        for(int i=0;i<bln;i++)
            a[i]=1ll*a[i]*INV2%MOD*inv[i]%MOD;
        NTT(a,bln,IDFT);
        for(int i=0;i<len;i++)
            b[i]=(1ll*b[i]*INV2%MOD+a[i])%MOD;
        return b;
    }
}

Poly Exp(const Poly& a){
    size_t len=1;
    Poly ans(1,1),one(1,1);
    while(len<(a.size()<<1)){
        len<<=1;
        Poly b=a;
        b.resize(len);
        ans=ans*(one-Ln(ans)+b);
        ans.resize(len);
    }
    ans.resize(a.size());
    return ans;
}

Poly Ln(const Poly& a){
    Poly ans=Integral(Derivative(a)*Inverse(a));
    ans.resize(a.size());
    return ans;
}

Poly Integral(const Poly& a){
    int len=a.size();
    Poly ans(len+1);
    for(int i=1;i<len;i++)
        ans[i]=1ll*a[i-1]*Pow(i,MOD-2,MOD)%MOD;
    return ans;
}

Poly Derivative(const Poly& a){
    int len=a.size();
    Poly ans(len-1);
    for(int i=1;i<len;i++)
        ans[i-1]=1ll*a[i]*i%MOD;
    return ans;
}

Poly operator/(Poly a,Poly b){
    int n=a.size()-1,m=b.size()-1;
    Poly ans(1);
    if(n>=m){
        std::reverse(a.begin(),a.end());
        std::reverse(b.begin(),b.end());
        b.resize(n-m+1);
        ans=Inverse(b)*a;
        ans.resize(n-m+1);
        std::reverse(ans.begin(),ans.end());
    }
    return ans;
}

Poly operator%(Poly a,Poly b){
    int n=a.size()-1,m=b.size()-1;
    Poly ans;
    if(n<m)
        ans=a;
    else
        ans=a-(a/b)*b;
    ans.resize(m);
    return ans;
}

Poly operator*(Poly a,Poly b){
    int len=a.size()+b.size()-1;
    int bln=NTTPre(len);
    NTT(a,bln,DFT);
    NTT(b,bln,DFT);
    for(int i=0;i<bln;i++)
        a[i]=1ll*a[i]*b[i]%MOD;
    NTT(a,bln,IDFT);
    a.resize(len);
    return a;
}

Poly operator+(const Poly& a,const Poly& b){
    Poly ans(std::max(a.size(),b.size()));
    std::copy(a.begin(),a.end(),ans.begin());
    for(size_t i=0;i<b.size();i++)
        ans[i]=(ans[i]+b[i])%MOD;
    return ans;
}

Poly operator-(const Poly& a,const Poly& b){
    Poly ans(std::max(a.size(),b.size()));
    std::copy(a.begin(),a.end(),ans.begin());
    for(size_t i=0;i<b.size();i++)
        ans[i]=(ans[i]+MOD-b[i])%MOD;
    return ans;
}

Poly Inverse(Poly a){
    int len=a.size();
    if(len==1)
        return Poly(1,Pow(a[0],MOD-2,MOD));
    else{
        Poly b(a);
        b.resize((len+1)>>1);
        b=Inverse(b);
        int bln=NTTPre(b.size()*2+a.size());
        NTT(a,bln,DFT);
        NTT(b,bln,DFT);
        for(int i=0;i<bln;i++)
            b[i]=(2ll*b[i]%MOD-1ll*b[i]*b[i]%MOD*a[i]%MOD+MOD)%MOD;
        NTT(b,bln,IDFT);
        b.resize(len);
        return b;
    }
}

void NTT(Poly& a,int len,int opt){
    a.resize(len);
    for(int i=0;i<len;i++)
        if(rev[i]>i)
            std::swap(a[i],a[rev[i]]);
    for(int i=1;i<len;i<<=1){
        int step=i<<1;
        int wn=Pow(G,(PHI+opt*PHI/step)%PHI,MOD);
        for(int j=0;j<len;j+=step){
            int w=1;
            for(int k=0;k<i;k++,w=1ll*w*wn%MOD){
                int x=a[j+k];
                int y=1ll*a[j+k+i]*w%MOD;
                a[j+k]=(x+y)%MOD;
                a[j+k+i]=(x-y+MOD)%MOD;
            }
        }
    }
    if(opt==IDFT){
        int inv=Pow(len,MOD-2,MOD);
        for(int i=0;i<len;i++)
            a[i]=1ll*a[i]*inv%MOD;
    }
}

int NTTPre(int n){
    int bln=1,bct=0;
    while(bln<n){
        bln<<=1;
        ++bct;
    }
    for(int i=0;i<bln;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bct-1));
    return bln;
}

inline int Pow(int a,int n,int p){
    int ans=1;
    while(n>0){
        if(n&1)
            ans=1ll*a*ans%p;
        a=1ll*a*a%p;
        n>>=1;
    }
    return ans;
}

int ExGCD(int a,int b,int& x,int& y){
    if(b==0){
        x=1,y=0;
        return a;
    }
    else{
        int gcd=ExGCD(b,a%b,y,x);
        y-=x*(a/b);
        return gcd;
    }
}

inline int Sqrt(int a,int p){
    int s=Pow(G,Log(G,a,p)>>1,p);
    return std::min(s,MOD-s);
}

inline int Log(int a,int x,int p){
    int s=sqrt(p+0.5);
    int inv=Pow(Pow(a,s,p),p-2,p);
    std::unordered_map<int,int> m;
    m[1]=0;
    int pow=1;
    for(int i=1;i<s;i++){
        pow=1ll*a*pow%p;
        if(!m.count(pow))
            m[pow]=i;
    }
    for(int i=0;i<s;i++){
        if(m.count(x))
            return i*s+m[x];
        x=1ll*x*inv%MOD;
    }
    return -1;
}