[BZOJ3992][SDOI2015]序列统计（dp+NTT+快速幂）

题目描述

传送门

题解

首先考虑，如果题目让求的是和的方案数，怎么dp以及优化
令$F(i,j)$表示一共选了i个数，和在模m意义下为j的方案数
那么$F(i,j)=\sum\limits_{k=1}^sF(i-1,j-S_k)$
这个dp可以写成卷积的形式并用NTT优化
令$b(i)$表示是否存在i这个数，0/1
令$a$表示$F(i-1)$，$F'$表示$F(i)$
那么$F'(i)=\sum\limits_{j=0}^{m-1}a(i-j)b(j)$，强行把背包卷了一下
可以看出做一次卷积就相当于一次$F(i-1)->F(i)$的转移，那么做n次就能求出$F(n)$
卷积是满足快速幂性质的，用快速幂加速一下就行了

但是这道题的转移方程是乘除而非加减
$F(i,j)=\sum\limits_{k=1}^sF(i-1,{j\over S_k})$
要将乘除转换为加减，需要用到原根
若P的原根为g，那么当且仅当x=P-1时，$g^x = 1\pmod P$成立
原根有一个非常有趣的性质，就是当$0<i<P$时，$g^i\space mod \space p$互异（不存在0）
那么，对于[1..m-1]的中的所有的数，都可以对应找到原根唯一的一个指数
可以将dp方程转化一下$F(i,g^x)=\sum\limits_{k=1}^s F(i-1,g^{x-y})$，其中$g^x=j,g^y=S_k$
同时dp的目标状态就变成了$F(n,g^z)$，其中$g^z=x$
将写得更清楚一些$F(i,x)=\sum\limits_{y=1}^{m-1}\sum\limits_{g^y\in S}F(i-1,x-y)$
这样我们就成功地将乘除变成了加减，和上面类似
令$b(i)$表示是否存在$g^i$这个数，0/1
令$a$表示$F(i-1)$，$F'$表示$F(i)$
那么$F'(i)=\sum\limits_{i=0}^{m-1}a(i-j)b(j)$，又是一个可以用NTT和快速幂加速的卷积形式
注意由于原根的性质$b(0)=0$，细节比较多

代码

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
#define LL long long
#define N 20005

const LL Mod=1004535809LL;
int cnt,m,s,x,n,y,L,R[N];
LL g,a[N],b[N],f[N],ans[N];
bool vis[N];

LL fast_pow(LL a,LL p,LL Mod)
{
    LL ans=1LL;
    for (;p;p>>=1LL,a=a*a%Mod)
        if (p&1LL)
            ans=ans*a%Mod;
    return ans;
}
LL get_g(int x)
{
    if (x==2) return 1;
    for (int i=2;;++i)
    {
        bool flag=1;
        for (int j=2;j*j<x;++j)
            if (fast_pow(i,(x-1)/j,x)==1)
            {
                flag=false;
                break;
            }
        if (flag) return (LL)i;
    }
}

void NTT(LL a[N],int n,int opt)
{
    for (int i=0;i<n;++i)
        if (i<R[i]) swap(a[i],a[R[i]]);
    for (int k=1;k<n;k<<=1)
    {
        LL wn=fast_pow(3LL,(Mod-1)/(k<<1),Mod);
        for (int i=0;i<n;i+=(k<<1))
        {
            LL w=1LL;
            for (int j=0;j<k;++j,w=w*wn%Mod)
            {
                LL x=a[i+j],y=w*a[i+j+k]%Mod;
                a[i+j]=(x+y)%Mod,a[i+j+k]=(x-y+Mod)%Mod;
            }
        }
    }
    if (opt==-1) reverse(a+1,a+n);
}
void cvlt(LL f[N],LL g[N])
{
    memset(a,0,sizeof(a));memset(b,0,sizeof(b));
    for (int i=0;i<n;++i) a[i]=f[i];
    for (int i=0;i<n;++i) b[i]=g[i];
    NTT(a,n,1);NTT(b,n,1);
    for (int i=0;i<n;++i) a[i]=a[i]*b[i]%Mod;
    NTT(a,n,-1);
    LL inv=fast_pow(n,Mod-2,Mod);
    for (int i=0;i<n;++i) a[i]=a[i]*inv%Mod;
    for (int i=0;i<m-1;++i) f[i]=(a[i]+a[i+m-1])%Mod;
}
void fast_cvlt(int p)
{
    ans[0]=1LL;
    for (;p;p>>=1,cvlt(f,f))
        if (p&1)
            cvlt(ans,f);
}

int main()
{
    scanf("%d%d%d%d",&cnt,&m,&x,&s);
    for (int i=1;i<=s;++i) scanf("%d",&y),vis[y]=1;
    g=get_g(m);
    LL now=1LL;int pos=-1;
    for (int i=0;i<m-1;++i,now=(now*g)%m)
    {
        if (vis[now]) f[i]=1;
        if (now==x) pos=i;
    }

    for (n=1;n<=(m-1)<<1;n<<=1) ++L;
    for (int i=0;i<n;++i)
        R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));

    fast_cvlt(cnt);
    if (pos!=-1) printf("%I64d\n",ans[pos]);
    else puts("0");
}