矩阵快速幂算法详细解析

在上一篇博客我谈到了大数的快速幂，而相对于矩阵的指数运算同样可以有方法，在此之前我们来看看矩阵的乘法：

矩阵的乘法是需要矩阵A的行数与矩阵B的列数相等的（A*B的前提条件）但矩阵快速幂一般只用到方阵（行数和列数相等的情况），从而避免了前提条件。

如：
$$\begin{gathered} \begin{array}{cc} \\ \left[\begin{array}{cccc}{A}_{11}& A_{12}& \cdots & A_{1n}\\ A_{21}& A_{22}& \cdots & A_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ A_{n1}& A_{n2}& \cdots & A_{nn}\end{array}\right]\left[\begin{array}{cccc}{B}_{11}& B_{12}& \cdots & B_{1n}\\ B_{21}& B_{22}& \cdots & B_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ B_{n1}& A_{n2}& \cdots & A_{nn}\end{array}\right]=& \left[\begin{array}{cccc}{C}_{11}& C_{12}& \cdots & C_{1n}\\ C_{21}& C_{22}& \cdots & C_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ C_{n1}& C_{n2}& \cdots & C_{nn}\end{array}\right]\end{array} \end{gathered}$$
根据线性代数的知识我们也容易知道，俩个矩阵相乘，前提是A的行数与B的列数相等，结果是其在新矩阵的行所对应的A的行与在新矩阵的列所对应的B的列对应相乘相加，所以我们可以先得出矩阵乘法的代码：

for(i=1;i<=n;i++)
{
    for(j=1;j<=n;j++)
    {
        for(k=1;k<=n;k++)
        {
            c[i][j] = c[i][j] + a[i][k] * b[k][j];
        }
    }
}

所以我们给出例子:

假如有一个n * n的方阵A,所谓的方阵就是行数与列数相等的矩阵，给出数M,求矩阵A的M次幂A^M

这个例子就很直接了，就是直接求矩阵的M次幂，结合快速幂算法和矩阵相乘可以得出代码：

#include <iostream>
#include <algorithm>
#include <cstdio>
using namespace std;

int M,n;

struct node//定义一个矩阵类型的结构体
{
    int m[100][100];
} ans,res;//ans是结果，res是最初的方阵
node mul(node A,node B)
{
    int i,j,k;
    node temp;//定义一个临时矩阵，存放A*B的结果
    for(i=0; i<n; i++)//先全部定义为0
    {
        for(j=0; j<n; j++)
        {
            temp.m[i][j] = 0;
        }
    }
    for(i=0; i<n; i++)//矩阵相乘的代码
    {
        for(j=0; j<n; j++)
        {
            for(k=0; k<n; k++)
            {
                temp.m[i][j] += A.m[i][k] * B.m[k][j];
            }
        }
    }
    return temp;
}
void quickpower(int M,int n)
{
    int i,j;
    for(i=0; i<n; i++)
    {
        for(j=0; j<n; j++)
        {
            if(i == j)
            {
                ans.m[i][j] = 1;
            }
            else
                ans.m[i][j] = 0;
        }
    }//这里是思想的转换，之前我们定义为1去计算，所以我们先初始化ans为
    //单位矩阵，我们知道单位矩阵与任何矩阵的乘积为其本身
    while(M)//快速幂的步骤
    {
        if(M & 1)
        {
            ans = mul(ans,res);
        }
        res = mul(res,res);
        M = M >> 1;
    }
}
int main()
{
    cin>>n;//方阵的阶数
    cin>>M;//指数
    int i,j;
    for(i=0; i<n; i++)
    {
        for(j=0; j<n; j++)
        {
            cin>>res.m[i][j];//初始化方阵res
        }
    }
    quickpower(M,n);//进行快速幂
    for(i=0; i<n; i++)//输出
    {
        for(j=0; j<n; j++)
        {
            printf("%d ",ans.m[i][j]);
        }
        printf("\n");
    }
    return 0;
}

