UVA 11149 Power of Matrix 快速幂

题目链接：

http://acm.hust.edu.cn/vjudge/contest/122094#problem/G

Power of Matrix

Time Limit:3000MSMemory Limit:0KB

问题描述

给你一个矩阵A，求A+A^2+A^3+...+A^k

输入

Input consists of no more than 20 test cases. The first line for each case contains two positive integers n
(≤ 40) and k (≤ 1000000). This is followed by n lines, each containing n non-negative integers, giving
the matrix A.
Input is terminated by a case where n = 0. This case need NOT be processed.

输出

For each case, your program should compute the matrix A + A2 + A3 + . . . + Ak
. Since the values may
be very large, you only need to print their last digit. Print a blank line after each case.

样例

sample input
3 2
0 2 0
0 0 2
0 0 0
0 0

sample output
0 2 4
0 0 2
0 0 0

题解

A+A^2+...A^k=(I+A^(n/2))(A+...+A^(n/2))=...
一直递推下去，深度只有logn，只要logn次快速幂求A^x。

代码

WA四次的代码：

#include<iostream>
#include<cstring>
#include<cstdio>
using namespace std;

const int maxn = 55;
const int mod = 10;
typedef int LL;


struct Matrix {
    LL mat[maxn][maxn];
    Matrix() { memset(mat, 0, sizeof(mat)); }
    friend Matrix operator *(const Matrix& A, const Matrix& B);
    friend Matrix operator +(const Matrix &A,const Matrix &B);
    friend Matrix operator ^(Matrix A, int n);
};

Matrix I;

Matrix operator +(const Matrix& A, const Matrix& B) {
    Matrix ret;
    for (int i = 0; i < maxn; i++) {
        for (int j = 0; j < maxn; j++) {
            ret.mat[i][j] = (A.mat[i][j] + B.mat[i][j])%mod;
        }
    }
    return ret;
}

Matrix operator *(const Matrix& A, const Matrix& B) {
    Matrix ret;
    for (int i = 0; i < maxn; i++) {
        for (int j = 0; j < maxn; j++) {
            for (int k = 0; k < maxn; k++) {
                ret.mat[i][j] = (ret.mat[i][j]+A.mat[i][k] * B.mat[k][j]) % mod;
            }
        }
    }
    return ret;
}

Matrix operator ^(Matrix A, int n) {
    Matrix ret=I;
    while (n) {
        if (n & 1) ret = ret*A;
        A = A*A;
        n /= 2;
    }
    return ret;
}

Matrix solve(Matrix A, int n) {
    if (!n) return I;
    if (n == 1) return A;
    //这里要加括号！！！  ret=I+(A^(n/2))!!!
    Matrix ret = I + A ^ (n / 2);
    ret = ret*solve(A, n / 2);
    //这里也要加！！！！ ret=ret+(A^n)!!!
    if (n % 2) ret = ret + A^n;
    return ret;
}

int n, k;

int main() {
    for (int i = 0; i < maxn; i++) I.mat[i][i] = 1;
    while (scanf("%d%d", &n, &k) == 2 && n) {
        Matrix A;
        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=solve(A, k);
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n-1; j++) {
                printf("%d ",ans.mat[i][j]);
            }
            printf("%d\n", ans.mat[i][n - 1]);
        }
        printf("\n");
    }
    return 0;
}

保险一些的写法：

#include<iostream>
#include<cstring>
#include<cstdio>
using namespace std;

const int maxn = 55;
const int mod = 10;
typedef int LL;


struct Matrix {
    LL mat[maxn][maxn];
    Matrix() { memset(mat, 0, sizeof(mat)); }
    friend Matrix operator *(const Matrix& A, const Matrix& B);
    friend Matrix operator +(const Matrix &A,const Matrix &B);
    friend Matrix pow(Matrix A, int n);
};

Matrix I;

Matrix operator +(const Matrix& A, const Matrix& B) {
    Matrix ret;
    for (int i = 0; i < maxn; i++) {
        for (int j = 0; j < maxn; j++) {
            ret.mat[i][j] = (A.mat[i][j] + B.mat[i][j])%mod;
        }
    }
    return ret;
}

Matrix operator *(const Matrix& A, const Matrix& B) {
    Matrix ret;
    for (int i = 0; i < maxn; i++) {
        for (int j = 0; j < maxn; j++) {
            for (int k = 0; k < maxn; k++) {
                ret.mat[i][j] = (ret.mat[i][j]+A.mat[i][k] * B.mat[k][j]) % mod;
            }
        }
    }
    return ret;
}

Matrix pow(Matrix A, int n) {
    Matrix ret=I;
    while (n) {
        if (n & 1) ret = ret*A;
        A = A*A;
        n /= 2;
    }
    return ret;
}

Matrix solve(Matrix A, int n) {
    if (!n) return I;
    if (n == 1) return A;
    Matrix ret = I + pow(A,n/2);
    ret = ret*solve(A, n / 2);
    if (n % 2) ret = ret + pow(A,n);
    return ret;
}

int n, k;

int main() {
    for (int i = 0; i < maxn; i++) I.mat[i][i] = 1;
    while (scanf("%d%d", &n, &k) == 2 && n) {
        Matrix A;
        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=solve(A, k);
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n-1; j++) {
                printf("%d ",ans.mat[i][j]);
            }
            printf("%d\n", ans.mat[i][n - 1]);
        }
        printf("\n");
    }
    return 0;
}

转载于:https://www.cnblogs.com/fenice/p/5679392.html