跳转至

斐波那契數列

斐波那契數列(The Fibonacci sequence,OEIS A000045)的定義如下:

\[ F_0 = 0, F_1 = 1, F_n = F_{n-1} + F_{n-2} \]

該數列的前幾項如下:

\[ 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, \dots \]

盧卡斯數列

盧卡斯數列(The Lucas sequence,OEIS A000032)的定義如下:

\[ L_0 = 2, L_1 = 1, L_n = L_{n-1} + L_{n-2} \]

該數列的前幾項如下:

\[ 2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, 199, \dots \]

研究斐波那契數列,很多時候需要藉助盧卡斯數列為工具。

斐波那契數列通項公式

\(n\) 個斐波那契數可以在 \(\Theta (n)\) 的時間內使用遞推公式計算。但我們仍有更快速的方法計算。

解析解

解析解即公式解。我們有斐波那契數列的通項公式(Binet's Formula):

\[ F_n = \frac{\left(\frac{1 + \sqrt{5}}{2}\right)^n - \left(\frac{1 - \sqrt{5}}{2}\right)^n}{\sqrt{5}} \]

這個公式可以很容易地用歸納法證明,當然也可以通過生成函數的概念推導,或者解一個方程得到。

當然你可能發現,這個公式分子的第二項總是小於 \(1\),並且它以指數級的速度減小。因此我們可以把這個公式寫成

\[ F_n = \left[\frac{\left(\frac{1 + \sqrt{5}}{2}\right)^n}{\sqrt{5}}\right] \]

這裏的中括號表示取離它最近的整數。

這兩個公式在計算的時候要求極高的精確度,因此在實踐中很少用到。但是請不要忽視!結合模意義下二次剩餘和逆元的概念,在 OI 中使用這個公式仍是有用的。

盧卡斯數列通項公式

我們有盧卡斯數列的通項公式:

\[ L_n = \left(\frac{1 + \sqrt{5}}{2}\right)^n + \left(\frac{1 - \sqrt{5}}{2}\right)^n \]

與斐波那契數列非常相似。事實上有:

\[ \frac{L_n + F_n\sqrt{5}}{2} = \left(\frac{1 + \sqrt{5}}{2}\right)^n \]

也就是説,\(L_n\)\(F_n\) 恰好構成 \(\left(\frac{1 + \sqrt{5}}{2}\right)^n\) 二項式展開再合併同類項後的分子係數。也就是説,Pell 方程

\[ x^2-5y^2=-4 \]

的全體解,恰好是

\[ \frac{x_n + y_n\sqrt{5}}{2} = \frac{L_n + F_n\sqrt{5}}{2} \]

恰好是盧卡斯數列和斐波那契數列。因此有

\[ {L_n}^2-5{F_n}^2=-4 \]

矩陣形式

斐波那契數列的遞推可以用矩陣乘法的形式表達:

\[ \begin{bmatrix}F_{n-1} & F_{n} \cr\end{bmatrix} = \begin{bmatrix}F_{n-2} & F_{n-1} \cr\end{bmatrix} \begin{bmatrix}0 & 1 \cr 1 & 1 \cr\end{bmatrix} \]

\(P = \begin{bmatrix}0 & 1 \cr 1 & 1 \cr\end{bmatrix}\),我們得到

\[ \begin{bmatrix}F_n & F_{n+1} \cr\end{bmatrix} = \begin{bmatrix}F_0 & F_1 \cr\end{bmatrix} P^n \]

於是我們可以用矩陣乘法在 \(\Theta(\log n)\) 的時間內計算斐波那契數列。此外,前一節講述的公式也可通過矩陣對角化的技巧來得到。

快速倍增法

使用上面的方法我們可以得到以下等式:

\[ \begin{aligned} F_{2k} &= F_k (2 F_{k+1} - F_{k}) \\ F_{2k+1} &= F_{k+1}^2 + F_{k}^2 \end{aligned} \]

於是可以通過這樣的方法快速計算兩個相鄰的斐波那契數(常數比矩乘小)。代碼如下,返回值是一個二元組 \((F_n,F_{n+1})\)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
pair<int, int> fib(int n) {
  if (n == 0) return {0, 1};
  auto p = fib(n >> 1);
  int c = p.first * (2 * p.second - p.first);
  int d = p.first * p.first + p.second * p.second;
  if (n & 1)
    return {d, c + d};
  else
    return {c, d};
}

性質

斐波那契數列擁有許多有趣的性質,這裏列舉出一部分簡單的性質:

  1. 卡西尼性質(Cassini's identity):\(F_{n-1} F_{n+1} - F_n^2 = (-1)^n\)
  2. 附加性質:\(F_{n+k} = F_k F_{n+1} + F_{k-1} F_n\)
  3. 取上一條性質中 \(k = n\),我們得到 \(F_{2n} = F_n (F_{n+1} + F_{n-1})\)
  4. 由上一條性質可以歸納證明,\(\forall k\in \mathbb{N},F_n|F_{nk}\)
  5. 上述性質可逆,即 \(\forall F_a|F_b,a|b\)
  6. GCD 性質:\((F_m, F_n) = F_{(m, n)}\)
  7. 以斐波那契數列相鄰兩項作為輸入會使歐幾里德算法達到最壞複雜度(具體參見 維基 - 拉梅)。

斐波那契數列與盧卡斯數列的關係

不難發現,關於盧卡斯數列與斐波那契數列的等式,與三角函數公式具有很高的相似性。比如:

\[ \frac{L_n + F_n\sqrt{5}}{2} = \left(\frac{1 + \sqrt{5}}{2}\right)^n \]

\[ \cos nx + i\sin nx = \left(\cos x + i\sin x\right)^n \]

很像。以及

\[ {L_n}^2-5{F_n}^2=-4 \]

\[ \cos^2 x + \sin^2 x = 1 \]

很像。因此,盧卡斯數列與餘弦函數很像,而斐波那契數列與正弦函數很像。比如,根據

\[ \left(\frac{1 + \sqrt{5}}{2}\right)^m\left(\frac{1 + \sqrt{5}}{2}\right)^n = \left(\frac{1 + \sqrt{5}}{2}\right)^{m+n} \]

可以得到兩下標之和的等式:

\[ 2L_{m+n}=5F_mF_n+L_mL_n \]
\[ 2F_{m+n}=F_mL_n+L_mF_n \]

於是推論就有二倍下標的等式:

\[ L_{2n}={L_n}^2-2{\left(-1\right)}^n \]
\[ F_{2n}=F_nL_n \]

這也是一種快速倍增下標的辦法。同樣地,也可以仿照三角函數的公式,比如奇偶性、和差化積、積化和差、半角、萬能代換等等,推理出更多有關盧卡斯數列與斐波那契數列的相應等式。

斐波那契編碼

我們可以利用斐波那契數列為正整數編碼。根據 齊肯多夫定理,任何自然數 \(n\) 可以被唯一地表示成一些斐波那契數的和:

\[ N = F_{k_1} + F_{k_2} + \ldots + F_{k_r} \]

並且 \(k_1 \ge k_2 + 2,\ k_2 \ge k_3 + 2,\ \ldots,\ k_r \ge 2\)(即不能使用兩個相鄰的斐波那契數)

於是我們可以用 \(d_0 d_1 d_2 \dots d_s 1\) 的編碼表示一個正整數,其中 \(d_i=1\) 則表示 \(F_{i+2}\) 被使用。編碼末位我們強制給它加一個 1(這樣會出現兩個相鄰的 1),表示這一串編碼結束。舉幾個例子:

\[ \begin{aligned} 1 &=& 1 &=& F_2 &=& (11)_F \\ 2 &=& 2 &=& F_3 &=& (011)_F \\ 6 &=& 5 + 1 &=& F_5 + F_2 &=& (10011)_F \\ 8 &=& 8 &=& F_6 &=& (000011)_F \\ 9 &=& 8 + 1 &=& F_6 + F_2 &=& (100011)_F \\ 19 &=& 13 + 5 + 1 &=& F_7 + F_5 + F_2 &=& (1001011)_F \end{aligned} \]

\(n\) 編碼的過程可以使用貪心算法解決:

  1. 從大到小枚舉斐波那契數 \(F_i\),直到 \(F_i\le n\)
  2. \(n\) 減掉 \(F_i\),在編碼的 \(i-2\) 的位置上放一個 1(編碼從左到右以 0 為起點)。
  3. 如果 \(n\) 為正,回到步驟 1。
  4. 最後在編碼末位添加一個 1,表示編碼的結束位置。

解碼過程同理,先刪掉末位的 1,對於編碼為 1 的位置 \(i\)(編碼從左到右以 0 為起點),累加一個 \(F_{i+2}\) 到答案。最後的答案就是原數字。

模意義下週期性

考慮模 \(p\) 意義下的斐波那契數列,可以容易地使用抽屜原理證明,該數列是有周期性的。考慮模意義下前 \(p^2+1\) 個斐波那契數對(兩個相鄰數配對):

\[ (F_1,\ F_2),\ (F_2,\ F_3),\ \ldots,\ (F_{p^2 + 1},\ F_{p^2 + 2}) \]

\(p\) 的剩餘系大小為 \(p\),意味着在前 \(p^2+1\) 個數對中必有兩個相同的數對,於是這兩個數對可以往後生成相同的斐波那契數列,那麼他們就是週期性的。

皮薩諾週期

\(m\) 意義下斐波那契數列的最小正週期被稱為 皮薩諾週期(Pisano periods,OEIS A001175)。

皮薩諾週期總是不超過 \(6m\),且只有在滿足 \(m=2\times 5^k\) 的形式時才取到等號。

當需要計算第 \(n\) 項斐波那契數模 \(m\) 的值的時候,如果 \(n\) 非常大,就需要計算斐波那契數模 \(m\) 的週期。當然,只需要計算週期,不一定是最小正週期。

容易驗證,斐波那契數模 \(2\) 的最小正週期是 \(3\),模 \(5\) 的最小正週期是 \(20\)

顯然,如果 \(a\)\(b\) 互素,\(ab\) 的皮薩諾週期就是 \(a\) 的皮薩諾週期與 \(b\) 的皮薩諾週期的最小公倍數。

計算週期還需要以下結論:

結論 1:對於奇素數 \(p\equiv 1,4 \pmod 5\)\(p-1\) 是斐波那契數模 \(p\) 的週期。即,奇素數 \(p\) 的皮薩諾週期整除 \(p-1\)

證明:

此時 \(5^\frac{p-1}{2} \equiv 1\pmod p\)

由二項式展開:

\[ F_p=\frac{2}{2^p\sqrt{5}}\left(\dbinom{p}{1}\sqrt{5}+\dbinom{p}{3}\sqrt{5}^3+\ldots+\dbinom{p}{p}\sqrt{5}^p\right)\equiv\sqrt{5}^{p-1}\equiv 1\pmod p \]
\[ F_{p+1}=\frac{2}{2^{p+1}\sqrt{5}}\left(\dbinom{p+1}{1}\sqrt{5}+\dbinom{p+1}{3}\sqrt{5}^3+\ldots+\dbinom{p+1}{p}\sqrt{5}^p\right)\equiv\frac{1}{2}\left(1+\sqrt{5}^{p-1}\right)\equiv 1\pmod p \]

因為 \(F_p\)\(F_{p+1}\) 兩項都同餘於 \(1\),與 \(F_1\)\(F_2\) 一致,所以 \(p-1\) 是週期。

結論 2:對於奇素數 \(p\equiv 2,3 \pmod 5\)\(2p+2\) 是斐波那契數模 \(p\) 的週期。即,奇素數 \(p\) 的皮薩諾週期整除 \(2p+2\)

證明:

此時 \(5^\frac{p-1}{2} \equiv -1\pmod p\)

由二項式展開:

\[ F_{2p}=\frac{2}{2^{2p}\sqrt{5}}\left(\dbinom{2p}{1}\sqrt{5}+\dbinom{2p}{3}\sqrt{5}^3+\ldots+\dbinom{2p}{2p-1}\sqrt{5}^{2p-1}\right) \]
\[ F_{2p+1}=\frac{2}{2^{2p+1}\sqrt{5}}\left(\dbinom{2p+1}{1}\sqrt{5}+\dbinom{2p+1}{3}\sqrt{5}^3+\ldots+\dbinom{2p+1}{2p+1}\sqrt{5}^{2p+1}\right) \]

\(p\) 之後,在 \(F_{2p}\) 式中,只有 \(\dbinom{2p}{p}\equiv 2 \pmod p\) 項留了下來;在 \(F_{2p+1}\) 式中,有 \(\dbinom{2p+1}{1}\equiv 1 \pmod p\)\(\dbinom{2p+1}{p}\equiv 2 \pmod p\)\(\dbinom{2p+1}{2p+1}\equiv 1 \pmod p\),三項留了下來。

\[ F_{2p}\equiv\frac{1}{2}\dbinom{2p}{p}\sqrt{5}^{p-1}\equiv -1 \pmod p \]
\[ F_{2p+1}\equiv\frac{1}{4}\left(\dbinom{2p+1}{1}+\dbinom{2p+1}{p}\sqrt{5}^{p-1}+\dbinom{2p+1}{2p+1}\sqrt{5}^{2p}\right)\equiv\frac{1}{4}\left(1-2+5\right)\equiv 1 \pmod p \]

於是 \(F_{2p}\)\(F_{2p+1}\) 兩項與 \(F_{-2}\)\(F_{-1}\) 一致,所以 \(2p+2\) 是週期。

結論 3:對於素數 \(p\)\(M\) 是斐波那契數模 \(p^{k-1}\) 的週期,等價於 \(Mp\) 是斐波那契數模 \(p^k\) 的週期。特別地,\(M\) 是模 \(p^{k-1}\) 的皮薩諾週期,等價於 \(Mp\) 是模 \(p^k\) 的皮薩諾週期。

證明:

這裏的證明需要把 \(\frac{1+\sqrt{5}}{2}\) 看作一個整體。

由於:

\[ F_M=\frac{1}{\sqrt{5}}\left(\left(\frac{1+\sqrt{5}}{2}\right)^M-\left(\frac{1-\sqrt{5}}{2}\right)^M\right)\equiv 0\pmod {p^{k-1}} \]
\[ F_{M+1}=\frac{1}{\sqrt{5}}\left(\left(\frac{1+\sqrt{5}}{2}\right)^{M+1}-\left(\frac{1-\sqrt{5}}{2}\right)^{M+1}\right)\equiv 1\pmod {p^{k-1}} \]

因此:

\[ \left(\frac{1+\sqrt{5}}{2}\right)^M \equiv \left(\frac{1-\sqrt{5}}{2}\right)^M\pmod {p^{k-1}} \]
\[ 1\equiv\frac{1}{\sqrt{5}}\left(\frac{1+\sqrt{5}}{2}\right)^M\left(\left(\frac{1+\sqrt{5}}{2}\right)-\left(\frac{1-\sqrt{5}}{2}\right)\right)=\left(\frac{1+\sqrt{5}}{2}\right)^M\pmod {p^{k-1}} \]

因為反方向也可以推導,所以 \(M\) 是斐波那契數模 \(p^{k-1}\) 的週期,等價於:

\[ \left(\frac{1+\sqrt{5}}{2}\right)^M \equiv \left(\frac{1-\sqrt{5}}{2}\right)^M\equiv 1\pmod {p^{k-1}} \]

\(p\) 是奇素數時,由 升冪引理,有:

\[ v_p\left(a^t-1\right)=v_p\left(a-1\right)+v_p(t) \]

\(p=2\) 時,由 升冪引理,有:

\[ v_2\left(a^t-1\right)=v_2\left(a-1\right)+v_2\left(a+1\right)+v_2(t)-1 \]

代入 \(a\)\(\left(\frac{1+\sqrt{5}}{2}\right)\)\(\left(\frac{1-\sqrt{5}}{2}\right)\)\(t\)\(M\)\(Mp\),上述條件也就等價於:

\[ \left(\frac{1+\sqrt{5}}{2}\right)^{Mp} \equiv \left(\frac{1-\sqrt{5}}{2}\right)^{Mp}\equiv 1\pmod {p^k} \]

因此也等價於 \(Mp\) 是斐波那契數模 \(p^k\) 的週期。

因為週期等價,所以最小正週期也等價。

三個結論證完。據此可以寫出代碼:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
struct prime {
  unsigned long long p;
  int times;
};

struct prime pp[2048];
int pptop;

unsigned long long get_cycle_from_mod(
    unsigned long long mod)  // 這裏求解的只是週期,不一定是最小正週期
{
  pptop = 0;
  srand(time(0));
  while (n != 1) {
    __int128_t factor = (__int128_t)10000000000 * 10000000000;
    min_factor(mod, &factor);  // 計算最小素因數
    struct prime temp;
    temp.p = factor;
    for (temp.times = 0; mod % factor == 0; temp.times++) {
      mod /= factor;
    }
    pp[pptop] = temp;
    pptop++;
  }
  unsigned long long m = 1;
  for (int i = 0; i < pptop; ++i) {
    int g;
    if (pp[i].p == 2) {
      g = 3;
    } else if (pp[i].p == 5) {
      g = 20;
    } else if (pp[i].p % 5 == 1 || pp[i].p % 5 == 4) {
      g = pp[i].p - 1;
    } else {
      g = (pp[i].p + 1) << 1;
    }
    m = lcm(m, g * qpow(pp[i].p, pp[i].times - 1));
  }
  return m;
}

習題