【POJ 3233】Matrix Power Series（矩阵快速幂+等比数列二分求和）

Matrix Power Series

Time Limit: 3000MS		Memory Limit: 131072K
Total Submissions: 29484		Accepted: 11946

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 109) and m (m < 104). Then follow n lines each containing n nonnegative integers below 32,768, giving A’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

Source

​​POJ Monthly--2007.06.03​​, Huang, Jinsong

 

题意：

给定矩阵A。以及k, p ,求 A+A^2+A^3+......A^k    的和对p取余

分析：

当k为偶数时：

比方 k=6,那么  A+A^2+A^3+A^4+A^5+A^6=    A+A^2+A^3+   A^3*（A+A^2+A^3）   

s(k)=s(k/2)+A^(n/2) * s(k/2) 即s(k)=(E+A^(n/2))*s(n/2)  （E为单位矩阵）

当k为奇数时：

s(k)=s(k-1)+A^k ,  那么k-1为偶数。能够依照上面的二分

还有种做法：

//矩阵快速幂注意特判
#include<iostream>
#include<stdio.h>
#include<algorithm>
#include<iomanip>
#include<cmath>
#include<cstring>
#include<vector>
#define N 10005  
using namespace std;
const int MAXN= 35;
 int p;
typedef long long ll;
int len=35;//构造矩阵的长度
struct Matrix//注意定义的是方阵，*法也是
{
   int M[MAXN][MAXN];
   Matrix(const bool I = 0)  ///初始化对角矩阵
   {
      memset(M, 0, sizeof(M));
      if (I)
       for(int i = 0; i < len; i++) 
             M[i][i] = 1;
   }
   Matrix operator *(const Matrix &y) ///矩阵乘法，对乘法重新定义，z=x*y
   {
      Matrix z;
      for (int i = 0; i < len; i++)//x矩阵的n
         for (int j = 0; j < len; j++)//y矩阵的n
            for (int k = 0; k < len; k++)//xy矩阵的m
               z.M[i][j] = (z.M[i][j]+M[i][k]*y.M[k][j])%p;
      return z;
   }
};

Matrix Pow(Matrix A, long long b)///矩阵的快速幂
{
       Matrix ans = Matrix(1);
       for (; b; b>>=1)
       {
        if (b&1) ans = ans*A;
        A = A*A;
       }
      return ans; 
}
Matrix Add(Matrix a,Matrix b)  //a+b
{
    Matrix c;
    for(int i = 0; i < len; ++i)
    {
        for(int j = 0; j < len; ++j)
        {
            c.M[i][j] =  (a.M[i][j] + b.M[i][j]) % p;
        }
    }
    return c;
}
Matrix MatrixSum(Matrix a,int k)//a + a^2 + a^3 + … + a^k
//原理：等比数列二分求和
{
    if(k == 1)
        return a;

    if(k&1)
    {
      return Add(MatrixSum(a,k-1),Pow(a,k));  //当n是奇数,f[n]=f[n-1]+A^(n);
    }
    else
    {
      Matrix E(1);
      return MatrixSum(a,k>>1)*Add(Pow(a,k>>1),E); //当n是偶数,f[n]=f[n/2]*(A^(n/2)+E);
    }
    
}
int main()  
{  
    int n,k;
   scanf("%d%d%d",&n,&k,&p);
   len=n;///构造矩阵的长度
       Matrix a;
      
      for(int i = 0; i < n; ++i)
        for(int j = 0; j < n; ++j)
        {
            scanf("%d",&a.M[i][j]);
            
        }
        //a.out();
        Matrix ans = MatrixSum(a,k);

        for(int i = 0; i < n; ++i)
        {
            for(int j = 0; j < n-1; ++j)
            {
                printf("%d ",ans.M[i][j]);
            }
            printf("%d\n",ans.M[i][n-1]);
        }
    return 0;  
}