对数级矩阵斐波那契数列

前两天11月23日，一打开电脑就说是斐波那契日。也许是巧合！

$f(n)=f(n-1)+f(n-2)$

这是个时间复杂度很大的递归调用算法。按照前面算法复杂度文章的计算方法。应该是$T(n)=2\ast T(n/1.01)+O(n^0)$，问题规模几乎没有缩小，递归外的运算都是常数级。大多书上也是把它作为算法复杂度的例子。

import math
math.log(2,1)
Traceback (most recent call last):
  File "<pyshell#7>", line 1, in <module>
    math.log(2,1)
ZeroDivisionError: float division by zero
math.log(2,1.01)
69.66071689357483
math.log(2,1.1)
7.272540897341713
math.log(2,1.2)
3.8017840169239308

$b=1$是不可能的，问题规模还是有缩小的。1.01的时候是$n^{69}$,1.2的时候是$n^{3.8}$,它一定是个指数级的时间复杂度算法。

代码0

import time
#import functools
#@functools.lru_cache(maxsize = 128, typed = False)
def f(n):
    assert n >= 0 #"n>0"
    if n <= 1:
        return n
    return f(n - 1) + f(n - 2)
s = time.time()
a=f(40)
print('递归调用所用时间',a,time.time()-s)

计算第40个斐波那契数用了11.7秒。

递归调用所用时间 102334155 11.706673383712769
import functools
@functools.lru_cache(maxsize = 128, typed = False)

当然后面两行代码打开的话，求400也是用时0秒，4000的话就报超过调用深度的错了。

递归的改造

兔子数列一节有它的迭代实现方法。里面还使用了存储的方法，提高运行速度。

STL

正文如下：

图1：对数级矩阵斐波那契数列

实现起来好像不难。

代码1

def fibonacci_iterative(n):
    if 0==n:return 0
    v=[(0,1)]
    #std:: pair<int ,int> v ={0,1};
    for i in range(1,n):
        v.append((v[i-1][1],v[i-1][0]+v[i-1][1]))
    return v.pop()

用元组代替向量，我们不改变元组的值，只是构造的时候用前面一个向量的两个值。接下来就要用到矩阵了。

矩阵

图2：矩阵

教课书上也有矩阵，主要讲的是它的表示方法。但也不确定老师没有讲过它可以用来实现对数级的斐波那契数列。还好numpy直接有矩阵的操作。

代码2

import numpy as np

def logntuzi(n):
    matrix2 = np.array([[1],[0]])
    matrix1 = np.array([[1,1],[1,0]])
    for _ in range(1,n):
        matrix1 = np.dot(matrix1,matrix1 )
    return np.dot(matrix1,matrix2)

这段代码是错的。
#matrix1 = np.matmul(matrix1,matrix1)
#matrix1 = np.multiply(matrix1,matrix1)
另外两种乘法也都试了。同样用memorizeFib验证前面的代码1是可以的。先放一放，是不是要再仔细看看书。

递归调用

书上讲递归调用的时候一般用汉诺塔，为什么不用斐波那契数列。图1中的代码并没有保存中间变量，只用了一个向量变量。Python的代码也可以的吧。

代码3

def fibonacci_iterative(n):
    if 0==n:return 0
    v=[0,1]
    #std:: pair<int ,int> v ={0,1};
    for i in range(1,n):
        v[0],v[1] = v[1],v[0]+v[1]
    return v

是的。只要改成数组，同时修改两个元素变量值就可以了，Python正好可以这样做。应该很多聪明的同学都会发现这样的方法，所以它不能做递归调用的例子。其实代码2只是一个很小的错误。

代码4

import numpy as np
def logntuzi(n):
    matrix2 = np.array([[1],[0]])
    matrix1 = np.array([[1,1],[1,0]])
    matrix3 = np.array([[1,1],[1,0]])
    for _ in range(1,n):
        matrix1 = np.dot(matrix1,matrix3)
        #matrix1 = np.matmul(matrix1,matrix1)
        #matrix1 = np.multiply(matrix1,matrix1)
    return np.dot(matrix1,matrix2)

