poj3233（矩阵快速幂）

poj3233 http://poj.org/problem?id=3233

给定n ，k，m

然后是n*n行，

我们先可以把式子转化为递推的，然后就可以用矩阵来加速计算了。  矩阵是加速递推计算的一个好工具

我们可以看到，矩阵的每个元素都是一个矩阵，其实这计算一个分块矩阵，我们可以把分块矩阵展开，它的乘法和普通矩阵的乘法是一样的。

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <algorithm>
#include <iostream>
#include <queue>
#include <stack>
#include <vector>
#include <map>
#include <set>
#include <string>
#include <math.h>
using namespace std;
#pragma warning(disable:4996)
typedef long long LL;
const int INF = 1 << 30;
int MOD;
/*
矩阵的元素是矩阵，即分块矩阵
分块矩阵的乘法和普通矩阵的乘法是一样的
*/
struct Matrix
{
    int mat[66][66];
    int size;
    Matrix(int n)
    {
        size = n;
    }
    void makeZero()
    {
        for (int i = 0; i < size; ++i)
        {
            for (int j = 0; j < size; ++j)
                mat[i][j] = 0;
        }
    }
    void makeUnit()
    {
        for (int i = 0; i < size; ++i)
        {
            for (int j = 0; j < size; ++j)
                mat[i][j] = (i == j);
        }
    }
};

Matrix operator*(const Matrix &lhs, const Matrix &rhs)
{
    Matrix ret(lhs.size);
    ret.makeZero();
    for (int i = 0; i < lhs.size; ++i)
    {
        for (int j = 0; j < lhs.size; ++j)
        {
            for (int k = 0; k < lhs.size; ++k)
            {
                ret.mat[i][j] += (lhs.mat[i][k] * rhs.mat[k][j]);
                if (ret.mat[i][j] >= MOD)
                    ret.mat[i][j] %= MOD;
            }
        }
    }
    return ret;
}
Matrix operator^(Matrix &a, int k)
{
    Matrix ret(a.size);//单位矩阵
    ret.makeUnit();
    while (k)
    {
        if (k & 1)
        {
            ret = ret * a;
        }
        a = a * a;
        k >>= 1;
    }
    return ret;
}
int main()
{
    int n, k, m, i, j;
    while (scanf("%d%d%d", &n, &k, &MOD) != EOF)
    {
        Matrix a(2*n);
        Matrix tmp(2*n);;
        a.makeZero();
        for (i = 0; i < n; ++i)
        {
            for (j = 0; j < n; ++j)
            {
                scanf("%d", &a.mat[i][j]);
                if (a.mat[i][j] >= MOD)
                    a.mat[i][j] %= MOD;
                tmp.mat[i][j] = a.mat[i][j];
                tmp.mat[i][n + j] = a.mat[i][j];
            }
            a.mat[n + i][i] = a.mat[n + i][n + i] = 1;
        }
        a = a ^ (k - 1);
        Matrix ans(2*n);
        ans.makeZero();
        m = 2 * n;
        for (i = 0; i < n; ++i)
        {
            for (j = 0; j < n; ++j)
            {
                for (k = 0; k < m; ++k)
                {
                    ans.mat[i][j] += tmp.mat[i][k] * a.mat[k][j];
                    if (ans.mat[i][j] >= MOD)
                        ans.mat[i][j] %= MOD;
                }
            }
        }
        for (i = 0; i < n; ++i)
        {
            for (j = 0; j < n; ++j)
            {
                j == n - 1 ? printf("%d
", ans.mat[i][j]) : printf("%d ", ans.mat[i][j]);
            }
        }
    }
    return 0;
}

当然了，我们还可以用二分的方法方法来计算

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <algorithm>
#include <iostream>
#include <queue>
#include <stack>
#include <vector>
#include <map>
#include <set>
#include <string>
#include <math.h>
using namespace std;
#pragma warning(disable:4996)
typedef long long LL;
const int INF = 1 << 30;
int MOD;
/*
矩阵的元素是矩阵，即分块矩阵
分块矩阵的乘法和普通矩阵的乘法是一样的
*/
struct Matrix
{
    int mat[66][66];
    int size;
    Matrix(int n)
    {
        size = n;
    }
    void makeZero()
    {
        for (int i = 0; i < size; ++i)
        {
            for (int j = 0; j < size; ++j)
                mat[i][j] = 0;
        }
    }
    void makeUnit()
    {
        for (int i = 0; i < size; ++i)
        {
            for (int j = 0; j < size; ++j)
                mat[i][j] = (i == j);
        }
    }
};

Matrix operator*(const Matrix &lhs, const Matrix &rhs)
{
    Matrix ret(lhs.size);
    ret.makeZero();
    for (int i = 0; i < lhs.size; ++i)
    {
        for (int j = 0; j < lhs.size; ++j)
        {
            for (int k = 0; k < lhs.size; ++k)
            {
                ret.mat[i][j] += (lhs.mat[i][k] * rhs.mat[k][j]);
                if (ret.mat[i][j] >= MOD)
                    ret.mat[i][j] %= MOD;
            }
        }
    }
    return ret;
}
Matrix operator^(Matrix &a, int k)
{
    Matrix ret(a.size);//单位矩阵
    ret.makeUnit();
    while (k)
    {
        if (k & 1)
        {
            ret = ret * a;
        }
        a = a * a;
        k >>= 1;
    }
    return ret;
}
Matrix operator+(const Matrix &lhs, const Matrix &rhs)
{
    Matrix ret(lhs.size);
    ret.makeZero();
    for (int i = 0; i < lhs.size; ++i)
    {
        for (int j = 0; j < lhs.size; ++j)
        {
            ret.mat[i][j] = lhs.mat[i][j] + rhs.mat[i][j];
            if (ret.mat[i][j] >= MOD)
                ret.mat[i][j] %= MOD;
        }
    }
    return ret;
}
Matrix matrixSum(Matrix a, int k)
{
    if (k == 1)
        return a;
    Matrix ret = matrixSum(a, k / 2);
    if (k & 1)
    {
        Matrix tmp = a ^ (k / 2 + 1);
        ret = ret * tmp + ret + tmp;
    }
    else
    {
        Matrix tmp = a ^ (k / 2);
        ret = ret*tmp + ret;
    }
    return ret;
}
int main()
{
    int n, k, m, i, j;
    while (scanf("%d%d%d", &n, &k, &MOD) != EOF)
    {
        Matrix a(n);
        for (i = 0; i < n; ++i)
        {
            for (j = 0; j < n; ++j)
            {
                scanf("%d", &a.mat[i][j]);
                if (a.mat[i][j] >= MOD)
                    a.mat[i][j] %= MOD;
            }
        }
        Matrix ans = matrixSum(a, k);
        for (i = 0; i < n; ++i)
        {
            for (j = 0; j < n; ++j)
                j == n - 1 ? printf("%d
", ans.mat[i][j]) : printf("%d ", ans.mat[i][j]);
        }
        return 0;
    }
}