POJ 3233 Matrix Power Series 矩阵快速幂+二分

最新推荐文章于 2020-01-21 15:40:01 发布

原创最新推荐文章于 2020-01-21 15:40:01 发布 · 671 阅读

0 ·

CC 4.0 BY-SA版权

ACM_数论专栏收录该内容

4 篇文章

订阅专栏

本文介绍了一种结合矩阵快速幂和二分法解决特定数学问题的方法，具体实现包括矩阵的基本运算如加法、乘法及求模，并通过一个示例程序详细展示了如何运用这些方法来高效求解矩阵的幂次。

对数字的处理:

快速幂模+二分求和求(∑A^b)%C

对矩阵的处理同理；

只是自己定义一些函数来表示 + * ^ %

/*
 * POJ 3233
 * fuqiang
 * 矩阵快速幂+二分
 * 2013/7/31
*/
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cstdlib>
using namespace std;
#define maxn 30+3
int n, k, m;
struct Matrix
{
    int val[maxn][maxn];
    void unit()  //单位矩阵
    {
        for(int i = 0; i < maxn; i++) val[i][i] = 1;
    }
    void zero()  //零矩阵
    {
        memset(val, 0, sizeof(val));
    }
} x;
Matrix operator %(const Matrix &a, const int &m)  //矩阵求模
{
    Matrix temp;
    for(int i = 1; i <= n; i++)
        for(int j = 1; j <= n; j++)
            temp.val[i][j] = a.val[i][j] % m;
    return temp;
}

Matrix operator *(const Matrix &a, const Matrix &b)  //矩阵乘法
{
    Matrix tmp;
    tmp.zero();
    for(int k = 1; k <= n; k++)
    {
        for(int i = 1; i <= n; i++)
            if(a.val[i][k])
                for(int j = 1; j <= n; j++)
                {
                    tmp.val[i][j] += a.val[i][k] * b.val[k][j];
                }
    }
    return tmp%m;
}

Matrix operator ^(Matrix x, int n)  //矩阵快速幂
{
    Matrix tmp;
    tmp.zero();
    tmp.unit();
    while(n)
    {
        if(n & 1) tmp = tmp * x;
        x = x * x;
        n >>= 1;
    }
    return tmp;
}

Matrix operator +(const Matrix &a, const Matrix &b)  //矩阵加法
{
    Matrix tmp;
    for(int i = 1; i <= n; i++)
        for(int j = 1; j <= n; j++)
            tmp.val[i][j] = (a.val[i][j] + b.val[i][j]) % m;
    return tmp;
}

Matrix sum(Matrix x, int k,int m)
{
    if(k == 1) return x;
    else
    {
        Matrix tmp = sum(x, k>>1, m)%m;
        if(k & 1)
        {
            Matrix tmp2 = x ^ ((k >> 1) + 1);
            return (tmp + tmp2 + tmp * tmp2)%m;
        }
        else
        {
            Matrix tmp2 = x ^ (k >> 1);
            return (tmp + tmp * tmp2)%m;
        }
    }
}

int main()
{
//#ifndef ONLINE_JUDGE
//    freopen("in","r",stdin);
//#endif
    while(~scanf("%d%d%d", &n, &k, &m))
    {
        for(int i = 1; i <= n; i++)
        {
            for(int j = 1; j <= n; j++)
            {
                scanf("%d", &x.val[i][j]);
                x.val[i][j] %= m;
            }
        }
        Matrix ans = sum(x, k, m);
        for(int i = 1; i <= n; i++)
        {
            for(int j = 1; j <= n; j++)
                printf("%d ", ans.val[i][j]);
            printf("\n");
        }
        break;
    }
    return 0;
}

思想来源:

Matrix67 大神, 他的博客: http://www.matrix67.com/blog/