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

问题描述

忙碌了一个学期的 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 i天为止，穿 j j j衣服的快乐值，因此我们可以得到下面的状态转移方程：
f [ i ] [ j ] = m a x ( f [ i − 1 ] [ k ] + H [ k ] [ j ] ) , 1 ≤ k ≤ M (1) 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 ( 1 0 9 ) O(10^9) O(109)，所以必须做到优化，观察这个状态转移方程的形式，联想到矩阵乘法公式：
f [ i ] [ j ] = ∑ k = 1 M ( f 2 [ i ] [ k ] × f 3 [ k ] [ j ] ) (2) 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 f_2 f2矩阵中的 f 2 [ i ] [ k ] f_2[i][k] f2[i][k]看作 f [ i − 1 ] [ k ] f[i-1][k] f[i−1][k]，将 f 3 f_3 f3矩阵中的 f 3 [ k ] [ j ] f_3[k][j] f3[k][j]看作 H [ k ] [ j ] H[k][j] H[k][j]。我们可以得到：
f [ i ] [ j ] = ∑ k = 1 M ( f [ i − 1 ] [ k ] × H [ k ] [ j ] ) (3) 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) (1)与公式 ( 3 ) (3) (3)， m a x max max与 ∑ \sum ∑相对应， + + +与 × \times ×相对应。想到或许可以进行矩阵快速幂优化，矩阵快速幂实现的条件是矩阵乘法具有结合律，而公式 ( 1 ) (1) (1)是具有结合律的。证明方法如下：

设A、B、C分别为 a × b , b × c , c × d a\times b, b\times c, c\times d a×b,b×c,c×d的矩阵
( ( A B ) C ) [ i , j ] = max ⁡ l = 1 c ( ( A B ) [ 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 ] + B C [ k , j ] ) = ( A ( B C ) ) [ i , j ] \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) (1)可以使用矩阵快速幂来进行优化，矩阵乘法的定义由公式 ( 3 ) (3) (3)更改为公式 ( 1 ) (1) (1)。转移矩阵如下：
a n s [ i ] = [ f [ i ] [ 1 ] ⋮ f [ i ] [ M ] ] = [ H [ 1 ] [ 1 ] ⋯ H [ M ] [ 1 ] ⋮ ⋱ ⋮ H [ 1 ] [ M ] ⋯ H [ M ] [ M ] ] ∗ [ f [ i − 1 ] [ 1 ] ⋮ f [ i − 1 ] [ M ] ] 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 ] ∼ f [ 1 ] [ M ] f[1][1]\sim f[1][M] f[1][1]∼f[1][M]均为0，最终结果从 n − 1 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;//注意，乘法含义改变，单位矩阵全0while(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;
}