POJ 3233 Matrix Power Series（矩阵求A+…+A^k 的和）

最新推荐文章于 2023-08-28 09:28:46 发布

最新推荐文章于 2023-08-28 09:28:46 发布 · 121 阅读

本文介绍了一种使用矩阵快速幂解决特定数学问题的方法。通过两次二分法处理矩阵的幂次求和，针对不同奇偶性的幂次进行优化。代码实现了矩阵运算的基本操作，并提供了完整的算法实现。

题目链接：Click here~~

题意：

给你一个n*n的矩阵 A，求出 A + A^2 + … + A^k 的值。

解题思路：

需用用到两次二分。

考虑 k 为奇数和偶数两种情况。

1、奇数：A + A^2 + … + A^k = 【A + A^2 + … + A^(k/2)】 + A^(k/2+1) + A^(k/2+1)*【A + A^2 + … + A^(k/2)】。

2、偶数：A + A^2 + … + A^k = 【A + A^2 + … + A^(k/2)】 + A^(k/2)*【A + A^2 + … + A^(k/2)】。

其中，A^(k/2) 和 A^(k/2+1) 与【A + A^2 + … + A^k 】的值可以通过二分算出。

#include <stdio.h>
#include <string.h>

typedef __int64 LL;

#define N 32

int K,mod;

struct Matrix
{
    int r,c;
    int m[N][N];
    Matrix(){}
    Matrix(int r,int c):r(r),c(c){}
    Matrix operator *(const Matrix& B)
    {
        Matrix T(r,B.c);
        for(int i=1;i<=T.r;i++)
        {
            for(int j=1;j<=T.c;j++)
            {
                int tt = 0;
                for(int k=1;k<=c;k++)
                    tt += m[i][k]*B.m[k][j] % mod;
                T.m[i][j] = tt % mod;
            }
        }
        return T;
    }
    Matrix operator +(const Matrix& B)
    {
        Matrix T(r,c);
        for(int i=1;i<=T.r;i++)
            for(int j=1;j<=T.c;j++)
                T.m[i][j] = (m[i][j]+B.m[i][j]) % mod;
        return T;
    }
    Matrix Unit(int h)
    {
        Matrix T(h,h);
        memset(T.m,0,sizeof(T.m));
        for(int i=1;i<=h;i++)
            T.m[i][i] = 1;
        return T;
    }
    Matrix Pow(int n)
    {
        Matrix P = *this,Res = Unit(r);
        while(n)
        {
            if(n&1)
                Res = Res*P;
            P = P*P;
            n >>= 1;
        }
        return Res;
    }
    void Init()
    {
        int h;
        scanf("%d%d%d",&h,&K,&mod);
        r = c = h;
        for(int i=1;i<=r;i++)
            for(int j=1;j<=c;j++)
                scanf("%d",&m[i][j]);
    }
    void Print()
    {
        for(int i=1;i<=r;i++)
        {
            for(int j=1;j<=c;j++)
                printf("%d ",m[i][j]);
            printf("\n");
        }
    }
}Single;

Matrix solve(int k)
{
    if(k == 1)
        return Single;
    Matrix Sum = solve(k/2);
    if(k&1)
    {
        Matrix Ak_1 = Single.Pow(k/2+1);
        return Sum + Ak_1 + Ak_1*Sum;
    }
    else
    {
        Matrix Ak = Single.Pow(k/2);
        return Sum + Ak*Sum;
    }
}

int main()
{
    Single.Init();
    Matrix Ans(Single.r,Single.c);
    Ans = solve(K);
    Ans.Print();
    return 0;
}