洛谷题单指南-状态压缩动态规划-P1357 花园

原题链接：https://www.luogu.com.cn/problem/P1357

题意解读：长度为n的0/1环形数组，连续m长度的子串中1的数量不超过k个，求一共有多少种可能的数组方案，结果对10^9+7取模。

解题思路：

1、可行解

由于m不超5，可以对连续m个子串的状态进行压缩，压缩成二进制整数。

设f[i][j]表示前i个元素，从i-m+1~i的子串状态是j的所有方案数。

不难发现，前i-1个元素的最后m个元素与前i个元素的最后m个元素，有m-1个元素是重叠的，

设前i-1个元素的最后m个元素状态是k，前i个元素的最后m个元素状态是j，k和j的关系有两种可能：

第一种，k等于j的前m-1个二进制位然后在最前面加一个0，k = j >> 1；

第二种，k等于j的前m-1个二进制位然后在最前面加一个1，k = j >> 1 | 1 << m - 1；

因此，可以得到递推关系：

f[i][j] = f[i-1][j >> 1] + f[i-1][j >> 1 | 1 << m - 1]

由于是环形数组，要处理一圈，从第0个到第n个即可，

对于所有状态s，进行初始化f[0][s] = 1，然后执行dp过程，结果就是所有f[n][s]的累加。

对结果的解释：即对某一个初始状态s，执行一圈dp，回到起始位置时状态仍然是s，这就是最后m个元素状态为s的方案数，枚举所有的s累加结果即为所有的方案数。

注意：考虑所有状态时，要判断是否合法，也就是其中1个个数是否不超过k，可以提前预处理。

80分代码：

#include <bits/stdc++.h>
using namespace std;

const int N = 100005, MOD = 1e9 + 7;
int f[N][1 << 5]; //f[i][j]表示前i个元素，从i-m+1~i的子串状态是j的所有方案数
int ones[1 << 5]; //二进制有多少个1
int n, m, k, ans;

void init()
{
    for(int i = 0; i < 1 << m; i++)
    {
        int cnt = 0;
        for(int j = 0; j < m; j++)
            if(i >> j & 1) 
                cnt++;
        ones[i] = cnt;
    }
}

int main()
{
    cin >> n >> m >> k;
    init();
    for(int s = 0; s < 1 << m; s++)
    {
        if(ones[s] > k) continue;
        memset(f, 0, sizeof(f));
        f[0][s] = 1;
        for(int i = 1; i <= n; i++)
        {
            for(int j = 0; j < 1 << m; j++)
            {
                int k1 = j >> 1;
                int k2 = j >> 1 | 1 << m - 1;
                if(ones[k1] <= k) f[i][j] = (f[i][j] + f[i - 1][k1]) % MOD;
                if(ones[k2] <= k) f[i][j] = (f[i][j] + f[i - 1][k2]) % MOD;
            }
        }
        ans = (ans + f[n][s]) % MOD;
    }
    cout << ans;
    return 0;
}

2、矩阵乘法优化

由于n最大到10^15，必须用logn的算法才可能通过，可以想到将递推优化成矩阵乘法。

来分析一下

f[i][j] = f[i-1][j >> 1] + f[i-1][j >> 1 | 1 << m - 1]

这里j >> 1，j >> 1 | 1 << m - 1都是可以转移到j的合法状态，更一般的，我们可以枚举所有状态k，看k是否能转移到j

设q[k][j]=1表示状态k可以转移到j，0表示不可以，那么上面递推式可以转换为：

对于这段代码：

由于要枚举所有的状态j，设size = 1<<m可以直接定义1*size矩阵：

p[i] = [ f[i][0], f[i][1]……f[i][(1<<m)-1]] ，而q是size*size矩阵

因此有p[i] = p[i-1] * q

再来看上一层代码：

要从1递归到n，因此有p[n] = p[0] * q^n

再看上一层代码，由于会枚举所有状态s，且初始化f[0][s]=1

这样，就是对于所有的状态s，依次计算：

p[n] = [ f[n][0], f[n][1]……f[n][(1<<m)-1]] = [f[0][0], f[0][1]…f[0][s]…f[0][(1<<m)-1]] * q^n

