Mortal Fibonacci Rabbits (ID: FIBD)

Problem

假設兔子性成熟需一個月，性成熟後每對兔子每個月必繁殖一對子代，若兔子的壽命為 m 月，求第 n 月兔子對數。

Given: Positive integers n≤100 and m≤20.

Return: The total number of pairs of rabbits that will remain after the n-th month if all rabbits live for m months.

# input
6 3
# output
4

Background

此處模仿先前克服 Rabbits and Recurrence Relations (ID: FIB) 的推理方式來解題 [1]。令 $F_{n}$ 為第 n 個月的兔子總數數，$R_{n}$ 與 $r_{n}$ 分別為第 n 個月成熟和未成熟的兔子數。則

$$ F_{n} = R_{n} + r_{n} \tag{1} $$ 

其中第 n 個月的子代數取決於第 n - 1 個月的成熟兔子數，所以

$$ r_{n} = R_{n-1} \tag{2} $$

假設 n > m + 1，則 $R_{n}$ 則是第 n - 1 個月的成熟兔子數加上未成熟兔子數，再扣掉陽壽已盡的兔子數。其中，第 n 個月時，陽壽已盡的兔子數即為第 n - m 個月的未成熟兔子數。

$$ R_{n} = R_{n-1} + r_{n-1} - r_{n - m} \tag{3}$$

將 (2) 和 (3) 代入 (1)，可得

$$\begin{aligned} 

F_{n} &= R_{n-1} + r_{n-1} - r_{n-m} +  R_{n-1} \\

&= F_{n-1} - r_{n-m} + R_{n-1} \tag{4} \\

\end{aligned}$$

由於 (3)，$R_{n-1}$ 可表示為

$$R_{n-1} = R_{n-2} + r_{n-2} - r_{n - m - 1} \tag{5} $$

因此，(4) 可轉換為

$$\begin{aligned} 

F_{n} &= F_{n-1} - r_{n-m} + R_{n-1}  \\

&= F_{n-1} - r_{n-m}  + R_{n-2} + r_{n-2} - r_{n - m - 1} \\

&= F_{n-1} + R_{n-2} + r_{n-2} - r_{n-m} - r_{n - m - 1} \\

&= F_{n-1} + F_{n-2} - (r_{n-m} + r_{n - m - 1}) \tag{6}

\end{aligned}$$

由於 (2)，$r_{n-m}$ 可表示為 $R_{n-m-1}$，因此

$$\begin{aligned} 

F_{n} &= F_{n-1} + F_{n-2} - (R_{n-m-1} + r_{n - m - 1}) \\

&= F_{n-1} + F_{n-2} - F_{n-m-1} \tag{7}

\end{aligned}$$

即第 n 月兔子對數為第 n - 1 和 n - 2 月兔子數之和，扣掉第 n - m - 1 個月的兔子數。

前述推導的是 n > m + 1 的狀況，而 n < m + 1 （沒有兔子死掉，等同一般的費氏數列）和 n = m + 1（最初的兔子死掉）的狀況如下：

$$F_{n}=

\begin{cases} 

F_{n-1} + F_{n-2}& \text{, if n < m + 1}\\ 

F_{n-1} + F_{n-2} - 1& \text{, if n = m + 1}

\end{cases}$$

Solution

處理實際案例時，隨著 $n$ 增加，$F_{n}$ 的數值也會增加。當輸入數值的位數超過儲存數值的位元限制，便無法保留所有位數的資訊。在此情況下，任何計算都會遺漏位數，導致浮點數誤差 [2]。

Python 會自動處理大數運算碰到的浮點數誤差，但 R 要額外動用 gmp package（用於大數運算的套件）[3]，將原生的數據類型轉換為gmp定義的數據類型，才能避免同樣的問題。

Python

由於 Python 計數從 0 開始，所以條件判斷的比較對象為 m 而非 m + 1。

def FibD (n, m):
    Fn = [1, 1]
    for i in range(2, n):
        if i < m:
            tmp = 0
        elif i == m:
            tmp = 1
        else:
            tmp = Fn[i - m - 1]
        Fn.append(Fn[i - 1] + Fn[i - 2] - tmp)
    return Fn[-1]

R

R 的計數則從 1 開始，所以依照前一節的推論寫條件判斷即可。為了避免浮點數誤差，所有參與計算的數值皆須先用 as.bigz() 傳換為 bigz 類型。

library(gmp)
FibD <- function (n, m) {
  Fn <- as.bigz(c(1, 1))
  for (i in 3:n) {
    if (i <= m) {
      tmp <- as.bigz(0)
    } else if (i == m + 1) {
      tmp <- as.bigz(1)
    } else {
      tmp <- Fn[i - m - 1]
    }
    Fn <- c(Fn, Fn[i - 1] + Fn[i - 2] - tmp)
  }
  return(Fn[n])
}


Discussion


