矩阵加速

· · 个人记录

updata 2024.8.27

好久没有这样整理过算法的笔记了，因为之前那个 cstdios 的号没了，所以打算重新写一篇这样的笔记。

先来放上前置知识，我的神仙同学写的。

相信各位 dalao 已经会了矩阵快速幂了！那我们来进一步了解这个算法的用途。

因为这个东西几年前学的，所以有些都忘记了(，如果有问题还请指出！

例题1

相信各位 dalao 都会拿到这个 60 分，就是暴力递推的那档部分分。

f_i=f_{i-1}+f_{i-3}

递推式大概长这样，我们发现这个东西空间的利用率极低，所以改成滚动数组：

f_{4}=f_{3}+f_{1} f_{3}=f_{2} f_{2}=f_{1}

注：这段东西依旧是过不了这道题目的。

但是这时我们会发现，这个东西它只需要维护一个类似于数列一样的东西。

\begin{bmatrix}f_3 \\ f_2 \\ f_1 \end{bmatrix}

需要变成

\begin{bmatrix}f_4 \\ f_3 \\ f_2 \end{bmatrix}

然后我们考虑对于第一个数怎么推出来，也就是之前矩阵的 f_1 \times 1+f_2 \times 0+f_3 \times 1 那么对于初始矩阵的第一行应该长这样：

\begin{bmatrix}1 & 0 & 1\end{bmatrix}

那么对于第二个数呢？是不是应该继承第一个数？

因此第二行应该是：

\begin{bmatrix}1 & 0 & 0\end{bmatrix}

同理可得第三行是：

\begin{bmatrix}0 & 1 & 0\end{bmatrix}

那么整一个初始矩阵就是：

\begin{bmatrix}1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0\end{bmatrix}

上面那个数列左乘这个矩阵，就可以进一步递推。

对于第 n 个数，就是不停的左乘 n-3 次。

然后对于这种不停左乘的操作，我们可以用矩阵快速幂去优化他！

如果你对于矩阵乘法不太熟练，建议去这里多看看。

#include <iostream>
#include <cstdio>
#include <cstring>
#define int long long 
using namespace std;
const int INF=5;
const int Mod=1e9+7;
struct _node_mul {
    int aa[INF][INF];
}ba,ans,c;
_node_mul mul(_node_mul x,_node_mul y,int lenx,int leny) {
    memset(c.aa,0,sizeof c.aa);
    for (int i=1;i<=lenx;i++)
        for (int j=1;j<=leny;j++)
            for (int k=1;k<=lenx;k++)
                c.aa[i][j]+=x.aa[i][k]*y.aa[k][j]%Mod,c.aa[i][j]%=Mod;
    return c;
}
int t,n;
void print(_node_mul kk) {
    for (int i=1;i<=3;i++) {
        for (int j=1;j<=3;j++)
            cout<<kk.aa[i][j]<<" ";
        cout<<"\n";
    }
}
int ksm(int x) {
    memset(ba.aa,0,sizeof ba.aa);
    memset(ans.aa,0,sizeof ans.aa);
    ba.aa[1][1]=1;ba.aa[1][3]=1;
    ba.aa[2][1]=1;ba.aa[3][2]=1;
    ans.aa[1][1]=1;ans.aa[2][1]=1;
    ans.aa[3][1]=1;
    while (x) {
        if (x&1) ans=mul(ba,ans,3,1);
        ba=mul(ba,ba,3,3);
        x>>=1;
//      print(ans);
    }
    return ans.aa[1][1];
}
signed main()
{
    ios::sync_with_stdio(false);
    cin>>t;
    while (t--) {
        cin>>n;
        if (n<=3) cout<<"1\n";
        else cout<<ksm(n-3)<<"\n";
    }
    return 0;
}

边界要稍微注意一下，还有就是左乘和右乘的区别。

例题2

对于这道题目，与上面那个的差不多，请读者自行思考。

例题3

稍微变形了一下的矩阵加速，对于那个数列应当修改为：

\begin{bmatrix}a_2 \\ a_1\end{bmatrix}

对于初始矩阵应该修改为：

\begin{bmatrix}p & q \\ 1 & 0\end{bmatrix}

因为当前的这个常数变了，因此在进行矩阵快速幂的时候需要将这个矩阵稍微进行修改。

当然，如果你认为矩阵快速幂仅仅是只能做到这些的，那就大错特错了！

例题4 例题4

双倍经验！

对于这道题目来说，其实本质上还是一个 DP 的想法。

设状态 f_{i,j} 表示到第 i 个点，时刻为 j 的方案数，对于原地不动，可以自己向自己连一条边，对于自爆，可以将边连到 0 号节点，并且 0 号节点没有出边。

#include <iostream>
#include <cstdio>
using namespace std;
const int INF=1e6+5;
const int INFN=35;
const int Mod=2017;
const int INFM=1e5+5;
struct _node_edge{
    int to_,next_;
}edge[INFM<<1];
int n,m,head[INFN],tot,t;
void add_edge(int x,int y) {
    edge[++tot]=(_node_edge){y,head[x]};
    head[x]=tot;return ;
}
short f[INFN][INF],ans;
signed main()
{
    ios::sync_with_stdio(false);
    cin>>n>>m;
    for (int i=1;i<=m;i++) {
        int x=0,y=0;cin>>x>>y;
        add_edge(x,y);add_edge(y,x);
    }
    for (int i=0;i<=n;i++) {
        add_edge(i,i);
        if (i) add_edge(i,0);
    }
    cin>>t;
    f[1][0]=1;
    for (int i=1;i<=t;i++) {
        for (int j=0;j<=n;j++) {
            for (int k=head[j];k;k=edge[k].next_) {
                int v=edge[k].to_;
                f[v][i]+=f[j][i-1];
                f[v][i]%=Mod;
            }
        }
    }
    for (int i=0;i<=n;i++) {
        ans+=f[i][t],ans%=Mod;
    }
    cout<<ans<<"\n";
    return 0;
}

上面这个代码开了 O2 可以过第一道例题 4。

但是我们发现这个 DP 非常的优美，这个时刻的影响只有前一层，这个性质非常有趣。

我们考虑将时刻那个给滚掉，就会发现本质上就是一个数列！

我们再将链式前向星修改成矩阵存图！

我们发现最核心的内容变成了：

f[1][0]=1;
    for (int i=1;i<=t;i++) {
        for (int j=0;j<=n;j++) {
            for (int k=0;k<=n;k++) {
//              if (!Map[j][k]) continue;
                int v=k;
                f[v][i]+=f[j][i-1]*Map[j][k];
                f[v][i]%=Mod;
            }
        }
    }

这个东西看起来就可以对 f 进行加速的操作！

类比一下矩阵快速幂，发现两个东西实际上是差不多的！

#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
const int INF=1e6+5;
const int INFN=305;
const int Mod=2017;
const int INFM=1e5+5;
int n,m,t,Map[INFN][INFN],Mapc[INFN][INFN];
short f[INFN],ans,c[INFN];
void mulf() {
    memset(c,0,sizeof c);
    for (int i=0;i<=n;i++)
        for (int j=1;j<=1;j++)
            for (int k=0;k<=n;k++)
                c[i]+=f[k]*Map[k][i]%Mod,c[i]%=Mod;
    for (int i=0;i<=n;i++)
        f[i]=c[i];
    return ;
}
// DP 的核心部分类比这里！
void print() {
    for (int i=0;i<=n;i++)
        cout<<f[i]<<" ";
    cout<<"\n";
}
void print1() {
    for (int i=0;i<=n;i++) {
        for (int j=0;j<=n;j++)
            cout<<Map[i][j]<<" ";
        cout<<"\n";
    }
}
void mulm() {
    memset(Mapc,0,sizeof Mapc);
    for (int i=0;i<=n;i++)
        for (int j=0;j<=n;j++)
            for (int k=0;k<=n;k++)
                Mapc[i][j]+=Map[i][k]*Map[k][j]%Mod,Mapc[i][j]%=Mod;
    for (int i=0;i<=n;i++)
        for (int j=0;j<=n;j++)
            Map[i][j]=Mapc[i][j];
    return ;
}
void ksm(int x) {
    while (x) {
//      print();print1();
        if (x&1) mulf();
        mulm();x>>=1;
    }
    for (int i=0;i<=n;i++)
        ans+=f[i],ans%=Mod;
    return ;
}
signed main()
{
    ios::sync_with_stdio(false);
    cin>>n>>m;
    for (int i=1;i<=m;i++) {
        int x=0,y=0;cin>>x>>y;
        Map[x][y]=Map[y][x]=1;
//      add_edge(x,y);add_edge(y,x);
    }
    for (int i=0;i<=n;i++) {
//      add_edge(i,i);
//      if (i) add_edge(i,0);
        Map[i][i]=1;
        if (i) Map[i][0]=1;
    }
    cin>>t;
    f[1]=1;
    ksm(t);
    cout<<ans<<"\n";
//  for (int i=1;i<=t;i++) {
//      for (int j=0;j<=n;j++) {
//          for (int k=0;k<=n;k++) {
////                if (!Map[j][k]) continue;
//              int v=k;
//              f[1][k]+=f[1][j]*Map[j][k];
//              f[v]%=Mod;
//          }
//      }
//  }
    return 0;
}

难不成矩阵快速幂就只能做计数类的题目？不！它甚至可以类似 floyd 跑最短路！

例题5

这种题目就比较的难了，考虑 k=0 就是朴素的 floyd，k=1 时，可以枚举每一条边，然后强制让其走过这条边。

上面的东西是 O(n^2m+n^3) 的，然后下面的这个东西就需要用到矩阵快速幂了！

考虑 k 较大的情况，其实我们完全可以分配成 k-1 和 1，更具体的，我们可以将中转点前面的视为 k-1 次，后面的视为 1 次，但这样子转移复杂度就会爆炸。

其实我们做了这么多矩阵快速幂的题目，已经可以发现，这个 min 其实具有结合律，而且他们长的很想矩阵乘法的形式。

于是我们就可以对他进一步优化了！

#include <iostream>
#include <cstdio>
#include <cstring>
#define int long long 
using namespace std;
const int INF=105;
const int INFN=5e3+5;
int Map[INF][INF],f[INF][INF],f1[INF][INF],c[INF][INF];
int n,m,kk,aa[INFN],bb[INFN],cc[INFN];
void mulans() {
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
            c[i][j]=1e18;
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
            for (int k=1;k<=n;k++) {
//              if (f[i][k]>1e17 || f1[k][j]>1e17) continue;
                c[i][j]=min(c[i][j],f1[i][k]+f[k][j]);
            }
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
            f[i][j]=c[i][j];
    return ;
}
void mulf() {
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
            c[i][j]=1e18;
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
            for (int k=1;k<=n;k++) {
//              if (f1[i][k]>1e17 || f1[k][j]>1e17) continue;
                c[i][j]=min(c[i][j],f1[i][k]+f1[k][j]);
            }
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
            f1[i][j]=c[i][j];
}
void print() {
    for (int i=1;i<=n;i++) {
        for (int j=1;j<=n;j++)
            cout<<f[i][j]<<" ";
        cout<<"\n";
    }
}
void ksm(int x) {
    while (x) {
        if (x&1) mulans();
        mulf();x>>=1;
//      print(); 
    }
    return ;
}
signed main()
{
    memset(Map,63,sizeof Map);
    memset(f,63,sizeof f);
    ios::sync_with_stdio(false);
    cin>>n>>m>>kk;
    for (int i=1;i<=n;i++)
        f[i][i]=Map[i][i]=0;
    for (int i=1;i<=m;i++) {
        int x=0,y=0,z=0;cin>>x>>y>>z;
        aa[i]=x;bb[i]=y;cc[i]=z;
        Map[x][y]=z;f[x][y]=z;
    }
    for (int k=1;k<=n;k++)
        for (int i=1;i<=n;i++)
            for (int j=1;j<=n;j++)
                f[i][j]=min(f[i][k]+f[k][j],f[i][j]);
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
            f1[i][j]=f[i][j];
    for (int i=1;i<=m;i++)
        for (int j=1;j<=n;j++)
            for (int k=1;k<=n;k++) {
//              f1[j][k]=min(f1[j][k],f[j][k]);
                f1[j][k]=min(f1[j][k],f[j][aa[i]]+f[bb[i]][k]-cc[i]);
            }
//  for (int i=1;i<=n;i++)
//      for (int j=1;j<=n;j++)
//          f1[i][j]=f[i][j]; 
//  for (int t1=2;t1<=kk;t1++)
//      for (int k=1;k<=n;k++)
//          for (int i=1;i<=n;i++)
//              for (int j=1;j<=n;j++)
//                  f[t1][i][j]=min(min(f[t1][i][j],f[t1-1][i][j]),f[t1-1][i][k]+f[1][k][j]);
//  cout<<f[kk][1][n]<<"\n";
//  cout<<f1[1][n]<<"\n";
    if (!kk) {cout<<f[1][n]<<"\n";return 0;}
    ksm(kk);
    cout<<f[1][n]<<"\n";
    return 0;
}

其实看看样子，然后再看看是不是有结合律，就可以看出是否是矩阵优化的题目了 (