[OI 向] 深入理解二阶线性递推

本文主要面向普及/提高组 OIer 和 ACMer。考虑大多数 OIer 的情况，本文默认读者只会矩阵乘法，不了解矩阵的行列式，矩阵的秩等内容。本文使用 C++ 编写代码示例。

什么是二阶线性递推

二阶线性递推数列在 OI 界还有个著名的名字：广义斐波那契数列。其所指为如下数列 ：

其中，，数列前两项为 。

由于该数列的递推关系是线性的，且每一项都和前两项有关，因此称为二阶线性递推数列。

斐波那契数列与朴素算法

当前文的数列 中 时，该数列就是大名鼎鼎的斐波那契数列 ：

基于这样的递推关系，我们可以写出线性复杂度 的朴素算法：

int f(int n) {
    if (n == 0 || n == 1) {
        return 1;
    }
    return f(n - 1) + f(n - 2);
}

对于较为简单的情况，这样的算法是可以接受的。

我们可以很容易的推广到更一般的情况，所以对于一般的二阶线性递推的朴素算法就留给读者练习吧。

通过矩阵乘法来优化递推

直接根据递推式计算虽然简单，但是实在是太慢了，我们需要优化。如果你还没有学过矩阵乘法，请移步至 OI Wiki（记得顺便把方阵的逆看看，后面要用）。如果你学过矩阵乘法，应当对下面的式子不陌生（如果你对这个式子有点陌生的话，请尝试手动计算左面的矩阵乘法）：

这个式子旨在将斐波那契递推关系用矩阵乘法表示，以进行进一步的操作。其中 被称为状态矩阵（在本例中，它也是个向量）， 被称为转移矩阵。这样的转移在线性代数中属于线性变换，其中转移矩阵被称为线性算子。我们可以很容易的推广到一般的二阶线性递推上：

据此可以得出使用矩阵乘法表示的二阶线性递推的通项公式：

通过矩阵快速幂算法就可以将算法优化到对数复杂度 ，这在一般情况下已经是最优的了（至于是不是真的是最优的，我也不知道，我也不会证）。

代码也是很好实现的，以斐波那契为例：

using vec = array<int, 2>;
struct matrix : public array<vec, 2> {
    using array<vec, 2>::array;
    matrix(const vec &a, const vec &b) : array{a, b} {
        (*this)[0] = a;
        (*this)[1] = b;
    }
    matrix &operator*=(const matrix &other) {
        matrix res{};
        for (int i{0}; i < 2; ++i)
            for (int j{0}; j < 2; ++j)
                for (int k{0}; k < 2; ++k)
                    res[i][j] += (*this)[i][k] * other[k][j];
        return *this = res;
    }
} e{{1, 0}, {0, 1}};

matrix qpow(matrix x, int n) {
    auto res{e};
    while (n) {
        if (n & 1)
            res *= x;
        x *= x, n >>= 1;
    }
    return res;
}

int fib(int n) {
    return qpow(matrix{{1, 1}, {1, 0}}, n)[0][0];
}

一种更有趣的优化算法

接下来的这部分内容会有点抽象，如果你学过抽象代数那看这部分内容就再好不过了。

我们可以找到另一种递推的矩阵表示：

如果你乍一眼看不出来，可以计算一下。为了找到这个式子，我思考了一个晚上。注意到在这个式子中所有的矩阵都形如：

每个矩阵都只与两个元素有关，因此我们可以使用一个二元组来表示这类矩阵，并定义二元组上的乘法：

事实上，在抽象代数领域，这是找到了一个特定类型矩阵的同构：

因此矩阵原有的代数性质对于二元组也是成立的，比如最重要的结合律，这是快速幂得以使用的基础。很容易得出使用这种二元组表示的通项公式：

很容易观察到二元组乘法的单位元：

而且不难找到 的逆元：

这样我们就可以方便的实现倒推了。

理论上讲由于这个方法相比矩阵会少算两三次乘法或加法，因此常数会小一点，实际上没什么差距。下面是针对斐波那契数列的代码：

struct phi_tuple : pair<int, int> {
    using pair<int, int>::pair;
    phi_tuple &operator*=(const phi_tuple &other) {
        // 这里计算的时候提取了公因式以减少一次乘法
        return *this = {first * (other.first + other.second) +
                            second * other.first,
                        first * other.first + second * other.second};
    }
} e{0, 1};

