蓝书（算法竞赛进阶指南）刷题记录——POJ3233 Matrix Power Series（分治+矩阵乘法）

题目：POJ3233.
题目大意：给定$k$和一个$n\ast n$的矩阵$A$，求${\sum }_{i=1}^k{A}^i$.
$1\le k\le 10^9,1\le n\le 30$.

一开始想到的是直接用等比数列求和公式，一看好像不会矩阵求逆…

然后想到${\sum }_{i=l}^r{A}^i=A^{l-1}{\sum }_{i=1}^{r-l+1}{A}^i$，顿时有了思路.

考虑分治解决这个问题，开始处理区间$[1,n]$的时候考虑处理出$lans=[1,\lfloor \frac{n}{2}\rfloor ]$的答案，然后分情况讨论：
1.当$2|n$时，直接计算出$[\lfloor \frac{n}{2}\rfloor ,n]$的答案为$rans=lans\ast A^{\frac{n}{2}}$.
2.否则，计算出$[\lfloor \frac{n}{2}\rfloor ,n-1]$的答案为$rans=lans\ast A^{\lfloor \frac{n}{2}\rfloor }$，然后再加上一个$A^n$即可.

考虑到矩阵的幂可以用快速幂直接求出，那么我们就可以做到$O(\log k)$次递归，递归到每一层都只用$O(\log k)$的时间计算，总复杂度$O(n^3{\log }^{2}k)$.

代码如下：

#include<iostream>
#include<cstdio>
#include<algorithm>
  using namespace std;

#define Abigail inline void
typedef long long LL;

const int N=30;

int n,k,mod;
struct matrix{
  int n,m,v[N+9][N+9];
  
  matrix operator + (const matrix &p)const{
  	matrix tmp=matrix();
  	tmp.n=n;tmp.m=m;
	for (int i=1;i<=n;++i)
	  for (int j=1;j<=m;++j)      //手残把j写成了i
	    tmp.v[i][j]=(v[i][j]+p.v[i][j])%mod;
	return tmp;
  }
  
  matrix operator * (const matrix &p)const{
    matrix tmp=matrix();      //然后没有初始化...
    tmp.n=n;tmp.m=p.m;
    for (int i=1;i<=tmp.n;++i)
      for (int k=1;k<=m;++k)
        for (int j=1;j<=tmp.m;++j)
          tmp.v[i][j]=(tmp.v[i][j]+v[i][k]*p.v[k][j])%mod;
    return tmp;
  }
  
}a;

matrix e(int n){
  matrix e=matrix();
  for (int i=1;i<=n;++i)
    e.v[i][i]=1;
  e.n=e.m=n;
  return e;
}

matrix power(matrix a,int k){
  matrix s=e(a.n);
  for (;k;k>>=1,a=a*a)
    if (k&1) s=s*a;
  return s;
}

matrix Calc(matrix a,int k){
  if (k==1) return a;
  matrix ans=Calc(a,k>>1);
  ans=ans+power(a,k>>1)*ans;
  if (k&1) ans=ans+power(a,k);
  return ans; 
}

Abigail into(){
  scanf("%d%d%d",&n,&k,&mod);
  a.n=a.m=n;
  for (int i=1;i<=n;++i)
    for (int j=1;j<=n;++j){
      scanf("%d",&a.v[i][j]);
      a.v[i][j]%=mod;
	}
}

Abigail work(){
  a=Calc(a,k);
}

Abigail outo(){
  for (int i=1;i<=n;++i){
  	for (int j=1;j<n;++j)
  	  printf("%d ",a.v[i][j]);
  	printf("%d\n",a.v[i][n]);
  }
}

int main(){
  into();
  work();
  outo();
  return 0;
}