POJ3233（矩阵快速幂+求矩阵的1~k次幂的和模m+输出结果矩阵）

Matrix Power Series

Time Limit: 3000MS
Memory Limit: 131072K

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 109) and m (m < 104). Then follow n lines each containing n nonnegative integers below 32,768, giving A’s elements in row-major order.

Output

Output the elements of S modulo m in the same way as A is given.

Sample Input

2 2 4
0 1
1 1

Sample Output

1 2
2 3

思路：

题意：

• 有一个 n * n 的矩阵 a，和一个正整数 k
• 求a + a² + a³ + … + a^k
• 最后模 m 输出结果矩阵

方法：

• 矩阵快速幂
• 二分求和
• 当k为偶数时，s(k) = (1 + A^k/2) * (A + A² + … + A^k/2)，1为单位矩阵
• 当k为奇数时，s(k) = A^(k+1)/2 + (1 + A^(k+1)/2) * (A + A² + … + A^k/2)，1为单位矩阵

AC代码：

#include <iostream>
#include <cstring>
#include <cstdio>
using namespace std;
#define maxn 33
struct matrix{
    int a[maxn][maxn];
}original,result,one;//初始,结果,单位矩阵 
int n,m;//n*n矩阵模m

//矩阵乘法 
matrix mul(matrix m1,matrix m2)
{
    matrix result;
    memset(result.a,0,sizeof(result.a));
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
        {
            for(int k=1;k<=n;k++)
            {
                result.a[i][j]+=m1.a[i][k]*m2.a[k][j];
                result.a[i][j]%=m;
            }
        }
    }
    return result;
}

//矩阵加法 
matrix add(matrix m1,matrix m2)
{
    matrix result;
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
        {
            result.a[i][j]=m1.a[i][j]+m2.a[i][j];
            result.a[i][j]%=m;
        }
    }
    return result;
}

//矩阵快速幂 
matrix mpow(matrix m1,int k)
{
    if(k==1) return m1;//1次幂为原矩阵 
    matrix result=mpow(m1,k/2);//result=m1^(k/2)
    result=mul(result,result);//result=m1^k
    if(k%2!=0) result=mul(result,m1);//若k为奇数再乘一次m1 
    return result;
} 

//将所有幂运算结果加一起 
matrix ans(matrix aa,int k)
{
    matrix m1,m2,m3,m4;
    if(k==1) return aa;
    if(k%2==0)//k为偶数 
    {
        m1=mpow(aa,k/2);//aa^(k/2) 
        m2=add(one,m1);//E+aa^(k/2)
        m3=ans(aa,k/2);//aa+aa^2+...+aa^(k/2)
        return mul(m2,m3);
    }
    else//k为奇数 
    {
        m1=mpow(aa,(k+1)/2);//aa^((k+1)/2)
        m2=add(one,m1);//E+aa^((k+1)/2)
        m3=ans(aa,k/2);//aa+aa^2+...+aa^(k/2)
        m4=mul(m2,m3);
        return add(m1,m4);
    }
} 

int main()
{
    int k;
    cin>>n>>k>>m;
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
        {
            cin>>original.a[i][j];
            original.a[i][j]%=m;
        }
    }
    for(int i=1;i<=n;i++)
    {
        one.a[i][i]=1;//单位矩阵 
    }
    result=ans(original,k);//结果矩阵 
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
        {
            if(j!=1) cout<<" ";//除了第一个数每个数后面有空格 
            cout<<result.a[i][j]%m;
        }
        cout<<endl;//输出一行后换行 
    }
    return 0;
}