问题描述

忙碌了一个学期的 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

3 2
0 1
1 0
4 3
1 2 3
1 2 3
1 2 3

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=1maxc​k=1maxb​(A[i,k]+B[k,l]+C[l,j])=k=1maxb​l=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;//注意，乘法含义改变，单位矩阵全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;
}