题解：洛谷P1357 花园
Description
小 L 有一座环形花园，沿花园的顺时针方向，他把各个花圃编号为 (1∼n)。花园 (1) 和 (n) 是相邻的。
他的环形花园每天都会换一个新花样，但他的花园都不外乎一个规则：任意相邻 (m) 个花圃中都只有不超过 (k) 个 C 形的花圃，其余花圃均为 P 形的花圃。
例如，若 (n=10 , m=5 , k=3) ，则
· CCPCPPPPCC 是一种不符合规则的花圃。
· CCPPPPCPCP 是一种符合规则的花圃。
请帮小 L 求出符合规则的花园种数对 (1e9 + 7) 取模的结果。
Hints
对于(40\%)的数据，(n leq 20)
对于(60\%)的数据，(m = 2)
对于(80\%)的数据，(n leq 10 ^5)
对于(100\%)的数据，(2 leq n leq 10^{15}, 2 leq m leq min(n,5), 1 leq k < m)
Algorithm
单看Description可能以为这是个纯数学题，然而上手就能发现这玩意不太可做。
其实，这题其实已经把算法写在脸上了。
首先可以观察到 (m) 的值很小，因此显然需要状态压缩。
接着容易发现 (n) 的量级对于状态压缩来说大得离谱，因此显然是矩阵加速递推做的。
不管计数DP有没有最优子结构，姑且也算DP吧，那么这题就是个矩阵加速的计数类环形状压DP。
接下来把递推式搞出来就完了。
先给出一些约定：
我们考虑把放 C 当做 1 ，放 P 当做 0 ，将相邻 (m) 位的摆花方法当成一个二进制数压缩。
因为转移时左移操作会带来麻烦，所以把二进数的位置与原串摆法相反。
所有位运算符号均采用 C++ 中的样式。
Step 1
先考虑普通的状压DP：
我们用 dp[i - 1][j] 表示考虑到第 (i - 1) 个位置，后 (m) 位顺次摆花的状态为 (j) 时的方案数。
如果把第 (i) 位设为 0 ，那么显然有新状态 j' = j >> 1 ，
如果把第 (i) 位设为 1 ，那么显然有新状态 j' = (j >> 1) | (1 << (m - 1)) ，
那么 dp[i][j] 显然就是dp[i - 1][j >> 1] 与 dp[i - 1][(j >> 1) | (1 << (m - 1))]求和得到的。
注意后者不一定合法，需要 __builtin_popcount() 判断一下。 //这是个统计二进制数中 1 的个数的函数，奇快无比。
最后，其实这题环的特征其实并不明显，
对于前 (m) 位状态为 t 的情况来说，我们只需要递推完事之后将 dp[n][t] 计入答案即可。
当然，容易注意到 dp[i][] 只与 dp[i - 1][] 有关，可以滚动优化一下：
细节的小技巧就不细说了，容易写出代码：
int ans = 0, ful = 1 << m, p = m - 1;
for(int t = 0; t != ful; ++t)
{
if(__builtin_popcount(t) > k) continue;
memset(dp, 0, sizeof(dp));
dp[0][t] = 1;
for(int i = 1; i <= n; ++i)
for(int j = 0; j <= ful; ++j)
{
dp[i & 1][j] = dp[i - 1 & 1][j >> 1];
if(__builtin_popcount((j >> 1) | (1 << p)) <= k)
dp[i & 1][j] += dp[i - 1 & 1][(j >> 1) | (1 << p)];
}
ans += dp[n & 1][t];
}
复杂度显然是 (O(ncdot 4^m)) 的。
Step 2
接下来就需要构造递推矩阵了。
以 (m = 2, k = 1) 为例，我们有
dp[n][00] = dp[n - 1][00] + dp[n - 1][01];
dp[n][01] = dp[n - 1][10];
dp[n][10] = dp[n - 1][00] + dp[n - 1][01];
dp[n][11] = 0;
从而有转移方程：
（(LaTeX) 的矩阵也太难写了吧……）
类似地，对于具体的 (m,k) ，
我们可以使用和上面的递推程序相同的方法构造出类似的矩阵。
最后同样来考虑环的问题。
根据 Step 1 中的分析，答案实际上就是(sum dp[ n ][ j ])，在我们这里又有
（特别应当注意，右侧列向量并非全 1 向量）
如果把右侧列向量表示出来乘上去再求和，自然能做，不过有些麻烦了。
这里实际上有个绝妙的性质：
考虑构造出来的矩阵 (A) , (A_{i,j}) 实际上表示了状态 (j) 可以向状态 (i) 转移。
而由于环首尾相连的特性，对于一个合法的状态（也就是右侧列向量中对应值为 1 的状态），
在 (n) 次转移后，该状态一定可以转移回原状态，也即是说 (A_{i,i}) 一定非零。
而在矩阵做了 (n) 次幂乘之后， (A_{i,i}) 也就是起始状态为 (i) 的方案数。
因此答案实际上即为 (sum A_{i,i})。
(sum A_{i,i} = sum dp[ n ][ j ]) 的代数证明留给读者。（这能证吗？不管了，我们读者肯定会）
最后写写代码：
#include<bits/stdc++.h>
using namespace std;
#define cnt __builtin_popcount
#define rep(i, n) for(int i = 0; i != n; ++i)
typedef long long ll;
const int P = 1e9 + 7;
ll n, m, k;
class Matrix {
private:
vector<vector<int> > num;
public:
Matrix() {}
Matrix(int n) {
num.resize(n);
rep(i, n) num[i].resize(n, 0);
}
int &operator () (int i, int j) {
return num[i][j];
}
inline void print()
{
int len = num.size();
for(int i = 0; i != len; ++i)
{
printf("%d", num[i][0]);
for(int j = 1; j != len; ++j)
printf(" %d", num[i][j]);
puts("");
}
}
inline int solve()
{
int ans = 0, len = num.size();
for(int i = 0; i != len; ++i)
(ans += num[i][i]) %= P;
return ans;
}
friend Matrix operator * (Matrix a, Matrix b)
{
int len = a.num.size();
Matrix ret(len);
rep(k, len) rep(i, len) rep(j, len) {
(ret(i, j) += (1LL * a(i, k) * b(k, j)) % P) %= P;
}
return ret;
}
friend Matrix operator ^ (Matrix a, ll n)
{
int len = a.num.size();
Matrix ret(len);
rep(i, len) ret(i, i) = 1;
while(n)
{
if(n & 1) ret = ret * a;
a = a * a, n >>= 1;
}
return ret;
}
};
int main()
{
cin >> n >> m >> k;
int ful = 1 << m, p = m - 1;
Matrix cal(ful);
for(int t = 0; t != ful; ++t)
{
if(cnt(t) > k) continue;
cal(t >> 1, t) = 1;
int sta = (t >> 1) | (1 << p);
if(cnt(sta) <= k) cal(sta, t) = 1;
}
cout << (cal ^ n).solve() << endl;
return 0;
}
这个代码非常之丑。最开始设计模板类翻车了，没有重构而是左改右改的……
总之不是好代码，不过反正是这个意思。
复杂度 (O(log n cdot 8 ^ m)) 很容易过。