[OI direction] Understand second-order linear recursion

This article is mainly aimed at popularization/improving groups OIer and ACMer. Considering most OIer cases, this article defaults that readers only know matrix multiplication and do not understand the determinant of the matrix, the rank of the matrix, etc. This article uses C++ to write code examples.

What is second-order linear recursion

Second-order linear recursive sequences have a famous name in the OI world:Generalized Fibonacci sequence. It refers to the following sequence\(\{a_n\}\)：

in,\(p,q\in \mathbb{R}\), the first two terms of the sequence are\(a_0, a_1\)。

Since the recursive relationship of this sequence is linear and each item is related to the first two terms, it is called a second-order linear recursive sequence.

Fibonacci sequence and naive algorithm

The current text sequence\(\{a_n\}\)middle\(p=q=1,a_0=a_1=1\)At that time, this sequence is the famous Fibonacci sequence\(\{f_n\}\)：

Based on such a recursive relationship, we can write linear complexity\(\mathcal{O}(n)\)The naive algorithm:

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

For moreSimple situation, such an algorithm is acceptable.

We can easily generalize to more general situations, so let the general simple algorithm of second-order linear recursion be left to readers for practice.

Optimize recursiveness by matrix multiplication

Although it is simple to directly calculate based on recursive form, it is really too slow and we need to optimize it. If you haven't learned matrix multiplication, please move toOI Wiki(Remember to check the reverse of the square array by the way, and you will use it later). If you have learned matrix multiplication, you should be familiar with the following equation (if you are a little unfamiliar with this equation, please try to manually calculate the matrix multiplication on the left):

This equation is intended to represent the Fibonacci recursive relationship by matrix multiplication for further operation. in\(\begin{bmatrix}f_{n-1} & f_{n-2}\end{bmatrix}\)It is called a state matrix (in this case, it is also a vector),\(\begin{bmatrix} 1 & 1\\ 1 & 0\end{bmatrix}\)It is called a transfer matrix. Such transfers are linear transformations in linear algebra, where the transfer matrix is ​​called a linear operator. We can easily generalize to general second-order linear recursion:

From this we can obtain the general formula for second-order linear recursion represented by matrix multiplication:

The algorithm can be optimized to logarithmic complexity through the matrix fast power algorithm\(\mathcal{O}(\log n)\), this is already the best in general (as for whether it is really the best, I don’t know, and I won’t prove it).

The code is also easy to implement, take Fibonacci as an example:

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];
}

A more interesting optimization algorithm

The next part of the content will be a bit abstract. If you have learned abstract algebra, it would be better to read this part of the content.

We can find another recursive matrix representation:

If you can't see it at first glance, you can calculate it. To find this formula, I thought about it all night. Note that in this formula all the matrices are as follows:

Each matrix is ​​related to only two elements, so we can use a binary to represent such matrices and define multiplication on the binary:

In fact, in the realm of abstract algebra, this isomorphism of a specific type matrix:

Therefore, the original algebraic properties of the matrix are also valid for quadratic groups, such as the most important law of bonding, which is the basis for the use of fast powers. It is easy to come up with a general formula that uses this binaries representation:

Unit elements of binary multiplication are easily observed:

And it's not difficult to find\((1,0)\)The inverse element:

In this way, we can easily implement the reverse.

Theoretically, since this method will calculate two or three times less multiplication or addition than a matrix, the constant will be smaller, and there is actually no difference. Here is the code for the Fibonacci sequence:

