POJ 3233 Matrix Power Series(矩阵快速幂+二分求和)

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A² + A³ + … + A^k.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 10⁹) and m (m < 10⁴). 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的矩阵，求矩阵k阶阶乘，即求A^1+A^2+......+A^k的和
思路：因为k的范围是1到10^9，所以不能暴力去做，于是应该想到比较快的时间复杂度：二分。在快速幂的基础上，通过二分来求和

具体代码如下：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <algorithm>
using namespace std;

const int m=33;

struct node
{
    int matrix[m][m];
    node(){}
    friend node operator * (node a,node b);
    friend node operator + (node a,node b);
};

int mod;
int n;

void init(node &a)
{
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
            a.matrix[i][j]=0;
    }

}

node operator * (node a,node b)
{
    node c;
    init(c);
    for(int k=0;k<n;k++)
    {
        for(int i=0;i<n;i++)
        {
            for(int j=0;j<n;j++)
            {
                c.matrix[i][j]=(c.matrix[i][j]+a.matrix[i][k]*b.matrix[k][j])%mod;
            }
        }
    }
    return c;
}

node operator + (node a,node b)
{
    node tmp;
    init(tmp);
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
            tmp.matrix[i][j]=(a.matrix[i][j]+b.matrix[i][j])%mod;
    }
    return tmp;
}

node fast_pow(node a,long long k)
{
    node b;
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            if(i==j) b.matrix[i][j]=1;
            else b.matrix[i][j]=0;
        }
    }
    while(k)
    {
        if(k&1) b=a*b;
        a=a*a;
        k>>=1;
    }
    return b;
}

node sum(node a,long long k)    //二分求和
{
    node tmp1,tmp2;
    if(k==1) return a;
    tmp1=sum(a,k/2);      //递归
    if(k&1)
    {
        tmp2=fast_pow(a,k/2+1);
        tmp1=tmp2*tmp1+tmp1;
        return tmp1+tmp2;
    }
    else
    {
        tmp2=fast_pow(a,k/2);
        return tmp2*tmp1+tmp1;
    }
}

int main()
{
    long long k;
    while(~scanf("%d%I64d%d",&n,&k,&mod))
    {
        node a;
        for(int i=0;i<n;i++)
        {
            for(int j=0;j<n;j++)
            {
                scanf("%d",&a.matrix[i][j]);
                a.matrix[i][j]%= mod;
            }
        }
        node b=sum(a,k);
        for(int i=0;i<n;i++)
        {
            for(int j=0;j<n;j++)
            {
                printf("%d%c",b.matrix[i][j],j==n-1?'
':' ');
            }
        }
    }
    return 0;
}