AtCoder Beginner Contest 204 F - Hanjo 2（状压DP + 矩阵快速幂加速）

F - Hanjo 2

题意

给出不限数量的 $1\ast 1$ 和 $1\ast 2$ 的方块，要求填满 $n\ast m$ 的空间的方案数。其中 $n<=6,m<=10^{12}$。

解析

本题和 Acwing 的状压DP题很像， 蒙德里安的梦想

我们先不考虑巨大的第二维，先考虑状态和转移，根据蒙德里安的梦想的dp状态定义：$f_{ij}:$ $j$ 代表的是二进制，而二进制数位上的 $1$ 代表第 $i$ 列的该行有一个方块是横着放置伸展到第 $i+1$ 列，蒙德里安没有 $1\ast 1$ 的方块所以要求剩余的连续空位必须是偶数，且方案唯一。但我们无论当前是怎样的形状都能用 $1\ast 1$ 的方块去填补剩余的空位，所以我们无须讨论放置方案是否合法而是需要计算放置 $1\ast 1$ 和 $1\ast 2$ 方案数。

我们方程的入口就是 $f_{00}=1$，即第 $0$ 列不存在一个 $1\ast 2$ 的方块伸展到 $1$ 列，方程的出口就是 $ans=f_{m0}$ 同样是没有一个 $1\ast 2$ 的方块伸到第 $m+1$ 列的方案。

那么这里先给出 DP 方程转移，以及方案计算，有详细的代码注释。

bool check(int x, int y){
    for(int i = 0; i < n; i ++)
        if((x >> i & 1) && (y >> i & 1)) return false; // 前一列的第i行有横着放的，那么这列的这一行就不能有横着放
    return true;
}

int c[10] = {1, 1, 2, 3, 5, 8, 13}; // c[i]:长度为i的连续空位，随意放1*1和1*2方块并填满的方案数
int get_res(int x, int y){ // 放竖着放1*2和1*1的方案数
    int sum = 0, res = 1;
    for(int i = 0; i < n; i ++){
        if(!(x >> i & 1) && !(y >> i & 1)) sum ++; // 前一列第i行没有横着放，以及本列第i行也没有横着放就是空位，可以考虑放1*1/1*2
        else res *= c[sum], sum = 0;
    }
    res *= c[sum];
    return res;
}

void dp(){ // 具体可以参考Acwing的蒙德里安的梦想，只是稍作修改
    ll f[110][110] = {0}; // 二进制上的 1 代表 第 i 列的 该行是否有一个横着放的2*1方块，剩余空位可以放1*1和1*2
    f[0][0] = 1;
    for(int i = 1; i <= m; i ++){
        for(int j = 0; j < (1 << n); j ++){
            for(int k = 0; k < (1 << n); k ++){
                if(check(j, k)){
                    f[i][k] = (f[i][k] + f[i - 1][j] * get_res(j, k) % mod) % mod;
                }
            }
        }
    }   
    cout << f[m][0] % mod;
    return ; 
}

我们现在已经写出DP状态和方程转移，但是这样转移是 $O(m\ast 2^{2n})$ 矩阵巨大的第二维 $m<=10^{12}$，告诉我们肯定需要优化。

观察每次的转移，发现每个第 $i$ 列到第 $i+1$ 列之间的转移都是相同的，唯一变化只是第 $i$ 列的值。 于是我们就想到用矩阵快速幂来加速这一过程。具体看代码以及注释。

完整代码

#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int SIZ = 100, mod = 998244353;

ll n, m;

bool check(int x, int y){
    for(int i = 0; i < n; i ++)
        if((x >> i & 1) && (y >> i & 1)) return false; // 前一列的第i行有横着放的，那么这一行就不能有横着放
    return true;
}

int c[10] = {1, 1, 2, 3, 5, 8, 13}; // c[i]:长度为i的连续空位，随意放1*1和1*2方块并填满的方案数
int get_res(int x, int y){ // 放竖着放1*2和1*1的方案数
    int sum = 0, res = 1;
    for(int i = 0; i < n; i ++){
        if(!(x >> i & 1) && !(y >> i & 1)) sum ++; // 前一列第i行没有横着放，以及本列第i行也没有横着放就是空位，可以考虑放1*1/1*2
        else res *= c[sum], sum = 0;
    }
    res *= c[sum];
    return res;
}

void dp(){ // 具体可以参考Acwing的蒙德里安的梦想，只是稍作修改
    ll f[110][110] = {0}; // 二进制上的 1 代表 第 i 列的 该行是否有一个横着放的2*1方块，剩余空位可以放1*1和1*2
    f[0][0] = 1;
    for(int i = 1; i <= m; i ++){
        for(int j = 0; j < (1 << n); j ++){
            for(int k = 0; k < (1 << n); k ++){
                if(check(j, k)){
                    f[i][k] = (f[i][k] + f[i - 1][j] * get_res(j, k) % mod) % mod;
                }
            }
        }
    }   
    cout << f[m][0] % mod;
    return ; 
}

struct Matrix{
    ll Mat[SIZ][SIZ];
    Matrix(){
        for(int i = 0; i < (1 << n); i ++)
            for(int j = 0; j < (1 << n); j ++) Mat[i][j] = 0;
    }
    inline void build(){
        for(int i = 0; i < (1 << n); i ++) Mat[i][i] = 1;
    }
};

Matrix operator *(const Matrix &A, const Matrix & B){ // 重载矩阵乘法
    Matrix a;
    for(int i = 0; i < (1 << n); i ++){
        for(int j = 0; j < (1 << n); j ++){
            for(int k = 0; k < (1 << n); k ++){
                a.Mat[i][j] = (a.Mat[i][j] + A.Mat[i][k] * B.Mat[k][j] % mod) % mod;
            }
        }
    }
    return a;
}

Matrix ksm(Matrix a, ll b){ // 矩阵快速幂加速递推
    Matrix res; res.build();
    while(b){
        if(b & 1) res = res * a;
        a = a * a;
        b >>= 1;
    }
    return res;
}

int main(){
    cin >> n >> m;

    Matrix A, B;
    for(int i = 0; i < (1 << n); i ++){ // 因为每次转移的都是相同的，根据dp方程式转移来构造矩阵 参考dp()
        for(int j = 0; j < (1 << n); j ++){
            if(check(i, j)){
                A.Mat[j][i] = get_res(i, j); // 可以从状态 i -> j 并且 fj + fi * res ,所以矩阵系数为res
            }
        }
    }

    B.Mat[0][0] = 1; // 初始矩阵 即 f[0][0] = 1

    Matrix ans = ksm(A, m); // 从 0 -> m 需要转移 m 次 m 次幂
    B = ans * B; // 最后乘上初始矩阵

    cout << B.Mat[0][0];
    return 0;
}

提交记录：Submission #40530870