Q老师度假（变形矩阵快速幂优化DP）

2024-04-02 21:41:00

问题描述

忙碌了一个学期的 Q老师决定奖励自己 N 天假期。

假期中不同的穿衣方式会有不同的快乐值。

已知 Q老师一共有 M 件衬衫，且如果昨天穿的是衬衫 A，今天穿的是衬衫 B，则 Q老师今天可以获得 f[A][B] 快乐值。

在 N 天假期结束后，Q老师最多可以获得多少快乐值？

Input

输入文件包含多组测试样例，每组测试样例格式描述如下：

第一行给出两个整数 N M，分别代表假期长度与 Q老师的衬衫总数。（2 ≤ N ≤ 100000, 1 ≤ M ≤ 100）

接下来 M 行，每行给出 M 个整数，其中第 i 行的第 j 个整数，表示 f[i][j]。（1 ≤ f[i][j] ≤ 1000000）

测试样例组数不会超过 10。

Output

每组测试样例输出一行，表示 Q老师可以获得的最大快乐值。

Sample input

Sample output

2
9

解题思路

由于天数是连续的，明显具有子结构的性质，所以这个题我们可以用DP来做。令f[i][j]表示第 $i$ i天为止，穿 $j$ j衣服的快乐值，因此我们可以得到下面的状态转移方程：
$f[i][j]=max(f[i-1][k]+H[k][j]), 1\le k\le M \tag{1}$ f[i][j]=max(f[i−1][k]+H[k][j]),1≤k≤M(1)
其中，k是枚举前一天穿的衣服。经过计算，我们可以发现时间复杂度达到了 $O(10^9)$ O(109)，所以必须做到优化，观察这个状态转移方程的形式，联想到矩阵乘法公式：
$f[i][j]=\sum_{k=1}^{M}(f_2[i][k]\times f_3[k][j]) \tag{2}$ f[i][j]=k=1∑M(f2[i][k]×f3[k][j])(2)
将 $f_2$ f2矩阵中的 $f_2[i][k]$ f2[i][k]看作 $f[i-1][k]$ f[i−1][k]，将 $f_3$ f3矩阵中的 $f_3[k][j]$ f3[k][j]看作 $H[k][j]$ H[k][j]。我们可以得到：
$f[i][j]=\sum_{k=1}^{M}(f[i-1][k]\times H[k][j])\tag{3}$ f[i][j]=k=1∑M(f[i−1][k]×H[k][j])(3)
观察公式 $(1)$ (1)与公式 $(3)$ (3)， $max$ max与 $\sum$ ∑相对应， $+$ +与 $\times$ ×相对应。想到或许可以进行矩阵快速幂优化，矩阵快速幂实现的条件是矩阵乘法具有结合律，而公式 $(1)$ (1)是具有结合律的。证明方法如下：

设A、B、C分别为 $a\times b, b\times c, c\times d$ a×b,b×c,c×d的矩阵
$\begin{aligned} ((AB)C)[i,j] &=\max_{l=1}^{c}((AB)[i,l]+C[l,j])\\ &=\max_{l=1}^{c}(\max_{k=1}^{b}(A[i,k]+B[k,l])+C[l,j])\\ &=\max_{l=1}^{c}\max_{k=1}^{b}(A[i,k]+B[k,l]+C[l,j])\\ &=\max_{k=1}^{b}\max_{l=1}^{c}(A[i,k]+B[k,l]+C[l,j])\\ &=\max_{k=1}^{b}(A[i,k]+\max_{l=1}^{c}(B[k,l]+C[l,j]))\\ &=\max_{k=1}^{b}(A[i,k]+BC[k,j])\\ &=(A(BC))[i,j] \end{aligned}$ ((AB)C)[i,j]=l=1maxc((AB)[i,l]+C[l,j])=l=1maxc(k=1maxb(A[i,k]+B[k,l])+C[l,j])=l=1maxck=1maxb(A[i,k]+B[k,l]+C[l,j])=k=1maxbl=1maxc(A[i,k]+B[k,l]+C[l,j])=k=1maxb(A[i,k]+l=1maxc(B[k,l]+C[l,j]))=k=1maxb(A[i,k]+BC[k,j])=(A(BC))[i,j]
因此公式 $(1)$ (1)可以使用矩阵快速幂来进行优化，矩阵乘法的定义由公式 $(3)$ (3)更改为公式 $(1)$ (1)。转移矩阵如下：
$ans[i]=\begin{bmatrix} f[i][1]\\ \vdots\\ f[i][M] \end{bmatrix}= \begin{bmatrix} H[1][1]&\cdots&H[M][1]\\ \vdots&\ddots&\vdots\\ H[1][M]&\cdots&H[M][M] \end{bmatrix}* \begin{bmatrix} f[i-1][1]\\ \vdots\\ f[i-1][M] \end{bmatrix}$ ans[i]=⎣⎢⎡f[i][1]⋮f[i][M]⎦⎥⎤=⎣⎢⎡H[1][1]⋮H[1][M]⋯⋱⋯H[M][1]⋮H[M][M]⎦⎥⎤∗⎣⎢⎡f[i−1][1]⋮f[i−1][M]⎦⎥⎤
初始状态 $f[1][1]\sim f[1][M]$ f[1][1]∼f[1][M]均为0，最终结果从 $n-1$ n−1次幂的矩阵中找最大值即可。

