矩阵与等比数列求和

矩阵乘法是一个很奇妙的东西，他经常被用于描述递推式然后运用矩阵快速幂解决问题

如果我们要求An=1+q+q*q+q*q*q...+q^n

那么可以想出递推式：

An=An-1*q+1

像其他递推式一样，我们也可以用矩阵描述它：

q    1  An-1 An-1*q+1=An

* =

0 1 1 1

这样我们就可以用矩阵快速幂来解决等比数列求和，不用再去写什么二分

但还有更精彩的：

q^n An-1   q 1 q^(n+1)q^n+An-1=An

*  =

0     1 0 1 0  1

这意味着

q             1

0             1

的（N+1）次方后，

右上角元素就是我们要求的An

对于矩阵的等比数列求和，也可以这么做

只要将1和0换为同大小的I，O矩阵

就可以算了，还不用递归（没有了爆栈的威胁）

还有一些方法也可以非递归求矩阵等比数列和：

比如倍增，用空间换常数（比矩阵套矩阵快1倍左右，但代码真恶心）：

#include<cstdio>
#include<cstring>
#include<cctype>
#include<algorithm>
#define maxn 105
#define lim 24
using namespace std;

int n,N,m,k;

//Mat type
int tmp[maxn][maxn],f[2][lim][maxn][maxn];//0 Pow , 1 Sum
int A[maxn][maxn],B[maxn][maxn];
char S[maxn];

inline void Clear(int ret[maxn][maxn])
{ 
	/*
	for(int i=0;i<n;i++) 
		for(int j=0;j<n;j++) 
			ret[i][j]=0;
	*/
	memset(ret,0,sizeof tmp);
}
inline void Unit(int ret[maxn][maxn])
{ 	
	Clear(ret);
	for(int i=0;i<n;i++) ret[i][i]=1;
}
inline void copy(int ret[maxn][maxn],int base[maxn][maxn])
{
	/*
	for(int i=0;i<n;i++) 
		for(int j=0;j<n;j++) 
			ret[i][j]=base[i][j]; 
	*/
	memcpy(ret,base,sizeof tmp);
}
inline void mul(int ret[maxn][maxn],int A[maxn][maxn],int B[maxn][maxn])
{
	//Clear(tmp);
	memset(tmp,0,sizeof tmp);
	for(int i=0;i<n;i++)
		for(int k=0;k<n;k++) if(A[i][k])
			for(int j=0;j<n;j++) if(B[k][j])
				tmp[i][j]=(tmp[i][j]+1ll*A[i][k]*B[k][j]%m)%m;
	//copy(ret,tmp);
	memcpy(ret,tmp,sizeof tmp);
}
inline void add(int ret[maxn][maxn],int A[maxn][maxn],int B[maxn][maxn])
{
	for(int i=0;i<n;i++) 
		for(int j=0;j<n;j++) 
			ret[i][j]=(A[i][j]+B[i][j])%m;
}
inline void Prep(int base[maxn][maxn],int k)
{
	//copy(f[0][1],base);copy(f[1][1],base);
	memcpy(f[0][1],base,sizeof f[0][1]);
	memcpy(f[1][1],base,sizeof f[1][1]);
	for(int i=2,j=2;i<=k;i<<=1,j++) 
		mul(f[0][j],f[0][j-1],f[0][j-1]),
		mul(f[1][j],f[1][j-1],f[0][j-1]),add(f[1][j],f[1][j],f[1][j-1]);
}

inline void Pow(int ret[maxn][maxn],int base[maxn][maxn],int k)
{
	for(Unit(ret);k;k>>=1,mul(base,base,base)) 
		if(k&1) 
			mul(ret,ret,base);
}
inline void SumPow(int ret[maxn][maxn],int base[maxn][maxn],int k)
{
	//Clear(ret);
	memset(ret,0,sizeof ret);
	Prep(base,k);
 	for(int i=1;k;k>>=1,i++) 
		if(k&1) 
			mul(ret,ret,f[0][i]),add(ret,ret,f[1][i]);
}

int main()
{
	
	freopen("tour.in","r",stdin);
	freopen("tour.out","w",stdout);
	
	scanf("%d",&n);N=n<<1;
	for(int i=0;i<n;i++) 
	{ 	
		scanf("%s",S);
		for(int j=0;j<n;j++) A[i][j]=(S[j]=='Y');
	}	
	scanf("%d%d",&k,&m);
	SumPow(B,A,k-1);
	int ans=0;
	for(int i=0;i<n;i++)
		ans=(ans+B[i][i])%m;
	printf("%d\n",ans);
}

但是数组法写矩阵看上去方便实际上一大堆坑，memset后面sizeof不能用传进来的数组，甚至还要用外面的数组辅助才能乘法，代码不清晰，反而又臭又长。。。。。。以后再也不用了。

但是公式法（q^n-1）/(q-1)

也有矩阵乘法所不能的功能

比如求 a^(b^c)%p

那么log（b^c）都是 O(c)级别的

那么可以用欧拉函数，a^(b^c)%p=a^(b^c % phi(p) ) %p

如果知道phi

复杂度是log(c)+log(b)

矩阵乘法不行

比如求公差为a，前(b^c)项和，有逆元就行，（没有逆元就。。。。）

矩阵显然没有啥欧拉函数吧。