poj 3233(等比矩阵的性质)

转载自：http://www.cnblogs.com/jiangjing/archive/2013/05/28/3103336.html

题意：求A^1+A^2+....+A^(n-1) + A^n.

分析：求a^1+..a^n这是矩阵乘法中关于等比矩阵的求法：

|A  E|

|0  E|

其中的A为m阶矩阵，E是单位矩阵，0是零矩阵。而我们要求的是：                                                                             

A^1+A^2+..A^L，由等比矩阵的性质

|A  ,  1|                 |A^n , 1+A^1+A^2+....+A^(n-1)|

|0  ,  1| 的n次方等于|0     ,                   1                   |

所以我们只需要将A矩阵扩大四倍，变成如上形式的矩阵B，然后开L+1次方就可以得到1+A^1+A^2+....+A^L。由于多了一个1，所以最后得到的答案我们还要减去1。同理我们把矩阵A变成B：

                                                          |A  E|

                                                          |0  E|

然后我们就是求B的n+1次幂之后得到的矩阵为|x1   x2|

                                                          |x3   x4|

右上角的矩阵x2减去单位矩阵E，得到就是要求的矩阵了！

Code:

#include<iostream>#include<cstdio>#include<algorithm>#include<cstring>using namespace std;const int maxn = 60 + 10;#define INF 0x3f3f3f3ftypedef long long ll;struct Matrix{    Matrix()    {        memset(a,0,sizeof(a));    }    int a[maxn][maxn];};int n,k,mod;Matrix Mul(Matrix x,Matrix y,int mod_val){    Matrix z;    for(int i = 0; i < 2 * n; i ++)    {        for(int k = 0; k < 2 * n; k ++)        {            if(x.a[i][k] == 0)                continue;            for(int j = 0;j < 2 * n; j ++)                z.a[i][j] = (z.a[i][j] + x.a[i][k] * y.a[k][j] % mod_val) % mod_val;        }    }    return z;}Matrix pow_mod(Matrix p,int t,int mod_val){    Matrix ret;    for(int i = 0; i < 2 * n ; i ++)        ret.a[i][i] = 1;    while(t)    {        if(t & 1)            ret = Mul(ret,p,mod_val);        p = Mul(p,p,mod_val);        t >>= 1;    }    return ret;}int temp[maxn][maxn];int E[maxn][maxn];int ans[maxn][maxn];Matrix p;int main(){    while( ~ scanf("%d%d%d",&n,&k,&mod))    {        for(int i = 0; i < n; i ++)            for(int j = 0; j < n; j ++)                scanf("%d",&temp[i][j]);        memset(E,0,sizeof(E));        for(int i = 0; i < n; i ++)            E[i][i] = 1;        for(int i = 0; i < n * 2; i ++)            for(int j = 0; j < n * 2; j ++)        {            if(i < n && j < n)                p.a[i][j] = temp[i][j];            if(i < n && j >= n)                p.a[i][j] = E[i][j - n];            if(i >= n && j >= n)                p.a[i][j] = E[i - n][j - n];            if(i >= n && j < n)                p.a[i][j] = 0;        }        Matrix t = pow_mod(p,k + 1,mod);        for(int i = 0; i < n; i ++)        {            for(int j = n; j < n * 2; j ++)            {                ans[i][j - n] = t.a[i][j];                if(i == j - n)                {                    ans[i][j - n] -= 1;                    if(ans[i][j - n] < 0)                        ans[i][j - n] += mod;                }            }        }        for(int i = 0; i < n; i ++)        {            for(int j = 0; j < n; j ++)            printf("%d%c",ans[i][j],j < n - 1 ? ' ' : '\n');        }    }    return 0;}