2018.01.02【CQOI2018】【BZOJ5298】【洛谷P4456】交错序列（矩阵快速幂）（轻微卡常）

BZOJ传送门

洛谷传送门

---

解析：

本来还可以继续把矩阵乘法一分为八卡卡常数的，今天就算了，CQOI2018真的做的我没心情卡了。。。

不知道出题人是哪根筋搭错了，5道版题还有两道卡常。。。

我只能说，这届CQOI出题人的水平估计都没有选手高（此话不负任何法律或人身责任（滑稽）

思路：

考虑原来要求的式子：$x^a{y}^b$怎么化简。

$$\begin{array}{rl}{x}^a{y}^b& =(n-x{)}^a{y}^b\\ & =\sum _{i=0}^a{C}_a^i(-1{)}^{a-i}{n}^i{y}^{a+b-i}\end{array}$$

显然前面这个东西${\sum }_{i=0}^a{C}_a^i(-1{)}^{a-i}{n}^i$可以很轻易的算出。

那么我们要求的就是所有长度为$n$的交错序列中，$y^{a+b-i}$的和。

考虑DP，设$f[k][i][1/0]$表示所有长度为$k$的，结尾是$0/1$的交错序列中，$1$的个数，即$y$的$i$次方之和。

显然，我们可以得到如下转移：
$$\begin{array}{rl}f[k][i][0]& =f[k-1][i][0]+f[k-1][i][1]\\ f[k][i][1]& =\sum _{j=0}^i{C}_i^{j}f[k-1][j][0]\end{array}$$

其中第二个式子由二项式定理得到，考虑所有$y$加一之后的二项式展开。

这个显然是矩阵乘法的形式，构造矩阵然后快速幂转移就好了。

卡常：

当然你要是认为这道题就完了也可以，慢慢卡常吧。

首先上矩阵快速幂卡常基础操作：除去$0$，优化循环结构。

然而这道题还有一个令人蛋疼的性质。

大矩阵是由四个小下三角矩阵拼起来的。。

所以乘法我们只需要对于每个下三角矩阵做一遍就行了。

当然，你也可以将矩阵一分为八，我已经不想写了。。。

---

代码：

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define re register
#define gc get_char
#define pc put_char
#define cs const

cs int N=190;
int mod;
inline int add(int a,int b){
	a+=b;return a>=mod?a-mod:a;
}

int siz1,siz2;
struct matrix{
	int a[N][N];
	matrix(){memset(a,0,sizeof a);}
	inline friend matrix operator*(cs matrix &A,cs matrix &B){
        matrix C;
		for(int re i=0;i<siz1;++i)
		for(int re j=0;j<siz2;++j)if(A.a[i][j])
		for(int re k=0;k<=i;++k)
		C.a[i][k]=add(C.a[i][k],1ll*A.a[i][j]*B.a[j][k]%mod);
		for(int re i=0;i<siz1;++i)
		for(int re j=0;j<siz2;++j)if(A.a[i][j])
		for(int re k=siz1;k<=siz1+i;++k)
		C.a[i][k]=add(C.a[i][k],1ll*A.a[i][j]*B.a[j][k]%mod);
		for(int re i=siz1;i<siz2;++i)
		for(int re j=0;j<siz2;++j)if(A.a[i][j])
		for(int re k=0;k<=i-siz1;++k)
		C.a[i][k]=add(C.a[i][k],1ll*A.a[i][j]*B.a[j][k]%mod);
		for(int re i=siz1;i<siz2;++i)
		for(int re j=0;j<siz2;++j)if(A.a[i][j])
		for(int re k=siz1;k<=i;++k)
		C.a[i][k]=add(C.a[i][k],1ll*A.a[i][j]*B.a[j][k]%mod);
		return C;
	}	
}mat;

inline matrix quickpow(matrix a,int b){
	matrix res;
	for(int re i=0;i<siz2;++i)res.a[i][i]=1;
	for(;b;b>>=1,a=a*a)if(b&1)res=res*a;
	return res;
}

int fac[100],inv[100],ifac[100];
inline int C(int n,int m){
	return (ll)fac[n]*ifac[n-m]%mod*ifac[m]%mod;
}
 
int n,a,b,ans;
signed main(){
	scanf("%d%d%d%d",&n,&a,&b,&mod);
	fac[0]=ifac[0]=fac[1]=ifac[1]=inv[1]=inv[0]=1;
	for(int re i=2;i<=a+b;++i){
		fac[i]=(ll)i*fac[i-1]%mod;
		inv[i]=(ll)(mod-mod/i)*inv[mod%i]%mod;
		ifac[i]=(ll)ifac[i-1]*inv[i]%mod;
	}
	siz1=a+b+1,siz2=siz1<<1;
	for(int re i=0;i<siz1;++i){
		mat.a[i][i]=mat.a[i][siz1+i]=1;
		for(int re j=0;j<=i;++j)mat.a[siz1+i][j]=C(i,j);
	}
	mat=quickpow(mat,n);
	int pown=1;
	for(int re i=0;i<=a;++i){
		ans+=(((a-i)&1)?-1ll:1ll)*C(a,i)*pown%mod*add(mat.a[a+b-i][0],mat.a[a+b-i+siz1][0])%mod;
		ans%=mod;
		pown=(ll)pown*n%mod;
	}
	cout<<(ans%mod+mod)%mod;
	return 0;
}