POJ3233 Matrix Power Series（快速幂求等比矩阵和）

题面

---

---

$ solution: $

首先，如果题目只要我们求 $ A^K $ 那这一题我们可以直接模版矩乘快速幂来做，但是它现在让我们求 $ sum_{i=1}^{k}{(A^i)} $ 所以我们思考一下这两者是否有什么关系。仔细一想，不难发现几个东西：

1. 一次矩阵乘法复杂度为 $ O(n^3) $ ，所以我们不能进行太多次矩阵乘法
2. 快速幂的复杂度为 $ O(logk) $ 再乘一下矩阵乘法的复杂度，我们现在只能再接受 $ O(log) $ 级别的处理了
3. 矩阵乘法满足交换律和结合律！！！！
4. 若我们已经知道了 $ A^1+A^2+A^3+A^4 $ 的值，我们 需要 $ A^5+A^6+A^7+A^8 $ 的值，我们可以直接将前者乘上一个 $ A^4 $ 就可以了！

根据以上发现，我们不妨再设一个矩阵B，来帮助我们理解！

我们设 $ B_x=A^1+A^2+.....+A^x $ 根据上面第四个的原理我们可以得到：

1. $ B_1=A^1 $
2. $ B_2=A^1+A^1 imes A^1=B_1+B_1 imes A^1 $
3. $ B_4=(A^1+A^2)+A^2 imes (A^1+A^2)=B_2+A^2 imes B_2 $
4. $ B_8=(A^1+A^2+A^3+A^4)+A^4 imes (A^1+A^2+A^3+A^4)=B_4+A^4 imes B_4 $
5. $ B_{16}=B_8+A^8 imes B_8 $
6. $ B_{2 imes x}==B_x+A_x imes B_x $
7. 在根据上面第四个规律可以得出： $ B_{x imes y}==B_x+A_x imes B_y $

而我们要得到的最终结果就是 $ B_K $ 嘛。如果上述的B矩阵的底数能相加，那我们不就可以仿照二进制来得到 $ B_K $ 了吗？就像K等于19，如果我们可以直接 $ B_{19}=B_{1+2+16}=B_1+B_2+B_{16} $ （假设！）那我们就能log预处理所有 $ B_{1<<i} $ 然后log求出 $ B_K $ 了！

我们再思考一下，发现上述 $ B $ 矩阵的底数是可以相加的，但是不像上面那个式子那么直接相加，我们应该这么求： $ B_{x+y}=B_x+A^x imes B_y $ 这样底数就可以相加了！！！

然后我们再仿照二进制，就能log求出 $ B_K $ 了！

---

---

$ code: $

#include<iostream>
#include<cstdio>
#include<iomanip>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<ctime>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>

#define ll long long
#define db double
#define inf 0x7fffffff
#define rg register int

using namespace std;

int n,m,k;

inline int qr(){ char ch;//快读
	while((ch=getchar())<'0'||ch>'9');
	int res=ch^48;
	while((ch=getchar())>='0'&&ch<='9')
		res=res*10+(ch^48);
	return res;
}

struct su{
	int s[30][30];
	inline void read(){//读入一个矩阵
		for(rg i=0;i<n;++i)
			for(rg j=0;j<n;++j)
				s[i][j]=qr();
	}
	inline void write(){//输出一个矩阵
		for(rg i=0;i<n;++i){
			for(rg j=0;j<n;++j)
				printf("%d ",s[i][j]);
			puts("");
		}
	}
	inline su operator *(su x){//矩阵乘法
		su y; memset(y.s,0,sizeof(y.s));
		for(rg i=0;i<n;++i)
			for(rg j=0;j<n;++j)
				for(rg o=0;o<n;++o)
					y.s[i][j]+=(ll)s[i][o]*x.s[o][j]%m,y.s[i][j]%=m;
		return y;
	}
	inline su operator +(su x){//矩阵加法
		for(rg i=0;i<n;++i)
			for(rg j=0;j<n;++j)
				x.s[i][j]+=s[i][j],x.s[i][j]%=m;
		return x;
	}
}ans,a[33],b[33];

int main(){
	//freopen("t1.in","r",stdin);
	//freopen("t1.out","w",stdout);
	n=qr();k=qr();m=qr();
	a[1].read(); b[1]=a[1];
	for(rg i=1;i<32;++i){//我们求出对应所有的a与b
		a[i+1]=a[i]*a[i]; //a数组表示A[1<<i-1]
		b[i+1]=b[i]*a[i]+b[i]; //b数组表示b[1<<i-1]
	}
	for(rg i=1;i<=32;++i){//根据二进制，不断累加，一直到b[k]
		if(k&1)ans=ans*a[i]+b[i]; //ans相当于b[ans]
		k>>=1; //k&1表示这一位上是一，不懂可以先学下快速幂的原理
	} ans.write();//累加完毕输出
	return 0;
}

我们发现上面main函数中的两个for循环可以换成一个，于是：

int main(){
	n=qr();k=qr();m=qr(); a.read(); b=a;
	while(k){
		if(k&1)ans=ans*a+b;
		b=b*a+b; a=a*a; k>>=1;
	} ans.write(); //不仅代码短，跑的还十分快！！！！！
	return 0;
}