Linear Recurrence

Recursive Function

遞迴函數就是同一個函數複合好幾次!

A recursive function f:
f(0) = 5
f(x) = 2 f(x-1) + 1

Let g(y) = 2y + 1
f(0) = 5
f(1) = g(5)
f(2) = g(g(5))
f(3) = g(g(g(5)))

Linear Recurrence

線性遞迴函數就是同一個線性函數複合好幾次。

{ f(n) = a1 ⋅ f(n-1) + a2 ⋅ f(n-2) + ... + ak ⋅ f(n-k)   when n > k
{ f(k) = ck
{   :     :
{ f(2) = c2
{ f(1) = c1

著名的例子是Fibonacci Sequence:僅兩項,係數皆是1。

{ f(n) = f(n-1) + f(n-2)   when n = 3, 4, 5, .....
{ f(2) = 1
{ f(1) = 1
{ f(1) = 1
{ f(2) = 1
{ f(3) = f(2) + f(1)
{ f(4) = f(3) + f(2)
{ f(5) = f(4) + f(3)
{   :      :      :

【註:recursion和recurrence,中文都翻譯為「遞迴」,然而兩者意義大不相同,讀者切莫混淆。】

演算法(Dynamic Programming)

由小到大一項一項慢慢算。採用「動態規劃」,有許多種實作方式。時間複雜度為O(N + (N-K)K),可以簡單寫成O(NK)。

演算法(Companion Matrix)

線性函數(線性變換)可以表示成矩陣。線性遞迴函數也可以表示成矩陣,該矩陣稱為Companion Matrix。莫名奇妙的稱呼。

1.
線性遞迴函數,由大到小的風格:
{ f(n) = c ⋅ f(n-1) + b ⋅ f(n-2) + a ⋅ f(n-3)
{ f(3) = 2
{ f(2) = 1
{ f(1) = 0

由小到大的風格:
{ f(1) = 0
{ f(2) = 1
{ f(3) = 2
{ a ⋅ f(n) + b ⋅ f(n+1) + c ⋅ f(n+2) = f(n+3)

2.
構造一個A矩陣來表示線性遞迴函數。
構造F1、F2、...、Fn來表示函數輸入。
F1經過A函數,就變成F2。

     [ f(1) ]   [ 0 ]        [ f(n)   ]       [ 0 1 0 ]
F1 = [ f(2) ] = [ 1 ]   Fn = [ f(n+1) ]   A = [ 0 0 1 ]
     [ f(3) ]   [ 2 ] ,      [ f(n+2) ] ,     [ a b c ]

         [ 0 1 0 ]   [ f(1) ]
A ⋅ F1 = [ 0 0 1 ] ⋅ [ f(2) ]
         [ a b c ]   [ f(3) ]

         [   0    + 1⋅f(2) +   0    ]   [ f(2) ]
       = [   0    +   0    + 1⋅f(3) ] = [ f(3) ] = F2
         [ a⋅f(1) + b⋅f(2) + c⋅f(3) ]   [ f(4) ]

重點在於最後一橫列。其他只是簡單位移。
A就是Companion Matrix。

3.
F2 = A ⋅ F1       = A¹ ⋅ F1
F3 = A ⋅ F2       = A² ⋅ F1
 ⋮     ⋮               ⋮
Fn = A ⋅ Fn-1     = Aⁿ⁻¹ ⋅ F1

4.
         [ f(f-2) ]
Fn-3+1 = [ f(f-1) ]     = Aⁿ⁻³ ⋅ F1
         [ f(n)   ]
           ^^^^ here is f(n)!

遞迴函數,即是同一個函數複合許多次。線性遞迴函數,即是同一個線性函數(線性變換)複合許多次──也就是矩陣次方。而矩陣次方有著效率不錯的Divide and Conquer演算法。

求得線性遞迴函數的第N項,僅需O(log(N-K))次矩陣乘法。K是線性遞迴函數的階數,也是Companion Matrix的大小。

時間複雜度O(K³log(N-K) + K²),可以簡單寫成O(K³logN)。此處假設矩陣相乘是O(K³)。

演算法(Polynomial Multiplication and Modulo)

2015年由競賽選手台灣大學李瑞珉《High Speed Linear Recursion》提出。我不清楚是否已有學術論文。

[shift]   Fn = w1 Fn-1 + ... + wk Fn-k
             = (w1 w1 Fn-2 + ... + w1 wk Fn-k-1) + ... + wk Fn-k
             = v2 Fn-2 + ... + vk-1 Fn-k-1

 RHS第一項再展開,然後同類項的係數相加。
 特色是,我們可以持續位移到我們喜歡的那k項(連續的)!
 位移的時間複雜度O(k),移到我們喜歡的那k項是 O(nk)。
 (但是以下不會用到。)

[expand]   Fn = w1 Fn-1 + ... + wk Fn-k
              = (w1 w1 Fn-2   + ... + w1 wk Fk-2  ) + ...... +
                (wk w1 Fn-k-1 + ... + wk wk Fn-k-k)
              = v2 Fn-2 + ............ + v2k Fn-2k

 RHS每個項各自再展開,然後同類項的係數相加。
 直接窮舉,時間複雜度 O(k²)。

 其實這就是多項式直式乘法(摺積),摺積有O(klogk)算法。
 (如果用快速數論轉換去算摺積,如何保證w1...wk不會爆表呢?)
 (所以只好請出中國餘數定理)

[expand]   Fn = w1 Fn-x-1   + ... + wk  Fn-x-k
              = v2 Fn-x-y-2 + ... + v2k Fn-x-y-2k   第二次展開,也可以用別的
              = uk Fn-y-k-1 + ... + u2k Fn-x-y-2k   縮成剛好k項

 因為有了[shift]的概念,所以我們取的那k項
 沒必要是 Fn-1  ... Fn-k
 也可以是 Fn-x-1 ... Fn-x-k
 最後再用 k-1次 [shrink] 縮成剛好k項

[shrink]   Fn = v1 Fn-1 + v2 Fn-2 + ... + vk Fn-k + ... + vK Fn-K
              =           u2 Fn-2 + ... + uk Fn-k + ... + uK Fn-K

 項數大於k的時候,RHS第一項再展開,可以縮減一項。
 縮一項的複雜度 O(k)。
 縮到剩下k項的時間複雜度是 O(k(K-k))。
 
 其實這就是多項式長除法(減法變加法。除數乘上一個負號)
 最後要取餘數。
 既然已經採取快速數論轉換,那就是要取模了。
 取模的多項式餘數,有O(klogk)算法。

令那k項是 Fn/2-1 ... Fn/2-k,總是取一半的位置。
如此一來,仿照整數次方(Exp)或餘數次方(ModExp)的演算法。
利用 [expand] [shrink] ,分割 logn 次。
     O(klogk) O(klogk)
求第n項只需要 O(klogk logn) 時間。

總結

K階的線性遞迴函數,計算第N項:

一、Dynamic Programming。逐步計算第一項到第N項,答案儲存於表格。時間複雜度為O(NK)。

二、Companion Matrix的次方。直接算第N項。時間複雜度為O(K³logN);假設矩陣相乘是O(K³)。

三、多項式乘法、多項式餘數。直接算第N項。時間複雜度為O(K²logN),可以加速為O(KlogKlogN)。

當K是常數,一變成O(N),二和三變成O(logN)。

四、Generating Function。找到一般式,直接算第N項。詳情請查閱離散數學書籍。

UVa 10229 10518 10655 10743 10754 10870 12653