P1357 花园（状压dp + 矩阵快速幂）

传送门

此题是一个很好的题目，可以加深对矩阵快速幂的理解。若对基本的矩阵快速幂不太理解，需要先了解一下简单的比较裸一点的题目后再看本题目，更有利于理解。可以先学习一下 HDU_2842。本人写过这道题目的题解，但是并不详细，仅供参考我的题解HUD2842。

Solution

参考题解1
参考题解2

假设：长度m的花圃所有的状态为 [0,S],S=2m−1[0, S], S = 2^m - 1[0,S],S=2m−1。状态1表示C形花圃，0表示P形花圃
我们假设花园为链状的，则：状态转移方程为

ll t1=j>>1,t2=(j>>1)|(1<<(m-1));
if(num_1(t2)<=k)//1的个数小于等于kf[i][j]=f[i-1][t1]+f[i-1][t2];
elsef[i][j]=f[i-1][t1];

由于是环形的，我只需要由某种长度为m的合法状态的花园进行状态转移n次，然后保证最后一次的状态与初始状态相同即可。
假设状态转移矩阵 X 中某个点a[i, j] = 1,表示状态 i 可由状态 j 转移过来。
现在假设 [0,S][0, S][0,S]中每一个合法状态的矩阵为 A
(000i0...0)\begin{pmatrix} 0\\0\\0\\i\\0\\...\\0 \end{pmatrix}⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛000i0...0⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞
此矩阵共S + 1 行，i∈[0,S]i\in[0, S]i∈[0,S]，表示某一种长度为m的状态为i的花园。
设 F=Xn×AF = X ^ n \times AF=Xn×A。显然，F 为 1×(S+1)1 \times (S + 1)1×(S+1) 的矩阵。
则此种情况下的方案数为 F[i]F[i]F[i]。
由此，我只需要分别求出所有状态的F矩阵然后求和即可。

但是，有一种更为简便的方法:
5. P=XnP = X ^ nP=Xn , ans=∑i=0SP[i][i]，i为合法状态ans = \sum^S_{i = 0} P[i][i]， i为合法状态ans=∑i=0SP[i][i]，i为合法状态。
6. 仔细思考后发现，其实就是把所有的合法方案放在一个 (S+1)×(S+1)(S + 1) \times (S + 1)(S+1)×(S+1)的矩阵M中，第i 行为状态为 i 的 1×(S+1)1 \times (S + 1)1×(S+1)的矩阵
(也可看做单位矩阵，只不过最后求和只加上合法的方案即可, 即 P=Xn×EP = X^n \times EP=Xn×E。)。
8. 然后对于每行中的状态 i 其结果自然是第 i 列的末状态（由于是原型花园嘛，得保证初始状态和最终的状态相同）。所以最终的答案就是 P 矩阵的斜对角线了。

Code

#define mt(a,b) memset(a, b, sizeof a)
#define ll long long int
const ll mod = 1000000007;
using namespace std;
const int N = (1 << 5) + 5;
ll n;
int m, k, S;
int ok[N];//合法状态//矩阵
struct matrix
{ll a[N][N];matrix() { mt(a, 0); }matrix operator * (const matrix x) const{matrix ans;for (int i = 0; i <= S; i++)for (int j = 0; j <= S; j++)for (int k = 0; k <= S; k++)ans.a[i][j] = (ans.a[i][j] + x.a[i][k] * a[k][j]) % mod;//?return ans;}
};ll power(ll n)
{matrix e, dp;for (int i = 0; i <= S; i++)//单位矩阵初始化e.a[i][i] = 1;for (int i = 0; i <= S; i++)// 求转移矩阵if (ok[i]){int t1 = i >> 1;int t2 = t1 | 1 << (m - 1);dp.a[i][t1] = 1;if (ok[t2])dp.a[i][t2] = 1;}for (; n; n >>= 1){if (n & 1) e = e * dp;dp = dp * dp;}ll ans = 0;for (int i = 0; i <= S; i++)if (ok[i])ans = (ans + e.a[i][i]) % mod;return ans;
}int main()
{cin >> n >> m >> k;S = (1 << m) - 1;//状态for (int i = 0; i <= S; i++){int cnt = 0, x = i;while (x) x -= x & -x, cnt++;if (cnt <= k) ok[i] = 1;}cout << power(n) << endl;return 0;
}