UVA 10870 Recurrences（矩阵乘法）

题意

求解递推式 (f(n)=a_1*f(n-1)+a_2*f(n-2)+....+a_d*f(n-d)) 的第 (n) 项模以 (m)。

(1 leq n leq 2^{31}-1)

(1 leq m leq 46340)

(1 leq d leq 15)

思路

矩阵乘法最经典的运用之一。先大致介绍一下矩阵乘法：

对于一个矩阵 (A_{np}) ，另一个矩阵 (B_{pm}) ，设它们的乘积为 (C_{n,m}) ，有 (C_{i,j}=displaystylesum_{k=1}^pA_{i,k}B_{k,j}) .

例如对于一个矩阵 (egin{pmatrix}a_{1,1}&a_{1,2}&a_{1,3}\a_{2,1}&a_{2,2}&a_{2,3}end{pmatrix}​) ，和另一个矩阵 (egin{pmatrix}b_{1,1}&b_{1,2}\b_{2,1}&b_{2,2}\b_{3,1}&b_{3,2}end{pmatrix}​) ，它们的积为：

[egin{pmatrix} a_{1,1}b_{1,1}+a_{1,2}b_{2,1}+a_{1,3}b_{3,1} & a_{1,1}b_{1,2}+a_{1,2}b_{2,2}+a_{1,3}b_{3,2}\ a_{2,1}b_{1,1}+a_{2,2}b_{2,1}+a_{2,3}b_{3,1} & a_{2,1}b_{1,2}+a_{2,2}b_{2,2}+a_{2,3}b_{3,2} end{pmatrix} ]

从定义式可以看出来，矩阵乘法不满足交换律，但满足结合律。满足结合律，就说明了可以快速幂。

矩阵乘法的题目的根本想法是构造矩阵。对于这道题，可以先构造出矩阵 (A_{1d}) ，分别表示数列 (f) 的前 (d) 项，那么只需要再构造出一个 (B_{dd}) ，使得 (A_{1d}B_{dd}) 得到 (f) 数列的第 (2) 项到第 (d+1) 项即可。具体构造见代码：

代码

#include<bits/stdc++.h>
#define FOR(i,x,y) for(int i=(x),i##END=(y);i<=i##END;++i)
#define DOR(i,x,y) for(int i=(x),i##END=(y);i>=i##END;--i)
typedef long long LL;
using namespace std;
const int N=20;
int P;
struct Matrix
{
	int n,m,a[N][N];
	int *operator [](const int x){return a[x];}
	void resize(int _n,int _m){n=_n,m=_m;}
	Matrix operator *(const Matrix &_)const
	{
		Matrix res;
		res.n=n,res.m=_.m;
		FOR(i,1,n)FOR(j,1,_.m)
		{
			res[i][j]=0;
			FOR(k,1,m)(res[i][j]+=(a[i][k]*_.a[k][j])%P)%=P;
		}
		return res;
	}
	Matrix operator *=(const Matrix &_){return (*this)=(*this)*_;}
};
int n,d;

Matrix Pow(Matrix a,int p)
{
	Matrix res;res.resize(a.n,a.n);
	FOR(i,1,res.n)FOR(j,1,res.m)res[i][j]=(i==j);	//res初始值是一个"单位1"的矩阵
	for(;p>0;p>>=1,a*=a)if(p&1)res*=a;
	return res;
}

int main()
{
	while(scanf("%d%d%d",&d,&n,&P),d|n|P)
	{
		Matrix A,B;A.resize(1,d),B.resize(d,d);
		FOR(i,1,d)FOR(j,1,d-1)B[i][j]=(i==j+1);
		FOR(i,1,d)scanf("%d",&B[d-i+1][d]),B[d-i+1][d]%=P;
		FOR(i,1,d)scanf("%d",&A[1][i]),A[1][i]%=P;
		if(n<=d)printf("%d
",A[1][n]);
		else
		{
			A*=Pow(B,n-d);
			printf("%d
",A[1][d]);
		}
	}
    return 0;
}