传送门
题解:
还行的一道题。其实我主要是来复习FFT怎么写的,突然发现自己不会FFT只会NTT
首先,如果注意到 f f f是一个二阶线性递推数列,这道题就已经做完了。
考虑利用特征方程解出 f f f的通项公式。
该数列特征方程为 λ 2 − 2 λ − 3 = 0 \lambda^2-2\lambda-3=0 λ2−2λ−3=0,两个根为 λ 1 = 3 , λ 2 = − 1 \lambda_1=3,\lambda_2=-1 λ1=3,λ2=−1,设 C 1 , C 2 C_1,C_2 C1,C2为待定系数。则其通项公式为 f n = C 1 λ 1 n + C 2 λ 2 n f_n=C_1\lambda_1^n+C_2\lambda_2^n fn=C1λ1n+C2λ2n
由于题目给出了 f 0 , f 1 f_0,f_1 f0,f1,我们可以把 C 1 , C 2 C_1,C_2 C1,C2解出来。
考虑题目要求和的东西,我们可以对于 λ 1 \lambda_1 λ1和 λ 2 \lambda_2 λ2两个部分分别计算。
其中由于
λ
2
=
−
1
\lambda_2=-1
λ2=−1,可以直接用组合数学算,注意由于模数的关系要用Lucas定理(结果数据最大的点
n
=
5
e
4
n=5e4
n=5e4,完全不用Lucas)
于是现在问题我们要求出选 k k k个元素,设 ∑ \sum ∑为其和,对于 λ Σ \lambda^\Sigma λΣ求和。
。。。。
非常显然直接构造多项式 ∏ i = 1 n ( λ a i x + 1 ) \prod_{i=1}^n(\lambda^{a_i}x+1) i=1∏n(λaix+1)
分治乘起来,然后看第 k k k项是啥就行了。
代码:
#include<bits/stdc++.h>
#define ll long long
#define re register
#define cs const
namespace IO{
static cs int Rlen=1<<22|1;
static char buf[Rlen],*p1,*p2;
inline char gc(){return (p1==p2)&&(p2=(p1=buf)+fread(buf,1,Rlen,stdin),p1==p2)?EOF:*p1++;}
template<typename T>T get(){
char c;T num;
while(!isdigit(c=gc()));num=c^48;
while(isdigit(c=gc()))num=(num+(num<<2)<<1)+(c^48);
return num;
}inline int gi(){return get<int>();}
}
using namespace IO;
using std::cerr;
using std::cout;
cs double PI=acos(-1);
struct cp{
double x,y;
cp(){}cp(double _x,double _y=0):x(_x),y(_y){}
friend cp operator+(cs cp &a,cs cp &b){return cp(a.x+b.x,a.y+b.y);}
friend cp operator-(cs cp &a,cs cp &b){return cp(a.x-b.x,a.y-b.y);}
friend cp operator*(cs cp &a,cs cp &b){return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
cp conj()cs{return cp(x,-y);}
cp& operator+=(cs cp &b){return *this=*this+b;}
cp& operator-=(cs cp &b){return *this=*this-b;}
cp& operator*=(cs cp &b){return *this=*this*b;}
cp& operator/=(int b){x/=b,y/=b;return *this;}
};
cs int mod=99991;
inline int add(int a,int b){return a+b>=mod?a+b-mod:a+b;}
inline int dec(int a,int b){return a-b<0?a-b+mod:a-b;}
inline int mul(int a,int b){ll r=(ll)a*b;return r>=mod?r%mod:r;}
inline int po(int a,int b){
int r=1;for(;b;b>>=1,a=mul(a,a))
if(b&1)r=mul(r,a);return r;
}
inline void Inc(int &a,int b){a+=b-mod;a+=a>>31&mod;}
inline void Dec(int &a,int b){a-=b;a+=a>>31&mod;}
inline void Mul(int &a,int b){a=mul(a,b);}
inline void ex_gcd(int a,int b,int &x,int &y){
if(!b){x=1,y=0;return ;}ex_gcd(b,a%b,y,x);y-=a/b*x;
}
inline int inv(int a){a%=mod;
int y,x;ex_gcd(mod,a,y,x);
return x+(x>>31&mod);
}
cs int bit=17,SIZE=1<<bit|7;
int r[SIZE];cp *w[bit+1];
int fac[mod],ifc[mod];
void init_FFT(){
for(int re i=1;i<=bit;++i)w[i]=new cp[1<<i-1];
for(int re d=1;d<=bit;++d){
for(int re i=0;i<(1<<d-1);++i)
w[d][i]=cp(cos(2*PI*i/(1<<d)),sin(2*PI*i/(1<<d)));
}
fac[0]=1;for(int re i=1;i<mod;++i)fac[i]=mul(fac[i-1],i);
ifc[mod-1]=inv(fac[mod-1]);
for(int re i=mod-1;i;--i)ifc[i-1]=mul(ifc[i],i);
}
void FFT(cp *A,int len,int typ){
for(int re i=1;i<len;++i)if(i<r[i])std::swap(A[i],A[r[i]]);
for(int re i=1,d=1;i<len;i<<=1,++d)
for(int re j=0;j<len;j+=i<<1)
for(int re k=0;k<i;++k){
cp &t1=A[j+k],&t2=A[i+j+k],t=t2*w[d][k];
t2=t1-t,t1+=t;
}
if(typ==-1){
std::reverse(A+1,A+len);
for(int re i=0;i<len;++i)A[i]/=len;
}
}
void init_rev(int l){for(int re i=1;i<l;++i)r[i]=r[i>>1]>>1|((i&1)?l>>1:0);}
void mul(int *a,int n,int *b,int m,int *c){
static cp A[SIZE],B[SIZE];int deg=n+m-1,l=1;
while(l<deg)l<<=1;init_rev(l);
for(int re i=0;i<l;++i){
A[i]=cp(i<n?a[i]:0,0);
B[i]=cp(i<m?b[i]:0,0);
}
FFT(A,l,1),FFT(B,l,1);
for(int re i=0;i<l;++i)
A[i]*=B[i];FFT(A,l,-1);
for(int re i=0;i<deg;++i)
c[i]=(ll)(A[i].x+0.5)%mod;
}
inline int CC(int n,int m){return mul(fac[n],mul(ifc[m],ifc[n-m]));}
inline int C(int n,int m){
int res=1;
do{
int nn=n%mod,mm=m%mod;
if(nn<mm)return 0;
Mul(res,CC(nn,mm));
n/=mod,m/=mod;
}while(n&&m);
return res;
}
cs int N=1e5+7;
int n,k,C1,C2,a[N],b0,b1,pw[mod];
int poly[35][N],ans[N],*st[35],tp;
void solve(int l,int r,int *f){
if(l==r){
f[0]=1;
f[1]=pw[a[l]%(mod-1)];
return ;
}
int mid=l+r>>1,*L,*R;
L=st[tp--];solve(l,mid,L);
R=st[tp--];solve(mid+1,r,R);
mul(L,mid-l+2,R,r-mid+1,f);
st[++tp]=L,st[++tp]=R;
}
signed main(){
#ifdef zxyoi
freopen("nothing_to_see.in","r",stdin);
#endif
n=gi(),k=gi();init_FFT();
for(int re i=1;i<=n;++i)a[i]=gi(),++(a[i]&1?b1:b0);
int f0=gi(),f1=gi();
C1=mul(f0+f1,inv(4));C2=dec(f0,C1);
pw[0]=1;for(int re i=1;i<mod;++i)pw[i]=mul(pw[i-1],3);
for(int re i=0;i<35;++i)st[++tp]=poly[i];
solve(1,n,ans);int t1=ans[k],t2=0;
for(int re i=std::max(0,k-b1),li=std::min(b0,k);i<=li;++i)
(k^i)&1?Dec(t2,mul(C(b0,i),C(b1,k-i))):Inc(t2,mul(C(b0,i),C(b1,k-i)));
cout<<add(mul(t1,C1),mul(t2,C2));
return 0;
}