添加一个matrix3，让matrix1连乘不变的matrix3，最后再乘以matrix2。为什么要添加一个matrix3？直觉告诉我这还是$n$数量级的算法，一个$for$循环一直乘。要去找个幂的算法。

泛型编程之乘法

或许网上有，辗转反侧几天。埃及乘法的对半再叠加给了我灵感。
从欧几里得算法到这一章有现代数论的梅森、业余数学家费马、欧拉、模运算、抽象群、拉格朗日，很多的理论知识。

代码5

def power7_2(n,matrix1 = np.array([[1,1],[1,0]]),matrix2 = np.array([[1,1],[1,0]])):
    while True:

这里应该留白，惊喜留给有缘人。只四五行代码，在前面的代码上改几个字母。

验证

在50的时候出现了负数。

power7_2(50)
array([[-1408458269, -1109825406],
       [-1109825406,  -298632863]])

四个数字都是负数。但是40的时候是可以的。

power7_2(40)
array([[267914296, 165580141],
       [165580141, 102334155]])

这几个数字是对的吗？可以用前面的方法验证一下。

fibonacci_iterative(40)
[63245986, 102334155]
fibonacci_iterative(41)
[102334155, 165580141]
fibonacci_iterative(42)
[165580141, 267914296]

没有问题，他们分别是第40，41，42个斐波那契数。用memorizeFib()我也测试过。在确定它是正确的之前我就写好了测试速度的代码。

代码6

def test(n):
    import time
    s = time.time()
    a=memorizeFib(n)
    print('memorizeFib',a,time.time()-s)
    s = time.time()
    a=fibonacci_iterative(n)
    print('fibonacci_iterative',a,time.time()-s)
    s = time.time()
    a=power7_2(n-2)
    print('power7_2',a,time.time()-s)

第一万个斐波那契数是一个有2091位数的很大的数字。但是只用了0.0022923946380615234秒，代码5应该用时更短的。但它总是负数或不对，改了好多次才得到代码5。
测试结果：

