Loj#6183. 看无可看

Loj#6183. 看无可看

题目描述
首先用特征根求出通项公式(A_n=pcdot 3^n+qcdot(-1)^n)。通过给定的(f_0,f_1)可以解出(p,q)。

然后我们要求的就是(sum_{|s'|=k}Pi_{xin s'}a_x)。这就是个背包。

考虑它的生成函数就是(Pi(1+a_ix))。用分治(FFT)求解。

代码：

#include<bits/stdc++.h>
#define ll long long
#define N 100005

using namespace std;
inline int Get() {int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch=getchar();}while('0'<=ch&&ch<='9') {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}return x*f;}

const ll mod=99991;
int Mod(int a) {return a<mod?a:a-mod;}
struct matrix {
	ll a[2][2];
	void Init() {memset(a,0,sizeof(a));}
}g,f;
matrix operator *(const matrix &x,const matrix &y) {
	static matrix tem;
	tem.Init();
	tem.a[0][0]=(x.a[0][0]*y.a[0][0]+x.a[0][1]*y.a[1][0])%mod;
	tem.a[0][1]=(x.a[0][0]*y.a[0][1]+x.a[0][1]*y.a[1][1])%mod;
	tem.a[1][0]=(x.a[1][0]*y.a[0][0]+x.a[1][1]*y.a[1][0])%mod;
	tem.a[1][1]=(x.a[1][0]*y.a[0][1]+x.a[1][1]*y.a[1][1])%mod;
	return tem;
}
matrix operator +(const matrix &x,const matrix &y) {
	static matrix tem;
	tem.Init();
	tem.a[0][0]=Mod(x.a[0][0]+y.a[0][0]);
	tem.a[0][1]=Mod(x.a[0][1]+y.a[0][1]);
	tem.a[1][0]=Mod(x.a[1][0]+y.a[1][0]);
	tem.a[1][1]=Mod(x.a[1][1]+y.a[1][1]);
	return tem;
}

matrix ksm(matrix g,int x) {
	static matrix ans;
	ans.Init();
	for(int i=0;i<2;i++) ans.a[i][i]=1;
	for(;x;x>>=1,g=g*g) if(x&1) ans=ans*g;
	return ans;
}

int n,k;
int a[N];
bool vis[N];
matrix G[N];
ll ans;
struct Com {
	double a,r;
	Com() {a=0,r=0;}
	Com(double x,double y) {a=x,r=y;}
};

Com operator +(const Com &x,const Com &y) {return Com(x.a+y.a,x.r+y.r);}
Com operator -(const Com &x,const Com &y) {return Com(x.a-y.a,x.r-y.r);}
Com operator *(const Com &x,const Com &y) {return Com(x.a*y.a-x.r*y.r,x.a*y.r+x.r*y.a);}
Com operator /(const Com &x,const double &y) {return Com(x.a/y,x.r/y);}
const double pi=acos(-1);
Com W[20][N<<2],M[20][N<<2];

void FFT(Com *a,int d,int flag) {
	static int rev[N<<2];
	int n=1<<d;
	for(int i=0;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<d-1);
	for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int s=1;s<=d;s++) {
		int len=1<<s,mid=len>>1;
		for(int i=0;i<n;i+=len) {
			for(int j=0;j<mid;j++) {
				Com u=a[i+j],v=(flag==1)?a[i+j+mid]*W[d][j*n/len]:a[i+j+mid]*M[d][j*n/len];
				a[i+j]=u+v;
				a[i+j+mid]=u-v;
			}
		}
	}
	if(flag==-1) for(int i=0;i<n;i++) a[i]=a[i]/n;
}

ll ksm(ll t,ll x) {
	ll ans=1;
	for(;x;x>>=1,t=t*t%mod)
		if(x&1) ans=ans*t%mod;
	return ans;
}

ll F[N<<2];
Com L[N<<2],R[N<<2];
Com A[N<<2];
Com tem[N<<2];
ll t;

void solve(int l,int r) {
	if(l==r) {
		F[l]=ksm(t,a[l]);
		return ;
	}
	int mid=l+r>>1;
	solve(l,mid),solve(mid+1,r);
	int d=ceil(log2(r-l+2));
	for(int k=0;k<1<<d;k++)
		L[k]=R[k]=Com(0,0);
	for(int i=l;i<=mid;i++) {
		L[i-l+1]=Com(F[i],0);
	}
	for(int i=mid+1;i<=r;i++) {
		R[i-mid]=Com(F[i],0);
	}
	L[0]=R[0]=Com(1,0);
	int len=1<<d;
	FFT(L,d,1);
	FFT(R,d,1);
	for(int k=0;k<1<<d;k++) {
		tem[k]=(L[k]*R[k]);
	}
	FFT(tem,d,-1);
	for(int k=l;k<=r;k++) {
		F[k]=((ll)(tem[k-l+1].a+0.5))%mod;
	}
}

int main() {
	
	n=Get(),k=Get();
	g.a[0][0]=0,g.a[0][1]=3;
	g.a[1][0]=1,g.a[1][1]=2;
	int d=ceil(log2(n*2));
	for(int i=1;i<=d;i++)
		for(int j=0;j<1<<i;j++)
			W[i][j]=Com(cos(2*j*pi/(1<<i)),sin(2*j*pi/(1<<i))),
			M[i][j]=Com(cos(-2*j*pi/(1<<i)),sin(-2*j*pi/(1<<i)));
	for(int i=1;i<=n;i++) {
		a[i]=Get();
		G[i]=ksm(g,a[i]);
	}
	f.a[0][0]=Get(),f.a[0][1]=Get();
	ll q=(f.a[0][0]+f.a[0][1])*ksm(4,mod-2)%mod,p=(3*f.a[0][0]-f.a[0][1]+mod)*ksm(4,mod-2)%mod;

	ll ans=0;
	t=3;
	solve(1,n);
	ans=q*F[k]%mod;
	t=mod-1;
	solve(1,n);
	(ans+=p*F[k])%=mod;
	cout<<ans;
	return 0;
}