Matrix Power Series【矩阵乘法】

>Link

POJ 3233

---

>Description

给出一个$n\ast n$的矩阵$A$，求出模$m$的情况下 $S_n=A+A^1+A^2+...+A^n$

---

>解题思路

我们先考虑不是矩阵的情况，很容易推出求$s_n$如下
$\left[\begin{array}{cc}{a}^{n-1}& s_{n-2}\end{array}\right]\to \left[\begin{array}{cc}{a}^n(a^{n-1}\ast a)& s_{n-1}(s_{n-2}+a^{n-1})\end{array}\right]$

所以$B$矩阵应为$\left[\begin{array}{cc}a& 1\\ 0& 1\end{array}\right]$

但是$a$和$s$实际上是矩阵，我们要求的是
$\left[\begin{array}{cc}{A}^{n-1}& S_{n-2}\end{array}\right]\to \left[\begin{array}{cc}{A}^n(A^{n-1}\ast A)& S_{n-1}(S_{n-2}+A^{n-1})\end{array}\right]$
其实可以把$B$矩阵扩大，也就是把$a$和$1$和$0$置换成等价的$n\ast n$大小的$X$矩阵

a：正常操作，$A\ast X=A^2$，所以$X$与$A$相同
1：一个矩阵乘上另一个矩阵$X$与原本的相同，可推出$X$矩阵的对角线为1，其余全为0
0：一个矩阵乘上一个矩阵为0，所以$X$矩阵为全为0的矩阵

最后加进去，$B$矩阵变成：
$\left[\begin{array}{cc}\left[\begin{array}{ccc}{a}_{1,1}& a_{1,2}& ...\\ ...& ...& ...\\ ...& ...& a_{n,n}\end{array}\right]& \left[\begin{array}{ccc}1& 0& ...\\ 0& 1& ...\\ 0& ...& 1\end{array}\right]\\ \left[\begin{array}{ccc}0& 0& 0\\ ...& ...& ...\\ 0& 0& 0\end{array}\right]& \left[\begin{array}{ccc}1& 0& ...\\ 0& 1& ...\\ 0& ...& 1\end{array}\right]\end{array}\right]$

---

>代码

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define LL long long
using namespace std;

struct matrix
{
	int x, y;
	LL a[70][70];
} A, B, ans;
LL m;

matrix operator *(matrix a, matrix b)
{
	matrix c;
	c.x = a.x, c.y = b.y;
	for (int i = 1; i <= c.x; i++)
	  for (int j = 1; j <= c.y; j++) c.a[i][j] = 0;
	for (int k = 1; k <= a.y; k++)
	  for (int i = 1; i <= c.x; i++)
	    for (int j = 1; j <= c.y; j++)
	      c.a[i][j] = (c.a[i][j] + a.a[i][k] * b.a[k][j] % m) % m;
	return c;
}
void power (LL k)
{
	if (k == 1) {B = A; return;}
	power (k >> 1);
	B = B * B;
	if (k & 1) B = B * A;
}

int main()
{
	int n;
	LL k;
	scanf ("%d%lld%lld", &n, &k, &m);
	ans.x = n, ans.y = 2 * n;
	for (int i = 1; i <= n; i++)
	  for (int j = 1; j <= n; j++)
	   scanf ("%lld", &ans.a[i][j]), ans.a[i][j] %= m;
	A.x = A.y = 2 * n;
	for (int i = 1; i <= n; i++)
	  for (int j = 1; j <= n; j++) A.a[i][j] = ans.a[i][j];
	for (int i = 1; i <= n; i++)
	  A.a[i][i + n] = A.a[i + n][i + n] = 1;
	power (k);
	ans = ans * B;
	for (int i = 1; i <= n; i++)
	{
		for (int j = n + 1; j <= 2 * n; j++)
		  printf ("%lld ", ans.a[i][j]);
		printf ("\n");
	}
	return 0;
}