docs/math/combinatorics/fibonacci.md
斐波那契数列（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})$．
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};
}
斐波那契数列拥有许多有趣的性质，这里列举出一部分简单的性质：
不难发现，关于卢卡斯数列与斐波那契数列的等式，与三角函数公式具有很高的相似性．比如：
$$ \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，对于编码为 1 的位置 $i$（编码从左到右以 0 为起点），累加一个 $F_{i+2}$ 到答案．最后的答案就是原数字．
对于模 $m$ 意义下的斐波那契数列，可以容易地使用抽屉原理证明，该数列是有周期性的．由于斐波那契数每一项的计算都依赖于前两项的取值，所以需要用相邻斐波那契数组成的数对描述数列当且所处的状态．考虑模意义下前 $m^2+1$ 个斐波那契数对：
$$ (F_0,\ F_1),\ (F_1,\ F_2),\ \ldots,\ (F_{m^2},\ F_{m^2 + 1}) $$
模 $m$ 的剩余系大小为 $m$，这意味着至多只可能有 $m^2$ 种互不相同的数对．因此，在前 $m^2+1$ 个数对中必有两个相同的数对，于是从这两个数对可以往后生成相同的斐波那契数列．那么，斐波那契数列就是周期性的，且（最小正）周期不会超过 $m^2$．
模 $m$ 意义下斐波那契数列的最小正周期被称为 Pisano 周期（Pisano period，皮萨诺周期，OEIS A001175）．本文中用 $\pi(m)$ 表示模 $m$ 的 Pisano 周期．
这一观察可以用于计算第 $n$ 项斐波那契数模 $m$ 的值．如果 $n$ 非常大，就需要计算斐波那契数模 $m$ 的周期．当然，只需要计算周期，不一定是最小正周期．
为此，有如下结论：
综合这些情形，可以说明：模 $m$ 的 Pisano 周期不会超过 $6m$，等号当且仅当 $m = 2\times 5^e~(e\in\mathbf N_+)$ 时取得．
利用上述结论，可以基于素因数分解算法，得到如下快速计算 Pisano 周期的方法：
??? example "参考代码"
cpp --8<-- "docs/math/code/combinatorics/fibonacci/pisano_estimate.cpp:pisano"
这样得到的周期可能只是 Pisano 周期的一个倍数．要得到精确的 Pisano 周期，可以进一步考察该周期的因数；或者，可以直接通过 BSGS 算法 以 $O(\sqrt{m})$ 的时间复杂度计算．
最后，本文简要证明上述关于 Pisano 周期的结论．值得说明的是，利用下文说明的方法，类似的结论可以推广到一般的二阶常系数线性齐次递推数列．尽管具体的常数有所差异，这些数列模 $m$ 的 Pisano 周期都是 $O(m)$ 的．
第一个观察是：利用 中国剩余定理，可以将讨论限制在素数幂模的情形．设 $m_1,m_2$ 是两个互素的模数．斐波那契数列在模 $m_1$ 下的周期是 $\pi(m_1)$ 及其倍数，在模 $m_2$ 下的周期是 $\pi(m_2)$ 及其倍数，所以它在模 $m_1m_2$ 下的最小正周期则恰为 $\pi(m_1)$ 和 $\pi(m_2)$ 的最小公倍数．这就是前文的结论 1．
另一个观察是：模 $m$ 下的 Pisano 周期，其实是最小的正整数 $k$，使得
$$ A^k = \begin{pmatrix} 1&1\1&0 \end{pmatrix}^k \equiv I \pmod{m}. $$
对于素数幂模 $m=p^e$ 的情形，可以通过经典的升幂论证联系到相应的素数模的情形．设 $k=\pi(p^e)$，就存在二阶方阵 $\Lambda$，使得
$$ A^k = p^e\Lambda + I $$
成立．故而，由 二项式定理 可知
$$ A^{kp} = (p^e\Lambda + I)^p = I + \sum_{i=1}^p\binom{p}{i}(p^e\Lambda)^i \equiv I\pmod{p^{e+1}}. $$
因此，由 阶的性质，有 $\pi(p^{e+1})\mid kp = p\pi(p^e)$．对 $e$ 归纳可知，$\pi(p^e)\mid p^{e-1}\pi(p)$ 总是成立．
对于素数模 $p$ 的情形，本文讨论两种证明方式．
=== "利用通项公式" 一种是利用斐波那契数列的通项公式：
$$
F_n = \dfrac{1}{\sqrt{5}}\left(\dfrac{1+\sqrt{5}}{2}\right)^n - \dfrac{1}{\sqrt{5}}\left(\dfrac{1-\sqrt{5}}{2}\right)^n.
$$
将它用二项式定理展开，并消去根式项：
$$
F_n = \dfrac{1}{2^{n-1}}\sum_{i=0}^{\lfloor(n-1)/2\rfloor}\binom{n}{2i+1}5^i.
$$
对于 $p=2$，这一表达式无法直接取模，但可以验证对应的 Pisano 周期为 $\pi(2)=3$．对于 $p=5$，有 $F_n\equiv n\cdot 3^{n-1}\pmod{p}$，可以直接验证对应的 Pisano 周期为 $\pi(5)=20$．对于剩余的奇素模数，可以分为两种情形：
- 如果 $p\equiv 1,4\pmod{5}$，就有
$$
\begin{aligned}
F_{p} &\equiv \dfrac{1}{2^{p-1}}\binom{p}{p}5^{(p-1)/2} \equiv 1 \pmod{p},\\
F_{p+1} &\equiv \dfrac{1}{2^p}\left(\binom{p+1}{1} + \binom{p+1}{p}5^{(p-1)/2}\right) \equiv 1 \pmod{p}.
\end{aligned}
$$
化简过程中，利用了如下结论：由 [Lucas 定理](../number-theory/lucas.md)，对于 $0 < k < p$ 都有 $\dbinom{p}{k}\equiv 0\pmod{p}$，而对于 $1 < k < p$ 都有 $\dbinom{p+1}{k}\equiv 0\pmod{p}$；由 [Fermat 小定理](../number-theory/fermat.md#费马小定理)，有 $2^{p-1}\equiv 5^{p-1}\equiv 1\pmod{p}$；对于 $p\equiv 1,4\pmod{5}$，都有 $p$ 是模 $5$ 的二次剩余，利用 [二次互反律](../number-theory/quad-residue.md#二次互反律)，也有 $5$ 是模 $p$ 的二次剩余，故而 $5^{(p-1)/2} \equiv 1\pmod{p}$．由此，有 $(F_p,F_{p+1}) \equiv (F_1,F_2) \pmod{p}$，所以 $(p-1)$ 是模 $p$ 的一个周期．所以，$\pi(p)\mid(p-1)$．
- 如果 $p\equiv 2,3\pmod{5}$，就有
$$
\begin{aligned}
F_{2p} &\equiv \dfrac{1}{2^{2p-1}}\binom{2p}{p}5^{(p-1)/2} \equiv -1 \pmod{p},\\
F_{2p+1} &\equiv \dfrac{1}{2^{2p}}\left(\binom{2p+1}{1} + \binom{2p+1}{p}5^{(p-1)/2} + \binom{2p+1}{2p+1}5^p\right) \equiv -1\pmod{p}.
\end{aligned}
$$
化简过程中，利用了如下结论：由 Lucas 定理，对于 $0 < k < p$ 和 $p < k < 2p$ 都有 $\dbinom{p}{k}\equiv 0\pmod{p}$，以及 $\dbinom{2p}{p}\equiv 2\pmod{p}$，而对于 $1 < k < p$ 和 $p + 1 < k < 2p$ 都有 $\dbinom{p}{k}\equiv 0\pmod{p}$，以及 $\dbinom{2p+1}{p}\equiv 2\pmod{p}$；由 Fermat 小定理，有 $2^{p-1}\equiv 5^{p-1}\equiv 1\pmod{p}$；对于 $p\equiv 2,3\pmod{5}$，都有 $p$ 是模 $5$ 的二次非剩余，利用二次互反律，也有 $5$ 是模 $p$ 的二次非剩余，故而 $5^{(p-1)/2} \equiv -1\pmod{p}$．由此，有 $(F_{2p},F_{2p+1}) \equiv (F_{-2},F_{-1}) \pmod{p}$，所以 $2(p+1)$ 是模 $p$ 的一个周期．所以，$\pi(p)\mid 2(p+1)$．
这就完成了证明．这一方法的局限性在于它高度依赖于斐波那契数列的通项公式，所以较难直接推广到一般的情形．
=== "利用扩域" 另一种证明方式则是试图直接计算矩阵 $A=\begin{pmatrix}1&1\1&0\end{pmatrix}$ 的阶．它的 特征多项式 是 $f(x) = x^2-x-1$，对应的判别式为 $\Delta = 5$．对于模 $p=5$，有 $\Delta\equiv 0\pmod{5}$，矩阵 $A$ 有两个相同特征值 $\lambda=3$，且不能对角化，需要单独计算．对于模 $p\equiv 1,4\pmod{5}$，由二次互反律可知，判别式 $\Delta=5$ 是模 $p$ 的二次剩余，矩阵 $A$ 在域 $\mathbf F_p$ 内有两个相异特征值 $\lambda_1\neq\lambda_2$，矩阵 $A$ 的阶就是 $\operatorname{lcm}(\operatorname{ord}(\lambda_1),\operatorname{ord}(\lambda_2))$，必然整除 $|\mathbf F_p^\times|=p-1$．对于模 $p\equiv 2,3\pmod{5}$，由二次互反律可知，判别式 $\Delta=5$ 是模 $p$ 的二次非剩余，矩阵 $A$ 在域 $\mathbf F_p$ 内没有特征值，而只有在 扩域 $\mathbf F_p[\sqrt{5}]$ 内才有两个相异特征值 $\lambda_1\neq\lambda_2$，由于 Frobenius 自同态 $x\mapsto x^p$ 将两根交换，有 $\lambda_2=\lambda_1^p$，故而 $\lambda_1^{p+1}=\lambda_2^{p+1}=\lambda_1\lambda_2=-1$，亦即 $\lambda_1^{2(p+1)}=\lambda_2^{2(p+1)}=1$，由此，矩阵 $A$ 的阶就是 $\operatorname{lcm}(\operatorname{ord}(\lambda_1),\operatorname{ord}(\lambda_2))$，必然整除 $2(p+1)$．这就得到了与前种方法一致的结论．
综上，对于不同的情形，相应地有：
所以，利用结论 1，对于一般的模数 $m=\prod_i p_i^{e_i}$，有
$$ \begin{aligned} \pi(m)&=\operatorname{lcm}{\pi(p_i^{e_i}):p_i\in\mathbf P} \ &\le \operatorname{lcm}{\pi(p_i^{e_i}):p_i=2\text{ or }p_i\equiv\pm1~(\operatorname{mod}{10})}\ &\quad \cdot 4\cdot\operatorname{lcm}{\pi(p_i^{e_i})/4:p_i=5\text{ or }p_i\equiv\pm3~(\operatorname{mod}{10})}\ &\le \prod{\pi(p_i^{e_i}):p_i=2\text{ or }p_i\equiv\pm1~(\operatorname{mod}{10})}\ &\quad \cdot 4\cdot\prod{\pi(p_i^{e_i})/4:p_i=5\text{ or }p_i\equiv\pm3~(\operatorname{mod}{10})}\ &\le \dfrac{3}{2}\cdot\prod{p_i^{e_i}:p_i=2\text{ or }p_i\equiv\pm1~(\operatorname{mod}{10})}\ &\quad \cdot 4\cdot\prod{p_i^{e_i}:p_i=5\text{ or }p_i\equiv\pm3~(\operatorname{mod}{10})}\ &= 6m. \end{aligned} $$
这就说明了斐波那契数列模 $m$ 的 Pisano 周期总是不超过 $6m$，而且等号当且仅当在 $m=2\cdot 5^e$ 处取得．
本页面主要译自博文 Числа Фибоначчи 与其英文翻译版 Fibonacci Numbers．其中俄文版版权协议为 Public Domain + Leave a Link；英文版版权协议为 CC-BY-SA 4.0．内容有改动．
严格来说，它是矩阵 $A$ 在一般线性群 $GL_2(\mathbf Z_m)$ 中的阶． ↩