题目背景

矩阵快速幂

题目描述

给定 n×nn\times n 的矩阵 A,求 AkA^k

输入格式

第一行两个整数 n,k 接下来 n 行,每行 n 个整数,第 i 行的第 j 的数表示 Ai,jA_{i,j}

输出格式

输出 AkA^k

共 n 行,每行 n 个数,第 i 行第 j 个数表示 (Ak)i,j(A^k)_{i,j},每个元素对 109+710^9+7 取模。

输入输出样例

输入 #1

2 1
1 1
1 1

输出 #1

1 1
1 1

说明/提示

【数据范围】
对于 100% 的数据:1n1000k1012,Ai,j10001\le n \le 100,0 \le k \le 10^{12}, |A_{i,j}| \le 1000∣

题目分析

先来了解一些矩阵相关的知识

矩阵

矩阵(matrix)。n×mn\times m 的矩阵指的是n行,m列的矩阵。

(123456)\begin{pmatrix} 1 & 2 & 3\\ 4 & 5 & 6\\ \end{pmatrix}

就是指的 2×32\times 3 的矩阵。

单位矩阵

单位矩阵指的是 对角线上为1,其他位置为0的矩阵。

(1001)\begin{pmatrix} 1 & \cdots & 0\\ \vdots & \ddots & \vdots\\ 0 & \cdots & 1\\ \end{pmatrix}

常用 I 来表示单位矩阵。

矩阵的幂次方

A0=IA^0=I

An=An1×A(n>0)A^n=A^{n-1}\times A (n>0)

性质

  1. 矩阵乘法满足分配率,结合律,不一定满足交换律
  2. 加法满足交换律和结合律

矩阵满足结合律,所以在求矩阵的幂的时候,可以使用 矩阵快速幂加速。

矩阵快速幂

分治的思路解决矩阵快速幂

An={An2×An2,n是偶数An2×An2×A,n是奇数A^n= \left\{\begin{aligned} &A^\frac{n}{2}\times A^\frac{n}{2}&,n是偶数 \\ &A^\frac{n}{2}\times A^\frac{n}{2} \times A &,n是奇数 \end{aligned} \right.

node matrixPow(node a,ll k){//矩阵的幂次方
	if(k==0){// 0次方
		return I;//矩阵的0次方是单位矩阵
	}
	node t=matrixPow(a,k/2);//求 a^{n/2} 次方
	if(k&1){//判断k是否是奇数
		return matrixMins(matrixMins(t,t),a);
	}else{//k是偶数
		return matrixMins(t,t);
	}
}

代码实现

#include <iostream>
#include <cstdio>
using namespace std;
typedef long long ll;
const int N=105;
const int M=1e9+7;
struct node{
	ll a[N][N]={0};
	int row,col;
};
node I;//单位矩阵

node matrixMins(node a,node b){//矩阵乘法
	node c;//答案矩阵
	c.row=a.row;
	c.col=b.col;
	int n=c.row,p=c.col,m=a.col;
	//计算矩阵乘法
	for(int i=1;i<=n;i++){
		for(int j=1;j<=p;j++){
			for(int k=1;k<=m;k++){
				c.a[i][j]+=a.a[i][k]*b.a[k][j]%M;
				c.a[i][j]%=M;
			}
		}
	}
	return c;
}
node matrixPow(node a,ll k){//矩阵的幂次方
	if(k==0){// 0次方
		return I;//矩阵的0次方是单位矩阵
	}
	node t=matrixPow(a,k/2);//求 a^{n/2} 次方
	if(k&1){//判断k是否是奇数
		return matrixMins(matrixMins(t,t),a);
	}else{//k是偶数
		return matrixMins(t,t);
	}
}
int main(){
	int n;
	node a;
	ll k;
	cin>>n>>k;
	a.col=a.row=n;
	I.col=I.row=n;
	for(int i=1;i<=n;i++){//输入矩阵A
		for(int j=1;j<=n;j++){
			cin>>a.a[i][j];
			//处理单元矩阵
			if(i==j) I.a[i][j]=1;
			else	I.a[i][j]=0;
		}
	}
	node ans=matrixPow(a,k);//计算a^k
	//输出结果矩阵
	for(int i=1;i<=ans.row;i++){
		for(int j=1;j<=ans.col;j++){
			cout<<ans.a[i][j]<<" ";
		}
		cout<<endl;
	}
	return 0;
}

Q.E.D.


( ノ^ω^)ノ゚ 稻 花 香 里 说 丰 年 , 听 取 WA 声 一 片 。(╥╯^╰╥)