矩阵乘幂优化k 阶常系数线性递推关系一、认识k 阶常系数线性递推关系我们熟悉的Fibonacci 数列：F[n]=F[n-1]+F[n-2]，就是一个2阶常系数线性递推关系，由此我们得出k 阶常系数线性递推关系的一般形式：k n k n n n F a F a F a F −−−+++=Λ2211其中：，是常数，有k 项，所以叫着k 阶常系数线性递推关系；k a a a 、、、Λ21∑=−×=k i in i n F a F 1二、对矩阵的认识矩阵就是一个数字阵列，一个n 行r 列的矩阵可以表示为： ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡nr n n r r a a a a a a a a a ΛΛΛΛ212222111211 我们称上面的矩阵为r n ×矩阵。

例如下面一个2×3矩阵； ⎥⎦⎤⎢⎣⎡662341如果一个行数和列数相等的矩阵，我们叫作方阵。

例如下面3×3方阵：⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡679762341 其实，矩阵对我们来说，并不陌生，因为它正好对应Pascal 中的一个二维数组。

Type matrix=array[1..n,1..r] of longint;三、矩阵的运算1、加法，减法⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡964687652342312345⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡−⎥⎦⎤⎢⎣⎡0100036523426623452、乘法：设A ，B 是两个矩阵，令C=A ×B ；那么：（1）A 的列数必须和B 的行数相等；设A 是n ×r 矩阵，B 是r ×m 矩阵；（2）A 和B 的乘积C 是一个n ×m 的矩阵；（3）∑=×=×++×+×+×=r k jk k i jr r i j i j i j i j i b a b a b a b a b a c 1,,,,,33,,22,,11,,Λ例如：已知： ，，求A ×B 。

⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=635241A ⎥⎦⎤⎢⎣⎡=654321B ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡+++++++++=×4536273629222722176*63*35*62*34*61*36*53*25*52*24*51*26*43*15*42*14*41*1B A 由矩阵乘法的运算法则，我们得出下列结论：矩阵乘法满足结合律，不满足交换律。

矩阵乘法的程序实现：function mul(a,b:matrix;n,r,m:integer):matrix; //矩阵乘法var c:matrix;i,j,k:longint;beginfillchar(c,sizeof(c),0);for i:=1 to n dofor j:=1 to m dofor k:=1 to r doc[i,j]:=c[i,j]+a[i,k]*b[k,j]mul:=c;end;练习：并观察有什么样的结论（和F[n]=F[n-1]+F[n-2]）已知： ，，求A ×B ；A ×B ×B ；A ×B ×B ×B 。

[]21−−=n n F F A ⎥⎦⎤⎢⎣⎡=0111B 3、方阵乘幂方阵乘幂是指，A 是一个方阵，将A 连乘n 次，即：。

n A C =【想一想】：如果A 不是方阵，能否进行n 次幂运算？【练习试题】： ⎥⎦⎤⎢⎣⎡=0111A ，求 3A C =由于矩阵乘法满足“乘法结合律”，因此我们可以用二分快速幂的思想来求方阵乘幂。

（1）求a^n 的二分快速幂程序//程序中的mm 为取余数function power(a,n:longint):longint;var x:longint;beginif n=1 then power:=a elsebeginx:=power(a,n shr 1) mod mm;x:=(x*x) mod mm;if odd(n) then power:=(x*a) mod mmelse power:=xend;end;（2）求矩阵A^xfunction mpower(a:matrix;n,r,m,x:longint):matrix;var c:matrix;beginif x=1 then mpower:=a elsebeginc:=mpower(a,n,r,m,x shr 1);c:=mul(c,c,n,r,m);if odd(x) then mpower:=mul(c,a,n,r,m)else mpower:=c;end;end;四、线性递推方程的优化算法【例题一】：Fibonacci 数列求Fibonacci 数列的第n 项，其中1<=n<=10^9，第一项是1，第2项是1。

输出Fibonacci 数列的第n 项Mod 2008 的值。

【问题分析】：考虑1×2的矩阵[f[n-2],f[n-1]]。

