POJ_3233 Matrix Power Series【矩阵递推方法】和 【二分+矩阵快速幂方法】

Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3+ … + Ak.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 109) and m (m < 104). Then follow n lines each containing n nonnegative integers below 32,768, giving A’s elements in row-major order.

Output

Output the elements of S modulo m in the same way as A is given.

Sample Input

2 2 4
0 1
1 1

Sample Output

1 2
2 3

题意理解：读题意可知，就是求矩阵的前k次方和。

解题思路：该题可以有两种解题思路，一是通过矩阵递推式，直接通过状态转移矩阵加快矩阵变化，推出最终结果。二是通过二分的思想解决。

一、关于递推，k最大可以达到1e9，那么考虑用矩阵快速幂求解。先要确定这初始矩阵是

，再根据递推式,得到状态转移矩阵是,这样我们就可以通过状态转移矩阵的n-1次方再乘以初始矩阵就是我们所要求的目标矩阵，因为我选用的初始矩阵是，所以状态转移矩阵就需要n-1次方，如果初始矩阵不同，那么转移矩阵所需要转移的次数也会不同。而通过这题，我们可以非常明显的看到在运用矩阵快速幂的时候，我们需要确定的有两个量，一个是初始矩阵，另一个是状态转移矩阵，而使用矩阵快速幂，使状态转移矩阵转移的更快。

二、关于二分主要是推出

                                                         

当k为偶数时不需要加 A^k，因为刚好二分完，当k为奇数时，因为k/2是向下取整，会掉了最后的A^k，所以要补回来。

关于代码，矩阵相乘非常好写，而矩阵快速幂就是在快速幂的基础上将数的相乘转换为了矩阵相乘，也是相当容易的。别忘了取模。

代码：

方法一：

//对状态转移矩阵用矩阵快速幂加速
#include <iostream>
#include <cstring>
#define ll long long
using namespace std;
ll N, Mod;
struct Matrix
{
    int m[61][61];
};

Matrix multi(Matrix a, Matrix b)
{
    Matrix ans;
    memset(ans.m, 0, sizeof(ans.m));
    for(int i = 1; i <= N; i++)
        for(int j = 1; j <= N; j++)
            for(int k = 1; k <= N; k++)
                ans.m[i][j] = (ans.m[i][j] + a.m[i][k]*b.m[k][j] % Mod) % Mod;
    return ans;
}

Matrix pow(Matrix a, ll n)
{
    Matrix ans;
    memset(ans.m, 0, sizeof(ans.m));
    for(int i = 1; i <= N; i++)
        ans.m[i][i] = 1;
    while(n)
    {
        if(n&1)
            ans = multi(ans, a);
        n >>= 1;
        a = multi(a, a);
    }
    return ans;
}

int main()
{
    Matrix ans, T;
    ll n, k, data;
    cin >> n >> k >> Mod;
    N = n*2;
    memset(ans.m, 0, sizeof(ans.m));
    memset(T.m, 0, sizeof(T.m));
    for(int i = 1; i <= n; i++)
    {
        ans.m[i][i] = 1;

    }
    for(int i = 1; i <= n; i++)
        for(int j = 1; j <= n; j++)
        {
            cin >> data;
            ans.m[i][j+n] = data % Mod;
            ans.m[i+n][j+n] = data % Mod;
            T.m[i][j] = data % Mod;
            T.m[i+n][j] = data % Mod;
        }
//    for(int i = 1; i <= N; i++)
//    {
//        for(int j = 1; j < N; j++)
//        {
//            cout << T.m[i][j]<< ' ';
//        }
//        cout << T.m[i][N] << endl;
//    }
    ans = pow(ans, k-1);
    ans = multi(ans, T);
    for(int i = 1; i <= n; i++)
    {
        for(int j = 1; j < n; j++)
        {
            cout << ans.m[i][j]<< ' ';
        }
        cout << ans.m[i][n] << endl;
    }
    return 0;
}

方法二：

因为该方法是重新做该题时重新码的代码，可能风格会有不同。

//矩阵快速幂 + 二分思想
#include <iostream>
#include <cstring>
#include <cstdio>

#define LL long long
using namespace std;
LL Mod;
int N, K;

struct Matrix
{
    LL m[35][35];
};
Matrix A, pr;

Matrix Multi(Matrix a, Matrix b)
{
    Matrix ans;
    memset(ans.m, 0, sizeof(ans.m));
    for(int i = 1; i <= N; i++)
        for(int j = 1; j <= N; j++)
        {
            for(int k = 1; k <= N; k++)
                ans.m[i][j] = (ans.m[i][j] + a.m[i][k]*b.m[k][j])%Mod;
        }
    return ans;
}

Matrix Add(Matrix a, Matrix b)
{
    Matrix ans;
    memset(ans.m, 0, sizeof(ans.m));
    for(int i = 1; i <= N; i++)
    {
        for(int j = 1; j <= N; j++)
        {
            ans.m[i][j] = (a.m[i][j] + b.m[i][j])% Mod;
        }
    }
    return ans;
}

Matrix Power(Matrix a, int k)
{
    Matrix ans = pr;
    while(k)
    {
        if(k&1)
        {
            ans = Multi(ans, a);
        }
        k >>= 1;
        a = Multi(a, a);
    }
    return ans;
}

Matrix solve(int k)
{
    if(k == 1)  return A;
    Matrix b, temp = Power(A, k>>1);
    temp = Add(temp, pr);
    b = Multi(temp, solve(k>>1));
    if(k&1)
    {
        b = Add(b, Power(A, k));
    }
    return b;
}

int main()
{
    scanf("%d %d %I64d", &N, &K, &Mod);
    memset(pr.m, 0, sizeof(pr.m));
    for(int i = 1; i <= N; i++)
    {
        pr.m[i][i] = 1;
        for(int j = 1; j <= N; j++)
        {
            scanf("%I64d", &A.m[i][j]);
        }
    }
    Matrix ans = solve(K);
    for(int i = 1; i <= N; i++)
    {
        printf("%I64d", ans.m[i][1]);
        for(int j = 2; j <= N; j++)
        {
            printf(" %I64d", ans.m[i][j]);
        }
        printf("\n");
    }
    return 0;
}