POJ—3233—Matrix power series—【矩阵快速幂】【二分，递归，分治】

2 2 4
0 1
1 1

1 2
2 3

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

typedef struct MATRIX
{
int mat[35][35];
}MATRIX;
MATRIX A;
int n,m;

MATRIX mat_multiply (MATRIX a,MATRIX b)     //矩阵乘法
{
MATRIX c;
memset(c.mat,0,sizeof(c.mat));
for(int k=0;k<n;k++)
for(int i=0;i<n;i++)
if(a.mat[i][k])
for(int j=0;j<n;j++)
if(b.mat[k][j])
{
c.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
if(c.mat[i][j]>=m)
c.mat[i][j]%=m;
}
return c;
}

MATRIX mat_add (MATRIX a,MATRIX b)          //矩阵加法
{
MATRIX c;
memset(c.mat,0,sizeof(c.mat));
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
{
if(a.mat[i][j] || b.mat[i][j])
{
c.mat[i][j]=a.mat[i][j] + b.mat[i][j];
if(c.mat[i][j]>=m)
c.mat[i][j]%=m;
}
}
return c;
}

MATRIX mat_quick_index (int N)           //矩阵求幂
{
MATRIX E,a;
a=A;
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
E.mat[i][j]=(i==j);
while(N)
{
if(N & 1)
E=mat_multiply(E,a);
N>>=1;
a=mat_multiply(a,a);
}
return E;
}

MATRIX mat_sum (int k)
{
if(k==1)
return A;
MATRIX temp,tnow;
temp=mat_sum(k/2);
if(k & 1)
{
tnow=mat_quick_index(k/2+1);
temp=mat_add(temp,mat_multiply(temp,tnow));
temp=mat_add(temp,tnow);
}
else
{
tnow=mat_quick_index(k/2);
temp=mat_add(temp,mat_multiply(temp,tnow));
}
return temp;
}

int main()
{
MATRIX B;
int k;
memset(B.mat,0,sizeof(B.mat));
scanf("%d%d%d",&n,&k,&m);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
{
scanf("%d",&A.mat[i][j]);
A.mat[i][j]%=m;     //初始化的时候就先取余。。。
}
B=mat_sum(k);
for(int i=0;i<n;i++)
{
for(int j=0;j<n;j++)
{
printf("%d",B.mat[i][j]);
if(j!=n-1)
printf(" ");
}
printf("\n");
}
return 0;
}

#include<iostream>
#include<cstdio>

#define MAXN 100

using namespace std;

struct Matrix
{
int size;
long modulo;
long element[MAXN][MAXN];
void setSize(int);
void setModulo(int);
Matrix operator* (Matrix);
Matrix power(int);

};

void Matrix::setSize(int a)
{
for (int i=0; i<a; i++)
for (int j=0; j<a; j++)
element[i][j]=0;
size = a;
}

void Matrix::setModulo(int a)
{
modulo = a;
}

Matrix Matrix::operator* (Matrix param)
{
Matrix product;
product.setSize(size);
product.setModulo(modulo);
for (int i=0; i<size; i++)
for (int j=0; j<size; j++)
for (int k=0; k<size; k++)
{
product.element[i][j]+=element[i][k]*param.element[k][j];
product.element[i][j]%=modulo;
}
return product;
}

Matrix Matrix::power(int exp)
{
Matrix tmp = (*this) * (*this);
if (exp==1)
return *this;
else if (exp & 1)
return tmp.power(exp/2) * (*this);
else
return tmp.power(exp/2);
}

int n,m,k;
Matrix m1,ans;
int main()
{
while(~scanf("%d%d%d",&n,&k,&m))
{
m1.setModulo(m);
m1.setSize(n*2);
for(int i=0;i<n;i++)
{
for(int j=0;j<n;j++)
{
scanf("%d",&m1.element[i][j]);
m1.element[i][j]%=m;
}
m1.element[i][i+n]=1;
m1.element[n+i][n+i]=1;
}
m1.element

=1;
ans=m1.power(k+1);
for(int i=0;i<n;i++)
for(int j=n;j<n*2;j++)
{
int out=ans.element[i][j];
if(j==i+n)
out=out==0?m-1:out-1;
printf("%d%s",out,j!=n*2-1?" ":"\n");
}
return 0;

}
return 0;
}