# [POJ 3233] Matrix Power Series【分治 / 矩乘优化】

• 2018-01-07
• 0
• 0

## Problem:

 Time Limit: 3000MS Memory Limit: 131072K

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

## Solution:

• 若 k 为偶数，则 S = A + A2 + A3 + … + Ak = (I + Ak/2) × (A + A2 + A3 + … + Ak/2)
• 若 k 为奇数，则 S = A + A2 + A3 + … + Ak = (I + A[k/2]) × (A + A2 + A3 + … + A[k/2]) + Ak

## Code #1: O(n3log2k) [1000K, 1735MS]

```#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<cassert>
#include<iostream>
#include<algorithm>
using namespace std;

int k, m;

struct Matrix{
int ele[32][32], size;

inline void setIdentity(int idensize){
memset(ele, 0, sizeof(ele));
size = idensize;
for(register int i = 1; i <= size; i++) ele[i][i] = 1;
}

inline Matrix operator + (const Matrix &mat2) const {
//	assert(size == mat2.size);
Matrix res; res.size = size;
for(register int i = 1; i <= size; i++)
for(register int j = 1; j <= size; j++)
res.ele[i][j] = (ele[i][j] + mat2.ele[i][j]) % m;
return res;
}

inline Matrix operator * (const Matrix &mat2) const {
//	assert(size == mat2.size);
Matrix res; res.size = size;
for(register int i = 1; i <= size; i++)
for(register int j = 1; j <= mat2.size; j++){
res.ele[i][j] = 0;
for(register int k = 1; k <= size; k++)
res.ele[i][j] = (res.ele[i][j] + ele[i][k] * mat2.ele[k][j]) % m;
}
return res;
}

inline Matrix pow(int ex){
Matrix bas = *this, res; res.setIdentity(size);
while(ex){
if(ex & 1) res = res * bas;
ex >>= 1, bas = bas * bas;
}
return res;
}
} A, I;

inline Matrix sum(int ex){  // Return the sum of the geometric sequence
if(ex == 1) return A;
if(ex & 1) return (I + A.pow(ex >> 1)) * sum(ex >> 1) + A.pow(ex);
return (I + A.pow(ex >> 1)) * sum(ex >> 1);
}

int main(){
scanf("%d%d%d", &A.size, &k, &m);
I.setIdentity(A.size);
for(register int i = 1; i <= A.size; i++)
for(register int j = 1; j <= A.size; j++)
scanf("%d", &A.ele[i][j]);
A = sum(k);
for(register int i = 1; i <= A.size; i++){
for(register int j = 1; j < A.size; j++)
printf("%d ", A.ele[i][j]);
printf("%d\n", A.ele[i][A.size]);
}
return 0;
}
```

## Code #2: O(n3logk) [372K, 719MS]

```#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<cassert>
#include<iostream>
#include<algorithm>
using namespace std;

int k, m;

struct Matrix{
int ele[32][32], size;

inline void setIdentity(int idsize, int x){
size = idsize;
for(register int i = 1; i <= size; i++)
for(register int j = 1; j <= size; j++)
if(i == j) ele[i][j] = x;
else ele[i][j] = 0;
}

inline Matrix operator + (const Matrix &mat2) const {
//	assert(size == mat2.size);
Matrix res; res.size = size;
for(register int i = 1; i <= size; i++)
for(register int j = 1; j <= size; j++)
res.ele[i][j] = (ele[i][j] + mat2.ele[i][j]) % m;
return res;
}

inline Matrix operator * (const Matrix &mat2) const {
//	assert(size == mat2.size);
Matrix res; res.size = size;
for(register int i = 1; i <= size; i++)
for(register int j = 1; j <= mat2.size; j++){
res.ele[i][j] = 0;
for(register int k = 1; k <= size; k++)
res.ele[i][j] = (res.ele[i][j] + ele[i][k] * mat2.ele[k][j]) % m;
}
return res;
}
} A, I, O;

struct matMatrix{
Matrix ele[3][3];
int size;

inline matMatrix operator * (const matMatrix &mat2) const {
//	assert(size == mat2.size);
matMatrix res; res.size = size;
for(register int i = 1; i <= size; i++)
for(register int j = 1; j <= mat2.size; j++){
res.ele[i][j].setIdentity(A.size, 0);
for(register int k = 1; k <= size; k++)
res.ele[i][j] = res.ele[i][j] + ele[i][k] * mat2.ele[k][j];
}
return res;
}

inline matMatrix pow(int ex){
matMatrix bas = *this, res = *this; ex--;
while(ex){
if(ex & 1) res = res * bas;
ex >>= 1, bas = bas * bas;
}
return res;
}
} B;

int main(){
scanf("%d%d%d", &A.size, &k, &m);
I.setIdentity(A.size, 1), O.setIdentity(A.size, 0);
for(register int i = 1; i <= A.size; i++)
for(register int j = 1; j <= A.size; j++)
scanf("%d", &A.ele[i][j]);
B.size = 2;
B.ele[1][1] = A, B.ele[1][2] = A;
B.ele[2][1] = O, B.ele[2][2] = I;
B = B.pow(k), A = B.ele[1][2];
for(register int i = 1; i <= A.size; i++){
for(register int j = 1; j < A.size; j++)
printf("%d ", A.ele[i][j]);
printf("%d\n", A.ele[i][A.size]);
}
return 0;
}
```

#### 评论

darkleafin.cf
(该域名已过期且被抢注。。)
darkleafin.github.io

49750

https://github.com/Darkleafin

OPEN AT 2017.12.10

Please refresh the page if the code cannot be displayed normally.

https://visualgo.net/en

- Theme by Qzhai