然后将每次计算之后的矩阵p[n]的第0行第s列的数值累加到答案即可。

对于矩阵q可以预计算，q^n可以用矩阵快速幂。

100分代码：

#include <bits/stdc++.h>
using namespace std;

const int M = 1 << 5, MOD = 1e9 + 7;
long long n, m, k, ans;
int ones[M]; //二进制有多少个1
struct Matrix
{
    int a[M][M];

    Matrix()
    {
        memset(a, 0 ,sizeof(a));
    }

    Matrix operator*(const Matrix &to) 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] + 1ll * a[i][k] * to.a[k][j] % MOD) % MOD;
        return res;
    }
} q;

void init_ones()
{
    for(int i = 0; i < 1 << m; i++)
    {
        int cnt = 0;
        for(int j = 0; j < m; j++)
            if(i >> j & 1) 
                cnt++;
        ones[i] = cnt;
    }
}

void init_q()
{
    for(int j = 0; j < 1 << m; j++)
    {
        int k1 = j >> 1;
        int k2 = j >> 1 | 1 << m - 1;
        if(ones[k1] <= k) q.a[k1][j] = 1;
        if(ones[k2] <= k) q.a[k2][j] = 1;
    }
}

Matrix qmi(Matrix &x, long long y)
{
    Matrix res;
    for(int i = 0; i < 1 << m; i++) res.a[i][i] = 1;
    while(y)
    {
        if(y & 1) res = res * x;
        x = x * x;
        y = y >> 1;
    }
    return res;
}

int main()
{
    cin >> n >> m >> k;
    init_ones();
    init_q();
    q = qmi(q, n);
    for(int s = 0; s < 1 << m; s++)
    {
        if(ones[s] > k) continue;
        Matrix p;
        p.a[0][s] = 1;
        p = p * q;
        ans = (ans + p.a[0][s]) % MOD;
    }
    cout << ans;
    return 0;
}

3、更进一步

在上面方法中，用p[i]表示1*size矩阵，每一行矩阵只有第s列初始值是1，且1的位置随着s的遍历而右移。

结果也是对第s列的值进行累加。

根据矩阵乘法的特性，何不直接定义p为size*size的矩阵，每一行初始值如下：

[1,0,0…0]

[0,1,0…0]

….

[0,0,0…1]

即是一个单位矩阵，用这个单位矩阵乘上q^n，既可以得到结果矩阵，而结果矩阵所有对角线的值之和即是答案。

因此结果矩阵就是q^n，将对角线值累加即可。

100分代码：

#include <bits/stdc++.h>
using namespace std;

const int M = 1 << 5, MOD = 1e9 + 7;
long long n, m, k, ans;
int ones[M]; //二进制有多少个1
struct Matrix
{
    int a[M][M];

    Matrix()
    {
        memset(a, 0 ,sizeof(a));
    }

    Matrix operator*(const Matrix &to) 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] + 1ll * a[i][k] * to.a[k][j] % MOD) % MOD;
        return res;
    }
} q;

void init_ones()
{
    for(int i = 0; i < 1 << m; i++)
    {
        int cnt = 0;
        for(int j = 0; j < m; j++)
            if(i >> j & 1) 
                cnt++;
        ones[i] = cnt;
    }
}

void init_q()
{
    for(int j = 0; j < 1 << m; j++)
    {
        int k1 = j >> 1;
        int k2 = j >> 1 | 1 << m - 1;
        if(ones[k1] <= k) q.a[k1][j] = 1;
        if(ones[k2] <= k) q.a[k2][j] = 1;
    }
}

Matrix qmi(Matrix &x, long long y)
{
    Matrix res;
    for(int i = 0; i < 1 << m; i++) res.a[i][i] = 1;
    while(y)
    {
        if(y & 1) res = res * x;
        x = x * x;
        y = y >> 1;
    }
    return res;
}

int main()
{
    cin >> n >> m >> k;
    init_ones();
    init_q();
    q = qmi(q, n);
    for(int i = 0; i < 1 << m; i++) 
        ans = (ans + q.a[i][i]) % MOD;
    cout << ans;
    return 0;
}