以上就是矩阵快速幂的思路与代码，但是一般的题并不会简单的让你去求矩阵的快速幂，其通常是与斐波那切数相结合，我们知道最简单的斐波那切的函数就是f(n)=f(n-1)+f(n-2)（当然不清楚的伙伴可以先去了解一下斐波那切数），倘若n很大的情况用快速矩阵可以得到：
$$\begin{gathered} \begin{array}{ccc} \\ \left[\begin{array}{c}f(n)\\ f(n-1)\end{array}\right]=& \left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]\left[\begin{array}{c}f(n-1)\\ f(n-2)\end{array}\right]=& {\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]}^{n-2}\left[\begin{array}{c}f(2)\\ f(1)\end{array}\right]\end{array} \end{gathered}$$
原式斐波那切可以转换为矩阵相乘，然而不断推演下去就是求这个矩阵的n-2次方，当然有时也有其他形式f(n)=af(n-1) +bf(n-2) + c:
$$\begin{gathered} \begin{array}{cc} \\ \left[\begin{array}{c}f(n)\\ f(n-1)\\ c\end{array}\right]=& \left[\begin{array}{ccc}a& b& 1\\ 1& 0& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}f(n-1)\\ f(n-2)\\ c\end{array}\right]\end{array} \end{gathered}$$
f(n) = cⁿ - f(n-1):
$$\begin{gathered} \begin{array}{cc} \\ \left[\begin{array}{c}f(n)\\ c^n\end{array}\right]=& \left[\begin{array}{cc}-1& c\\ 0& c\end{array}\right]\left[\begin{array}{c}f(n-1)\\ c^{n-1}\end{array}\right]\end{array} \end{gathered}$$
下面给出一个例子：

勇者龙宫院圣哉与n名冒险者一起去讨伐神秘魔物，龙宫院圣哉十分谨慎，他只会在最后一刻出手，每名冒险者轮流攻击魔物，冒险者的攻击有着某种规律，目前造成的总伤害是上一名冒险者攻击后造成的总伤害的4倍与上上名冒险者攻击后造成的总伤害的3倍之和，即当前总伤害f(n)=4f(n-1)+3f(n-2)（魔物的奇怪设定使总伤害忽高忽低），又由于异世界的奇异设定，冒险者们的总伤害不会超过666666，即对666666取模，龙宫院圣哉清楚的知道这个魔物的血量为m（m>666666），他想知道在所有的冒险者攻击完了以后，自己需要造成多少点伤害才能杀死魔物？目前第一名冒险者攻击后总共造成了4点伤害，第二名冒险者攻击后总共造成了233点伤害。2<n<1e9.输入勇者n和魔王血量m,
输入 3
666667
输出 665723

这个很明显是矩阵快速幂，f(n)=4f(n-1)+3f(n-2)

#include <iostream>
#include <algorithm>
#include <cstdio>
using namespace std;
struct node 
{
    long long int m[10][10];
}ans,res;
node mul(node A,node B)
{
    int i,j,k;
    node temp;//定义一个临时矩阵，存放A*B的结果
    for(i=0; i<n; i++)//先全部定义为0
    {
        for(j=0; j<n; j++)
        {
            temp.m[i][j] = 0;
        }
    }
    for(i=0; i<n; i++)//矩阵相乘的代码
    {
        for(j=0; j<n; j++)
        {
            for(k=0; k<n; k++)
            {
                temp.m[i][j] += (A.m[i][k] * B.m[k][j])%666666;//取模
            }
        }
    }
    return temp;
}
node quickpower(node a,long long n)
{
    node c;
    memset(c.m,0,sizeof(c.m));
    int i;
    for(i=0;i<2;i++)  c.m[i][i]=1;//定义一个单位矩阵
    while(n)
    {
        if(n & 1)
        {
            c=mul(c,a);
        }
        a=mul(a,a);
        n=n>>1;
    }
    return c;
}
int main()
{
    long long n,k;
    cin>>n>>k;
    memset(ans.m,0,sizeof(ans.m));
    ans.m[0][0]=4;
    ans.m[0][1]=3;
    ans.m[1][0]=1;
    ans.m[1][1]=0;//这是构建常数矩阵，上面的公式可以参考，主要是其参数
    n=n-2;//前俩个冒险者减去
    ans=quickpower(ans,n);//求出常数矩阵的n-2次方
    memset(res.m,0,sizeof(res.m));
    res.m[0][0]=233;
    res.m[1][0]=4;
    res=mul(ans,res);//求出最后的结果矩阵
    k=k-res.m[0][0];
    cout<<k<<endl;//减去最后的造成的总伤害
    return 0;
}

以上就是我对矩阵快速幂的整理和理解，希望对大家有所帮助！