POJ 3233 Matrix Power Series（矩陣快速冪）

POJ 3233 Matrix Power Series（矩陣快速冪）

Description

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

Source

題意：給出一個n*n的矩陣A，求A+A^2+A^3+……+A^k mod m的結果是多少？
方法一：兩次二分
Sk = A + A2 + A3 + … + Ak   
    =(1+A^(k/2))*(A + A^2 + A^3 + … + A^(k/2)  )+{A^k}
    =(1+A^(k/2))*(S(k/2) )+ {Ak} 當k為偶數時沒有{Ak}
即
k%2==0: S[k]=F[k/2] (1+A[k/2]);
k%2==1: S[k]=F[k-1]+A[k];
#include 
#include 
#include 
using namespace std;

const int N = 32;
struct Matrix {
    int mat[N][N];
    Matrix() {
        memset(mat, 0, sizeof(mat));
        for(int i = 0; i < N; i++)
            mat[i][i] = 1;
    }
} E;
int n, k, mod;

Matrix Multi(Matrix a, Matrix b) {
    Matrix res;
    for(int i = 0; i < n; i++) {
        for(int j = 0; j < n; j++) {
            res.mat[i][j] = 0;
            for(int k = 0; k < n; k++)
                res.mat[i][j] = (res.mat[i][j] + a.mat[i][k] * b.mat[k][j]) % mod;
        }
    }
    return res;
}

Matrix Add(Matrix a, Matrix b) {
    Matrix res;
    for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++)
            res.mat[i][j] = (a.mat[i][j] + b.mat[i][j]) % mod;
    return res;
}

Matrix Pow(Matrix a, int n) {
    Matrix res;
    while(n) {
        if(n&1) res = Multi(res, a);
        a = Multi(a, a);
        n >>= 1;
    }
    return res;
}

Matrix Get_Ans(Matrix a, int k) {
    if(k == 1) return a;
    if(k&1) return Add(Pow(a, k), Get_Ans(a, k-1));
    if(k % 2 == 0) {
        Matrix A = Get_Ans(a, k/2);
        Matrix B = Pow(a, k/2);
        Matrix C = Multi(A, B);
        return Add(C, A);
    }
}

int main() {
    Matrix A;
    while(~scanf("%d%d%d", &n, &k, &mod)) {
        for(int i = 0; i < n; i++)
            for(int j = 0; j < n; j++) {
                scanf("%d", &A.mat[i][j]);
                A.mat[i][j] %= mod;
            }
        Matrix ans = Get_Ans(A, k);
        for(int i = 0; i < n; i++) {
            for(int j = 0; j < n; j++) {
                if(j) printf(" ");
                printf("%d", ans.mat[i][j]);
            }
            printf("\n");
        }
    }
    return 0;
}


方法二：
構造一個新矩陣B = | A   E |
                  | 0   E |

則 B^(k+1) = | A^(k+1)      A^k+A^(k-1)+……+A^2 + A + 1 |
             |     0                 E                   |
所以可以先直接用矩陣快速冪求出B^(k+1),然後用左上角的那一部分減去單位矩陣就是最後要求的矩陣。
#include 
#include 
#include 
using namespace std;
const int N = 61;
int n, mod, k;
struct Matrix {
    int mat[N][N];
    Matrix() {
        memset(mat, 0, sizeof(mat));
        for(int i = 0; i < N; i++)
            mat[i][i] = 1;
    }
};

Matrix Multi(Matrix a, Matrix b) {
    Matrix res;
    memset(res.mat, 0, sizeof(res.mat));
    for(int i = 0; i < n * 2; i++) {
        for(int k = 0; k < n * 2; k++) {
            if(a.mat[i][k]) {
                for(int j = 0; j < n * 2; j++) {
                    res.mat[i][j] += a.mat[i][k] * b.mat[k][j];
                    res.mat[i][j] %= mod;
                }
            }
        }
    }
    return res;
}

Matrix Pow(Matrix x, int m) {
    Matrix res;
    while(m) {
        if(m&1) res = Multi(res, x);
        x = Multi(x, x);
        m >>= 1;
    }
    return res;
}

int main() {
    Matrix A;
    while(~scanf("%d%d%d",&n, &k, &mod)) {
        memset(A.mat, 0, sizeof(A.mat));
        for(int i = 0; i < n; i++) {
            for(int j = 0; j < n; j++) {
                scanf("%d", &A.mat[i][j]);
                A.mat[i][j] %= mod;
            }
            A.mat[i][i+n] = 1;
        }
        for(int i = n; i < 2 * n; i++)
            A.mat[i][i] = 1;
        Matrix ans = Pow(A, k + 1);
        for(int i = 0; i < n; i++) {
            for(int j = n; j < n * 2; j++) {
                if(j > n) printf(" ");
                if(j - i == n)
                    printf("%d", (ans.mat[i][j] - 1 + mod) % mod);
                else
                    printf("%d", ans.mat[i][j]);
            }
            printf("\n");
        }
    }
    return 0;
}