Page 201 : 專欄 更快地計算遞迴式

事實上，要求m項遞推式的第n項的值可以不使用矩陣，而是使用初項的線性表示，通過快速冪在O(m ^ 2 * log n)的時間內求出答案。有興趣的讀者可以試著思考看看。

看到這裡，我首先想到的是《SICP》裡有類似的內容，在Page31，練習1.19。

下文第一段是一次“失敗”的思考，我最後發現自己並沒有給出O(m ^ 2 * log n)的演算法，而是自己重新建立了對矩陣冪的理解。

經過一天的思考（在缺覺的頭疼中），我終於可以給出問題的一個解答，見第二段

從問題本身推演到矩陣乘積（O(m ^ 3 * log n)）

以斐波那契為例，首先有遞推式

Fn = Fn-1 + Fn-2

注意，這裡所討論的遞推式的一般結構不僅僅侷限於

fn = a1 * fn-1 + a2 * fn-2 ...

可以是更廣義的

Ai = a1 * Ai-1 + b1 * Bi-1 ...
Bi = a2 * Ai-1 + b2 * Bi-1 ...
...

回來，討論斐波那契。

為了便於理解，變換問題的形式，對於Fn = Fn-1 + Fn-2，對應有變換f1接受一個元組(a, b)得到(c, d)，使得c = a + b, d = a。

例如，輸入是(Fn-1, Fn-2)，結果就是(Fn, Fn-1)。

我們有了給出第i項值給出第i+1項值的變換f1。

f1跟原遞推式不同的地方在於它是一個變換，有統一的輸入和輸出。

現在，我們可以定義類同於函式複合的操作，例如

f3(x) = f1(f1(f1(x)))

即

f3 = f1 · f1 · f1 = f1 ^ 3

於是，可以理解

fn = f1 ^ n

我們現在定義了變換f1的冪乘，但是還不知道如何實現快速冪。

關鍵是如何定義fi(i > 1)上的·（乘法）操作，即

fi+j = fi · fj

回頭觀察斐波那契的例子，定義f1的，實際上是

c = a + b
d = a

更明白一點，是

c = 1 * a + 1 * b
d = 1 * a + 0 * b

可以觀察到一個固定的結構，而所有的fi都基於這一結構，構成一族變換。

對於任意fi，有唯一一組p,q,s,t，使

c = p * a + q * b
d = s * a + t * b

給定fi{pi, qi, si, ti}和fj{pj, qj, sj, tj}，類比函式複合，定義fi+j = fi · fj

c = pi * (pj * a + qj * b) + qi * (sj * a + tj * b)
d = si * (pj * a + qj * b) + ti * (sj * a + tj * b)

調整一下，

c = a * (pi * pj + qi * sj) + b * (pi * qj + qi * tj)
d = a * (si * pj + ti * sj) + b * (si * qj + ti * tj)

即

pi+j = pi * pj + qi * sj
qi+j = pi * qj + qi * tj
si+j = si * pj + ti * sj
ti+j = si * qj + ti * tj

這就是m=2時變換的一般結構。

在這裡我注意到這個部分的複雜度是O(m ^ 3)（簡單地說，此處m取2，有8個乘）。

實際上，上面的等式就是矩陣相乘的展開。

[pi+j, qi+j]        [pi, qi]        [pj, qj]
               =                *   
[si+j, ti+j]        [si, ti]        [sj, tj]

這樣，可以理解矩陣其實是變換的一種非常自然的表現形式。

注意，因為我們的例子“斐波那契”的遞推式中沒有表現常數項，所以這裡的結構是不包括常數項的。

O(m ^ 2)

經過思考，我發現自己擴大了問題域，即上段一開始給出的廣義結構。

如果我們僅著眼於形如

Fn = a * Fn-1 + b * Fn-2 ...

這樣的遞推式，問題可以進一步簡化。

上段所提出的模型的問題在於它維護了多條“軌道”，具體說，是m條軌道。

當m=2時，p，q維護了一條軌道，s，t維護一條。

對於我們現在所關注的遞推式，這種維護其實是不必要的。

m條軌道的意義在於我們必須維護m個值以支援變換操作。

但是

Fn = a * Fn-1 + b * Fn-2 ...

可以共用一條軌道：初始有m項初值和f1

Fm+1 = f1(F1 ... Fm)

繼續直到

F2m-1 = f1(Fm-1 ... F2m-2)

現在，我們有F1 ... F2m-1。

對於任意fi，都可以通過帶入這m個m連續區間，得到

Fx+j = fi(Fj ... Fj+m)

即fi變換後又有m個值，可以支援下一步變換（fj），因為我們只關心變換的結構，我們這次不再擴充套件，只做一個值的變換

Fy = fj(Fx ... Fx+m)

接下來就可以像上一段一樣計算引數了。

以斐波那契為例

對於斐波那契，有

//任意的初值組
a
b

//初值擴充所使用的常數，就是f1的結構
p1 = 1
q1 = 1

做一次（m-1次)擴充

c = p1 * a + q1 * b

對fi{pi, qi}變換

x = pi * a + qi * b
y = pi * b + qi * c

對fj{pj, qj}變換

z = pj * x + qj * y

化簡

pj * x  = pj * (pi * a + qi * b)
        = pi * pj * a + qi * pj * b

qi * c  = qi * (p1 * a + q1 * b)
        = p1 * qi * a + q1 * qi * b

qj * y  = qj * (pi * b + qi * c)
        = qj * (pi * b + p1 * qi * a + q1 * qi * b)
        = p1 * qi * qj * a + pi * qj * b + q1 * qi * qj * b

z       = pi * pj * a + qi * pj * b
            + p1 * qi * qj * a + pi * qj * b + q1 * qi * qj * b
        = (pi * pj + p1 * qi * qj) * a
            + (2 * qi * pj + q1 * qi * qj) * b

即有

pi+j = pi * pj + p1 * qi * qj
qi+j = 2 * qi * pj + q1 * qi * qj

這就是m=2時

Fn = a * Fn-1 + b * Fn-2 ...

的一般變換組合結構

其中，p1和q1是常數，唯一對應於f1{p1, q1}。

複雜度分析

這個變換顯然比之前的矩陣形式有更少的乘積項，不過每一個乘積項乘的次數增多了（上例為7）

經過觀察可以發現，2, p1, q1都是常數項，共有m ^ 2個乘積項，即只有不超過m ^ 2個常數項，不影響複雜度。

而每個乘積項中之多有2個變數出現（觀察上式，計算x，y時引入一個，計算z時引入另一個，這裡也解釋為什麼只有m ^ 2個乘積項），

所以單次變換複合的複雜度為O(m ^ 2)。

注意，這裡的變換都是以m的倍數為跨度的，比如f1通過1...m個初值得到第m+1項，f2則得到第2m + 1項