poj 3233 Matrix Power Series

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

Input

　　The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 10⁹) and m (m < 10⁴). 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

---

　　这道题直接暴力呢。。$Oleft(n^{3}k ight)$的时间复杂度，必死无疑。所以只能另寻出路，注意到矩阵满足分配律

矩阵乘法的右分配律 若$A$和$B$皆为$n$行$p$列的矩阵，$C$为$p$行$m$列的矩阵，则有

$(A + B)C = AB + AC$

　　左分配律差不多（因为矩阵乘法不满足交换律，所以矩阵乘法的分配律分为左分配律和右分配律）

　　证明就水水的写一下

　　对于$A + B$，我们有$left(A + B ight )_{ij} = A_{ij} + B_{ij}$

　　那么

$left[left(A + B ight )C ight ]_{ij}\=sum_{i = 1}^{p}left(A + B ight )_{ip}C_{pj}\=sum_{i = 1}^{p}A_{ip}C_{pj} + sum_{i = 1}^{p}B_{ip}C_{pj}\=left(AC + BC ight )_{ij}$

　　所以$(A + B)C = AB + AC$

　　因为有了分配律这个神奇东西，现在计算$A + A^{2} + cdots + A^{6}$的和，就等价于$left(A + A^{2} + A^{3} ight ) + left(A + A^{2} + A^{3} ight )A^{3}$

　　现在设$fleft(n ight ) = sum_{i = 1}^{n}A^{i}$

　　那么有

$fleft ( n ight )=left{egin{matrix}A left(n = 1 ight )\fleft(n - 1 ight ) + A^{n} left(2 mid n ight ) \ fleft(frac{n}{2} ight) left(I + A^{frac{n}{2}} ight ) left(2 mid n ight )  end{matrix} ight.$

　　时间复杂度是$Oleft(n^{3}log^{2}k ight )$，但是如果代码写得比较丑就有TLE的风险，可以考虑一下以下优化

1. 不要没事就作无用的初始化，很耗时间]
2. 降低取模次数
3. 适当地使用取模优化
4. 减少无用的内存拷贝的次数

Code

/**
 * poj
 * Problem#3233
 * Accepted
 * Time:2266ms
 * Memory:3892k 
 */
#include<iostream>
#include<cstdio>
#include<cctype>
#include<cstring>
#include<cstdlib>
#include<fstream>
#include<sstream>
#include<algorithm>
#include<map>
#include<set>
#include<queue>
#include<vector>
#include<stack>
using namespace std;
typedef bool boolean;
#define INF 0xfffffff
#define smin(a, b) a = min(a, b)
#define smax(a, b) a = max(a, b)
template<typename T>
inline void readInteger(T& u){
    char x;
    int aFlag = 1;
    while(!isdigit((x = getchar())) && x != '-');
    if(x == '-'){
        x = getchar();
        aFlag = -1;
    }
    for(u = x - '0'; isdigit((x = getchar())); u = (u << 1) + (u << 3) + x - '0');
    ungetc(x, stdin);
    u *= aFlag;
}

typedef class Matrix {
    public:
        int* p;
        int lines, cols;
        int moder;
        Matrix():p(NULL), lines(0), cols(0), moder(1)    {    } 
        Matrix(int lines, int cols, int moder):lines(lines), cols(cols), moder(moder) {
            p = new int[(lines * cols)];
        }
        
        int* operator [](int pos) {
            return p + pos * cols;
        }
        
        Matrix operator *(Matrix b) {
            Matrix res(lines, b.cols, moder);
            for(int i = 0; i < lines; i++) {
                for(int j = 0; j < b.cols; j++) {
                    res[i][j] = 0;
                    for(int k = 0; k < cols; k++) {
                        (res[i][j] += (*this)[i][k] * b[k][j] % moder) %= moder;
                    }
                }
            }
            return res;
        }
        
        Matrix operator +(Matrix b) {
            Matrix res(lines, cols, moder);
            for(int i = 0; i < lines; i++)
                for(int j = 0; j < cols; j++)
                    res[i][j] = ((*this)[i][j] + b[i][j]) % moder;
            return res;
        }
}Matrix;

template<typename T>
T pow(T a, int pos) {
    if(pos == 1)    return a;
    T temp = pow(a, pos / 2);
    if(pos & 1)    return temp * temp * a;
    return temp * temp;
}

int n, k, m;
Matrix a;

inline void init() {
    readInteger(n);
    readInteger(k);
    readInteger(m);
    a = Matrix(n, n, m);
    for(int i = 0; i < n; i++) {
        for(int j = 0; j < n; j++) {
            readInteger(a[i][j]);
            a[i][j] %= m;
        }
    }
}

template<typename T>
T pow_sum(T a, int pos) {
    if(pos == 1)    return a;
    T temp = pow_sum(a, pos / 2);
    temp = temp + temp * pow(a, pos / 2);
    if(pos & 1)    return temp + pow(a, pos);
    return temp;
}

Matrix res;
inline void solve() {
    res = pow_sum(a, k);
    for(int i = 0; i < n; i++) {
        for(int j = 0; j < n; j++) {
            printf("%d ", res[i][j]);
        }
        putchar('
');
    }
}

int main() {
    init();
    solve();
    return 0;
}