POJ 3233 Matrix Power Series——快速幂&&等比&&分治

题目

给定一个 $n imes n$  的矩阵 $A$ 和正整数 $k$ 和 $m$。求矩阵 $A$ 的幂的和。

$$S = A + A^2 + ... + A^k$$

输出 $S$ 的各个元素对 $M$ 取余后的结果（$1 leq n leq 30, 1 leq k leq  10^9, 1 leq M leq 10^4$）。

分析

数据范围 $n$ 很小，$k$ 很大，不肯能逐一求得。

由于具有等比性质，

设  $S_k = I + A + ... + A^{k-1}$

则有

$$egin{pmatrix} A^k\  S_k end{pmatrix} = egin{pmatrix} A & 0\  I & I end{pmatrix} egin{pmatrix} A^{k-1}\  S_{k-1} end{pmatrix} = egin{pmatrix} A & 0\  I & I end{pmatrix}^kegin{pmatrix} I\  0 end{pmatrix}$$

对这个新矩阵使用快速幂即可。

代码中的输出，使用了分块矩阵乘法的性质进行了简化。

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

typedef long long ll;
struct matrix
{
    int r, c;
    int mat[61][61];
    matrix(){
        memset(mat, 0, sizeof(mat));
    }
};
int n, k, p;

matrix mul(matrix A, matrix B)   //矩阵相乘
{
    matrix ret;
    ret.r = A.r; ret.c = B.c;
    for(int i = 0;i < A.r;i++)
        for(int k = 0;k < A.c;k++)
            for(int j = 0;j < B.c;j++)
            {
                ret.mat[i][j] = (ret.mat[i][j] + A.mat[i][k] * B.mat[k][j]) % p;
            }
    return ret;
}

matrix mpow(matrix A, int n)
{
    matrix ret;
    ret.r = A.r; ret.c = A.c;
    for(int i = 0;i < ret.r;i++)  ret.mat[i][i] = 1;
    while(n)
    {
        if(n & 1)  ret = mul(ret, A);
        A = mul(A, A);
        n >>= 1;
    }
    return ret;
}

int  main()
{
    scanf("%d%d%d", &n, &k, &p);
    matrix a, b;
    a.r = a.c = n;
    for(int i = 0;i < n;i++)  for(int j = 0;j < n;j++)  scanf("%d", &a.mat[i][j]);
    b.r = b.c = 2*n;
    for(int i = 0;i < n;i++)
    {
        for(int j = 0;j < n;j++)  b.mat[i][j] = a.mat[i][j];
        b.mat[n+i][i] = b.mat[n+i][n+i] = 1;
    }
    b = mpow(b, k+1);
    for(int i = 0;i < n;i++)
        for(int j = 0;j < n;j++)
        {
            int tmp = b.mat[n+i][j] % p;
            if(i == j)  tmp = (tmp + p - 1) % p;
            printf("%d%c", tmp, j == n-1 ? '
' : ' ');
        }
}

还有一种经典的分治方法，

例如，

$A+A^2+A^3+A^4 = (A+A^2) + A^2(A + A^2)$，

$A+A^2+A^3+A^4+A^5 = (A+A^2) +A^3 +  A^3(A + A^2)$.

因此，分k的奇偶递归一下就可以了。

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

typedef long long ll;
struct matrix
{
    int r, c;
    int mat[61][61];
    matrix(){
        memset(mat, 0, sizeof(mat));
    }
};
int n, k, p;

matrix mul(matrix A, matrix B)   //矩阵相乘
{
    matrix ret;
    ret.r = A.r; ret.c = B.c;
    for(int i = 0;i < A.r;i++)
        for(int k = 0;k < A.c;k++)
            for(int j = 0;j < B.c;j++)
            {
                ret.mat[i][j] = (ret.mat[i][j] + A.mat[i][k] * B.mat[k][j]) % p;
            }
    return ret;
}

matrix mpow(matrix A, int n)
{
    matrix ret;
    ret.r = A.r; ret.c = A.c;
    for(int i = 0;i < ret.r;i++)  ret.mat[i][i] = 1;
    while(n)
    {
        if(n & 1)  ret = mul(ret, A);
        A = mul(A, A);
        n >>= 1;
    }
    return ret;
}

matrix add(matrix A, matrix B)
{
    matrix ret;
    ret.r = A.r; ret.c = A.c;
    for(int i = 0;i < A.r;i++)
        for(int j = 0;j < A.c;j++)
            ret.mat[i][j] = (A.mat[i][j] + B.mat[i][j]) % p;
    return ret;
}

matrix sum(matrix x, int k)  //A+A^2+..+A^k
{
    if(k == 1)  return x;
    matrix s = sum(x, k/2);
    if(k & 1)
    {
        matrix tmp = mpow(x, k/2+1);
        return add(s, add(tmp, mul(tmp, s)));
    }
    else
    {
        matrix tmp = mpow(x, k/2);
        return add(s, mul(tmp, s));
    }
}

int  main()
{
    scanf("%d%d%d", &n, &k, &p);
    matrix a, ans;
    a.r = a.c = n;
    for(int i = 0;i < n;i++)  for(int j = 0;j < n;j++)  scanf("%d", &a.mat[i][j]);
    ans = sum(a, k);
    for(int i = 0;i < n;i++)  for(int j = 0;j < n;j++)  printf("%d%c", ans.mat[i][j], j == n-1 ? '
' : ' ');
}

这个时间复杂度咋算啊？知道的大犇请留言。

参考链接：

1. https://blog.csdn.net/rowanhaoa/article/details/21023599

2. https://blog.csdn.net/scf0920/article/details/39345197