Poj 3233 Matrix Power Series

题目连接：http://poj.org/problem?id=3233

本题练习矩阵运算。

需要两次用到二分法，求和一次，求幂一次。

二分法的知识参考我以前写过的一篇关于二分法的思考：http://blog.youkuaiyun.com/niuox/article/details/7226471 。美其名曰：快速幂。

代码:

#include <iostream>
#include <stdio.h>
#include <string.h>
using namespace std;
#define max 32

struct Matrix
{
    int data[max][max];
}matrix;

int n,mod;

//矩阵相加
Matrix add(Matrix u,Matrix v)
{
    Matrix t;
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            t.data[i][j] = (u.data[i][j] + v.data[i][j]) % mod;
        }
    }
    return t;
}

//矩阵相乘
Matrix mul(Matrix u,Matrix v)
{
    Matrix t;
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            t.data[i][j] = 0;
            for(int k=0;k<n;k++)
            {
                t.data[i][j] += (u.data[i][k] * v.data[k][j])%mod;
                t.data[i][j] %= mod;
            }
        }
    }
    return t;
}

//矩阵的幂
Matrix power(int k)
{
    int k_temp = k;
    Matrix result,a;

    memset(result.data,0,sizeof(result.data));

    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            a.data[i][j] = matrix.data[i][j];
            if(i == j)
            {
                result.data[i][j] = 1;
            }
        }
    }
    while(k_temp)
    {
        if(k_temp&1)
        {
            result = mul(result,a);
        }
        a = mul(a,a);
        k_temp = k_temp>>1;
    }
    return result;

}
//二分法求和
Matrix halfDiv(int k)
{
    if(k == 1)
    {
        return matrix;
    }
    if(k&1)
    {
        return add(halfDiv(k-1),power(k));
    }
    else
    {
        Matrix t = halfDiv(k>>1);
        return add(t,mul(t,power(k>>1)));
    }
}
int main()
{
    #ifndef ONLINE_JUDGE
        freopen("in.txt","r",stdin);
    #endif
    int k;
    scanf("%d %d %d",&n,&k,&mod);
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            scanf("%d",&matrix.data[i][j]);
        }
    }

    Matrix result = halfDiv(k);
    //Matrix result = power(2);
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n-1;j++)
        {
            printf("%d ",result.data[i][j]);
        }
        printf("%d\n",result.data[i][n-1]);
    }

    return 0;
}

另一种方法，构造法，能够构造出求和结果的矩阵。

思想是：

B= 

{A,A

0,E

},那么 B^k = 

{A,A+A^2+A^3...+A^k

0,E

}

代码：

#include <iostream>
#include <stdio.h>
#include <string.h>
using namespace std;
#define max 65

struct Matrix
{
    int data[max][max];
}matrix;

int n,mod;

//矩阵相乘
Matrix mul(Matrix u,Matrix v)
{
    Matrix t;
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            t.data[i][j] = 0;
            for(int k=0;k<n;k++)
            {
                t.data[i][j] += (u.data[i][k] * v.data[k][j])%mod;
                t.data[i][j] %= mod;
            }
        }
    }
    return t;
}

//矩阵的幂
Matrix power(int k)
{
    int k_temp = k;
    Matrix result,a;

    memset(result.data,0,sizeof(result.data));

    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            a.data[i][j] = matrix.data[i][j];
            if(i == j)
            {
                result.data[i][j] = 1;
            }
        }
    }
    while(k_temp)
    {
        if(k_temp&1)
        {
            result = mul(result,a);
        }
        a = mul(a,a);
        k_temp = k_temp>>1;
    }
    return result;

}

int main()
{
    #ifndef ONLINE_JUDGE
        freopen("in.txt","r",stdin);
    #endif
    int k;
    scanf("%d %d %d",&n,&k,&mod);
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            scanf("%d",&matrix.data[i][j]);
            matrix.data[i][j+n] = matrix.data[i][j];
            matrix.data[i+n][j] = 0;
            if(i == j)
            {
                matrix.data[i+n][i+n] = 1;
            }
        }
    }
    n = n<<1;
    //Matrix result = halfDiv(k);
    Matrix result = power(k);
    n = n>>1;
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n-1;j++)
        {
            printf("%d ",result.data[i][j+n]);
        }
        printf("%d\n",result.data[i][n+n-1]);
    }

    return 0;
}

推荐一个学习矩阵的好资料： http://www.matrix67.com/blog/archives/276