struct phi_tuple : pair<int, int> {
     using pair<int, int>::pair;
     phi_tuple &operator*=(const phi_tuple &other) {
         // When calculating here, the common factor is extracted to reduce the multiplication once
         return *this = {first * ( + ) +
                             second * ,
                         first * + 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;
 }

Linear algebra understanding of general term formulas by feature root method

For second-order linear recursion, a well-known method of finding general term formulas is the feature root method. For second-order linear recursion:

Featured equations:

Two complex roots of this equation can be solved:

when\(\phi\neq\psi\)When (i.e.\(\Delta\neq0\)),have:

when\(\phi=\psi\)When (i.e.\(\Delta=0\)),have:

Next, just need to\(a_0,a_1\)Substitute in and out\(s,t\)Just do it.

Let's take the Fibonacci sequence as an example. The characteristic equations of the Fibonacci sequence are:

Solution:

then:

Will\(f_0=1,f_1=1\)Substitute it and you will get a very complicated piece of stuff, so you can get it for convenience\(f_0=0,f_1=1\)：

This way we can get\(f_0=0,f_1=1\)The Fibonacci formula is now:

Please continue to read the program implementation...

Similar diagonalization and matrix eigenvalues ​​and eigenvectors

For convenience, if there are no special instructions, the contents of this section shall be\(2\times 2\)Matrix as an example

There are many proofs in the above method, and here we consider it from the perspective of linear algebra. A very natural idea is to find a way to expand the general formula of matrix representation:

Obviously the key to the problem is to expand the square matrix power on the right. At first glance, it will be a little unthinkable, so consider it from a simpler diagonal matrix. A diagonal matrix refers to a square matrix shaped like the following:

The blank area indicates omitted\(0\). Observe the square and cube of the focus matrix:

Guess the power of the diagonal matrix satisfies:

It is easy to prove by mathematical induction. So we found a method to calculate the power of the diagonal matrix, and then consider generalizing to a general square matrix. The most direct idea is to convert the square matrix into a diagonal matrix. Consider the following formula:

Note:

in\(\boldsymbol{P}\)and\(\boldsymbol{P}^{-1}\)Appears in pairs, so:

Just find it next\(\boldsymbol{A}\)Corresponding\(\boldsymbol{P}\)and\(\boldsymbol{\Lambda}\)That is, this process is called similar diagonalization of the matrix.

Before looking for similar diagonal methods, let’s add some knowledge. You may already know that linear system of equations can be expressed equivalently by matrix multiplication:

Based on this, the system of linear equations can be regarded as\(\boldsymbol{Ax}=\boldsymbol{b}\)form. when\(\boldsymbol{b}=\boldsymbol{0}\)At that time, it is called a homogeneous linear system of equations, and each square matrix has a unique corresponding to a homogeneous linear system of equations.

Also, you may have heard of determinants:

phalanx\(\boldsymbol{A}\)The determinant can show some properties of a square matrix, such as the situation of the solutions of the homogeneous system of equations corresponding to the square matrix and whether it is reversible:

1. \(|\boldsymbol{A}| = 0\)At the time, the square matrix is ​​irreversible, there are infinitely many solutions, and any solutions are linearly correlated, which means that each solution can be obtained by multiplying the other solution by number.
2. \(|\boldsymbol{A}| \neq 0\)When the square matrix is ​​reversible, there is only zero solution

Due to limited space, no proof will be made here.

Similar diagonalization requires the use of the eigenvalues ​​of the matrix\(\lambda\)and feature vectors\(\boldsymbol{x}\), which satisfies:

The right side of the equal sign can be multiplied by the unit matrix\(\boldsymbol{E}\)To transform:

because\(\boldsymbol{x}\neq\boldsymbol{0}\)Therefore

This is actually about\(\lambda\)The polynomial equation, in fact, to the left of the equal sign is called the feature polynomial of the matrix. Since this article takes a second-order square matrix as an example, the feature polynomial is mono-quadratic, so there are two roots\(\lambda_1,\lambda_2\), both roots are eigenvalues ​​of the matrix. Will\(\lambda\)You can find the eigenvector\(\boldsymbol{x}\), each\(\lambda_i\)All correspond to infinite sets of solutions\(k\boldsymbol{x_i},k\in\mathbb{R}\)。

Going back to similar diagonalization, it can be proved that if a matrix can similar diagonalize, then there must be the following similar diagonalization scheme:

in\(\boldsymbol{x}_i\)is any feature value\(\lambda_i\)The corresponding non-zero eigenvector. in\(\boldsymbol{P}\)It is called a feature matrix. If and only\(\lambda_i\)The two are different, the matrix\(\boldsymbol{P}\)It is reversible, and it is not difficult for readers to prove it. Solving linear equation systems and inverse matrices can both be done by the Gaussian elimination method. Due to space limitations, I will not expand it here.

Next, we will still take the Fibonacci sequence as an example. by\(f_0=1\)First item. There are matrix formulas (where\(\odot\)Indicates unimportant elements, the same below):

Eigenpolynomials, eigenvalues, and eigenvectors of transfer matrix:

Similar diagonal matrix and feature matrix:

Inverse of feature matrix:

Similar diagonalization:

So we get the general formula:

Next is to write the code. For convenience, the following text will be\(f_0=0\)As the first item... (To be honest, I don't know why I have to go around here, but I'm too lazy to change it) We have the general formula:

If you just want the value of a certain item, then you can just calculate it in a floating point environment. In fact, we can find a formula for cutting constants in half:

in\(\mathrm{round}(x)\)Indicates rounding to an integer. This is due to\(\left|\tfrac{1}{\sqrt{5}}\left(\tfrac{1-\sqrt{5}}{2}\right)^{n}\right|<\tfrac{1}{2}\), so rounding can ignore this part of the impact. More radical, based on engineering experience (actually I am too lazy to analyze errors), because\(\tfrac{1}{\sqrt{5}}\approx 0.447,\tfrac{1-\sqrt{5}}{2}\approx 1.618\), the following formula can be perfectly used as general formula:

The code is also very easy to implement and is suitable forNo need to take a mold：

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));
}

