bzoj 3992 [SDOI2015] 序列统计 —— NTT (循环卷积+快速幂)

时间:2023-03-09 13:25:33
bzoj 3992 [SDOI2015] 序列统计 —— NTT (循环卷积+快速幂)

题目:https://www.lydsy.com/JudgeOnline/problem.php?id=3992

(学习NTT:https://riteme.github.io/blog/2016-8-22/ntt.html

https://www.cnblogs.com/Mychael/p/9297652.html

http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform#i-15 )

首先,如果把方案数和乘积分别放在系数和次数上,就可以用多项式做了;

方案数放在系数上好说,但次数是相加的,如何表示乘积?

考虑乘积与加法的关系 —— 幂的相乘就是指数相加;

所以可以找出乘积的模数 m 的原根,用其次数相加代表乘积,这个次数好像被称为“指标”;

构造出多项式,由于要取模,所以用 NTT 做;

也就是要把初始的多项式做 n 次幂,可以用快速幂,但注意累乘起来的是系数而不是点值;

指标从0开始或从1开始都可以,也就是把 0 次方作为 1 和把 m-1 次方作为 1 的区别,对应系数的时候要根据这个注意一下(代码中注释里的方案也可);

初始化一个多项式并不是把每个系数都赋值1!而只有第0项是1,这样别的多项式乘过来还是那个多项式;

然后要特别注意读入时去掉0!因为原根系列中没有模出0的,所以以原根为基础的 NTT 算的时候不能考虑0,而反正最后要求的方案中,x >= 1,一旦有0,乘积就是0了,所以0对答案没有影响,就当没给这个数算即可;

一下午的心血...

代码如下:

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
int const xn=(<<),xm=,mod=;
int n,m,rev[xn],g,a[xn],b[xn],lim,r[xm],cnt,pri[xm],inv;
int rd()
{
int ret=,f=; char ch=getchar();
while(ch<''||ch>''){if(ch=='-')f=; ch=getchar();}
while(ch>=''&&ch<='')ret=(ret<<)+(ret<<)+ch-'',ch=getchar();
return f?ret:-ret;
}
int pw(ll a,int b,int md)
{
ll ret=;
for(;b;b>>=,a=(a*a)%md)if(b&)ret=(ret*a)%md;
return ret;
}
void div(int x)
{
for(int i=;i*i<=x;i++)
{
if(x%i)continue;
pri[++cnt]=i; while(x%i==)x/=i;
}
if(x>)pri[++cnt]=x;
}
void init()
{
lim=; int l=;
while(lim<=m+m)lim<<=,l++;
//while(lim<=2*(m-1))lim<<=1,l++;
for(int i=;i<lim;i++)
rev[i]=((rev[i>>]>>)|((i&)<<(l-)));
inv=pw(lim,mod-,mod); if(m==){g=; return;}
div(m-);
for(g=;;g++)
{
bool f=;
for(int j=;j<=cnt;j++)
if(pw(g,(m-)/pri[j],m)==){f=; break;}
if(!f)break;
}
for(int i=,k=g;i<m;i++,k=(ll)k*g%m)r[k]=i;
//for(int i=0,k=1;i<m-1;i++,k=(ll)k*g%m)r[k]=i;//k=1
}
int upt(int x){while(x>=mod)x-=mod; while(x<)x+=mod; return x;}
void ntt(int *a,int tp)
{
for(int i=;i<lim;i++)
if(i<rev[i])swap(a[i],a[rev[i]]);
for(int mid=;mid<lim;mid<<=)
{
int wn=pw(,(mod-)/(mid<<),mod);
if(tp==-)wn=pw(wn,mod-,mod);//
for(int j=,len=(mid<<);j<lim;j+=len)
{
int w=;
for(int k=;k<mid;k++,w=(ll)w*wn%mod)
{
int x=a[j+k],y=(ll)w*a[j+mid+k]%mod;
a[j+k]=upt(x+y); a[j+mid+k]=upt(x-y);
}
}
}
}
void pww()
{
ntt(a,);
for(int i=;i<lim;i++)a[i]=(ll)a[i]*a[i]%mod;
ntt(a,-);
for(int i=;i<lim;i++)a[i]=(ll)a[i]*inv%mod;
for(int i=m;i<lim;i++)a[i%m+]=upt(a[i%m+]+a[i]),a[i]=;//%m+1
//for(int i=m-1;i<lim;i++)a[i%(m-1)]=upt(a[i%(m-1)]+a[i]),a[i]=0;
}
void mul()
{
ntt(a,); ntt(b,);//
for(int i=;i<lim;i++)b[i]=(ll)b[i]*a[i]%mod;
ntt(a,-); ntt(b,-);//
for(int i=;i<lim;i++)
a[i]=(ll)a[i]*inv%mod,b[i]=(ll)b[i]*inv%mod; for(int i=m;i<lim;i++)
a[i%m+]=upt(a[i%m+]+a[i]),a[i]=,
b[i%m+]=upt(b[i%m+]+b[i]),b[i]=;
/*
for(int i=m-1;i<lim;i++)
a[i%(m-1)]=upt(a[i%(m-1)]+a[i]),a[i]=0,
b[i%(m-1)]=upt(b[i%(m-1)]+b[i]),b[i]=0;
*/
}
int main()
{
n=rd(); m=rd(); init();
int p=rd(),num=rd();
for(int i=,x;i<=num;i++)
{
x=rd();
if(x)a[r[x]]=;//x!=0 !!
}
int t=n;
//for(int i=0;i<lim;i++)b[i]=1;
b[]=;//!
for(;t;t>>=,pww())if(t&)mul();
printf("%d\n",b[r[p]]);
return ;
}