1. 其他的解題方式

除了以代數運算來推論各個月兔子對數的關係以外，另一種策略是先追蹤全數兔子的年齡，再統計當前仍存活的兔子數量。 在以下的範例中，abeliangrape 以 list 儲存各年齡的兔子對數。其中，list 的序號對應年齡，而 list 的長度則反映最大壽命。即 list 的首項儲存 0 歲的兔子對數（子代對數）、次項儲存 1 歲的兔子對數、……以此類推。

因為兔子在成熟後每個月都會繁殖一對子代，所以下個月的子代對數即為這個月所有成熟兔子的年齡和。這個概念對應到迴圈中的 [sum(ages[1:])]，即透過[1:] 排除未成熟的兔子。 而下個月的成熟兔子對數則為這個月的兔子對數，扣掉下個月壽命已盡的兔子對數。這個概念對應到迴圈中的 ages[:-1]，即透過 [:-1] 排除下個月將死亡的兔子。

最後加總各年齡的兔子對數，即取得當前的兔子總數。

# by Rosalind user: abeliangrape (https://rosalind.info/users/abeliangrape/)
def fib(n,k=1):
  ages = [1] + [0]*(k-1)
  for i in xrange(n-1):
    ages = [sum(ages[1:])] + ages[:-1]
  return sum(ages)

2. 大數計算時遭遇的浮點數誤差

浮點數誤差 (floating point number) 是指數值超出浮點數位元限制所致的資訊喪失。在 R 語言中，numeric 類別以雙精度浮點數 (double, double precision floating point numbers) 儲存數值，因此無論儲存或運算，numberic 類別皆無法保證第 16 位以後的數值精確。

11111111111111111
[1] 11111111111111112
> 11111111111111111 + 11111111111111111
[1] 22222222222222224

然而，浮點數誤差是欲以有限位元表示實數的必然限制，並非 R 或 Python 特有的狀況。目前，通用的 IEEE 754 分別定義了以 32 bit (單精度浮點數) 和 64 bit (雙精度浮點數) 儲存實數的規範，詳細介紹可參考從 IEEE 754 標準來看為什麼浮點誤差是無法避免的 及 浮點數誤差IEEE-754。

以單精度浮點數為例，在儲存數值前，實數要先轉換為二進位制，再正規化為$\pm 1.mmmmm\times 2^e$形式。IEEE 754 將 32 bit 劃分為符號 (sign)、指數 (exponent) 和尾數 (mantissa) 區段，每個區段都儲存正規化形式的一部分。

  • 符號 (sign)：佔 1 bit，對應數值的正負號($\pm$)
  • 指數 (exponent)：佔 8 bit，對應正規化後的指數 ($e$)，此指數須以超127格式轉換為二進位制。
  • 尾數 (mantissa)：佔 23 bit，對應正規化後的小數 ($mmmmm$)

由於位元數有限，凡位數大於浮點數儲存限制的數值皆無法完整保存。此時，除了與誤差共存或完全使用十進位計算以外，也可使用位元數較多的雙精度浮點數。雙精度浮點數有 64 bit，其中指數佔 11 bit，尾數佔 52 bit，能更精確地表達數值。

不過，儘管雙精度浮點數原則上能表示 $2^{64}$ 種數值，但為了劃分符號、指數和尾數，實際上的有效位數僅有 52 個。 52 位數有 $2^{52}$ 種可能，換算為十進位得 4,503,599,627,370,496，共 16 位數。這解釋了為什麼 R 語言的 numeric 只能確保數值前 16 位的精確。

除了大數值以外，具有無限小數的其他數值，浮點數也只能以近似值表示，例如，

  • 無限循環小數：$\frac{1}{7}$ $=0.\overline{142857}$
  • 無理數（無限不循環小數）：$\pi$、$\sqrt{3}$、$e$
  • 轉換為二進位制後成為無限小數的數值：0.1、2.3、8.9

關於浮點數誤差，還可參考以下文章：

3. R: gmp package

Python 內建大數運算的功能，但 R 需要載入 gmp 套件來克服大數運算引發的浮點數誤差。R 原生的數據類型經轉換為 gmp 物件後，即可照常以二元運算子或function操作。

as.bigz(num)# 整數
as.bigq(num)# 有理數

須留意的是，此套件只保證 gmp 物件的精確性。若轉換前數值已超出浮點數儲存範圍，需要先將數字以 character 形式輸入，以免在轉換格式前便損失精確性。

# input with numeric object
> as.bigz(11111111111111111)
Big Integer ('bigz') :
[1] 11111111111111112

# input with character object
> as.bigz("11111111111111111")
Big Integer ('bigz') :
[1] 11111111111111111

關於 gmp 的介紹和操作，可參考 really large numbers in r 和 R: wrong arithmetic for big numbers

於

沒有留言:

張貼留言

訂閱： 張貼留言 (Atom)
Back to top