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}=\displaystyle\sum_{k=1}^pA_{i,k}B_{k,j}\) .

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

\[\begin{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\n",A[1][n]);
		else
		{
			A*=Pow(B,n-d);
			printf("%d\n",A[1][d]);
		}
	}
    return 0;
}