poj 3233（矩陣的快速冪+等比矩陣的性質）

題意：容易理解.

分析：求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，得到就是要求的矩陣了！

代碼實現：

 

#include<stdio.h>
#include<string.h>
struct node{
    int p[65][65];
};
int mod,len;
struct node suan(struct node a,struct node b)//矩陣a乘以矩陣b
{
    int i,j,k;
    struct node c;
    for(i=1;i<=len;i++)
    {
        for(j=1;j<=len;j++)
        {
            c.p[i][j]=0;
            for(k=1;k<=len;k++)
                c.p[i][j]=(a.p[i][k]*b.p[k][j]+c.p[i][j])%mod;
        }
    }
    return c;
}
struct node haha(struct node a,struct node b,int n)
{
    while(n)//矩陣的快速冪
    {
        if(n%2==1)
            b=suan(b,a);
        n=n/2;
        a=suan(a,a);
    }
    return b;
}
int main()
{
    int i,j,n,k;
    struct node a,b;
    while(scanf("%d%d%d",&n,&k,&mod)!=EOF)
    {
        len=n*2;
        for(i=1;i<=n;i++)
            for(j=1;j<=n;j++)
                scanf("%d",&a.p[i][j]);
        for(i=1;i<=n;i++)//右上部分
            for(j=n+1;j<=n*2;j++)
                if(i+n==j)
                    a.p[i][j]=1;
                else
                    a.p[i][j]=0;
        for(i=n+1;i<=n*2;i++)//左下部分
            for(j=1;j<=n;j++)
                a.p[i][j]=0;
        for(i=n+1;i<=2*n;i++)//右下部分
            for(j=n+1;j<=n*2;j++)
                if(i==j)
                    a.p[i][j]=1;
                else
                    a.p[i][j]=0;
        for(i=1;i<=n*2;i++)//把b變成單位矩陣
            for(j=1;j<=n*2;j++)
                if(i==j)
                    b.p[i][j]=1;
                else
                    b.p[i][j]=0;
        a=haha(a,b,k+1);
        for(i=1;i<=n;i++)//減去單位矩陣
            for(j=n+1;j<=len;j++)
            {
                if(i+n==j)
                    a.p[i][j]--;
                while(a.p[i][j]<0)//為了防止溢出
                    a.p[i][j]+=mod;
            }
        for(i=1;i<=n;i++)
        {
            for(j=n+1;j<len;j++)
                printf("%d ",a.p[i][j]);
            printf("%d\n",a.p[i][len]);
        }
    }
    return 0;
}