But if your requirement is Fibonacci sequence modulo, because it involves\(\sqrt5\), that's more complicated. The first idea is to have a secondary residual, but\(5\)There is no quadratic residual under many common modulo (if you don’t know what quadratic residual is, you can simply understand it as a square root in the modulo sense, which is the content of number theory). If you have learned about abstract algebra or algebraic number theory, you will probably think of something - quadratic integer rings.

Quadratic integer ring

A quadratic integer ring is a set of all numbers composed of integers added to a special irrational number, and obtained by addition, subtraction and multiplication. These numbers not only retain the properties of integer rings, but also allow operations such as residual division in a wider range. If you have learned complex numbers, just like complex numbers are composed of real numbers and imaginary parts, a quadratic integer ring is composed of integers and a specific irrational number (e.g.\(\sqrt5\)) is composed of integer multiples of , which is usually expressed as\(\mathbb{Z}[\sqrt{d}]\). Quadratic integer ring\(\mathbb{Z}[\sqrt{d}]\)The operations on the same page are consistent with the quadratic root operation:

1. Addition and subtraction:\((a_1+a_2\sqrt{d})+(b_1+b_2\sqrt{d}) = (a_1+b_1)+(a_2+b_2)\sqrt{d}\)
2. multiplication:\((a_1+a_2\sqrt{d})(b_1+b_2\sqrt{d})=(a_1b_1+a_2b_2d)+(a_1b_2+a_2b_1)\sqrt{d}\)

Since division is more complicated and we can't use it now, we won't discuss it. In the modulus sense, each time the calculation is completed, you can modulate the coefficients of the integer part and the root part, which actually constitute the remaining ring of the quadratic integer ring, which is generally recorded as\(\mathbb{Z}[\sqrt{d}]/\{p\}\)，\(p\)It is modulus. In the modulo sense, the quadratic integer divided by a pure integer is still converted into a number theory reciprocal (inverse) multiplied by the integer. The situation of dividing by quadratic integers is too complicated and cannot be used in this article for the time being, so I won't discuss it.

Next is the pleasant code implementation link:

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 * ) % 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 + , s + );
    }
    constexpr sr5m64 operator-(const sr5m64 &o) const {
        return *this + (-o);
    }
    constexpr sr5m64 operator*(const sr5m64 &o) const {
        return sr5m64(x *  + 5 * s * , x *  + s * );
    }
    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;
}

Some properties of second-order linear recursion

Periodicity in the model sense

Any second-order linear recursive in the mode\(m\)All of them are periodic. Consider the state modulus in the matrix representation\(m\)

Obviously, each item in the state has only at most\(m\)Therefore, all possible states have at most\(m^2\)Therefore, after the second-order linear recursion is sufficiently sufficiently repeated, two identical states will inevitably occur, and therefore periodicity will appear. This conclusion can also be generalized to arbitrary order linear recursion.

Since binomial theorems are required to solve the minimum positive cycle, we will not expand it here. Please refer to the details for details.Pisano Cycle. OneA widely circulated topicThis part of knowledge may be required.

The ratio of two adjacent terms of the Fibonacci sequence tends to be golden ratio

The golden ratio refers to

It is obviously a characteristic root of the Fibonacci sequence. The subtitle of this section means:

Just substitute the formula directly to force the limit:

Similar conclusions can be drawn for a part of the second-order linear recursion.

Conclusion

If you are tired, just write this one first...