poj 3233 Matrix Power Series

H - Matrix Power Series

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⁹) and m (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

 

 

题意：已知一个n*n的矩阵A，和一个正整数k，求S = A + A2 + A3 + … + Ak。

 

思路：矩阵快速幂。首先我们知道 A^x 可以用矩阵快速幂求出来（具体可见poj 3070）。其次可以对k进行二分，每次将规模减半，分k为奇偶两种情况，如当k = 6和k = 7时有：

      k = 10 有： S(10) = ( A^1+A^2+A^3+A^4+ A^5 ) + A^5 * ( A^1+A^2+A^3+A^4+A^5 ) = S(5) + A^5 * S(5)

      k = 5 有： S(5) = ( A^1+A^2 ) + A^3 + A^3 * ( A^1+A^2 ) = S(2) + A^3 + A^3 * S(2)

  k = 2 有 :  S(2) = A^1 + A^2 = S(1) + A^1 * S(1)

 

二分 递归 矩阵 快速幂

 

#include <cstdio>

#include <cstring>

using namespace std ;

struct node

{

    int a[50][50] ;

};

int n,k,m;

node mul(node p,node q)//矩阵乘法

{

    node s;

    int i,j,k;

    for(i=0;i<n;i++)

    {

        for(j=0;j<n;j++)

        {

            s.a[i][j]=0 ;

            for(k=0;k<n;k++)

            {

                s.a[i][j]=(s.a[i][j]+p.a[i][k]*q.a[k][j])%m;

            }

        }

    }

    return s ;

}

node add(node p,node q)//矩阵加法

{

    int i,j;

    node s;

    for(i=0;i<n;i++)

        for(j=0;j<n;j++)

            s.a[i][j]=(p.a[i][j]+q.a[i][j])%m ;

    return s ;

}

node pow(node p,int k)//矩阵快速幂

{

    if(k==1)

        return p;

    node s=pow(p,k/2);

    s=mul(s,s);

    if(k%2)

        s=mul(s,p);

    return s;

}

node f(node p,int k)//二分递归

{

    if(k==1)

        return p;

    node s=f(p,k/2),q,temp;

    int i,j;

    for(i=0;i<n;i++)

    {

        for(j=0;j<n;j++)

        {

            if(i==j)

                temp.a[i][j]=1;

            else

                temp.a[i][j]=0;

        }

    }////标准矩阵*任何矩阵都为1

    if(k%2)//奇数

    {

        q=pow(p,k/2+1) ;

        s=add(q,mul(s,add(q,temp))) ;

    }

    else

    {

        q=pow(p,k/2);

        s=mul(s,add(q,temp));

    }

    return s ;

}

int main()

{

    int i,j;

    node p,s;

    while(scanf("%d%d%d",&n,&k,&m)!=EOF)

    {

        for(i=0;i<n;i++)

            for(j=0;j<n;j++)

                scanf("%d",&p.a[i][j]);

        s=f(p,k);

        for(i=0;i<n;i++)

        {

            for(j=0;j<n;j++)

            {

                if(j==n-1)

                    printf("%d\n",s.a[i][j]);

                else

                    printf("%d ",s.a[i][j]);

            }

        }

    }

    return 0;

}