phi_tuple qpow(phi_tuple x, int n) {
    auto res{e};
    while (n) {
        if (n & 1)
            res *= x;
        x *= x, n >>= 1;
    }
    return res;
}

int fib(int n) {
    return qpow(phi_tuple{1, 0}, n + 1).first;
}

特征根法求通项公式的线性代数理解

对于二阶线性递推，一种广为人知的求通项公式的方法是特征根法。对于二阶线性递推：

有特征方程：

可以解出该方程的两个复根：

当 时（即 ），有：

当 时（即 ），有：

接下来只需要将 代入求出 即可。

---

下面我们以斐波那契数列为例。斐波那契数列的特征方程为：

解得：

于是：

将 代入，会得到一坨很复杂的东西，于是方便起见取 ：

这样我们就可以得到取 时的斐波那契通项公式了：

程序实现请继续往后看……

相似对角化和矩阵特征值与特征向量

方便起见，如无特殊说明，本节内容均以 矩阵为例

上面的方法有多种证明，在这里我们从线性代数的角度考虑。一种很自然的想法是想办法展开矩阵表示的通项公式：

显然问题的关键在于展开右侧的方阵幂。乍一看会有点没思路，不妨从更简单的对角矩阵考虑。对角矩阵指的是形如下面的方阵：

其中空白的区域表示省略的 。观察对焦方阵的平方以及立方：

猜测对角矩阵的幂满足：

很容易通过数学归纳法证明。于是我们找到了计算对角矩阵的幂次的方法，接下来考虑推广到一般的方阵。最直接的想法是将方阵转化为对角矩阵。考虑下面的式子：

注意到：

其中 和 成对出现，于是：

接下来只需找到 对应的 和 即可，这个过程被称为矩阵的相似对角化。

在寻找相似对角化的方法之前，我们先补充一些知识。你也许已经知道线性方程组可以用矩阵乘法等价表示：

基于此，线性方程组均可视为 的形式。当 时，称其为齐次线性方程组，每个方阵都唯一对应一个齐次线性方程组。

另外，你还可能听说过行列式：

方阵 的行列式可以显示方阵的一些性质，比如方阵对应的齐次方程组的解的情况和是否可逆：

1. 时，方阵不可逆，有无穷多解，且任意解之间都线性相关，也就是说每个解都可以通过数乘另一个解得到
2. 时，方阵可逆，只有零解

因篇幅有限，在此不做证明。

相似对角化需要利用矩阵的特征值 和特征向量 ，其满足：

等号右侧可以乘上单位矩阵 来变形：

由于 故

这其实就是一个关于 的多项式方程，实际上等号左边被称为矩阵的特征多项式。由于本文以二阶方阵为例，因此特征多项式是一元二次的，那么就有两个根 ，这两个根都是矩阵的特征值。将 回代即可得求出特征向量 ，每个 均对应无穷多组解 。

回到相似对角化，可以证明，如果一个矩阵是可以相似对角化的，那么必然有下面的相似对角化方案：

其中 是任意一个特征值 对应的非零特征向量。其中 被称为特征矩阵。当且仅当 两两不同时，矩阵 可逆，读者自证不难。解线性方程组和求逆矩阵都可以通过高斯消元法，限于篇幅，这里就不展开了。

---

接下来仍然以斐波那契数列为例。以 为首项。有矩阵通项公式（其中 表示不重要的元素，下同）：

转移矩阵的特征多项式、特征值和特征向量：

相似对角矩阵和特征矩阵：

特征矩阵的逆：

相似对角化求幂次：

于是我们就得到了通项公式：

---

接下来就是写代码了。方便起见，后文就以 作为首项……（说实话我也不知道我为啥要在这里绕来绕去，但是懒得改了）我们有通项公式：

如果只是要某一项的值的话，那直接在浮点数环境下计算就可以了。事实上，我们可以发现一个常数对半砍的公式：

其中 表示四舍五入到整数。这是由于 ，因此四舍五入可以忽略这部分影响。更加激进的，根据工程经验（其实是懒得误差分析），由于 ，下面的公式可以完美作为通项公式：

代码也很好实现，适用于不需要取模的情况：

float qpow(float x, int n) {
    float r = 1;
    while (n) {
        if (n & 1)
            r *= x;
        x *= x, n >>= 1;
    }
    return r;
}

