矩阵快速幂优化状压DP--骨牌

题目：
有一个$H\times W$ 的棋盘。 她想知道，用 $1\times 2$ 的骨牌覆盖这个棋盘，有多少种不同的方案。 一个合法的方案满足所有骨牌的边都与棋盘的边平行，骨牌之间没有重叠，并且骨牌的所有部分都在棋盘内。 两种方案不同当且仅当在一种方案中，一个格子被骨牌覆盖，在另一种方案中，这个格子没有被骨牌覆盖。
$H\le 5,W\le 10^{18}$

solution:
（这是一道巨巨巨难的题，虽然代码不长但是非常难想

第一眼看到这道题的时候想到了以前的一道类似的题，也是矩阵优化$dp$，但不同的是那道必须填满，而这道不需要填满而且两种方案不同仅当有格子的覆盖情况不同，这样的话，就会出现多种摆放骨牌方法对应同一种覆盖情况，就不能用普通的状压做了

首先，对这种多种对应同一种的情况，考虑状态再次压缩，设$f[i][S]$表示前$i$列，$S$表示一个凸出的插头形状的集合，这个集合中的所有插头形状都可以对应相同的一种覆盖方式，比如说$10110$和$10000$就可以放到一个集合里

这样设状态看起来是$2^{2^H}$的，只能通过$H\le 3$的点，但通过搜索可以发现合法状态很少，最多只有$90$多种，对于这种看起来很多其实很少的状态，可以先搜索预处理出所有合法的集合记到一个数组$g[i]$里面，然后通过排序二分等操作就可以预处理出矩阵，然后的步骤就很简单了

要特别说一下这里的$dfs$过程，比较神奇，首先$dfs(intx)$表示当前已经搜到了$x$集合，要继续搜和$x$覆盖情况不一样的，这里采用的方法是，$2^H$循环一次所有的插头形状假设为$i$，然后看$i$和$x$集合内的插头的关系，$if((i\&j)==j)$，就说明$ixorj$是一种不同的情况，然后就把可和$ixorj$覆盖情况相同的都处理出来，这里的方法是写一个递归函数，看有没有连续的两个$1$，有的话说明可以用竖块替代，也是一个合法的。

具体看代码吧，挺难理解的$qwq$

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define maxn 100
#define int long long
using namespace std;
int n,m,mod,ed,g[maxn],tot,ans;

inline int rd(){
	int x=0,f=1;char c=' ';
	while(c<'0' || c>'9') f=c=='-'?-1:1,c=getchar();
	while(c<='9' && c>='0') x=x*10+c-'0',c=getchar();
	return x*f;
}

struct Mat{
	int a[maxn][maxn];
	Mat(){memset(a,0,sizeof a);}
	Mat operator *(const Mat &x) const{
		Mat ret;
		for(int k=0;k<tot;k++)
			for(int i=0;i<tot;i++)
				for(int j=0;j<tot;j++)
				(ret.a[i][j]+=a[i][k]*x.a[k][j]%mod)%=mod;
		return ret;
	}
}f,s;

inline Mat qpow(Mat x,int k){
	Mat ret;
	ret.a[0][0]=1;
	while(k){
		if(k&1) ret=ret*x;
		x=x*x; k>>=1;
	} return ret;
}

inline void change(int &k,int x){
	k|=(1<<x);
	for(int i=0;i<n-1;i++)
		if((x>>i&3)==3) change(k,x^(3<<i));
}

inline void dfs(int x){
	int i,j,y,l,r,mid;
	l=0,r=tot;
	while(l<r){
		mid=(l+r)>>1;
		if(g[mid]<x) l=mid+1;
		else r=mid-1;
	}r++;
	if(g[r]==x) return;
	g[tot++]=x; sort(g,g+tot);
	for(i=0;i<ed;i++){
		for(y=j=0;j<ed;j++)
			if((x>>j&1) && (i&j)==j)
				change(y,i^j);
		if(y) dfs(y);
	}
}

signed main(){
	n=rd(); m=rd(); mod=rd(); ed=1<<n;
	dfs(1);
	for(int cur=0;cur<tot;cur++){
		int x=g[cur];
		for(int i=0;i<ed;i++){
			int y=0;
			for(int j=0;j<ed;j++)
				if((x>>j&1) && (i|j)==i) change(y,i^j);
			if(y){
				int l=0,r=tot,mid;
				while(l<r){
					mid=(l+r)>>1;
					if(g[mid]<y) l=mid+1;
					else r=mid-1;
				}r++;
				f.a[cur][r]++;
			}
		}
	}
	s=qpow(f,m);
	for(int i=0;i<tot;i++)
		if(g[i]&1) (ans+=s.a[0][i])%=mod;
	printf("%lld\n",ans);
	return 0;
}