POJ3233 Matrix Power Series 矩阵乘法

最新推荐文章于 2020-12-19 11:24:01 发布

Albafica

最新推荐文章于 2020-12-19 11:24:01 发布

阅读量485

点赞数

分类专栏：矩阵乘法数学

本文链接：https://blog.csdn.net/lenleaves/article/details/8865906

版权

数学同时被 2 个专栏收录

38 篇文章 0 订阅

订阅专栏

矩阵乘法

13 篇文章 0 订阅

订阅专栏

这道题目是第三种矩阵乘法的应用。

然我们求 S = A + A ² + A ³ + … + A^k.的结果矩阵

M67大牛的博客上面讲的是一种分治的方法。

设f[n]=A^1+A^2+....A^n;
当n是偶数,f[n]=f[n/2]+f[n/2]*A^(n/2);
但n是奇数,f[n]=f[n-1]+A^(n);

而看了DISCUSS之后发现一种编程复杂度更加小的方法。再次涨了姿势。。

Let B=[ [A I];
[0 I] ]

B^(k+1) = [ [A^k I+A+...+A^k];
[0 I ] ]

只能说矩阵真是太神奇了。。。。。。。

我比较懒。所以。。。实现了第二种。。。。

Matrix Power Series

Time Limit: 3000MS		Memory Limit: 131072K
Total Submissions: 11648		Accepted: 4977

Description

Given a n × n matrix A and a positive integer k, find the sumS = A + A² +A³ + … + A^k.

Input

The input contains exactly one test case. The first line of input contains three positive integersn (n ≤ 30), k (k ≤ 10⁹) andm (m < 10⁴). Then follown lines each containing n nonnegative integers below 32,768, givingA’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

#include<iostream>
#include<cstdio>

#define MAXN 100

using namespace std;

struct Matrix
{
    int size;
    long modulo;
    long element[MAXN][MAXN];
    void setSize(int);
    void setModulo(int);
    Matrix operator* (Matrix);
    Matrix power(int);

};

void Matrix::setSize(int a)
{
    for (int i=0; i<a; i++)
        for (int j=0; j<a; j++)
            element[i][j]=0;
    size = a;
}

void Matrix::setModulo(int a)
{
    modulo = a;
}

Matrix Matrix::operator* (Matrix param)
{
    Matrix product;
    product.setSize(size);
    product.setModulo(modulo);
    for (int i=0; i<size; i++)
        for (int j=0; j<size; j++)
            for (int k=0; k<size; k++)
            {
                product.element[i][j]+=element[i][k]*param.element[k][j];
                product.element[i][j]%=modulo;
            }
    return product;
}

Matrix Matrix::power(int exp)
{
    Matrix tmp = (*this) * (*this);
    if (exp==1) return *this;
    else if (exp & 1) return tmp.power(exp/2) * (*this);
    else return tmp.power(exp/2);
}

int n,m,k;
Matrix m1,ans;
int main()
{
    while(~scanf("%d%d%d",&n,&k,&m))
    {
        m1.setModulo(m);
        m1.setSize(n*2);
        for(int i=0;i<n;i++)
        {
            for(int j=0;j<n;j++)
            {
                scanf("%d",&m1.element[i][j]);
                m1.element[i][j]%=m;
            }
            m1.element[i][i+n]=1;
            m1.element[n+i][n+i]=1;
        }
        m1.element[n][n]=1;
        ans=m1.power(k+1);
        for(int i=0;i<n;i++)
            for(int j=n;j<n*2;j++)
            {
                int out=ans.element[i][j];
                if(j==i+n)
                    out=out==0?m-1:out-1;
                printf("%d%s",out,j!=n*2-1?" ":"\n");
            }
        return 0;

    }
    return 0;
}