P3321 [SDOI2015]序列统计

题意

首先考虑暴力DP怎么做：

设\(f_{i,j}\)表示选了\(i\)个数，乘积为\(j\)的方案数。
\(f_{i,j}=\sum\limits_{a*b\%m= j}f_{i-1,a}*f_{1,b}\)
这是\(O(nm^2)\)的，我们考虑优化：

1.首先注意到\(n\)很大，于是考虑倍增地转移：
\(f_{2*i,j}=\sum\limits_{a*b\%m= j}f_{i,a}*f_{i,b}\)
于是就变为\(O(m^2logn)\)的了。

2.发现下标很有意思，似乎是卷积的形式，但是这是乘积的形式，我们考虑将它化为加法的形式，即取\(log\)。由于是模意义下的，我们考虑指标。

我们先求出\(m\)的一个原根\(g\)。根据原根的性质，我们对于一个\(x\in [0,m-1]\)，只有唯一一个\(k\)，满足\(g^k\equiv x\pmod{p}\)，因此此我们能将上式变为：
\(f_{2*i,j}=\sum\limits_{g^a*g^b\%m= g^j}f_{i,a}*f_{i,b}\)
\(f_{2*i,j}=\sum\limits_{a+b\%(m-1)= j}f_{i,a}*f_{i,b}\)

于是就可以\(NTT\)加速了。

注意我们对下标取了模，但实际上\(NTT\)时是不会对下标取模的，也就是说\(f_{i,j}\)的答案会存在第\(i\)位和第\(i+m-1\)位，我们要将它们加起来。

code:

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxm=32010;
const ll mod=1004535809;
const ll G=3;
const ll invG=334845270;
int n,m,K,S,lim=1,len;
int pos[maxm],pri[maxm];
ll f[maxm],g[maxm],a[maxm],b[maxm],c[maxm];
unordered_map<int,int>mp;
inline ll power(ll x,ll k,ll mod)
{
	ll res=1;
	while(k)
	{
		if(k&1)res=res*x%mod;
		x=x*x%mod;k>>=1;
	}
	return res;
}
inline int getroot(int mod)
{
	int tmp=mod-1,tot=0;
	for(int i=2;i*i<=tmp;i++)
	{
		if(tmp%i)continue;
		pri[++tot]=i;
		while(tmp%i==0)tmp/=i;
	}
	if(tmp>1)pri[++tot]=tmp;
	for(int i=2;i<=mod-1;i++)
	{
		bool flag=1;
		for(int j=1;j<=tot;j++)
			if(power(i,(mod-1)/pri[j],mod)==1){flag=0;break;}
		if(flag)return i;
	}
	return -1;
}
inline void NTT(ll* a,int op)
{
    for(int i=0;i<lim;i++)if(i<pos[i])swap(a[i],a[pos[i]]);
    for(int l=1;l<lim;l<<=1)
    {
        ll wn=power(op==1?G:invG,(mod-1)/(l<<1),mod);
        for(int i=0;i<lim;i+=l<<1)
        {
           	ll w=1;
            for(int j=0;j<l;j++,w=w*wn%mod)
            {
                int x=a[i+j],y=w*a[i+l+j]%mod;
                a[i+j]=(x+y)%mod;a[i+l+j]=(x-y+mod)%mod;
            }
        }
    }
    if(op==1)return;
    ll inv=power(lim,mod-2,mod);
    for(int i=0;i<lim;i++)a[i]=a[i]*inv%mod;
}
inline void mul(ll* f,ll* g,ll* h)
{
	for(int i=0;i<lim;i++)a[i]=f[i],b[i]=g[i];
	NTT(a,1);NTT(b,1);
	for(int i=0;i<lim;i++)a[i]=a[i]*b[i]%mod;
	NTT(a,-1);
	for(int i=0;i<m-1;i++)c[i]=(a[i]+a[i+m-1])%mod;
	for(int i=0;i<m-1;i++)h[i]=c[i];
}
int main()
{
	scanf("%d%d%d%d",&n,&m,&K,&S);
	int gi=getroot(m);
	for(int i=0;i<m-1;i++)mp[power(gi,i,m)]=i;
	for(int i=1;i<=S;i++)
	{
		int x;scanf("%d",&x);x%=m;
		if(x)f[mp[x]]++;
	}
	while(lim<=2*m)lim<<=1,len++;
	for(int i=0;i<lim;i++)pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
	g[mp[1]]=1;
	while(n)
	{
		if(n&1)mul(g,f,g);
		mul(f,f,f);n>>=1;
	}
	printf("%lld\n",g[mp[K]]);
	return 0;
}