jobdu1443矩阵快速幂

---

本文主要根据之前讲过的二分求幂，通过一道题目引入矩阵快速求幂的方法。

问题

先看下面的问题： [ jobdu-1443]

题目描述：
题目描述：
A为一个方阵，则Tr A表示A的迹（就是主对角线上各项的和），现要求Tr(A^k)%9973。
输入：
数据的第一行是一个T，表示有T组数据。
每组数据的第一行有n(2 <= n <= 10)和k(2 <= k < 10^9)两个数据。接下来有n行，每行有n个数据，每个数据的范围是[0,9]，表示方阵A的内容。
输出：
对应每组数据，输出Tr(A^k)%9973。
样例输入：
2
2 2
1 0
0 1
3 99999999
1 2 3
4 5 6
7 8 9
样例输出：
2
2686

思路

题目的意思非常明确。可以立马想到的是将矩阵的K次幂求出来，然后计算迹，最后求模即可。考虑到矩阵乘法的复杂度是O(N^3)，K次幂要进行K-1次乘法。计算量较大。

更好的办法是借鉴整数二分求幂的思路，将实数替换为矩阵。就得到了矩阵快速幂的办法。当然，这其中有几个注意点：

• 单位矩阵类似于整数中1的作用
• 矩阵求模的位置要小心。

下面回顾下二分求幂的算法：

POWER_INTEGER(x, n)
1 pow ← 1
2 while (n > 0)
3 do if (n mod 2 = 1)
4 then pow ← pow * x
5 x ← x * x
6 n ← n / 2
7 return pow

错误代码片段

    Matrix& operator%(int mod)
    {
        for( int i = 0; i < size_; ++i )
        {
            for( int j = 0; j < size_; ++j )
            {
                mat_[i][j] = mat_[i][j]%mod;
            }
        }
        return *this;
    }

Matrix fast_pow(Matrix& a, int b)
{
    Matrix ans(a.size_);
    ans.set_diag(1);

    while(b)
    {
        if(b%2)
        {
            ans = (ans * a)%MOD;
        }

        a = (a%MOD)*(a%MOD);
        b /= 2;
    }
    return ans;
}

这一块问题出在了，上面说道注意的第二点。矩阵做乘法的时候，中间的临时结果可能会溢出。所以，应该及时求模，避免溢出。所以，改进的办法是把求模换到乘法的地方。

正确代码

#include <iostream>
#include <cstring>
#include <fstream>
#define LOCAL
const int maxn = 10;
const int MOD = 9973;

struct Matrix
{

    int size_;
    int mat_[maxn+1][maxn+1];


    Matrix( int s = 0 ) : size_(s)
    {
        std::memset(mat_,0,sizeof(mat_));
    }
    Matrix( const Matrix& rhs )
    {
        size_ = rhs.size_;
        *this = rhs;
    }

    void init()
    {
        for( int i = 0; i < size_; ++i )
        {
            for( int j = 0; j < size_; ++j )
            {
                std::cin >> mat_[i][j];
            }
        }
    }
    void set_diag( int val )
    {
        for( int i = 0; i < size_; ++i )
        {
            for( int j = 0; j < size_; ++j )
            {
                mat_[i][j] = (i==j)?val:0;
            }
        }
    }
    void print()
    {
        for( int i = 0; i < size_; ++i )
        {
            for( int j = 0; j < size_; ++j )
            {
                std::cout << mat_[i][j] << " ";
            }
            std::cout << std::endl;
        }
    }

    Matrix operator*( const Matrix& b )
    {
        Matrix ret(b.size_);
        for( int row = 0; row < size_; ++row )
        {
            for( int col = 0; col < size_; ++col )
            {
                for( int i = 0; i < size_; ++i )
                {
                    ret.mat_[row][col] = (ret.mat_[row][col] +  (mat_[row][i]%MOD*b.mat_[i][col]%MOD)%MOD)%MOD;
                }
            }
        }
        return ret;
    }

    Matrix& operator=( const Matrix& b )
    {
        for( int i = 0; i < size_; ++i )
        {
            for( int j = 0; j < size_; ++j )
            {
                mat_[i][j] = b.mat_[i][j];
            }
        }
        return *this;
    }
    Matrix operator^(int b)
    {
        Matrix w(*this);

        Matrix ans(this->size_);
        ans.set_diag(1);

        while(b)
        {
            if(b%2)
            {
                ans = ans*w;
            }

            w = w*w;
            b /= 2;
        }
        return ans;
    }

    int get_tr()
    {
        int ans = 0;
        for( int i = 0; i < size_; ++i )
        {
            for( int j = 0; j < size_; ++j )
            {
                if(i==j)
                {
                    ans += mat_[i][j];
                    ans = ans%MOD;
                }
            }
        }
        return ans;
    }
};

int main()
{

    int t = 0;
    std::cin >> t;
    while( t-- )
    {
        int s,e;
        std::cin >> s >> e;
        Matrix a(s);
        a.init();

        Matrix ret = a^e;
        int ans = ret.get_tr();
        std::cout << ans << std::endl;

    }
}

总结

还是改bug的事，现在真的发现了。改bug才是真正提高能力的时候。因为它同样提高你分析能力的能力。找bug不是件容易的事。不过现在的心得还是，见到bug又能提高能力了。所以要积极。先圈定一个小范围，动手去试。

---