矩阵的运算 --- 矩阵快速幂（UVA10689 - Yet another Number Sequence）

        我们向大家介绍了快速幂取模(点我)，大家对快速幂都有了一定的了解。由此我们产生了一个想法，既然数值能够使用快速幂来提高幂运算的效率，那么同理矩阵也行。于是引出我们今天要介绍的“矩阵快速幂”。

       我们在学习数值的快速幂的时候，代码是这样的：

typedef long long ll; 
ll mod_pow(ll x, ll n, ll mod){
    ll res = 1;
    while( n > 0 ){ 
        if( n & 1 ) res = res * x % mod;  
        x = x * x % mod;
        n >>= 1;            
    }
    return res;
}

        我们思考，如果将上述的变量类型换为Matrix那么我们是不是就完成了矩阵的快速幂呢？答案是肯定的。

代码如下：

//矩阵快速幂
Matrix pow(Matrix a, int n){
    Matrix res(row, col);
    //将矩阵res初始化为单位矩阵
    for(int i = 0; i < row; i++){
        res.m[i][i] = 1;
    }
    //这时候直接使用快速幂的代码
    while(n > 0){
        if(n & 1) res = res * a;
        a = a * a;
        n >>= 1;
    }
    return res;
}

这时候我们就完成了矩阵的快速幂，完整代码如下：

#include <limits.h>
class Matrix{
public:
    int **m; //保存矩阵值的指针
    int row, col, mod; //分别保存矩阵的行、列以及取模的值
    Matrix(int r = 0, int c = 0, int d = INT_MAX);
    Matrix(const Matrix &value);
    ~Matrix();
    Matrix operator + (const Matrix &rht) const;
    Matrix operator * (const Matrix &rht) const;
    Matrix& operator = (const Matrix &rht);
    Matrix pow(int n) const;
};

Matrix::Matrix(int r, int c, int d):row(r), col(c), mod(d){
    m = new int*[row];
    for(int i = 0; i < row; i++){
        m[i] = new int[col];
    }
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            m[i][j] = 0;
        }
    }
}
Matrix::Matrix(const Matrix &value){
    row = value.row;
    col = value.col;
    mod = value.mod;
    m = new int*[row];
    for(int i = 0; i < row; i++){
        m[i] = new int[col];
    }
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            m[i][j] = value.m[i][j];
        }
    }
}
Matrix::~Matrix(){
    for(int i = 0; i < row; i++){
        delete[] m[i];
    }
    delete[] m;
}
Matrix Matrix::operator + (const Matrix &rht) const{
    Matrix temp(row, col, mod);
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            temp.m[i][j] = (m[i][j] + rht.m[i][j])%mod;
        }
    }
    return temp;
}
Matrix Matrix::operator * (const Matrix &rht) const{
    Matrix temp(row, rht.col, mod);
    for(int i = 0; i < row; i++){
        for(int k = 0; k < rht.row; k++){
            for(int j = 0; j < rht.col; j++){
                temp.m[i][j] = (temp.m[i][j] + m[i][k]*rht.m[k][j])%mod;
            }
        }
    }
    return temp;
}
Matrix& Matrix::operator = (const Matrix &rht){
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            m[i][j] = rht.m[i][j];
        }
    }
    return *this;
}
//矩阵快速幂
Matrix Matrix::pow(int n) const{
    Matrix a(*this), res(row, col, mod);
    //将矩阵res初始化为单位矩阵
    for(int i = 0; i < row; i++){
        res.m[i][i] = 1;
    }
    //这时候直接使用快速幂的代码
    while(n > 0){
        if(n & 1) res = res * a;
        a = a * a;
        n >>= 1;
    }
    return res;
}

        我们可以来一道题目练练手：UVA10689 - Yet another Number Sequence

乍一看，好像和这个矩阵快速幂没什么关系。我们仔细读题就会发现竟然需要数据竟然到了10^9,正常方法完全无解啊。对啊，我们现在不是正在学矩阵吗，于是我们换一种角度思考。很快我们就发现，斐波那数列还可以这样表示：

       通过递推求解方程我们可以得到以下方程：

        现在我们就成功的将问题转化为矩阵的幂了，很多问题也是如此。有时换一个角度，我们就可以看到一片光明。言归正传，既然我们已经将问题变为矩阵的幂的计算，那么此时我们关心的问题应当就在幂的大小，我们发现幂的大小在10^9的数量级，如果按照正常的方式，根本就无法再给定的时间计算完成。这时就是快速幂大显身手的时候了， O(log n)的时间复杂度，我们惊奇的发现快速幂只要20几次就完成了10^9的幂运算。

代码如下：

#include <cstdio>
#include <iostream>
#include <algorithm>
#include "Matrix.h" //上述的矩阵类
using namespace std;

int main()
{
    int a, b, n, m, T, mod;
    scanf("%d", &T);
    while(T-- && scanf("%d%d%d%d", &a, &b, &n, &m)){
        mod = 1;
        while(m--) mod *= 10;
        Matrix mat(2, 2, mod);
        mat.m[0][0] = mat.m[1][0] = mat.m[0][1] = 1;
        mat.m[1][1] = 0;
        Matrix ans(mat.pow(n-1));
        printf("%d\n", (b*ans.m[0][0] + a*ans.m[0][1])%mod);
    }
    return 0;
}

(如有错误，欢迎指正，转载请注明出处)