test(10000)
memorizeFib 33644764876431783266621612005107543310302148460680063906564769974680081442166662368155595513633734025582065332680836159373734790483865268263040892463056431887354544369559827491606602099884183933864652731300088830269235673613135117579297437854413752130520504347701602264758318906527890855154366159582987279682987510631200575428783453215515103870818298969791613127856265033195487140214287532698187962046936097879900350962302291026368131493195275630227837628441540360584402572114334961180023091208287046088923962328835461505776583271252546093591128203925285393434620904245248929403901706233888991085841065183173360437470737908552631764325733993712871937587746897479926305837065742830161637408969178426378624212835258112820516370298089332099905707920064367426202389783111470054074998459250360633560933883831923386783056136435351892133279732908133732642652633989763922723407882928177953580570993691049175470808931841056146322338217465637321248226383092103297701648054726243842374862411453093812206564914032751086643394517512161526545361333111314042436854805106765843493523836959653428071768775328348234345557366719731392746273629108210679280784718035329131176778924659089938635459327894523777674406192240337638674004021330343297496902028328145933418826817683893072003634795623117103101291953169794607632737589253530772552375943788434504067715555779056450443016640119462580972216729758615026968443146952034614932291105970676243268515992834709891284706740862008587135016260312071903172086094081298321581077282076353186624611278245537208532365305775956430072517744315051539600905168603220349163222640885248852433158051534849622434848299380905070483482449327453732624567755879089187190803662058009594743150052402532709746995318770724376825907419939632265984147498193609285223945039707165443156421328157688908058783183404917434556270520223564846495196112460268313970975069382648706613264507665074611512677522748621598642530711298441182622661057163515069260029861704945425047491378115154139941550671256271197133252763631939606902895650288268608362241082050562430701794976171121233066073310059947366875 0.0022923946380615234
fibonacci_iterative [20793608237133498072112648988642836825087036094015903119682945866528501423455686648927456034305226515591757343297190158010624794267250973176133810179902738038231789748346235556483191431591924532394420028067810320408724414693462849062668387083308048250920654493340878733226377580847446324873797603734794648258113858631550404081017260381202919943892370942852601647398213554479081823593715429566945149312993664846779090437799284773675379284270660175134664833266377698642012106891355791141872776934080803504956794094648292880566056364718187662668970758537383352677420835574155945658542003634765324541006121012446785689171494803262408602693091211601973938229446636049901531963286159699077880427720289235539329671877182915643419079186525118678856821600897520171070499437657067342400871083908811800976259727431820539554256869460815355918458253398234382360435762759823179896116748424269545924633204614137992850814352018738480923581553988990897151469406131695614497783720743461373756218685106856826090696339815490921253714537241866911604250597353747823733268178182198509240226955826416016690084749816072843582488613184829905383150180047844353751554201573833105521980998123833253261228689824051777846588461079790807828367132384798451794011076569057522158680378961532160858387223882974380483931929541222100800313580688585002598879566463221427820448492565073106595808837401648996423563386109782045634122467872921845606409174360635618216883812562321664442822952537577492715365321134204530686742435454505103269768144370118494906390254934942358904031509877369722437053383165360388595116980245927935225901537634925654872380877183008301074569444002426436414756905094535072804764684492105680024739914490555904391369218696387092918189246157103450387050229300603241611410707453960080170928277951834763216705242485820801423866526633816082921442883095463259080471819329201710147828025221385656340207489796317663278872207607791034431700112753558813478888727503825389066823098683355695718137867882982111710796422706778536913192342733364556727928018953989153106047379741280794091639429908796650294603536651238230626, 33644764876431783266621612005107543310302148460680063906564769974680081442166662368155595513633734025582065332680836159373734790483865268263040892463056431887354544369559827491606602099884183933864652731300088830269235673613135117579297437854413752130520504347701602264758318906527890855154366159582987279682987510631200575428783453215515103870818298969791613127856265033195487140214287532698187962046936097879900350962302291026368131493195275630227837628441540360584402572114334961180023091208287046088923962328835461505776583271252546093591128203925285393434620904245248929403901706233888991085841065183173360437470737908552631764325733993712871937587746897479926305837065742830161637408969178426378624212835258112820516370298089332099905707920064367426202389783111470054074998459250360633560933883831923386783056136435351892133279732908133732642652633989763922723407882928177953580570993691049175470808931841056146322338217465637321248226383092103297701648054726243842374862411453093812206564914032751086643394517512161526545361333111314042436854805106765843493523836959653428071768775328348234345557366719731392746273629108210679280784718035329131176778924659089938635459327894523777674406192240337638674004021330343297496902028328145933418826817683893072003634795623117103101291953169794607632737589253530772552375943788434504067715555779056450443016640119462580972216729758615026968443146952034614932291105970676243268515992834709891284706740862008587135016260312071903172086094081298321581077282076353186624611278245537208532365305775956430072517744315051539600905168603220349163222640885248852433158051534849622434848299380905070483482449327453732624567755879089187190803662058009594743150052402532709746995318770724376825907419939632265984147498193609285223945039707165443156421328157688908058783183404917434556270520223564846495196112460268313970975069382648706613264507665074611512677522748621598642530711298441182622661057163515069260029861704945425047491378115154139941550671256271197133252763631939606902895650288268608362241082050562430701794976171121233066073310059947366875] 0.0
power7_2 [[1242044891  890489442]
 [ 890489442  351555449]] 0.0

也许numpy的矩阵乘法就只能计算十位数的数，所以这个测试也就没有意义了。可以自己写个power方法吗？