poj 3233 Matrix Power Series---矩阵快速幂

要求矩阵A的1到k次幂的和，矩阵快速幂加上二分求和

其中，矩阵相乘二分：A^2k=A^k*A^k,
                                        A^(2k+1)=A^k*A^k*A.
求和二分：A+A^2+A...+A^(2k+1)=   A+A^2+...+A^k+A^(k+1)+A^(k+1)*(A+A^2+...+A^k).
                   A+A^2+...+A^2k           =   A+A^2+...+A^k+A^k*(A+A^2+...+A^k).

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

int N,k,m;

struct matrix
{
       int a[35][35];
}origin,res;

matrix multiply(matrix x,matrix y)
{
    matrix temp;
    memset(temp.a,0,sizeof(temp.a));
    for(int i=0;i<N;i++)
    {
        for(int j=0;j<N;j++)
        {
            for(int k=0;k<N;k++)
            {
                temp.a[i][j]+=(x.a[i][k]*y.a[k][j]);
            }
            temp.a[i][j]%=m;
        }
    }
    return temp;
}

/*matrix calc(int n)
{
    while(n)
    {
        if(n&1)
        {
            res=multiply(res,origin);
            n--;
        }
        else
        {
            n>>=1;
            origin=multiply(origin,origin);
        }
    }
    return multiply(res,origin);
}*/

matrix calc(int exp)
{
    matrix p,q;
    p=origin;
    q=res;
    while (exp!=1)
    {
        if (exp&1)
        {
            exp--;
            q=multiply(p,q);
        }
        else
        {
            exp>>=1;
            p=multiply(p,p);
        }
    }
    return multiply(p,q);
}

matrix add(matrix x,matrix y)
{
    matrix tmp;
    int i,j;
    for(i=0;i<N;i++)
        for(j=0;j<N;j++)
        {
            tmp.a[i][j]=x.a[i][j]+y.a[i][j];
            tmp.a[i][j]%=m;
        }
    return tmp;
}


void init()
{
    int i,j;
    memset(res.a,0,sizeof res.a);
    for(i=0;i<N;i++)
        for(j=0;j<N;j++)
        {
            scanf("%d",&origin.a[i][j]);
            origin.a[i][j]%=m;
        }
    for(i=0;i<N;i++)
    {
        res.a[i][i]=1;
    }
}

matrix sum(int n)
{
    matrix tmp,now;
    if(n==1) return origin;
    tmp=sum(n/2);
    if(n&1)
    {
        now=calc(n/2+1);
        tmp=add(tmp,multiply(tmp,now));
        tmp=add(tmp,now);
    }
    else
    {
        now=calc(n/2);
        tmp=add(tmp,multiply(tmp,now));
    }
    return tmp;
}


int main()
{
    int i,j;
    while(~scanf("%d%d%d",&N,&k,&m))
    {
        init();
        matrix tmp=sum(k);
        for(i=0;i<N;i++)
        {
            printf("%d",tmp.a[i][0]);
            for(j=1;j<N;j++)
            {
                printf(" %d",tmp.a[i][j]);
            }
            putchar('\n');
        }
    }
    return 0;
}