BZOJ3992 SDOI2015序列统计

Problem

BZOJ

Solution

离散对数+循环卷积NTT+快速幂
我做这道题颓了一个下午，我要报警了

求$\prod_{i=1}^n a_i≡x\pmod m$的方案数
我们从dp出发考虑这道题
设f[i][j]为前i项之积模m的值为j的方案数，则有

$$f[i][j]=\sum_{k*r\%m=j}f[i-1][k]*f[i-1][r]$$
这个k和r的条件很类似于卷积，只不过卷积是加法，那么我们会希望把乘法运算转化为加法运算，当然由于对数运算比较麻烦，我们先不考虑。由于m是一个奇素数，那么我们一定可以求出它的一个原根g，我们考虑用离散对数的乘法性质。这个性质在 上篇文章中已经介绍。把f的第二维改做离散对数，即得
$$f[i][j]=\sum_{k+r\%m=j}f[i-1][k]*f[i-1][r]$$
可以发现这其实是一个循环卷积的形式，NTT一下，注意在做循环卷积的时候记得把后面的清零即可
其次这个操作每次都是一样的，用快速幂优化

时间复杂度$O(mlogmlogn)$

Code

#include <algorithm>
#include <cstdio>
using namespace std;
typedef long long ll;
const int maxn=100010,mod=1004535809,G=3;
int n,m,x,sz,cnt,tot,q,fn=1,len,a[maxn],b[maxn],pri[maxn],vis[maxn],stk[maxn];
int ind[maxn],r[maxn];
template <typename Tp> inline void read(Tp &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
    if(ch=='-') f=1,ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
    if(f) x=-x;
}
inline int pls(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
void init()
{
    for(int i=2;i<100000;i++)
    {
        if(!vis[i]) pri[++cnt]=i;
        for(int j=1;j<=cnt&&(ll)i*pri[j]<100000;j++)
        {
            vis[i*pri[j]]=1;
            if(i%pri[j]) break;
        }
    }
}
ll power(ll x,ll y,ll mod)
{
    ll res=1;
    for(;y;y>>=1,x=x*x%mod)
      if(y&1)
        res=res*x%mod;
    return res;
}
int getg(int n)
{
    int tmp=--n;
    for(int i=1;i<=cnt;i++)
      if(tmp%pri[i]==0)
      {
        stk[++tot]=pri[i];
        while(tmp%pri[i]==0) tmp/=pri[i];
      }
    if(tmp^1) stk[++tot]=tmp;
    for(int i=2,f;i<=n;i++)
    {
        f=1;
        for(int j=1;j<=tot&&f;j++)
          if(power(i,n/stk[j],n+1)==1)
            f=0;
        if(f) return i;
    }
    return -1;
}
void input()
{
    read(n);read(m);read(x);read(sz);
    init();q=getg(m);
}
void ntt(int *a,int n,int f)
{
    for(int i=0;i<n;i++)
      if(i<r[i]) swap(a[i],a[r[i]]);
    for(int i=1;i<n;i<<=1)
    {
        int gn;
        gn=power(G,(mod-1)/(i<<1),mod);
        for(int j=0;j<n;j+=(i<<1))
        {
            int x,y,g=1;
            for(int k=0;k<i;k++,g=(ll)g*gn%mod)
            {
                x=a[j+k],y=(ll)g*a[j+k+i]%mod;
                a[j+k]=pls(x,y);a[j+k+i]=dec(x,y);
            }
        }
    }
    if(f==-1)
    {
        int inv=power(n,mod-2,mod);reverse(a+1,a+n);
        for(int i=0;i<n;i++) a[i]=(ll)a[i]*inv%mod;
    }
}
void sqr(int *a)
{
    ntt(a,fn,1);
    for(int i=0;i<fn;i++) a[i]=(ll)a[i]*a[i]%mod;
    ntt(a,fn,-1);
    for(int i=0;i<m-1;i++) a[i]=pls(a[i],a[i+m-1]),a[i+m-1]=0;
}
void mul(int *b,int *a)
{
    ntt(b,fn,1);ntt(a,fn,1);
    for(int i=0;i<fn;i++) b[i]=(ll)b[i]*a[i]%mod;
    ntt(b,fn,-1);ntt(a,fn,-1);
    for(int i=0;i<m-1;i++) b[i]=pls(b[i],b[i+m-1]),b[i+m-1]=0;
}
int main()
{
    #ifndef ONLINE_JUDGE
    freopen("in.txt","r",stdin);
    #endif
    input();
    while(fn<m) fn<<=1,len++;
    fn<<=1;len++;
    for(int i=1;i<fn;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(len-1));
    for(int i=1,now=1;i<m-1;i++)
    {
        now=(ll)now*q%m;
        ind[now]=i;
    }
    for(int i=1,now;i<=sz;i++)
    {
        read(now);
        if(now) a[ind[now]]=1;
    }
    b[0]=1;
    for(;n;n>>=1,sqr(a))
      if(n&1)
        mul(b,a);
    printf("%d\n",b[ind[x]]);
    return 0;
}