题目地址：洛谷：【P1357】花园 – 洛谷

题目描述

小L有一座环形花园，沿花园的顺时针方向，他把各个花圃编号为1~N(2<=N<=10^15)。他的环形花园每天都会换一个新花样，但他的花园都不外乎一个规则，任意相邻M(2<=M<=5,M<=N)个花圃中有不超过K(1<=K<M)个C形的花圃，其余花圃均为P形的花圃。
例如，N=10,M=5,K=3。则
CCPCPPPPCC 是一种不符合规则的花圃；
CCPPPPCPCP 是一种符合规则的花圃。
请帮小L求出符合规则的花园种数Mod 1000000007

由于请您编写一个程序解决此题。

输入输出格式

输入格式：
一行，三个数N,M,K。

输出格式：
花园种数Mod 1000000007

输入输出样例

输入样例#1：

10 5 3

输出样例#1：

458

输入样例#2：

6 2 1

输出样例#2：

18

说明

【数据规模】
40%的数据中，N<=20；
60%的数据中，M=2；
80%的数据中，N<=10^5。
100%的数据中，N<=10^15。

题解

我们看一下M=2怎么做。我们枚举后面2个花圃的状态，一共有4种：CC、CP、PC、PP，设计状态dp[i][S]表示到第i个且末尾2个花圃的状态为S的方案数，然后枚举后面接上去哪个花圃来转移，如CC可以转移到CC和CP。如果M更大，到5了，就可以考虑状压这个状态来转移。转移同上，只是要把操作换为位运算。这样直接做的复杂度是[eq]O(n \cdot 2^{m})[/eq]的。
当n变得无法承受的时候，我们需要一个[eq]O(\log n)[/eq]的算法。这个时候我们想到了快速幂。我们知道，从一个状态向下一个状态的转移是一个固定不变的线性变换，既然如此，我们构建转移矩阵，进行矩阵快速幂计算转移矩阵的n次方即可。复杂度是[eq]O(\log n)[/eq]的（忽略矩阵乘法带来的开销）。
对于构造转移矩阵的方法，实际上矩阵的第i行表示下一层的i这个元素由这一层的元素乘上怎么样的系数加和得来的，因此对于i \rightarrow i'这个转移，我们直接令trans_{i', i} = 1即可。由于初始的时候是一个单位矩阵，也可以不乘这个初始矩阵。
为什么不是n-m次幂？注意我们DP状态的定义，答案在dp[n][S]中，至于当i小于m的时候，实际上S多出来那部分的意义是在枚举这个序列末尾的情况，因为花园是个环嘛。答案计算的是[eq]\sum_S dp[n][S] = trans^n_{S, S}[/eq]，是因为初始的时候是个列向量，是一排1，扩展成方阵以后就是单位矩阵了，最后肯定要按照列向量的意义计算，就取主对角线上的数值。

代码

// Code by KSkun, 2018/4
#include <cstdio>
#include <cstring>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1;
    register char c = fgc();
    while(c < '0' || c > '9') {
        if(c == '-') neg = -1;
        c = fgc();
    }
    while(c >= '0' && c <= '9') {
        res = res * 10 + c - '0';
        c = fgc();
    }
    return res * neg;
}

const int MAXN = 50, MO = 1e9 + 7;

LL n, m, k, cnt[MAXN];

struct Matrix {
    LL a[MAXN][MAXN];
    Matrix() {
        memset(a, 0, sizeof(a));
    }
    inline Matrix operator*(const Matrix &rhs) const {
        Matrix res;
        for(int i = 0; i < 1 << m; i++) {
            for(int j = 0; j < 1 << m; j++) {
                for(int k = 0; k < 1 << m; k++) {
                    res.a[i][j] = (res.a[i][j] + a[i][k] * rhs.a[k][j] % MO) % MO;
                }
            }
        }
        return res;
    }
    inline Matrix& operator*=(const Matrix &x) {
        return *this = *this * x;
    }
};

inline Matrix fpow(Matrix mat, LL k) {
    Matrix res;
    for(int i = 0; i < 1 << m; i++) {
        res.a[i][i] = 1;
    }
    while(k) {
        if(k & 1) res *= mat;
        mat *= mat;
        k >>= 1;
    }
    return res;
}

int main() {
    n = readint(); m = readint(); k = readint();
    for(int i = 0; i < 1 << m; i++) {
        cnt[i] = cnt[i >> 1] + (i & 1);
    }
    Matrix trans;
    for(int i = 0; i < 1 << m; i++) {
        if(cnt[i] > k) continue;
        trans.a[i >> 1][i] = 1;
        trans.a[i >> 1 | (1 << (m - 1))][i] = 1;
    }
    trans = fpow(trans, n);
    LL res = 0;
    for(int i = 0; i < 1 << m; i++) {
        res = (res + trans.a[i][i]) % MO;
    }
    printf("%lld", res);
    return 0;
}