POJ 3233 && NYOJ 298 Matrix Power Series（矩阵快速幂）

Default

Matrix Power Series

Time Limit: 3000MS
Memory Limit: 131072KTotal Submissions: 15553
Accepted: 6658

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 40 11 1

Sample Output

1 22 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 <cstdio>#include <cstring>#include <algorithm>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 <cstdio>#include <cstring>#include <algorithm>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;}