【LOJ6183】看无可看（特征方程）（Lucas定理）（分治FFT）

传送门

---

题解：

还行的一道题。其实我主要是来复习FFT怎么写的，突然发现自己不会FFT只会NTT

首先，如果注意到$f$是一个二阶线性递推数列，这道题就已经做完了。

考虑利用特征方程解出$f$的通项公式。

该数列特征方程为${\lambda }^2-2\lambda -3=0$，两个根为${\lambda }_1=3,{\lambda }_2=-1$，设$C_1,C_2$为待定系数。则其通项公式为$f_n=C_1{\lambda }_1^n+C_2{\lambda }_2^n$

由于题目给出了$f_0,f_1$，我们可以把$C_1,C_2$解出来。

考虑题目要求和的东西，我们可以对于${\lambda }_1$和${\lambda }_2$两个部分分别计算。

其中由于${\lambda }_2=-1$，可以直接用组合数学算，注意由于模数的关系要用Lucas定理（结果数据最大的点$n=5e4$，完全不用Lucas）

于是现在问题我们要求出选$k$个元素，设$\sum$为其和，对于${\lambda }^{\Sigma }$求和。

。。。。

非常显然直接构造多项式$$\prod _{i=1}^n({\lambda }^{a_i}x+1)$$

分治乘起来，然后看第$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;
}