斐波那契數的O(logn)尾遞迴推導

尾遞迴 · 發表 2018-10-21 20:19:15

摘要：最近剛看《SICP》，第1章的練習1.19上有一題挺有意思的，用對數複雜度來實現斐波那契演算法，題目中不僅給了一些提示還給了程式碼的大部分，但第一次看到還是挺納悶的。我思考了一下，嘗試推導題目中的演算法，解釋一下它能夠運作的原理。首先介紹一下分解的方法，假設存在...

最近剛看《SICP》，第1章的練習1.19上有一題挺有意思的，用對數複雜度來實現斐波那契演算法，題目中不僅給了一些提示還給了程式碼的大部分，但第一次看到還是挺納悶的。

我思考了一下，嘗試推導題目中的演算法，解釋一下它能夠運作的原理。

首先介紹一下分解的方法，假設存在，那麼:

關於它的推導有多種方法，比如說利用 Line"/> 這個定義式也可以推出來，所以就不再贅述。

當n是偶數的時候，設k為n/2，且定義那麼

當n為奇數，設 ,且定義那麼

直接照著這樣的形式，不斷減半就就可以找出把fib(n)降低到的方法：

function fib(n: number): number {
if (n < 2) return n;
if (n % 2 === 0) {
const l = fib(n / 2 - 1);
const r = fib(n / 2);
return r * (r + 2 * l)
} else {
const l = fib((n - 1) / 2);
const r = fib((n + 1) / 2);
return l ** 2 + r ** 2;
}
}

雖然遞迴的方式可以很容易寫出來，但如果想寫出尾遞迴或者迭代的形式就比較困難。

《SICP》裡特別提到，如果想做成迭代，找出一個 不變數 是一個關鍵，也就是說在每次迭代的時候都能維持一個不變。

在迭代形式下的對數演算法，當n為奇數的時候，一般會通過遞減1的形式降低到偶數，當n為偶數的時候再通過減半的形式降低n，直到n降到0.

按照這個思路，我們可以分析出

n為奇數時， ,不變數的其他部分必須得能同時增加1，我們假定這個用來吸收奇數部分的變數為k，此時k+1。
n為偶數時，，不變數的其他部分必須能夠吸收到原本的n/2，假定這個變數為x。

所以的形式可以滿足上面的需求

n為奇數的時候
n為偶數的時候

所以我們可以得到不變數的形式，在上面任何步驟它的值都不會變，且：

初始時設，所以
在n降到0的時候，此時 ,所以可以通過返回fib(k)來得到結果

分解可以得到

設，那麼上面式子可以化成

的形式不好分解，但是由於

我們可以把它轉換成只包含fib(x)、fib(x-1)和n的形式，同理也可以轉換成fib(n)、fib(n-1)和x的形式，所以我們可以認定存在一個，使得

設，可得

其中變數定義為

下面開始實現：

偶數的時候，，基於的形式可得

所以當，結果應當轉移到上，所以

通過把代入到可得

所以可以得到

奇數的時候， ,結果應當轉移到上面

由可得

所以可得

所以可以得到

下面開始用程式碼實現：

function fib_iter(a: number, b: number, p: number, q: number, n: number): number {
if (n === 0) return b;
if (n % 2 === 0) {
const _p = p ** 2 + q ** 2;
const _q = 2 * p * q + q ** 2;
return fib_iter(a, b, _p, _q, n / 2)
} else {
const _a = (b * q) + (a * q) + (a * p);
const _b = (b * p) + (a * q);
return fib_iter(_a, _b, p, q, n - 1)
}
}

這個形式就已經和習題裡給出來的程式碼如同一轍了，而且可以很容易看清楚它的演算法複雜度只有。

這裡面的變數定義如前所述為

且初始化值如前所述為 ,所以計算初始值可得：

的函式實現就是

function fib(n: number): number {
return fib_iter(1, 0, 0, 1, n)
}

當然也可以把尾遞迴的形式變成迴圈：

function fib(n: number): number {
let a = 1;
let b = 0;
let p = 0;
let q = 1;
while (n > 0) {
if (n % 2 === 0) {
const _p = p ** 2 + q ** 2;
const _q = 2 * p * q + q ** 2;
p = _p;
q = _q;
n = n / 2;
} else {
const _a = (b * q) + (a * q) + (a * p);
const _b = (b * p) + (a * q);
a = _a;
b = _b;
n = n - 1;
}
}
return b;
}

順便一提，習題裡有這麼一句話：

Now consider T to be the special case of