需要注意的是，由于矩阵乘法定义改变，不要忘记在矩阵快速幂中的单位矩阵设置为全0矩阵。

完整代码

//#pragma GCC optimize(2)
//#pragma G++ optimize(2)
//#include <bits/stdc++.h>
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
#include <string>
#include <climits>
#include <algorithm>
#include <queue>
#include <vector>
using namespace std;

struct Matrix{
    Matrix(int _size=0) {
        Size=_size; mat=new long long*[Size];
        for (int i=0; i<Size; i++) mat[i]=new long long[Size];
        for (int i=0; i<Size; i++)
            for (int j=0; j<Size; j++)
                mat[i][j]=0;
    }
    Matrix(const Matrix& t){
        Size=t.Size; mat=new long long*[Size];
        for (int i=0; i<Size; i++) mat[i]=new long long[Size];
        memcpy(mat,t.mat,sizeof(mat));
    }
    ~Matrix(){
        for (int i=0; i<Size; i++) delete[] mat[i];
        delete[] mat;
    }
    Matrix operator*(const Matrix& t) const{//伪.乘法
        Matrix ret(Size);
        for (int i=0; i<Size; ++i)
            for (int j=0; j<Size; ++j)
                for (int k=0; k<Size; ++k)
                    ret.mat[i][j]=max(ret.mat[i][j],mat[i][k]+t.mat[k][j]);
        return ret;
    }
    Matrix operator%(const int m) const {
        Matrix ret(Size);
        for (int i=0; i<Size; i++)
            for (int j=0; j<Size; j++)
                ret.mat[i][j]=mat[i][j]%m;
        return ret;
    }
    Matrix& operator=(const Matrix& t){
        for (int i=0; i<Size; i++) delete[] mat[i];
        delete[] mat;
        Size=t.Size;
        mat=new long long*[Size];
        for (int i=0; i<Size; i++) mat[i]=new long long[Size];
        for (int i=0; i<Size; i++)
            for (int j=0; j<Size; j++)
                mat[i][j]=t.mat[i][j];
        return *this;
    }
    void quick_pow(int x){
        Matrix ret(Size);
        //for (int i=0; i<Size; i++) ret.mat[i][i]=1;//注意，乘法含义改变，单位矩阵全0
        while(x){
            if(x&1) {
                ret=ret*(*this);
            }
            *this=(*this)*(*this);
            x>>=1;
        }
        *this=ret;
    }
    void output(){
        printf("Size: %d\n",Size);
        for (int i=0; i<Size; i++){
            for (int j=0; j<Size; j++){
                printf("%lld ",mat[i][j]);
            }
            printf("\n");
        }
    }
    int Size;
    long long** mat;
};
int getint(){
    int x=0,s=1; char ch=' ';
    while(ch<'0' || ch>'9'){ ch=getchar(); if(ch=='-') s=-1;}
    while(ch>='0' && ch<='9'){ x=x*10+ch-'0'; ch=getchar();}
    return x*s;
}
int n,m;
int main(){
    //ios::sync_with_stdio(false);
    //cin.tie(0);
    while(scanf("%d %d",&n,&m)!=EOF){
        Matrix mat1(m);
        for (int i=0; i<m; i++)
            for (int j=0; j<m; j++)
                scanf("%lld",&mat1.mat[i][j]);

        mat1.quick_pow(n-1);
        long long ans=0;
        for (int i=0; i<m; i++)
            for (int j=0; j<m; j++)
                ans=max(ans,mat1.mat[i][j]);
        //mat1.output();
        printf("%lld\n",ans);
    }
    return 0;
}

码农公寓