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

题目

给定一个 $n \times 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}$

则有

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

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

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

#include<cstdio>

#include<cstring>

using namespace std;

typedef long long ll;

struct matrix

{

    int r, c;

    int mat[][];

    matrix(){

        memset(mat, , 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 = ;i < A.r;i++)

        for(int k = ;k < A.c;k++)

            for(int j = ;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 = ;i < ret.r;i++)  ret.mat[i][i] = ;

    while(n)

    {

        if(n & )  ret = mul(ret, A);

        A = mul(A, A);

        n >>= ;

    }

    return ret;

}

int  main()

{

    scanf("%d%d%d", &n, &k, &p);

    matrix a, b;

    a.r = a.c = n;

    for(int i = ;i < n;i++)  for(int j = ;j < n;j++)  scanf("%d", &a.mat[i][j]);

    b.r = b.c = *n;

    for(int i = ;i < n;i++)

    {

        for(int j = ;j < n;j++)  b.mat[i][j] = a.mat[i][j];

        b.mat[n+i][i] = b.mat[n+i][n+i] = ;

    }

    b = mpow(b, k+);

    for(int i = ;i < n;i++)

        for(int j = ;j < n;j++)

        {

            int tmp = b.mat[n+i][j] % p;

            if(i == j)  tmp = (tmp + p - ) % p;

            printf("%d%c", tmp, j == n- ? '\n' : ' ');

        }

}

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

例如，

$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[][];

    matrix(){

        memset(mat, , 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 = ;i < A.r;i++)

        for(int k = ;k < A.c;k++)

            for(int j = ;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 = ;i < ret.r;i++)  ret.mat[i][i] = ;

    while(n)

    {

        if(n & )  ret = mul(ret, A);

        A = mul(A, A);

        n >>= ;

    }

    return ret;

}

matrix add(matrix A, matrix B)

{

    matrix ret;

    ret.r = A.r; ret.c = A.c;

    for(int i = ;i < A.r;i++)

        for(int j = ;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 == )  return x;

    matrix s = sum(x, k/);

    if(k & )

    {

        matrix tmp = mpow(x, k/+);

        return add(s, add(tmp, mul(tmp, s)));

    }

    else

    {

        matrix tmp = mpow(x, k/);

        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 = ;i < n;i++)  for(int j = ;j < n;j++)  scanf("%d", &a.mat[i][j]);

    ans = sum(a, k);

    for(int i = ;i < n;i++)  for(int j = ;j < n;j++)  printf("%d%c", ans.mat[i][j], j == n- ? '\n' : ' ');

}

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

参考链接：

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

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

巴特西

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

题目

分析

最新文章

热门文章