Matrix Power Series

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak.

Analysis

一道矩阵背景下的简单sumdiv，同样是两种算法。

• 等比数列

ans=A^ k+1 / A-1

好像很简单，但是要用到矩阵除法，即矩阵逆元，我目前不会..

• 拆分子问题

当k为偶数时，原式就等于(A^ k/2 +1)(A1+A2+...+A^k/2)，这样就拆分出了子问题，可以用快速幂解决，至于奇数..把最后一个拎出来就行。

Code

#include <iostream>
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <algorithm>
const int N=110;
int n,MOD;
class Matrix{
	private:
		int h,s[N][N];
		void clear(){
			h=0;
			memset(s,0,sizeof(s));
		}
	public:
		Matrix operator + (Matrix x){
			Matrix y;
			y.clear();
			y.h=h;
			for(int i=0;i<h;i++)
				for(int j=0;j<h;j++)
					y.s[i][j]=(s[i][j]+x.s[i][j])%MOD;
			return y;
		}
		Matrix operator * (Matrix x){
			Matrix y;
			y.clear();
			y.h=h;
			for(int i=0;i<h;i++)
				for(int j=0;j<h;j++)
					for(int k=0;k<h;k++)
						y.s[i][j]=(y.s[i][j]+s[i][k]*x.s[k][j])%MOD;
			return y;
		}
		Matrix operator *= (Matrix x){
			return *this=*this*x;
		}
		void Init(int x){
			clear();
			h=x;
			for(int i=0;i<h;i++)
				for(int j=0;j<h;j++)
					scanf("%d",&s[i][j]);
		}
		void Enit(int x){
			clear();
			h=x;
			for(int i=0;i<x;i++)
				s[i][i]=1;
		}
		void Output(){
			for(int i=0;i<h;i++){
				for(int j=0;j<h;j++)
					printf("%d ",s[i][j]);
				puts("");
			}
		}
}A,E,ans;
Matrix Fp(int k){
	Matrix x,y;
	x.Enit(n);
	y=A;
	while(k){
		if(k&1)x*=y;
		k>>=1;
		y*=y;
	}
	return x;
}
void Dc(int k){
	if(k==1)ans=A;
	else{
		Dc(k/2);
		ans*=Fp(k/2)+E;
		ans=k&1?ans+Fp(k):ans;
	}
}
int main(){
	int k;scanf("%d%d%d",&n,&k,&MOD);
	A.Init(n),E.Enit(n);
	Dc(k);
	ans.Output();
	return 0;
}