矩阵乘幂优化k阶常系数线性递推关系

矩阵乘幂优化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 i

n 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

⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡−⎥⎦⎤⎢⎣

⎡010003652342662345

2、乘法：

设A ，B 是两个矩阵，令C=A ×B ；那么：

（1）A 的列数必须和B 的行数相等；设A 是n ×r 矩阵，B 是r ×m 矩阵；

（2）A 和B 的乘积C 是一个n ×m 的矩阵；

（3）

∑=×=×++×+×+×=r k j

k k i j

r 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;

begin

fillchar(c,sizeof(c),0);

for i:=1 to n do

for j:=1 to m do

for k:=1 to r do

c[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;

begin

if n=1 then power:=a else

begin

x:=power(a,n shr 1) mod mm;

x:=(x*x) mod mm;

if odd(n) then power:=(x*a) mod mm

else power:=x

end;

end;

（2）求矩阵A^x

function mpower(a:matrix;n,r,m,x:longint):matrix;

var c:matrix;

begin

if x=1 then mpower:=a else

begin

c:=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]]

又因为矩阵乘法满足结合律，故有： []

[][][][]11

11110111110]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 ，即： ⎥

⎥⎥⎦⎤⎢⎢⎢⎣⎡=110011010A