Matrix Power Series(行列高速べき乗)

行列の高速べき乗:http://www.cnblogs.com/kuangbin/archive/2012/08/17/2643347.html
Matrix Power Series
Time Limit：3000 MS
 
メモリLimit:13131313722 K
Total Submissions：16341
 
Acceepted:6966
Description
Given a n× n matix A and a positive integer k，find the sumS=A 2+A 3+…+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）.The n follown line eas each containing n nonegative integers belogint belogint
Output
Output the elemens 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

簡単で分かりやすいです。
二次二次マトリックス、具体的にコードを見ましょう。
ACコード:

#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>
#include <vector>
#include <cmath>
#include <queue>
#include <map>
#include <set>
#define eps 1e-9

using namespace std;
typedef long long ll;
typedef pair<int,int>P;
const int M = 1e5 + 100;
const int INF = 0x3f3f3f3f;
int n,mod,k;

struct Mac {
    int c[33][33]; //    
    void unit1() { //    
        for(int i = 0; i < n; i++) {
            for(int j = 0; j < n; j++) {
                scanf("%d",&c[i][j]);
            }
        }
    }
    void unit2() { //    
        c[0][0] = c[1][1] = 1;
        c[1][0] = c[0][1] = 0;
    }
};
Mac tmp,sum; //    ,     

Mac mul(Mac a,Mac b) { //    
    Mac p;
    memset(p.c,0,sizeof(p.c));
    for(int i = 0; i < n; i++) {
        for(int j = 0; j < n; j++) {
            if(a.c[i][j]) {
                for(int k = 0; k < n; k++) {
                    p.c[i][k] += (a.c[i][j] * b.c[j][k]) % mod;
                    p.c[i][k] %= mod;
                }
            }
        }
    }
    return p;
}

Mac add(Mac a,Mac b) { //    
    for(int i = 0; i < n; i++) {
        for(int j = 0; j < n; j++) {
            a.c[i][j] += b.c[i][j] % mod;
            a.c[i][j] %= mod;
        }
    }
    return a;
}

Mac quickmod(Mac a,int n) { //     
    Mac p;
    p.unit2();
    while(n) {
        if(n & 1) p = mul(p,a);
        a = mul(a,a);
        n >>= 1;
    }
    return p;
}

Mac binary_matrix1(int k) { //       
    Mac res;
    if(k == 1) {
        return tmp;
    }
    res = binary_matrix1(k / 2);
    res = mul(res,res);
    if(k & 1) res = mul(res,tmp);
    return res;
}

Mac binary_matrix2(int k){ //           A1+A2+A3...;
    if(k == 1){
        return tmp;
    }
    Mac res = binary_matrix2(k / 2);
    Mac tp = binary_matrix1(k / 2);
    tp = mul(tp,res);
    sum = add(sum,tp);
    if(k & 1) sum = add(sum,binary_matrix1(k));
    return sum;
}


int main() {
    cin>>n>>k>>mod;
    tmp.unit1();
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            sum.c[i][j] = tmp.c[i][j];
        }
    }
    Mac res = binary_matrix2(k);
    for(int i = 0; i < n; i++) {
        for(int j = 0; j < n; j++) {
            printf(j == n - 1 ? "%d
" : "%d ",res.c[i][j]);
        }
    }
    return 0;
}