POJ 3233 Matrix Power Series(矩阵+二分)

题目大意：求由矩阵 A构成的矩阵 S = A + A^2 + A^3 + … + A^k。k的取值范围是：10^9数据很大，应该二分。

对于一个k来说，s(k) = (1+A^(k/2)) *( A+A^2+……+A^(k/2))。如果k为奇数的话需要加上A^(k/2 + 1)。

所以二分求和，复杂度就降下来了，当然还得用到矩阵快速幂。

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 nnonnegative
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

#include <algorithm>
#include <iostream>
#include <stdlib.h>
#include <string.h>
#include <iomanip>
#include <stdio.h>
#include <string>
#include <queue>
#include <cmath>
#include <stack>
#include <map>
#include <set>
#define eps 1e-10
///#define M 1000100
#define LL __int64
///#define LL long long
///#define INF 0x7ffffff
#define INF 0x3f3f3f3f
#define PI 3.1415926535898
#define zero(x) ((fabs(x)<eps)?0:x)

///#define mod 9973

int mod;
const int maxn = 2010;

using namespace std;

struct matrix
{
    int f[40][40];
};

matrix mul(matrix a, matrix b, int n)///矩阵乘法
{
    matrix c;
    memset(c.f, 0, sizeof(c.f));
    for(int i = 0; i < n; i++)
    {
        for(int j = 0; j < n; j++)
        {
            for(int k = 0; k < n; k++) c.f[i][j] += a.f[i][k]*b.f[k][j];
            c.f[i][j] %= mod;
        }
    }
    return c;
}

matrix pow_mod(matrix a, int b, int n)///矩阵快速幂
{
    matrix s;
    memset(s.f, 0 , sizeof(s.f));
    for(int i = 0; i < n; i++) s.f[i][i] = 1;
    while(b)
    {
        if(b&1) s = mul(s, a, n);
        a = mul(a, a, n);
        b >>= 1;
    }
    return s;
}

matrix Add(matrix a,matrix b, int n)  ///矩阵加法
{
    matrix c;
    for(int i = 0; i < n; i++)
    {
        for(int j = 0; j < n; j++)
        {
            c.f[i][j] = a.f[i][j]+b.f[i][j];
            c.f[i][j] %= mod;
        }
    }
    return c;
}

matrix Matrix_Sum(matrix a, int k, int n)
{
    if(k == 1) return a;
    matrix dx,dy;
    dx = Matrix_Sum(a, k/2, n);///二分，递归；
    if(k&1)
    {
        dy = pow_mod(a, k/2+1, n);
        dx = Add(dx, mul(dx, dy, n), n);
        dx = Add(dy,dx, n);
    }
    else
    {
        dy = pow_mod(a, k/2, n);
        dx = Add(dx,mul(dx, dy, n), n);
    }
    return dx;
}

int main()
{
    int n, k;
    while(cin >>n>>k>>mod)
    {
        matrix c;
        for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n; j++)
            {
                scanf("%d",&c.f[i][j]);
                c.f[i][j] %= mod;
            }
        }
        matrix d = Matrix_Sum(c, k, n);
        for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n-1; j++) cout<<d.f[i][j]<<" ";
            cout<<d.f[i][n-1]<<endl;
        }
    }
    return 0;
}