int fib(int n) {
    return round(0.447 * qpow(1.618, n));
}

但是如果你的需求是斐波那契数列取模，因为涉及到 ，那就比较复杂了。第一想法是二次剩余，但是 在很多常见模数下没有二次剩余（如果你不知道什么是二次剩余，可以简单的理解为模意义下的平方根，这是数论的内容）。如果你了解过抽象代数或者代数数论，你肯能会想到一个东西——二次整数环。

二次整数环

二次整数环就是所有由整数加上一个特殊无理数后，通过加减乘得到的数所组成的集合。这些数不仅保留了整数环的性质，还可以在更广的范围内进行带余除法等操作。如果你学过复数，那就像复数是由实数和虚数部分组成一样，二次整数环是由整数和一个特定的无理数（如 ）的整数倍组合而成的数，这通常表示为 。二次整数环 上的运算和二次根式运算是一致的：

1. 加减法：
2. 乘法：

由于除法较为复杂而且我们现在用不到，就不做讨论了。在模意义下，每次计算完分别对整数部分和根式部分的系数取模即可，这实际上构成了二次整数环的剩余类环，一般记作 ， 是模数。在模意义下二次整数除以一个纯整数依旧是转化为乘上该整数的数论倒数（逆元）。除以二次整数的情况过于复杂且暂时本文用不到，就不做讨论了。

接下来就是愉快的代码实现环节：

using i64 = int64_t;

constexpr i64 m = 1e9 + 7;

template <typename num> constexpr num qpow(num a, i64 n) {
    num r = 1;
    while (n) {
        if (n & 1)
            r *= a;
        a *= a, n >>= 1;
    }
    return r;
}

struct m64 {
    i64 x;
    constexpr m64(const i64 x) : x(((x % m) + m) % m) {
    }
    constexpr operator i64() const {
        return x;
    }
    constexpr m64 operator-() const {
        return m64(-x);
    }
    constexpr m64 operator*(const m64 &o) const {
        return m64((x * o.x) % m);
    }
    constexpr m64 &operator*=(const m64 &o) {
        *this = *this * o;
        return *this;
    }
    constexpr m64 operator/(const m64 &o) const {
        return m64(*this * qpow(o, m - 2));
    }
};

struct sr5m64 {
    m64 x, s;
    constexpr sr5m64(m64 x, m64 s) : x(x), s(s) {
    }
    constexpr sr5m64(i64 x) : x(x), s(0) {
    }
    constexpr sr5m64 operator-() const {
        return sr5m64(-x, -s);
    }
    constexpr sr5m64 operator+(const sr5m64 &o) const {
        return sr5m64(x + o.x, s + o.s);
    }
    constexpr sr5m64 operator-(const sr5m64 &o) const {
        return *this + (-o);
    }
    constexpr sr5m64 operator*(const sr5m64 &o) const {
        return sr5m64(x * o.x + 5 * s * o.s, x * o.s + s * o.x);
    }
    constexpr sr5m64 &operator*=(const sr5m64 &o) {
        *this = *this * o;
        return *this;
    }
    constexpr sr5m64 operator/(const m64 &o) const {
        return sr5m64(x / o, s / o);
    }
};

constexpr m64 fib(i64 n) {
    return ((qpow(sr5m64(1, 1) / 2, n) - qpow(sr5m64(1, -1) / 2, n)) * sr5m64(0, 1) / 5).x;
}

二阶线性递推的一些性质

模意义下的周期性

任何的二阶线性递推在模 下都是具有周期性的。考虑矩阵表示中的状态模

显然状态的每一项最多只有 种，因此所有可能的状态最多有 种，所以在二阶线性递推充分多次后必然产生两个相同的状态，因此会呈现出周期性。这个结论也可以推广到任意阶线性递推。

由于求解最小正周期需要用到二项式定理等内容，在此不做展开，详情参考皮萨诺周期。一道广为流传的题目可能需要这部分知识。

斐波那契数列相邻两项之比趋于黄金分割比

黄金分割比指的是

很显然其是斐波那契数列的一个特征根。本节小标题的意思是：

直接代入通项公式暴力求极限即可：

对于一部分二阶线性递推，也能得出相似的结论。

结语

累了，就先写到这里吧…