MATLAB实现矩阵快速幂

武汉大学

MATLAB及其应用课程

课程论文

学生姓名：

学号：

电话：

班级：

指导教师：

2016年4月4日

目录

一、实验题目 (3)

二、实验目的 (3)

三、实验设备和条件 (3)

四、实验原理 (3)

五、实验内容与过程 (6)

六、程序代码截图 (8)

七、实验结果 (10)

八、实验收获与总结 (12)

一、实验题目

利用矩阵快速幂算法与符号运算求解大型斐波那契数。

二、实验目的

<1>实现在较短的时间内完成大规模运算。

<2>充分利用MATLAB方便的矩阵运算以及超高精度的符号运算的特色，进行对大规模数值的求解运算。

三、实验设备和条件

<1>设备：MATLAB R2015a软件

<2>条件：

①MATLAB拥有着极其方便的矩阵运算，每一个数组都可以表示为相应维度的矩阵，可以通过两个数组的直接相乘来获得“矩阵1×矩阵2”运算的结果矩阵。

②MATLAB拥有着精度超高的符号运算，而且使用起来非常方便。

四、实验原理

根据题目，“利用矩阵快速幂算法与符号运算求解大型斐波那契数”，需要将该实验的题目拆成“矩阵快速幂算法”、“大型”和“斐波那契数”几个关键点进行原理的阐述。

为了方便描述，稍排了个序进行原理讲解：

<1>“斐波那契数”

斐波那契数列（Fibonacci sequence），又称黄金分割数列、因数学家列昂纳多·斐波那契（Leonardoda Fibonacci ）以兔子繁殖为例子而引入，故又称为“兔子数列”，指的是这样一个数列：0、1、1、2、3、5、8、13、21、34、……在数学上，斐波纳契数列以如下被以递推的方法定义：F（0）=0，F（1）=1，F（n）=F(n-1)+F(n-2)（n≥2，n∈N*）在现代物理、准晶体结构、化学等领域，斐波纳契数列都有直接的应用。

而斐波那契数，就是指在已经定义了递推初始值的前提下，持续递推斐波那契数列所得的其中一个数。

<2>“大型”

对于诸如C，C++，JAVA等语言，进行高精度数字的求解总是非常具有挑战性的，即使只是进行大规模整数运算，C++也需要铺写大量的代码定义大数类和重载操作符。如果再加入矩阵运算，使用起来极其不方便。甚至在大数类定义不够严谨的情况下，轻则造成运行效率极低，数组空间溢出等问题，重则造成数据溢出还找半天找不到bug的“惨案”。

而MATLAB的符号类精度之高非常恐怖，且操作异常方便，可以像处理普通整型一样地

处理符号类。因此，使用MATLAB可以让“大型”运算变得具有非常强的可行性。

所谓“大型”，当然不会是求第10个、第20个、第100个斐波那契数这等简单循环便可以轻易得到结果的数据（当然，其实第100个斐波那契数已经达到了3.54×1020的整型存放不下的大数）。

所谓“大型”，应该是指第一万个、第十万、百万，乃至第千万、第一亿个斐波那契数。根据下图所示，斐波那契数列有着跟指数函数类似的曲线，到达一定程度后会有类似“指数爆炸”的曲线斜率。

<3>“矩阵快速幂算法”

“矩阵快速幂算法”又可细分为“矩阵”和“快速幂”：

1．“矩阵”

斐波那契数的递推式可以用一个2×2的矩阵A来表示。

A=a b

c d =01

11

初始值状态为

10 01×

fib(1)

fib(2)

根据线性代数知识，容易得以下表达式

fib(2) fib(3)=

fib2+0

fib1+fib(2)

=A×

fib(1)

fib(2)

fib(3) fib(4)=

fib1+fib(2)

fib2+fib(3)

=

fib1+fib(2)

fib1+2×fib(2)

=A2×

fib(1)

fib(2)

以此类推，第n个斐波那契数的值为1×2的ans0矩阵的第二个值，在MATLAB中表示的方式为ans0(2)：

ans0=A n−2×fib(1) fib(2)

2．“快速幂”

在讲述快速幂原理之前上几张图：

图一、图二、图三分别计算了求解第103个、第104个、第105个斐波那契数时候所用的时间。分别用了0.393秒，3.517秒以及34.5秒。不难推测，第106个斐波那契数可能需要花费350秒左右的时间。因为该算法的运算次数随着横坐标n的变化是线性的。

而我们要进行对第106、第107、乃至第108个斐波那契数的求解，用这种简单粗暴的运算方法是行不通的。因此，利用快速幂算法将O（N）的运算优化到O（log2N）来节省运算时间是非常有必要的。

那么，快速幂究竟是什么算法呢。

首先，我们来看看幂运算：

X n=X1+1+ ……+1 （n个1的和）=X∙X∙……∙X (n个X的积)

通常我们都会使用计算机的循环语句进行计算，这样，当我们计算一个n次方的幂值时，所花费的时间复杂度就是O（N）。因此，如果用每次求下一个斐波那契数的时候，都要乘上一个2×2的矩阵A。用矩阵代替递推公式求解斐波那契数需要进行n次的矩阵乘法。所花费的时间与上例三张例图相比好不了多少，甚至可能会花费更长时间。

A n=A1+1+ ……+1 （n个1的和）=A×A×……×A (n个A的积,A表示矩阵)

这时候就要探讨如何减少循环次数进行大规模的幂运算了。

我们考虑一个问题，当n为偶数时，设n=2k，我们是否可以把每次循环时候乘上的单元定为A2呢？这样可以表示为

A n=A2+2+ ……+2 （k个2的和）=A2×A2×……×A2 (k个A2的积,A表示矩阵)

这样使用循环计算，就可以节省n/2次的运算，虽然时间复杂度仍然为O（N），对于大规模运算几乎没有多少贡献，但也为我们提供了一个思考方向——改变累乘的单位（因为我讲的是矩阵快速幂，因此在之后的文章中，累乘单位称为unitMx），以节省循环次数——快速幂的核心想法因此诞生。

之前我们令unitMx=A与unitMx=A2称为静态地改变，即初始化之后，在循环中就不再变化。这样做的收效是十分微小的。即使当n%10==0的时候，令unitMx=A10，也无法动摇大规模运算的代价。当然，当n为平方数的时候，即n=t2的时候，令unitMx=A t，可以获得