[C#] Fibonacci sequence 費氏數列

-

最近把 Codewars : The Millionth Fibonacci Kata 解出來，想說都出到 BigInteger 還考慮了負整數 n，大概也夠完整了吧，就想把這個主題整理一下。

費氏數列定義

費波那契數列 (Fibonacci sequence) 或常簡稱費氏數列，是中學數學常見、程式設計界也常見的主題，他有的一些大自然現象（像兔子、蜜蜂數量與黃金比例）我就不討論了，純粹討論數學與程式。費氏數列是如下的數列：

0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, ...

規律是「每一項是前兩項的和」，例如 8 = 3 + 5，或 144 = 55 + 89，而「最初兩項是 0 和 1」（數學版本因為足標正整數從 1 開始，故最初兩項為 1 和 1，我這裡討論程式版本 index 從 0 開始，故最初兩項為 0 和 1），其遞迴關係式為：

費氏數列中的數稱為費氏數 (Fibonacci number)，若有一函數 F(n) 表費氏數列中 index = n 的數，例如 n = 6 的 a6 即 F(6) = 8，則此 F(n) 函數定義為：

這個 F(n) 的取法就成為常見問題。

遞迴函數 Recursive Function

顯然上面已經是個遞迴函數的 base case 和 recursive case 了，所以程式碼：

  public int Fib(int n)
  {
    if (n == 0)                         // base case
    {
      return 0;
    }
    if(n == 1)                          // base case
    {
      return 1;
    }
    return Fib(n - 2) + Fib(n - 1);     // recursive case
  }

這是個很好訓練 Recursion 的練習，但這個作法效能太差，例如我求 Fib(40) 花了近 5 秒才得到 102334155 ，更不用說 50、100，如果目標想很到更大的結果使用 long 甚至 BigInteger，但效能追不上，所以再思考別的方法。

費氏數的通解 Binet’s Formula

在高中數學的數列章節，有個費氏數列一般項公式與證明，是用數學歸納法求證的，這裡我不討論他的推導與證明，只使用他：

可參考 Mathematics LibreTexts : 10.4: Fibonacci Numbers and the Golden Ratio 。在 1843 年由數學家 雅可比內 Jacques Philippe Marie Binet 提出，故稱 Binet's Formula，可由 n 直接得到費氏數，若將那兩個無理數分式令為 Φ 和 φ 可得簡潔版本：

在使用這個方法前先看一下需要用到的兩個命名空間 System 下的方法：

  • Math.Sqrt(n) 回傳 n 的平方根，資料型態為 double，例如要得到 √5 就是 Math.Sqrt(5)
  • Math.Pow(a,b) 回傳 a 的 b 次方，資料型態為 double，例如要得到 35 就是 Math.Pow(3,5)

程式碼：

  using System;
  public int Fib(int n)
  {
    double Phi = (1 + Math.Sqrt(5)) / 2;        // Φ
    double phi = (1 - Math.Sqrt(5)) / 2;        // φ
    return (int)((Math.Pow(Phi, n) - Math.Pow(phi, n)) / Math.Sqrt(5));
  }

雖然一次取到費氏數是非常快的方法，但他有準確度的問題，如果取 n = 90 (上面要用 long 回傳）得到的是 2880067194370824704，可是實際上第 90 個費氏數 為 2880067194370816120，如果 n 不太大（例如 50 以內）是既快速且準確的方法，但 n 很大時會因為 double 的準確度開始有誤差，雖然雙精度浮點算很精確的，但畢竟還是逼近，差之毫釐失之千里。

排列組合等價問題

回想高中排列組合有一種問題：10 階的階梯，如果每次有兩種走法：一次走 1 階和一次走 2 階，那有幾種走完階梯的方法？看下列的基本例子：

  • 共 1 階，有 (1) 共 1 種走法。
  • 共 2 階，有 (2), (1,1) 共 2 種走法。
  • 共 3 階，有 (1,2), (2,1), (1,1,1) 共 3 種走法。
  • 共 4 階，有 (2,2), (1,1,2), (1,2,1), (2,1,1), (1,1,1,1) 共 5 種走法。
  • 共 5 階，有 (1,2,2), (2,1,2), (1,1,1,2), (2,2,1), (1,1,2,1), (1,2,1,1), (2,1,1,1), (1,1,1,1,1) 共 8 種走法。
  • 共 6 階，有 (2,2,2), (1,1,2,2), (1,2,1,2), (2,1,1,2), (1,1,1,1,2), (1,2,2,1), (2,1,2,1), (1,1,1,2,1), (2,2,1,1), (1,1,2,1,1), (1,2,1,1,1), (2,1,1,1,1), (1,1,1,1,1,1) 共 13 種走法。

會發現走法數1, 2, 3, 5, 8, 13 ... 正是費氏數列，原因是用最後一步來分類，最後一步有兩種走法，就是一次走 2 階的走法，和一次走 1 階的走法，

若總共要走 n 階，走法分成下面兩種：

  • 走法之中最後一步是走 2 階的話代表前面有 ( n - 2 ) 階要走，這些走法就是「共 ( n - 2 ) 階的走法，最後再一次走 2 階」
  • 走法之中最後一步是走 1 階的話代表前面有 ( n - 1 ) 階要走, 這些走法就是「共 ( n - 1 ) 階的走法，最後再一次走 1 階」

以上面「共 5 階的 8 種走法」：(1,2,2), (2,1,2), (1,1,1,2), (2,2,1), (1,1,2,1), (1,2,1,1), (2,1,1,1), (1,1,1,1,1)，可以分類成兩種：

  • 「最後一步是走 2 階」：(1,2,2), (2,1,2), (1,1,1,2)，這些走法也就是「共 3 階的走法，再多一步走 2 階」
  • 「最後一步是走 1 階」：(2,2,1), (1,1,2,1), (1,2,1,1), (2,1,1,1), (1,1,1,1,1) ，這些走法也就是「共 4 階的走法，再多一步走 1 階」

得到結論「共 5 階的走法數」即「共 4 階的走法數」+「共 3 階的走法數」：

  • 共 3 階，有 (1,2), (2,1), (1,1,1) 共 3 種走法。
  • 共 4 階，有 (2,2), (1,1,2), (1,2,1), (2,1,1), (1,1,1,1) 共 5 種走法。
  • 共 5 階
    • 最後一步走 2 階，有 (1,2,2), (2,1,2), (1,1,1,2)
    • 最後一步走 1 階，有 (2,2,1), (1,1,2,1), (1,2,1,1), (2,1,1,1), (1,1,1,1,1)

以此類推，即得到費氏數列的遞迴關係式 F(n) = F(n-1) + F(n-2)，所以我們的費氏數問題可以轉化成這種走樓梯問題，尋求 ( n - 1 ) 階的走法有多少。

而這種問題還有一個排列組合的解法，如果要解共 5 階的走法數，假設走 1 階的 x 次，走 2 階的 y 次，即 1 ‧ x + 2 ‧ y = 5，求 x + 2y = 5 的非負整數解，共有 (x,y) = (1,2), (3,1), (5,0) 三組，分別再求「1,2,2」的排列數 3!/2! = 3，求「1,1,1,2」的排列數 4!/3! = 4，求「1,1,1,1,1」的排列數 5!/5! = 1，得走法數共 3 + 4 + 1 = 8 種。

故要解 F(n) 第 n 個費氏數，可考慮成「共走 ( n - 1 ) 階的問題」轉化成等價的排列組合解法，去求「x + 2y = n - 1 的 x , y 非負整數解，再分別求出排列數 C(x + y, x)」，取總和」。

這個等價解法速度快也精確，就是前置作業稍微多了點，要有存放 x, y 非負整數解數對的物件：

  class Pair
  {
    public int x { get; set; }
    public int y { get; set; }
  }

然後是計算排列數的 Comibination() 方法，雖然是 x 個 1 和 y 個 2 的不盡相異物之排列數，但可換成 C ( x + y, x ) 的組合數，若選擇求出所有階乘的作法，數字偏大且再一次遞迴，又再被約分，我個人認為是蠻不必要的浪費，所以選擇求組合數的作法，只要分母分子依序列出等約分即可。程式碼：

  private int Combination(int n, int r)
  {
    r = (r > n - r) ? (n - r) : r;    // 餘組合定理，C(n,r) = C(n,n-r) 取小的算
    int result = 1;
    for(int i = 0; i < r; i++)
    {
      result = result * (n - i) / (i + 1);
    }
    return result;
  }

最後是 Fib() 的內容，程式碼：

  public int Fib(int n) 
  {
    int target = n - 1;
    int count = target / 2 + 1;        // 非負整數解個數
    Pair[] solutions = new Pair[count];

    for (int i = 0; i < count; i++)
    {
      solutions[i] = new Pair
      {
        x = 2 * i + target % 2,        // target 為奇數時 x = 1,3,5 ...
        y = target / 2 - i             // target 為偶數時 x = 0,2,4,...
      };
    }

    int result = 0;
    for(int i = 0; i < count; i++)
    {
      Pair pair = solutions[i];
      result += Combination(pair.x + pair.y, pair.x);  // 計算C( x + y, y )
    }
    return result;
  }

其實這個解法是意外中想到的，我本來這篇目標是寫下面要談的矩陣解法，想起之前教高中數學時有這個等價問題，嘗試將他應用在程式上發現效率頗佳，故附上。以上都可以改成回傳 BigInteger 更進一步。

矩陣與log(N)求費氏數

這個解法是在 Codewars : The Millionth Fibonacci Kata 中的提示 1.2 Procedures and the Processes They Generate 1.2.4 Exponentiation ，其下的 Exercise 1.19. ，利用 Divide-and-conquer algorithm ，其實就是我之前寫的 Pow(x,n)實作 想到的對半步進求法，後來看演算法的書才知道這是一個固定的解問題思維。

不過 Codewars 的提示，最開始我很疑惑這和費氏數有什麼關係，因為費氏數是用累加得出來的，並不像指數 Pow() 的實作中是利用連乘出來的，去看前人的 Discussion 發現線性組合化為矩陣的作法（我竟然沒聯想到！），如果是矩陣的連乘，確實可以將上面的對半步進作法用上。考慮每個費氏數是由前面兩個費氏數之和求得：

以 F(8) 為例是前兩個 F(7) 和 F(6) 之和，而 F(7) 當然是前兩個 F(6) 和 F(5) 之和，但這裡不這樣寫，F(7) 還是 F(7)，然後將這兩個式子用 F(7) 與 F(6) 來表示：

組合數對分別為 (1,1) 和 (1,0)，用這兩個列向量做成一個二階方陣，則 F(8), F(7) 組成行向量就是用這個方陣與 F(7), F(6) 組成的行向量做矩陣乘法得來的：

同理，F(7), F(6) 組成的行向量，可用此二階方陣乘上 F(6), F(5) 組成的行向量得到：

將上式代入上上式，就可逐一化簡：

反覆操作，將相同的二階方陣 [ [ 1 , 1 ] , [ 1 , 0 ] ] 合併成次方 ：

最後就可得到「 F(8), F(7) 所成的行向量」用「 F(1), F(0) 所成的行向量」表示的結果：

令 M = [ [ 1 , 1 ] , [ 1 , 0 ] ] ，則「 F(n), F(n-1) 所成的行向量」可利用 M 的 (n-1) 次方再乘上「 F(1), F(0) 所成的行向量」來求得：

由以上推論，就有了方向：

  • 想求 F(n) 第 n 個費氏數
  • 先找出 M 矩陣的 n-1 次方之結果
  • 要求 Mn-1，利用對半步進，如果 n-1 是偶數，找 M(n-1)/2，再平方
  • 如果 n-1 是奇數，找 M(n-2)/2，平方後再多乘一個 M
  • Mn-1 再乘上 [ [ F(1) ] , [ F(0) ] ] 後的第一列即 F(n)，也就是 Mn-1 的第一列第一行元素

首先先做出矩陣乘法，這裡我不考慮二維陣列，反正二階方陣也才四個元素，我將 [ [ a, b ] , [ c , d ] ] 視為 { a , b , c , d } ，而二階方陣的乘法為：

故程式碼：

  public int[] MatrixMultiply(int[] A, int[] B)
  {
    int a = A[0], b = A[1], c = A[2], d = A[3];
    int p = B[0], q = B[1], r = B[2], s = B[3];

    return new[] { a * p + b * r, a * q + b * s,
                   c * p + d * r, c * q + d * s };
  }

而求 M 的 n 次方，就利用對半步進來遞迴，base case 是 0 次與 1 次，recursive case 要判斷 n 是偶數還是奇數，

  • n 是偶數，n = 2 * k，則先求 k 次，用 k 次的結果來平方。
  • n 是奇數，n = 2 * k + 1，則先求 k 次，用 k 次的結果平方並多乘 1 次

求 M 的 n 次方的程式碼：

  public int[] MatrixExp(int n)
  {
    if (n == 0)
    {
      return new int[] { 1, 0, 0, 1 };       // base case
    }
    if (n == 1)
    {
      return new int[] { 1, 1, 1, 0 };       // base case
    }
    int[] matrixHalf = MatrixExp(n / 2);     // 求 k 次
    int[] matrixN = MatrixMultiply(matrixHalf, matrixHalf);

    if (n % 2 == 0)               // n 為偶數時 n/2 + n/2 = n
    {
      return matrixN;
    }
    else                          // n 為奇數時 n/2 + n/2 + 1 = n
    {
      return MatrixMultiply(matrixN, new[] { 1, 1, 1, 0 });
    }
  }

重要的遞迴完成後，求費氏數的函數本身就只要呼叫他即可，這邊要注意的是 F(0) 要獨自給他，還有本來還有一步 Mn-1 要乘以「 F(1), F(0) 成的行向量」，但是剛好 F(1)=1, F(0)=0 ，故乘開結果的第一列，其實就是 Mn-1 的第一列第一行，也就是 [0] 元素，故我直接回傳即完成。

Fib() 程式碼：

  public int Fib(int n)
  {
    if (n == 0)
    {
      return 0;
    }
    return MatrixExp(n - 1)[0];
  }

這個是該 Codewars kata 的提示作法，通關之後看大家作法都差不多，有興趣的人可以改成 BigInteger 來看他的強大。

留言

張貼留言