根据fibonacci 数列的递推关系，我们希望通过乘以一个2×2的矩阵，得到矩阵[f[n-1],f[n]]=[ f[n-1]，f[n-2]+f[n-1]]，很容易构造出这个2×2矩阵A ，即：[][]⎥⎦⎤⎢⎣⎡×−−=−1110]1[]2[][]1[n F n F n F n F 所以，有[f[1],f[2]]×A ＝[f[2],f[3]]又因为矩阵乘法满足结合律，故有： [][][][][]1111110111110]2[]1[]2[]1[1110]1[]2[][]1[−−−⎥⎦⎤⎢⎣⎡×=⎥⎦⎤⎢⎣⎡×=×=⎥⎦⎤⎢⎣⎡×−−=−n n n F F A F F n F n F n F n F这个矩阵的第一个元素即为所求。

至于如何快速求出A n-1，相信大家都会，即递归地：n 为偶数时，A n =(A n/2)2；n 为奇数时，A n =(A n/2)2*A 。

【例题二】：数列f[n]=f[n-1]+f[n-2]+1,f[1]=f[2]=1的第n 项的快速求法, 输出该数列的第n 项Mod 2008 的值。

【试题分析】：仿照前例，考虑1×3的矩阵 [f[n-2],f[n-1],1]，希望求得某3×3的矩阵A ，使得此1×3的矩阵乘以A 得到矩阵：[f[n-1],f[n],1]＝[f[n-1],f[n-1]+f[n-2]+1,1]容易构造出这个3×3的矩阵A ，即： ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=110011010A2110011010 f[2],1][f[1],110011010 1],1]-f[n 2],-[f[n A 1],1]-f[n 2],-[f[n 1,1]2]-f[n 1]-f[n 1],-[f[n f[n],1]1],-[f[n −⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡×==⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡×=×=++=n 五．常见例题【典型例题1】：线性递推式【问题描述】动态规划的实现形式之一是递推，因此递推在OI 中十分重要。

在某信息学的分支学科中，LC 学会了如何求一阶线性递推数列。

由于他想在正在学习主干学科，因此希望知道求出N 阶线性递推数列。

为此，他了解到了以下的内容：一个N 阶线性递推式是这样的式子：也就是说，这个数列的每一项都是由他之前连续N 项相加所得。

其中还包括一个常数An 。

例如，当N=2，A0=A1=1，A2=0时，这个式子就是我们熟悉的斐波那契数列。

当然，作为边界条件，F0、F1……Fn-1都是已知的。

LC 对如何去求这个式子一筹莫展，因此请你来帮助他。

你的任务就是对于一个给定的N 阶线性递推式，求出它的第K 项是多少。

【文件输入】第一行两个整数：N ，K ，其中N 表示这个式子是N 阶线性递推式，K 表示你需要求的那一项。

第二行有N+1个整数：A0,A1,…，An ，表示这个递推式的系数。

第三行有N 个整数：F0，F1，….，Fn-1，表示数列的初始值。

【文件输出】只有一行，其中只有一个整数，表示这个数列第K 项的值。

由于数据较大，你只需输出结果 mod 9973的值。

【样例输入】2 101 1 00 1【样例输出】55【数据规模】对于50%的数据：N<=K<=10^6对于100%的数据：1<N<=10N<=K<=101^81<=Ai ，Fi<=104【试题分析】线性递推式求第k 项可转换成矩阵连乘如：已知起始项F1,F2,F3，F n =2F n-1+3*f n-2+4F n-3+6，求第K 项值列矩阵方法1：]6,1*6322314,3,2[1100021003010400*]6,3,2,1[]6,4,3,2[+++=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=F F F F F F F F F F F ，，，，，，，，，，，，或]1,1*6322314,3,2[1600021003010400*]1,3,2,1[]6,4,3,2[+++=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=F F F F F F F F F F F ，，，，，，，，，，，，因此第K 项的值是[F1,F2,F3,6]* 矩阵的k-1次方得到的矩阵第一项的值。

构造方法，若不考虑常数项，则n 阶递推式构造n*n 的矩阵，从第二行第一列开始的(n-1)*(n-1)的子矩阵为对角线为1的矩阵，第n 列从上到下为n 阶递推式由前到后的n 个系数,其它项为0。

若有常数项，则列在最后，并增加矩阵维数。

注意：列项时由前到后列如，[F n-2 F n-1 F n ]列矩阵方法2：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+++=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡66322314326321*10001234010000106432F F F F F F F F F F F ，，，，，，，，，，，，同样，第K 项的值是矩阵的k-1次方得到的矩阵第一项的值。

构造方法，若不考虑常数项，则n 阶递推式构造n*n 的矩阵，从第一行第二列开始的(n-1)*(n-1)的子矩阵为对角线为1的矩阵，第n 行为n 阶递推式由前到后的n 个系数,其它为项为0。

注意：列项时由前到后列如，[F n-2 F n-1 F n ]因此，本题的方法核心是构造矩阵。

本题由于项是从第0项开始的，所以第k 项即为求矩阵的k 次方。

const m=9973;mxn=9;type matrix=array[1..mxn+1,1..mxn+1]of longint;var a,b:matrix;f,ans,n,i,j,c:longint;k:int64;function multiply(a,b:matrix):matrix;var c:matrix;i,j,k:longint;beginfillchar(c,sizeof(c),0);for i:=1 to n+1 dofor j:=1 to n+1 dofor k:=1 to n+1 doc[i,j]:=(c[i,j]+a[i,k]*b[k,j])mod m;exit(c);end;function pow(i:int64):matrix;var t:matrix;beginif i=1 then exit(a);t:=pow(i div 2);t:=multiply(t,t);if odd(i) then t:=multiply(t,a);exit(t);end;beginreadln(n,k);for i:=2 to n do a[i,i-1]:=1;for i:=1 to n+1 dobeginread(c);a[i,n]:=c;end;a[n+1,n+1]:=1;b:=pow(k);ans:=b[n+1,1];for i:=1 to n dobeginread(f);ans:=(ans+f*b[i,1])mod m;end;writeln(ans);end.【典型例题1】：病毒入侵【问题描述】H5N1型高致病性禽流感击了bzbz国，不可避免的，bzbz 国的大量鸡死于流感。