[状压 + 矩阵乘法] HDU 4332 Constructing Chimney

题目大意：
我们计划建一个大烟囱。烟囱的每一节是一个 3*3 的正方形，中间是空的，如图所示。
现在我们只有许多 1*1*2 的块。我们建造的烟囱的高度为 N。现在请计算有多少种不同
的建造烟囱的方法。答案模 10^9 +7。

据说是很经典的一道状压矩乘 dp 。

首先可以对前三层裸状压 dp 求出状态到状态的转移。可以发现，有用的状态只有 70 个。然后写到矩乘里面即可。

记得开 long long，还有 Case 的首字母要大写！！！

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>

#define mo 1000000007LL

int t, n, m;
int num[71];
int a[256][256], b[256][256], c[256][256], (* f)[256];
bool v1[256], v2[256], v3[256], * vst;
struct ma
{
	long long a[71][71];
	long long * operator [](const int & b) {return a[b];}
}	p, q, x, y;

ma operator *(ma & a, ma & b)
{
	ma c = ma();
	for (int i = 1; i <= 70; ++ i)
		for (int j = 1; j <= 70; ++ j)
			for (int k = 1; k <= 70; ++ k)
				c[i][j] = (c[i][j] + a[i][k] * b[k][j]) % mo;
	return c;
}

ma operator /(ma & a, ma & b)
{
	ma c = ma();
	for (int j = 1; j <= 70; ++ j)
		for (int k = 1; k <= 70; ++ k)
			c[1][j] = (c[1][j] + a[1][k] * b[k][j]) % mo;
	return c;
}

void show(int t, int k)
{
	if (t) show(t >> 1, k);
	fprintf(stderr, "%d", t & 1);
	if (t == k) fprintf(stderr, "\n");
}

void dfs(int r, int s, int t, int l)
{
	if (r == (1 << 8) - 1) return (void)(++ f[s][t], vst[t] = 1);
	for (int i = l; i < 7; ++ i)
		if (! (r >> i & 1))
		{
			dfs(r | 1 << i, s, t | 1 << i, i + 1);
			if (! (r >> i + 1 & 1)) dfs(r | 1 << i | 1 << i + 1, s, t, i + 2);
		}
	if (! (r >> 7 & 1))
	{
		dfs(r | 1 << 7, s, t | 1 << 7, 8);
		if (! (r & 1)) dfs(r | 1 << 7 | 1, s, t, 9);
	}
}

void prepare()
{
	f = a, vst = v1, dfs(0, 0, 0, 0);
	f = b, vst = v2;
	for (int i = 0; i < 1 << 8; ++ i)
		if (v1[i]) dfs(i, i, 0, 0);
	f = c, vst = v3;
	for (int i = 0; i < 1 << 8; ++ i)
		if (v2[i]) dfs(i, i, 0, 0), num[++ m] = i;
	for (int i = 1; i <= 70; ++ i)
	{
		p[1][i] = a[0][num[i]];
		for (int j = 1; j <= 70; ++ j) q[j][i] = c[num[i]][num[j]];
	}
}

int main()
{
	freopen("a.in", "r", stdin);
	freopen("a.out", "w", stdout);
	
	prepare(), scanf("%d", & t);
	for (int i = 1; i <= t; ++ i)
	{
		scanf("%d", & n), x = p, y = q;
		for (-- n; n; n >>= 1, y = y * y)
			if (n & 1) x = x / y;
		printf("Case %d: %I64d\n", i, x[1][1]);
	}
	
	return 0;
}