
## Catalan Numbers

Catalan numbers $H_n$ can be applied to the following problems:

1. There are $2n$ people queuing up to enter a theater. The entrance fee is 5 yuan. Among them, only $n$ people have a 5 yuan bill, while the other $n$ people only have a 10 yuan bill. The theater has no other bills. How many ways are there to ensure that whenever a person with a 10 yuan bill buys a ticket, the ticket office has a 5 yuan bill for change?
2. There is a grid of size $n\times n$ with the lower left corner at $(0, 0)$ and the upper right corner at $(n, n)$. Starting from the lower left corner, you can only move one unit to the right or up each time, without going above the diagonal line $y=x$. How many possible paths are there to reach the upper right corner without crossing the line $y=x$?
3. Choose $2n$ points on a circle and connect these points pairwise to obtain $n$ line segments. How many ways are there to ensure that the $n$ line segments obtained do not intersect?
4. Without intersecting diagonals, how many ways are there to divide a convex polygon area into triangular regions?
5. If the stack (infinite) has a push sequence of $1,2,3, \cdots ,n$, how many different pop sequences are there?
6. How many different binary trees can be constructed with $n$ nodes?
7. Given $2n$ numbers consisting of $n$ $+1$s and $n$ $-1$s, denoted as $a_1,a_2, \cdots ,a_{2n}$, such that the partial sum satisfies $a_1+a_2+ \cdots +a_k \geq 0~(k=1,2,3, \cdots ,2n)$, how many sequences satisfy this condition?

The corresponding sequence is:

| $H_0$ | $H_1$ | $H_2$ | $H_3$ | $H_4$ | $H_5$ | $H_6$ | ... |
| :---: | :---: | :---: | :---: | :---: | :---: | :---: | :-: |
|   1   |   1   |   2   |   5   |  14   |  42   |  132  | ... |

## Recurrence Relation

The solution to this recurrence relation is:

$$
H_n = \frac{\binom{2n}{n}}{n+1}(n \geq 2, n \in \mathbf{N_{+}})
$$

Common formulas for Catalan numbers are:

$$
H_n = \begin{cases}
    \sum_{i=1}^{n} H_{i-1} H_{n-i} & n \geq 2, n \in \mathbf{N_{+}}\\
    1 & n = 0, 1
\end{cases}
$$

$$
H_n = \frac{H_{n-1} (4n-2)}{n+1}
$$

$$
H_n = \binom{2n}{n} - \binom{2n}{n-1}
$$

> A common application: for a stack, the push sequence is $1,2, \ldots ,n$, find the total number of possible pop sequences.

```cpp
#include <iostream>
using namespace std;
int n;
long long f[25];

int main() {
  f[0] = 1;
  cin >> n;
  for (int i = 1; i <= n; i++) f[i] = f[i - 1] * (4 * i - 2) / (i + 1);
  // Here, formula 2 is used
  cout << f[n] << endl;
  return 0;
}
```

## Closed Form

The recurrence relation of Catalan numbers is

$$
H_n=\sum_{i=0}^{n-1}H_{i}H_{n-i-1} \quad (n\ge 2)
$$

where $H_0=1$ and $H_1=1$. Let its ordinary generating function be $H(x)$.

We find that the recurrence relation of Catalan numbers is similar to the form of convolution. Therefore, we construct an equation about $H(x)$ using convolution:

$$
\begin{aligned}
H(x)&=\sum_{n\ge 0}H_nx^n\\
&=1+\sum_{n\ge 1}\sum_{i=0}^{n-1}H_ix^iH_{n-i-1}x^{n-i-1}x\\
&=1+x\sum_{i\ge 0}H_{i}x^i\sum_{n\ge 0}H_{n}x^{n}\\
&=1+xH^2(x)
\end{aligned}
$$

Solving this equation, we get

$$
H(x)=\frac{1\pm \sqrt{1-4x}}{2x}
$$

Now, we have a question: which root should we take? Rationalizing the numerator, we get

$$
H(x)=\frac{2}{1\mp \sqrt{1-4x}}
$$

Substituting $x=0$, we obtain the constant term of $H(x)$, which is $H_0$. When $H(x)=\dfrac{2}{1+\sqrt{1-4x}}$, we have $H(0)=1$, satisfying the requirement. The other root leads to the denominator being $0$ (non-convergence), so it is discarded.

Therefore, we obtain the closed form of the generating function of Catalan numbers:

$$
H(x)=\frac{1- \sqrt{1-4x}}{2x}
$$

Next, we need to expand it. But note that its denominator is not in the form of a polynomial like the Fibonacci sequence, so it is not convenient to use the expansion form of geometric sequences. Here we need to use Newton's binomial theorem. Let's first expand $\sqrt{1-4x}$:

$$
\begin{aligned}
(1-4x)^{\frac{1}{2}}
&=\sum_{n\ge 0}\binom{\frac{1}{2}}{n}(-4x)^n\\
&=1+\sum_{n\ge 1}\frac{\left(\frac{1}{2}\right)^{\underline{n}}}{n!}(-4x)^n
\end{aligned} \tag{1}
$$

Note that

$$
\begin{aligned}
\left(\frac{1}{2}\right)^{\underline{n}}
&=\frac{1}{2}\frac{-1}{2}\frac{-3}{2}\cdots\frac{-(2n-3)}{2}\\
&=\frac{(-1)^{n-1}(2n-3

)!!}{2^n}\\
&=\frac{(-1)^{n-1}(2n-2)!}{2^n(2n-2)!!}\\
&=\frac{(-1)^{n-1}(2n-2)!}{2^{2n-1}(n-1)!}
\end{aligned}
$$

Here, we use the simplification technique of double factorials. Then, substitute into $(1)$, we get

$$
\begin{aligned}
(1-4x)^{\frac{1}{2}}
&=1+\sum_{n\ge 1}\frac{(-1)^{n-1}(2n-2)!}{2^{2n-1}(n-1)!n!}(-4x)^n\\
&=1-\sum_{n\ge 1}\frac{(2n-2)!}{(n-1)!n!}2x^n\\
&=1-\sum_{n\ge 1}\binom{2n-1}{n}\frac{1}{(2n-1)}2x^n
\end{aligned}
$$

Substitute back into the original expression, we get

$$
\begin{aligned}
H(x)&=\frac{1- \sqrt{1-4x}}{2x}\\
&=\frac{1}{2x}\sum_{n\ge 1}\binom{2n-1}{n}\frac{1}{(2n-1)}2x^n\\
&=\sum_{n\ge 1}\binom{2n-1}{n}\frac{1}{(2n-1)}x^{n-1}\\
&=\sum_{n\ge 0}\binom{2n+1}{n+1}\frac{1}{(2n+1)}x^{n}\\
&=\sum_{n\ge 0}\binom{2n}{n}\frac{1}{n+1}x^{n}\\
\end{aligned}
$$

Thus, we obtain the general formula for Catalan numbers.

## Path Counting Problem

Non-decreasing paths refer to paths that can only move upwards or to the right.

1. The number of non-decreasing paths from $(0,0)$ to $(m,n)$ is equal to the number of permutations of $m$ $x$ and $n$ $y$, i.e., $\dbinom{n + m}{m}$.

2. The number of non-decreasing paths from $(0,0)$ to $(n,n)$ that do not touch the line $y=x$ except at endpoints:

   First, consider the paths below the line $y=x$, which all start from $(0, 0)$, pass through $(1, 0)$ and $(n, n-1)$, and reach $(n,n)$. This can be seen as the number of non-decreasing paths from $(1,0)$ to $(n,n-1)$ without touching $y=x$.

   All non-decreasing paths total $\dbinom{2n-2}{n-1}$. For any path that touches $y=x$, we can symmetrically transform the part from the point where it leaves this line to $(1,0)$, obtaining a non-decreasing path from $(0,1)$ to $(n,n-1)$. The reverse is also true. Hence, the number of non-decreasing paths below $y=x$ is $\dbinom{2n-2}{n-1} - \dbinom{2n-2}{n}$. Due to symmetry, the answer is $2\dbinom{2n-2}{n-1} - 2\dbinom{2n-2}{n}$.

3. The number of non-decreasing paths from $(0,0)$ to $(n,n)$ that do not cross the line $y=x$: Similarly, we can obtain $\dfrac{2}{n+1}\dbinom{2n}{n}$ using similar methods.

> NOTE: This article is rewritten by [OI WIKI](https://oi-wiki.org/math/combinatorics/catalan/).


The article discusses Catalan numbers, a sequence used in combinatorial problems. It covers their recurrence relation, closed form, and applications in path counting. It provides insights into problem-solving using Catalan numbers.

The Catalan Series


When we navigate through the ocean of divide and conquer, there's a question that must be addressed—how is the time complexity of divide and conquer algorithms calculated?

## From Divide and Conquer to Recurrence Relation

Divide and conquer involves breaking down a large problem into smaller subproblems and then merging the solutions of these subproblems to obtain the solution to the original problem. Recurrence relation is a mathematical way of modeling the time cost of this process.

Starting with a concrete example makes it easier to understand. Taking merge sort as an example, let $T(n)$ represent the time complexity of the original problem, where $n$ is the size of the input data. In the first split, the original data is divided into two equal parts, each of size $n/2$. Then, recursive calls to merge sort are made on these two parts, and the time complexity of this part is $2 T(n/2)$. Finally, the two sorted subarrays are merged, requiring a time complexity of $\Theta(n)$.

Based on the analysis above, we can write the recurrence relation for merge sort as follows:

$$
T(n) = \left\{ \begin{aligned} &\Theta(1) & \text{if } n=1\\ &2 T(n/2) + \Theta(n) & \text{if } n>1 \end{aligned} \right.
$$

It can be imagined that for any divide and conquer algorithm, its time complexity can be expressed as a recurrence relation. The right-hand side of the recurrence relation usually consists of two parts: the total cost of the divided subproblems and the total cost of dividing and merging subproblems. The first part remains a function of $T$, while the second part represents the actual asymptotic complexity.

Having a recurrence relation is not sufficient to determine the exact time complexity of the algorithm. We need to further calculate the actual solution $T(n)$.

## Substitution Method

The first method for solving recurrence relations is called the substitution method. We need to guess the form of the solution first and then use mathematical induction to prove that the guess is correct.

For Equation (1), let's assume that our guess for the solution is $T(n) = O(n \log n)$. If we remove the asymptotic notation, this equation is equivalent to $T(n) \le c_1 n \log n$, where $c_1$ is any positive constant.

According to mathematical induction, we first need to establish that the guess holds for small values of $n$. When $n=1$, according to Equation (1), $T(1) = \Theta(1) = d_1$, where $d_1$ is some constant greater than 0. According to the guess, we hope $T(1) \le c_1 \log 1 = 0$, but no matter how we choose $c_1$, this equation cannot hold because $T(1) = d_1$ is always greater than 0. The mathematical induction fails before it even starts.

However, there is no need to worry; this is just a special case caused by $ \log 1 = 0$. We can place the initial state of mathematical induction at $n>1$ without affecting the results of mathematical induction. We are only interested in the asymptotic properties of $T(n)$ when $n$ is sufficiently large, not in the initial stage.

Now, let $n = 2$, then $T(2) = 2T(1) + d_2 = 2d_1 + d_2$. We hope that $T(2) \le c_1 \log 2$ holds. Simplifying, we get $c_1 \ge d_1 + d_2/2$. Since $d_1$ and $d_2$ are constants, such $c_1$ exists, and the initial condition holds.

Next, mathematical induction assumes that the solution holds at $n/2$, i.e., $T(n/2) \le c_1 \frac{n}{2} \log\frac{n}{2} = \frac{1}{2}c_1n(\log n - 1)$. Then we prove that $T(n) \le c_1 n \log n$ also holds.

Substituting $T(n/2)$ into Equation (1), we get

$$
\begin{equation} \begin{aligned} T(n) & = 2 T\left(\frac{n}{2}\right) + \Theta(n) \\ & \le 2 \frac{1}{2}c_1n(\log n - 1) + d_2n \\ & = c_1 n \log n + (d_2 - c_1) n \\ & \le c_1 n \log n \end{aligned} \end{equation}
$$

Where the last line's inequality requires $d_2 - c_1 < 0$, obviously such $c_1$ exists, so the proof is complete.

Thus, we have proved that $T(n) = O(n \log n)$ is the solution to the recurrence relation (1). However, this is not the end; this asymptotic notation is not a tight bound. Can we further guess that $T(n) = \Theta(n \log n)$ holds? This is indeed a good idea. So, as we proceed, as long as we can prove that $T(n) = \Omega(n \log n)$ holds, we can obtain the tight bound mentioned above. Removing the asymptotic notation, this equation is equivalent to $T(n) \ge c_2 n \log n$. Using mathematical induction again, confirming that this holds in the initial stage, and then assuming that the solution holds at $n/2$, we can derive that the solution holds at $n$. Therefore, the proof is complete. Since the proof steps are almost identical to those just mentioned, we will not elaborate on them, and you can try it yourself.

The substitution method relies on accurate guesses, which can be impractical. Therefore, in practice, a more common approach is to find a guess using other methods and then use the substitution method to verify the guess. The recursive tree method we will introduce next provides this possibility.

## Recursive Tree

Expanding the recurrence relation layer by layer directly yields a tree structure. For example, let's take $T(n) = 3 T(n/4) + \Theta(n^2)$ as an example and see how to draw its recursive tree.

Expanding the first layer, we draw three child nodes under the root node.

Here, $cn^2$ represents the cost of decomposition and merging at this layer, which corresponds to the part represented by $\Theta(n^2)$ in the recurrence relation.

Next, further expand to the second layer.

![Recursive Tree](https://image-assets.dreams.plus/20231218224727.webp)

Here, the cost of each node in the first layer is still $c(\frac{n}{4})^2$, still corresponding to the part represented by $\Theta(n^2)$ in the recurrence relation, but the data size has changed from $n$ to $n/4$.

Continue expanding in this way until the leaf nodes.

![Recursive Tree](https://image-assets.dreams.plus/20231218224626.webp)

This completes a complete recursive tree. The data size of each layer is reduced to one-fourth of the previous one, until the data size of the leaf node is 1, so the total number of layers is $\log_4 n$. The rightmost side of the figure calculates the total cost of all nodes in each layer, and the calculation is simple: multiply the cost of a single node by the number of nodes in that layer. In the last layer, the cost of a single node is a constant, and the number of nodes is $3^{\log_4 n} = n^{\log_4 3}$, and we denote the total cost as $\Theta(n^{\log_4 3})$.

Next, we just need to add up the costs of all nodes to obtain the solution to the recurrence relation.

$$
\begin{aligned} T(n) &= cn^2 + \frac{3}{16} cn^2 + \left(\frac{3}{16}\right)^2 cn^2 + \cdots + \left( \frac{3}{16} \right)^{\log_4 n-1} cn^2 + \Theta(n^{\log_4 3}) \\ &=\sum_{i=0}^{\log_4 n-1}{\left(\frac{3}{16}\right)^i cn^2} + \Theta(n^{\log_4 3})\\ &\lt \sum_{i=0}^{\infty}{\left(\frac{3}{16}\right)^i cn^2} + \Theta(n^{\log_4 3}) \\ &= \frac{1}{1 - 3/16}cn^2 + \Theta(n^{\log_4 3}) \\ &= \frac{16}{13} cn^2 + \Theta(n^{\log_4 3}) \\ &= O(n^2) \end{aligned}
$$

Here, in the third line, we made a simplification by replacing a finite geometric series with an infinite geometric series. Although the total size has increased, this infinite decreasing geometric series has a limit, so this replacement can simplify the derivation process. The second term in the fifth line is a lower-order term relative to the first term, so it is discarded, resulting in the final result.

Similar to before, we have only proven the upper bound of the recurrence relation. Naturally, we want to try whether it can also be proven as a lower bound. In the above equation, each term in the summation contains $cn^2$, and each term is positive, so $T(n) = \Omega(n^2)$ obviously holds. Therefore, the solution to the recurrence relation is $T(n) = \Theta(n^2)$.

In this example, we directly calculated the solution to the recurrence relation using the recursive tree method. Of course, we can use the substitution method to verify whether this solution is correct. In some cases, it may be slightly difficult to calculate the exact solution using the recursive tree method, but it does not prevent us from making a reasonable guess and then further verifying this guess using the substitution method.

## Master Theorem

Both the substitution method and the recursion tree method have their drawbacks. The former requires an accurate guess, while the latter involves drawing a recursion tree and deriving results. If you derive too much, you may find that the forms of most recursive solutions seem to be the same. From the derivation process above, it can be observed that the final solution depends entirely on the sum of costs in the middle layers and has nothing to do with the leaf layer because the latter is a lower-order term. In other recursion trees, the final solution may depend entirely on the cost of the leaf layer, unrelated to the middle layers. Therefore, we may be able to classify recursive formulas and directly write out the final solution according to some rules. This is the idea behind the master theorem.

When the divide-and-conquer method evenly splits the data each time, the recursive formula generally has the following form.

$$
T(n) = aT\left(\frac{n}{b}\right) + f(n) \tag{2}
$$

where $a$ is a positive integer representing the number of subproblems when dividing each time, $b$ is an integer greater than 1 representing the reduction factor in problem size each time, and $f(n)$ is an asymptotically positive function representing the cost of dividing and merging.

To find a unified solution applicable to any $a$, $b$, and $f(n)$, we still need to use the recursion tree to calculate the costs of all intermediate and leaf nodes and try to derive a simple result.

![](https://image-assets.dreams.plus/20231218224657.webp)

From the graph, we can see that the total cost of all nodes is $T(n) = \Theta(n^{\log_b a}) + \sum_{j=0}^{\log_b n-1}{a^j f(n/b^j)}$. Further derivation requires a breakthrough idea, and I don't know how the original researchers discovered the following pattern, but it is indeed magical.

Researchers found that the relationship between $f(n)$ and $\Theta(n^{\log_b a})$ determines the final simplified form of $\sum_{j=0}^{\log_b n-1}{a^j f(n/b^j)}$. Define

$$
g(n) = \sum_{j=0}^{\log_b n-1}{a^j f(n/b^j)} \tag{3}
$$

Now let's analyze three cases.

1. If $\Theta(n^{\log_b a})$ is larger, we can assume that there exists a constant $\epsilon > 0$ such that $f(n) = O(n^{\log_b a-\epsilon})$. Substituting this into the summation formula, it can be simplified to $g(n) = O(n^{\log_b a})$.

2. If both are equally large, directly substitute into equation (3) and simplify to get $g(n) = \Theta(n^{\log_b a} \log n)$.

3. If $f(n)$ is larger and, for some constant $c < 1$, for all sufficiently large $n$, $a f(n/b) \le c f(n)$ holds, then it can be deduced that $g(n) = \Theta(f(n))$.

Here, the complete proof process is not provided, partly because I don't want to make it too long, and secondly, I feel that these proofs are somewhat cumbersome and not very meaningful. Interested students can refer to the original book.

Substituting the functions $g(n)$ for these three cases into the complete solution, we get:

- For case 1: $T(n) = \Theta(n^{\log_b a}) + O(n^{\log_b a}) = \Theta(n^{\log_b a})$

- For case 2: $T(n) = \Theta(n^{\log_b a}) + \Theta(n^{\log_b a} \log n) = \Theta(n^{\log_b a} \log n)$

- For case 3: $T(n) = \Theta(n^{\log_b a}) + \Theta(f(n)) = \Theta(f(n))$

Now, for any recursive formula, we only need to determine which case it belongs to and directly substitute into the corresponding formula. To facilitate understanding, let's go through a few examples.

**Example 1:**

$T(n) = 9 T(n/3) + n$

For this recursive formula, we have $a = 9$, $b = 3$, $f(n) = n$. Since $f(n) = n$ is asymptotically smaller than $\Theta(n^2)$, it can be applied to the first case of the master theorem, and the solution is $T(n) = \Theta(n^2)$.

**Example 2:**

$T(n) = T(2n/3) + 1$

For this formula, $a = 1$, $b = 3/2$, $f(n) = 1$. Since $n^{\log_b a} = n^{\log_{3/2} 1} = n^0 = 1$, which is exactly equal to $f(n)$, it can be applied to the second case of the master theorem, and the solution is $T(n) = \Theta(\log n)$.

**Example 3:**

$T(n) = 3T(n/4) + n\log n$

For this formula, $a = 3$, $b = 4$, $f(n) = n\log n$. Since $n^{\log_b a} = n^{\log_4 3} = n^{0.793}$, and $f(n) = n\log n$ is asymptotically larger than $n^{0.793}$, it should consider the third case of the master theorem. However, we need to check if the additional condition holds: for some constant $c < 1$ and all sufficiently large $n$, $a f(n/b) \le c f(n)$ holds. Substituting $a$, $b$, and $f(n)$ into the inequality, we get

$$
3\frac{n}{4}\log\frac{n}{4} \le cn\log n \implies \frac{3}{4}( \log n - 2) \le c\log n \implies \left(\frac{3}{4} - c\right) \log n \le \frac{3}{2}
$$

When $c \ge \frac{3}{4}$, this inequality always holds, and the additional condition is satisfied. Therefore, it can be applied to the third case of the master theorem, and the solution is $T(n) = \Theta(n\log n)$.

**Example 4:**

$T(n) = 2T(n/2) + n\log n$

For this formula, $a = 2$, $b = 2$, $f(n) = n\log n$. Since $n^{\log_b a} = n^{\log_2 2} = n$, and $f(n) = n\log n$ is asymptotically larger than $n$, it should consider the third case of the master theorem. Substituting $a$, $b$, and $f(n)$ into the inequality, we get

$$
2\frac{n}{2}\log\frac{n}{2} \le cn\log n \implies \log n - 1 \le c\log n \implies (1-c)\log n \le 1
$$

When $c \ge 1$, this inequality always holds. However, the additional condition requires $c < 1$, so this case cannot be applied to the third case of the master theorem. In other words, not all recursive formulas of the form $T(n) = a T(n/b) + f(n)$ can be solved using the master theorem; there is a gap between cases two and three. In such cases, we can only use the recursion tree method. Students can try it, and the result will be $T(n) = \Theta(n \log^2 n)$.

After going through these four examples, have you mastered the master theorem? If you have any questions, please leave a comment, and I will reply as soon as possible.

The study of algorithms is endless, and the calculation of time complexity is fundamental and should not be underestimated.

## Details

Careful readers may notice that in the fourth case mentioned above, although $nlgn$ grows asymptotically faster than $n$, it is not asymptotically greater in the polynomial sense because there is no $\epsilon > 0$ such that $nlgn = \Omega \left( n^{1+\epsilon} \right)$ holds. On the contrary, $nlgn = o \left( n^{1+\epsilon} \right)$ holds for any $\epsilon > 0$. So, in fact, this case does not need to be judged for regular conditions to draw conclusions in advance. However, in the third case in the previous text, I did not emphasize asymptotic greater in the polynomial sense, not because I forgot to write it, but because this condition is not needed at all.

For any recurrence relation that fits the master theorem format, if it satisfies the regular conditions in case three, then it must satisfy $f\left(n\right) = \Omega \left( n^{log_b a + \epsilon} \right)$." The specific proof details can be found in [this answer](https://cs.stackexchange.com/a/121749/128497) on StackExchange, and are not detailed here.

Although regular conditions can imply $f\left(n\right) = \Omega \left( n^{log_b a + \epsilon} \right)$, the reverse is not true. We can construct the following function:

$T\left(n\right) = T\left(n/2\right) + n\left( 2 - \cos n \right)$

Here, $a=1$, $b=2$, and $f\left(n\right) = n \left( 2 - \cos n \right) = \Omega \left( n^{log_2 1+ \epsilon} \right) = \Omega \left( n^\epsilon \right)$ holds for any $\epsilon \le1$. However, when applying the regular condition:

$\begin{aligned} f\left(n / 2\right) &\le c f\left(n\right) \\ \frac{n}{2} \left(2 - \cos \frac{n}{2} \right) & \le c \left( 2 - \cos n \right) \end{aligned}$

For any $c \lt 1$, the right side must be $c \left( 2 - \cos n \right) \lt 3$. However, the left side, when $n$ is large enough, must be greater than $3$. Therefore, this equation cannot hold, and the regular condition is not satisfied.

## Reference

These cases are from the Wikipedia entry on [Master theorem](<https://en.wikipedia.org/w/index.php?title=Master_theorem_(analysis_of_algorithms)&section=7#Inadmissible_equations>).


In this article, we will model the time complexity of divide and conquer using mathematical methods, analyze its asymptotic properties, and provide three methods of calculation.

The Master Theorem For Time Complexity


## Prologue

_Fast Inverse Square Root_ is a rather famous algorithm, and the specific author is no longer traceable. Its prominence comes from its extensive use in the game Quake III, where it became well-known.

Why do we need the inverse square root instead of the square root itself? In computer graphics, many calculations involve the inverse square root rather than the square root. For example, the `normalize()` operation for a vector $\frac{1}{\sqrt{x^2 + y^2 + z^2}}(x,y,z)$ includes the inverse square root. In computer graphics, each frame's rendering requires at least tens of thousands of calls to normalize. If we can accelerate the calculation of the inverse square root, the rendering speed will significantly improve.

The algorithm originated more than 20 years ago. Today, the algorithm has research value but not practical use since modern CPUs have faster and more accurate instructions. We can simply use `1/sqrt(y)`.

However, this doesn't hinder us from understanding the elegance of this algorithm. In a brief overview, it involves:

1. IEEE 754 Standard
2. Undefined behavior in the C language
3. Newton's iteration method.

The focus here is on the second step—how to derive the magic number `0x5f3759df`. For more detailed information, you can refer to the source.

## Algorithm

```c
float Q_rsqrt(float number)
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;
    x2 = number * 0.5F;
    y = number;

    /* Core Code */
    // 1. reinterpret_cast from float to int
    i = *(long *)&y; // evil floating-point bit-level hacking
    // 2. Estimate the inverse square root
    i = 0x5f3759df - (i >> 1); // what the heck? (....why not 0x69696969?)
    // reinterpret_cast from int to float
    y = *(float *)&i;
    // 3. Newton's method
    y = y * (threehalfs - (x2 * y * y)); // 1st iteration
    // y = y * (threehalfs - (x2 * y * y)); // 2nd iteration, this can be removed

    return y;
}
```

The algorithm roughly consists of three parts:

- `reinterpret_cast` used for conversion between float and int because any compiler would otherwise prohibit bitwise operations on float types.
- Estimation using the magic number `0x5f3759df` for the inverse square root (explained later).
- One iteration of Newton's method (explained later).

## Estimating the Inverse Square Root

In the IEEE 754 standard, a 32-bit floating-point number contains 1 sign bit, 8 exponent bits $E$, and 23 mantissa bits $M$. Any floating-point number can be represented as: $(1.M) \times 2^{E - 127} = (1 + \frac{M}{2^{23}}) \times 2^{E - 127}, M \in[0, 2^{23} ), E \in [0,256)$.

To calculate $y$'s inverse square root $\frac{1}{\sqrt{y}} = y^{-\frac{1}{2}}$, the IEEE 754 form of $y$ is $y = (1 + \frac{M_y}{2^{23}}) \times 2^{E_y - 127}$.

For simplicity, let's take the logarithm of $y$ to the base 2:

$$
\begin{align}
\log y &= \log{((1 + M_y) \times 2^{E_y - 127})} \\
&= \log{(1+\frac{M_y}{2^{23}})} + E_y - 127
\end{align}
$$

It's clear that $\frac{M_y}{2^{23}}$ is in the interval $[0,1]$. In this interval, $\log(1+x) \approx x$. If we add a constant $\mu$, the approximation improves:

Using this approximation, we can further simplify:

$$
\begin{align}
\log y
&= \log{(1+\frac{M_y}{2^{23}})} + E_y - 127 \\
&\approx \frac{M_y}{2^{23}} + \mu + E_y - 127 \\
&= \frac{1}{2^{23}}(M_y + E_y \times 2^{23}) + \mu - 127
\end{align}
$$

Now, let's take the logarithm of $y^{-\frac{1}{2}}$:

$$
\log y^{-\frac{1}{2}} = -\frac{1}{2} \log y \approx -\frac{1}{2}\left( \frac{1}{2^{23}}(M_y + E_y \times 2^{23}) + \mu - 127 \right)
$$

Assuming the accurate value of $y^{-\frac{1}{2}}$ is $A$, with exponent and mantissa $E_A$ and $M_A$ respectively, we can take the logarithm of $A$:

$$
\log A \approx \frac{1}{2^{23}}(M_A + E_A \times 2^{23}) + \mu - 127
$$

By equating the two logarithmic expressions, we get:

$$
\begin{align}
M_A + E_A \times 2^{23} &= 2^{23} \times \frac{3}{2}(127 - \mu) - \frac{1}{2}(M_y + E_y \times 2^{23})
\end{align}
$$

When $\mu$ is set to $0.0450461875791687011756$, the first term becomes `1597463011` in hexadecimal (`0x5F3759e3`). If we replace the second term's $\frac{1}{2}$ with a right shift operation, the result is remarkably close to the line in the code:

```c
i = 0x5f3759df - (i >> 1);
```

The magic number is slightly different, but with an appropriate $\mu$ value, it should yield exactly the same number.

## Newton's Iteration Method

The update formula for Newton's iteration method is:

$$
x_{t+1} = x_t - \frac{f(x_t)}{f^\prime(x_t)}
$$

The precondition for Newton's iteration method is that the initial point should be near the function's zero. Since the estimation from `i = 0x5f3759df - (i >> 1);` should have a small error, we can assume the precondition is satisfied.

Now, our problem is, given $y$, we need to find $x$ such that $\frac{1}{\sqrt{y}} = x$, i.e., $\frac{1}{x^2} = y$. Rearranging, we get $f(x) = \frac{1}{x^2} - y = 0$. Applying Newton's iteration to $f(x)$ gives the iteration formula:

$$
x_{t+1} = x_t - \frac{y}{2}x_t^3
$$

This is translated into code as:

```c
y = y * (threehalfs - (x2 * y * y));
```

## References

[1] [Is fast inverse square root still used?](https://www.quora.com/Is-fast-inverse-square-root-still-used)

[2] [Fast Inverse Square Root — A Quake III Algorithm](https://www.youtube.com/watch?v=p8u_k2LIZyo)

[3] [Relevant Thesis](https://www.lomont.org/papers/2003/InvSqrt.pdf)


Magic Number - 0x5f3759df

The article introduces an algorithm called Fast Inverse Square Root, which is used to quickly calculate the inverse square root of a floating point number.

Fast Inverse Square Root


Brief Introduction to the Implementation Approach in glibc:

For the logarithm of any base, we can utilize the change of base formula:

$$
\log_k{x} = \frac{\ln{x}}{\ln k} \tag{1}
$$

The term $\ln k$ in equation (1) is a constant, so for logarithms of any base, we only need to calculate $\ln x$.

### Computation of $\ln x$

In advanced calculus, we learned the following Taylor series:

$$
\ln(1+x) = x - \frac{1}{2}x^2 + \frac{1}{3}x^3 - \frac{1}{4}x^4 + \dots = \sum_{k=1}^{\infty}\frac{(-1)^k}{k}x^k \tag{2}
$$

Suppose we want to find $\ln x$; naturally, we can set $t=x-1$, then:

$$
\ln(x) = \ln(1+t) = \sum_{k=1}^{\infty}\frac{(-1)^k}{k}t^k \tag{2}
$$

However, expanding the series to infinity in practical calculations is often not feasible. If we expand it to the $n$th term, the Taylor series with the Peano remainder term can be expressed as:

$$
\ln(1+t) = t - \frac{1}{2}t^2 + \frac{1}{3}t^3 - \frac{1}{4}t^4 + \dots = \sum_{k=1}^{n}\frac{(-1)^k}{k}t^k + o(t^n) \tag{3}
$$

Note that the Peano remainder term here is $o(t^n)$, representing the error term. If $t>1$ (or $x > 2$), the error in the above formula may become significantly large with the increase of $t$ or $x$.

However, another formula is available:

$$
\ln x = 2\left(\frac{x-1}{x+1} + \frac{1}{3}\left(\frac{x-1}{x+1}\right)^3 + \frac{1}{5}\left(\frac{x-1}{x+1}\right)^5 + \dots\right) = 2\sum_{k=0}^{\infty}\frac{1}{2k+1}\left(\frac{x-1}{x+1}\right)^{2k+1} \tag{4}
$$

Here's an informal proof for equation $(4)$:

Replace $x$ in equation $(2)$ with $\frac{1}{x}$ and $-\frac{1}{x}$, respectively:

$\begin{aligned}
\ln\left(1+\frac{1}{x}\right) &= \ln\left(\frac{x+1}{x}\right) = \sum_{k=1}^{\infty}\frac{(-1)^k}{k}\frac{1}{x^k} \\
\ln\left(1-\frac{1}{x}\right) &= -\sum_{k=1}^{\infty}\frac{1}{kx^k}
\end{aligned}$

Note:

$\begin{aligned}
\ln\left(\frac{x-1}{x}\right) &= -\sum_{k=1}^{\infty}\frac{1}{kx^k} \\
\Longrightarrow -\ln\left(\frac{x-1}{x}\right) &= \sum_{k=1}^{\infty}\frac{1}{kx^k} \\
\Longrightarrow \ln\left(\frac{x}{x-1}\right) &= \sum_{k=1}^{\infty}\frac{1}{kx^k}
\end{aligned}$

Then:

$\begin{aligned}
&\ln\left(\frac{x+1}{x}\right) + \ln\left(\frac{x}{x-1}\right) = \ln\left(\frac{x+1}{x-1}\right) \\
&\sum_{k=1}^n\frac{(-1)^k}{kx^k} + \sum_{k=1}^n\frac{1}{kx^k} = 2\sum_{k=0}^n{\frac{1}{2k+1}\frac{1}{x^{2k+1}}}
\end{aligned}$

This yields:

$$
\ln\left(\frac{x+1}{x-1}\right) = 2\sum_{k=0}^n{\frac{1}{2k+1}\frac{1}{x^{2k+1}}}
$$

Let $z = \frac{(x+1)}{(x-1)}$, then $x = \frac{(z-1)}{(z+1)}$. Substituting this into the above formula gives equation (4).

The advantage of this formula is that if we expand it to the $n$th term:

$$
\ln x = 2\sum_{k=0}^{n}\frac{1}{2k+1}\left(\frac{x-1}{x+1}\right)^{2k+1} + o\left(\left(\frac{x-1}{x+1}\right)^{2k+1}\right) \tag{5}
$$

The remainder term is $o\left(\frac{(x-1)}{(x+1)^{2k+1}}\right)$, and as $\lim_{x\rightarrow+\infty}\left(\frac{x-1}{x+1}\right)^{2k+1} = 1$, the remainder term ensures that the error remains within the range of $1$.

However, this is still not sufficient, as the error of $1$ is still too large. According to the analysis of the remainder term above, the value of $x$ should be closer to $1$ for higher precision.

Considering the IEEE 754 floating-point standard, we know that any floating-point number is represented in the computer as:

$$
x = 1.xxx \times 2^{exp}
$$

This means that we can directly obtain the exponent ($\text{exp}$) and the fractional part ($1.xxx$) using bitwise operations.

Since the floating-point number is already in this form, we can directly use the fractional part $1.xxx$. If we calculate $\ln(1.xxx)$ using the earlier formula, it is equivalent to calculating $\ln\left(\frac{x}{2^{exp}}\right)$. According to the logarithmic formula:

$$
\ln\left(\frac{x}{2^{exp}}\right) = \ln x - \text{exp}\times \ln 2
$$

Therefore:

$$
\ln x = \ln\left(\frac{x}{2^{exp}}\right) + \text{exp} \times \ln 2
$$

This completes the calculation of $\ln x$, equivalent to computing the logarithm of any base.

However, this is a general overview, and glibc employs various precision improvement techniques in its implementation.


Take the Glibc's Mindset for Example

How Do Computer Calculate the Log-Function


In Fibonacci's work _The Book of Calculation_ the Fibonacci sequence is defined as follows:

$$
F_n = \begin{cases}
0 & \text{if } n = 0 \\
1 & \text{if } n = 1 \\
F_{n-1} + F_{n-2} & \text{if } n \geq 2
\end{cases}
$$

It can be proven that its closed-form formula is:

$$
F_n = \frac{1}{\sqrt{5}}\left(\left(\frac{1+\sqrt{5}}{2}\right)^n - \left(\frac{1-\sqrt{5}}{2}\right)^n\right)
$$

With my current knowledge, I can only comprehend two proof methods as follows:

## Matrix Method

Let's first discuss the case when $n \geq 2$. Our goal now is to transform the recurrence formula of the Fibonacci sequence into matrix form. How can we do that? We can approach it from the perspective of a system of linear equations.

First, here is what we know:

$$
F_{n-1} + F_{n-2} = F_{n}
$$

We can add the equation $F_{n-1} + 0 \cdot F_{n-2} = F_{n-1}$ to form a system of linear equations:

$$
\begin{cases}
F_{n-1} + F_{n-2} = F_{n} \\
F_{n-1} + 0 \cdot F_{n-2} = F_{n-1}
\end{cases}
$$

This can be transformed into matrix form as follows:

$$
\begin{bmatrix}
1 & 1 \\
1 & 0
\end{bmatrix}
\begin{bmatrix}
F_{n-1} \\
F_{n-2}
\end{bmatrix}
=
\begin{bmatrix}
F_{n} \\
F_{n-1}
\end{bmatrix}
$$

Now, we can iterate this process:

$$
\begin{align}
\begin{bmatrix}
F_{n} \\
F_{n-1}
\end{bmatrix}
&=
\begin{bmatrix}
1 & 1 \\
1 & 0
\end{bmatrix}
\begin{bmatrix}
F_{n-1} \\
F_{n-2}
\end{bmatrix} \\
&=
\begin{bmatrix}
1 & 1 \\
1 & 0
\end{bmatrix}^2
\begin{bmatrix}
F_{n-2} \\
F_{n-3}
\end{bmatrix} \\
&=\cdots \\
&=
\begin{bmatrix}
1 & 1 \\
1 & 0
\end{bmatrix}^{n-1}
\begin{bmatrix}
F_{1} \\
F_{0}
\end{bmatrix}
\end{align}
$$

We denote the matrix $\boldsymbol{A} = \begin{bmatrix}
1 & 1 \\
1 & 0
\end{bmatrix}$. So, the problem becomes finding $\boldsymbol{A}^{n-1}$, and then we can calculate $\boldsymbol{A}^{n}$ and replace all instances of $n$ with $n-1$.

Notice that matrix $\boldsymbol{A}$ is a square matrix, and we can utilize the eigenvalues and eigenvectors of matrices.

Eigenvectors can be understood as vectors that, when right-multiplied by the matrix, result in a vector parallel to themselves. Eigenvalues are the scaling factors by which the eigenvectors are scaled when right-multiplied by the matrix. In other words:

$$
\boldsymbol{A}\boldsymbol{x} = \lambda\boldsymbol{x}\tag{1}
$$

Here, $\lambda$ is the eigenvalue, and the non-zero vector $\boldsymbol{x} \in \mathbb{R}^n$ (where $n$ is the order of the square matrix) is the eigenvector corresponding to the eigenvalue $\lambda$ of matrix $\boldsymbol{A}$.

So, the specific approach is to first find the eigenvalues $\lambda_1$ and $\lambda_2$ (usually, an $n$-order matrix has $n$ eigenvalues) and obtain the diagonal matrix $\rm{diag}\{\lambda_1, \lambda_2\}$. Then, we find an invertible matrix $\boldsymbol{P}$ such that:

$$
\boldsymbol{P}^{-1}\boldsymbol{A}\boldsymbol{P} = \rm{diag}\{\lambda_1, \lambda_2\}
$$

Using matrix multiplication properties:

$$
(\boldsymbol{P}^{-1}\boldsymbol{A}\boldsymbol{P})^n = \boldsymbol{P}^{-1}\boldsymbol{A}(\boldsymbol{P}\boldsymbol{P}^{-1})\boldsymbol{A}(\boldsymbol{P}\cdots\boldsymbol{P}^{-1})\boldsymbol{A}\boldsymbol{P} = \boldsymbol{P}^{-1}\boldsymbol{A}^n\boldsymbol{P}\tag{2}
$$

This allows us to calculate $\boldsymbol{A}^{n}$.

Let's first find the eigenvalues. We can rewrite equation $(1)$ as:

$$
(\boldsymbol{A}-\lambda\boldsymbol{E})\boldsymbol{x} = \boldsymbol{0}\tag{3}
$$

Here, $\boldsymbol{E}$ is the identity matrix, and we can compute that $\boldsymbol{A}-\lambda\boldsymbol{E} = \begin{bmatrix}
1-\lambda & 1 \\
1 & -\lambda
\end{bmatrix}$. To ensure that there is a non-zero solution, we solve for:

$$
\left| \boldsymbol{A}-\lambda\boldsymbol{E} \right| = \begin{vmatrix}
1-\lambda & 1 \\
1 & -\lambda
\end{vmatrix} = \lambda^2 - \lambda - 1 = 0
$$

Solving this equation, we obtain $\lambda_1 = \frac{1+\sqrt{5}}{2}$ and $\lambda_2 = \frac{1-\sqrt{5}}{2}$.

Therefore, the diagonal matrix is:

$$
\rm{diag}\{\lambda_1, \lambda_2\} = \begin{bmatrix}
\frac{1+\sqrt{5}}{2} & 0 \\
0 & \frac{1-\sqrt{5}}{2}
\end{bmatrix}
$$

Assuming the eigenvector is $\boldsymbol{x} = \begin{bmatrix}
x & y
\end{bmatrix}^T$, we can substitute $\lambda_1$ and $\lambda_2$ into equation $(2)$ to obtain two systems of equations:

$$
\begin{cases}
(1-\lambda_1)x + y = 0 \\
x - \lambda_1y = 0
\end{cases},
\begin{cases}
(1-\lambda_2)x + y = 0 \\
x - \lambda_2y = 0
\end{cases}
$$

Setting $y = 1$ in both systems of equations gives us the two eigenvectors of the matrix:

$$
\boldsymbol{x}_1 = \begin{bmatrix}
\frac{1+\sqrt{5}}{2} & 1
\end{bmatrix}^T,
\boldsymbol{x}_2 = \begin{bmatrix}
\frac{1-\sqrt{5}}{2} & 1
\end{bmatrix}^T
$$

So, the invertible matrix $\boldsymbol{P}$ is formed by these two eigenvectors:

$$
\boldsymbol{P} = \begin{bmatrix}
\frac{1+\sqrt{5}}{2} & \frac{1-\sqrt{5}}{2} \\
1 & 1
\end{bmatrix}
$$

Why is it like this? Let $\boldsymbol{P} = \begin{bmatrix}
x_1 & x_2 \\
y_1 & y_2
\end{bmatrix}$ (where $x_1, y_1, x_2, y_2$ are the components of eigenvectors $\boldsymbol{x}_1, \boldsymbol{x}_2$ respectively).

Now, if we calculate $\boldsymbol{P}\rm{diag}\{\lambda_1, \lambda_2\}$, it exactly equals $\begin{bmatrix}
\lambda_1x_1 & \lambda_2x_2 \\
\lambda_1y_1 & \lambda_2y_2
\end{bmatrix}$. This means that $\boldsymbol{A}\boldsymbol{P} = \boldsymbol{P}\rm{diag}\{\lambda_1, \lambda_2\}$ is inevitable. Left-multiplying by $\boldsymbol{P}^{-1}$, we get $\boldsymbol{P}^{-1}\boldsymbol{A}\boldsymbol{P} = \rm{diag}\{\lambda_1, \lambda_2\}$. So, $\boldsymbol{P} = \begin{bmatrix}
x_1 & x_2 \\
y_1 & y_2
\end{bmatrix}$ is reasonable.

Its inverse matrix is easy to calculate:

$$
\boldsymbol{P}^{-1} = \frac{\boldsymbol{P}^*}{\left|\boldsymbol{P}\right|} = \frac{1}{\sqrt{5}}\begin{bmatrix}
1 & \frac{1-\sqrt{5}}{2} \\
1 & \frac{1+\sqrt{5}}{2}
\end{bmatrix}
$$

Substituting into equation $(2)$:

$$
\boldsymbol{A}^n = \boldsymbol{P}(\boldsymbol{P}^{-1}\boldsymbol{A}\boldsymbol{P})^n\boldsymbol{P}^{-1} = \boldsymbol{P}\rm{diag}^n\{\lambda_1, \lambda_2\}\boldsymbol{P}^{-1}
$$

Where the diagonal matrix is:

$$
\rm{diag}^n\{\lambda_1, \lambda_2\} = \begin{bmatrix}
\left(\frac{1+\sqrt{5}}{2}\right)^n & 0 \\
0 & \left(\frac{1-\sqrt{5}}{2}\right)^n
\end{bmatrix}
$$

Therefore,

$$
\boldsymbol{A}^n = \frac{1}{\sqrt{5}}\begin{bmatrix}
\left(\frac{1+\sqrt{5}}{2}\right)^{n+1} + \left(\frac{1-\sqrt{5}}{2}\right)^{n+1} & \left(\frac{1+\sqrt{5}}{2}\right)^{n+1}\left(\frac{1-\sqrt{5}}{2}\right) + \left(\frac{1-\sqrt{5}}{2}\right)^{n+1}\left(\frac{1+\sqrt{5}}{2}\right) \\
\left(\frac{1+\sqrt{5}}{2}\right)^n + \left(\frac{1-\sqrt{5}}{2}\right)^n & \left(\frac{1+\sqrt{5}}{2}\right)^{n}\left(\frac{1-\sqrt{5}}{2}\right) + \left(\frac{1-\sqrt{5}}{2}\right)^{n}\left(\frac{1+\sqrt{5}}{2}\right)
\end{bmatrix}
$$

Then,

$$
\begin{bmatrix}
F_{n} \\
F_{n-1}
\end{bmatrix}
= \boldsymbol{A}^{n-1}
\begin{bmatrix}
F_1 \\
F_0
\end{bmatrix}
= \frac{1}{\sqrt{5}}\begin{bmatrix}
\left(\frac{1+\sqrt{5}}{2}\right)^{n} + \left(\frac{1-\sqrt{5}}{2}\right)^{n} \\
\left(\frac{1+\sqrt{5}}{2}\right)^{n-1} + \left(\frac{1-\sqrt{5}}{2}\right)^{n-1}
\end{bmatrix}
$$

Considering only the first row of the matrices on both sides of the equation, we obtain the closed-form formula for the Fibonacci sequence:

$$
F_n = \frac{1}{\sqrt{5}}\left(\left(\frac{1+\sqrt{5}}{2}\right)^n - \left(\frac{1-\sqrt{5}}{2}\right)^n\right)
$$

Substituting $n=0$ and $n=1$ to verify, we find that they both satisfy the equation.

---

## Difference Equation Method

Defining the difference of a sequence $\{a_n\}$ as $\Delta a_n = a_{n+1} - a_n$ (using backward difference here), we can define the second-order difference as:

$$
\Delta^2 a_n = \Delta a_{n+1} - \Delta a_n = a_{n+2} - 2a_{n+1} + a_n
$$

Further, we can define the $m$-th order difference:

$$
\Delta^m a_n = \Delta^{m-1} a_{n+1} - \Delta^{m-1} a_n = \sum_{i=0}^{m} (-1)^i C_m^i a_{n+m-i}
$$

And

$$
F(a_n, \Delta a_n,
\Delta a_n, \Delta^2 a_n, \Delta^3 a_n, \ldots)
$$

The Fibonacci sequence, defined as $F_n = F_{n-1} + F_{n-2}$, has a second-order difference that is constant:

$$
\Delta^2 F_n = F_{n+2} - 2F_{n+1} + F_n = F_{n+1} + F_{n+1} - 2F_{n+1} + F_n = F_{n+1} - F_n = F_n
$$

Now, we can write the second-order linear homogeneous difference equation for $F_n$:

$$
\Delta^2 F_n - F_n = 0
$$

This is a characteristic equation, and we can solve it by assuming $F_n = r^n$:

$$
r^2 - 1 = 0
$$

Solving this equation gives us two solutions: $r = 1$ and $r = -1$. Therefore, the general solution for the homogeneous difference equation is:

$$
F_n = c_1 \cdot 1^n + c_2 \cdot (-1)^n
$$

Now, we need initial conditions to find the particular solution. We know that $F_0 = 0$ and $F_1 = 1$. Substituting these into the general solution:

$$
\begin{align*}
F_0 &= c_1 \cdot 1^0 + c_2 \cdot (-1)^0 = c_1 + c_2 = 0 \\
F_1 &= c_1 \cdot 1^1 + c_2 \cdot (-1)^1 = c_1 - c_2 = 1
\end{align*}
$$

Solving this system of equations, we find $c_1 = \frac{1}{2}$ and $c_2 = -\frac{1}{2}$. Therefore, the particular solution for the Fibonacci sequence is:

$$
F_n = \frac{1}{2} \cdot 1^n - \frac{1}{2} \cdot (-1)^n
$$

This can be simplified further by noting that $(-1)^n$ is equal to $(-1)^{n-1} \cdot (-1) = -(-1)^{n-1}$:

$$
F_n = \frac{1}{2} - \frac{1}{2} \cdot (-1)^{n-1}
$$

This is indeed the closed-form formula for the Fibonacci sequence:

$$
F_n = \frac{1}{\sqrt{5}}\left(\left(\frac{1+\sqrt{5}}{2}\right)^n - \left(\frac{1-\sqrt{5}}{2}\right)^n\right)
$$

So, we have successfully derived the same result using the difference equation method.

In summary, both the matrix method and the difference equation method lead to the same closed-form expression for the Fibonacci sequence, demonstrating the beauty of mathematics in providing multiple ways to arrive at a solution.


This article gives two methods to derive Fibonacci sequence: matrix method and difference equation method

Use Differential Equation Method and Matrix Method to Find Fibonacci Sequence General Formula


## Proposition

Generally speaking, a proposition is typically composed of a declarative sentence expressing a judgment. Propositions are classified into true and false propositions. Propositions can often be assigned to a symbol for expression, such as $p, q, r$, where these symbols are called propositional variables. A simple proposition is the most basic form of a proposition that cannot be divided into sub-propositions.

For convenience in operations, the truth or falsity of a proposition is represented using binary values $p=0$ and $p=1$, where they respectively indicate false and true propositions.

### Connectives

These connectives are basic tools in formal logic, used to construct and analyze complex logical expressions that describe the combination of propositions.

- Negation (Not), denoted as $\neg$
- Conjunction (And), denoted as $\land$
- Disjunction (Or), denoted as $\lor$
- Conditional (Implies, ImpliedBy), denoted as $\implies$ and $\impliedby$
- Biconditional (If and only if), denoted as $\iff$

Multiple simple propositions are combined into compound propositions using these connectives. The precedence of these connectives decreases from top to bottom.

Their truth table is as follows:

$$
\begin{array}{|c|c|c|c|c|c|c|}
\hline
p & q & \neg p & p \land q & p \lor q & p \implies q & p \iff q \\
\hline
0 & 0 & 1 & 0 & 0 & 1 & 1 \\
\hline
0 & 1 & 1 & 0 & 1 & 1 & 0 \\
\hline
1 & 0 & 0 & 0 & 1 & 0 & 0 \\
\hline
1 & 1 & 0 & 1 & 1 & 1 & 1 \\
\hline
\end{array}
$$

### Well-formed Formulas (WFF)

1. A single propositional variable is a well-formed formula (WFF), also called an atomic proposition.
2. If $A$ is a WFF, then $(\neg A)$ is a WFF.
3. If $A$ and $B$ are WFFs, then $(A \land B)$, $(A \lor B)$, $(A \implies B)$, and $(A \iff B)$ are WFFs.
4. A string of symbols formed by applying (1)~(3) a finite number of times is a WFF.

WFFs are also called propositional formulas, or simply formulas. They are often represented as $f(p_{1}, p_{2}, ..., p_{n})$.

### Formula Levels

1. If the formula $A = f(p)$ is a single propositional variable, then $A$ is called a level-0 formula.
2. A formula $A$ is called a level-$(n+1)$ formula (for $n \geq 0$) if one of the following holds:

   - (a) $A = \neg B$, where $B$ is a level-$n$ formula;
   - (b) $A = B \land C$, where $B$ and $C$ are formulas of levels $i$ and $j$, and $n = \max(i,j)$;
   - (c) $A = B \lor C$, where the levels of $B$ and $C$ and $n$ are as in (b);
   - (d) $A = B \implies C$, where the levels of $B$ and $C$ and $n$ are as in (b);
   - (e) $A = B \iff C$, where the levels of $B$ and $C$ are as in (b).

3. If the level of formula $A$ is $k$, then $A$ is called a level-$k$ formula.
4. The formula $f(p_{1}, p_{2}, ..., p_{n})$ is used to express a formula of level $n$.

### Tautologies and Contradictions

If $\forall p_{i} \in \{0,1\}$, $A = f(p_{1}, p_{2}, ..., p_{n}) \equiv 1$, then $A$ is called a tautology, or a logically valid formula. Conversely, if $\forall p_{i} \in \{0,1\}$, $A = f(p_{1}, p_{2}, ..., p_{n}) \equiv 0$, then $A$ is called a contradiction.

If a formula is not a contradiction, it is also called a satisfiable formula.

### Dummy Variables

A dummy variable (also known as a neutral element) refers to a special element that, under a certain operation, does not change the value of other elements.

If the value of formula $A = f(p_{1}, p_{2}, ..., p_{n})$ depends only on certain variables $p_{i1}, p_{i2}, ..., p_{in}$, then the remaining variables $\{i \in Z^{*}, i \in [1,n] \setminus \{p_{i_1}, p_{i_2}, \dots, p_{in}\}\}$ are called dummy variables. Changing the dummy variables does not affect the value of the formula.

For example, in $p \land q \equiv 1$, where $q = 1$, the formula is always true regardless of the value of $p$, making $p$ a dummy variable.

### Formula Equivalence

If two formulas have the same number of variables and the same truth table (or logical behavior), then these two formulas are equivalent. In other words, if $\forall i \in \{0,1,2,...,n\}$, $p_{i} = q_{i}$, and $f(p_{1}, p_{2}, ..., p_{n}) = g(q_{1}, q_{2}, ..., q_{n})$, then $f \equiv g$, and we can say that $f$ and $g$ are logically equivalent.

## Logical Equivalence Calculus

Let formulas $A$ and $B$ share the same propositional variables, and $A$ or $B$ may contain dummy variables. If $A$ and $B$ have the same truth table, then $A \iff B$ is a tautology. This indicates that $A$ and $B$ are logically equivalent, denoted as $A \iff B$.

Below are some common identities, which can be used to perform more complex logical operations.

### Equivalence Identities

#### 1.1 Double Negation

$$
A \iff \neg \neg A
$$

#### 1.2 Idempotent Laws

$$
A \lor A \iff A, \quad A \land A \iff A
$$

#### 1.3 Commutative Laws

$$
A \lor B \iff B \lor A, \quad A \land B \iff B \land A
$$

#### 1.4 Associative Laws

$$
(A \lor B) \lor C \iff A \lor (B \lor C),\\ (A \land B) \land C \iff A \land (B \land C)
$$

#### 1.5 Distributive Laws

$$
A \lor (B \land C) \iff (A \lor B) \land (A \lor C)
$$

(Distributivity of $\lor$ over $\land$)

$$
A \land (B \lor C) \iff (A \land B) \lor (A \land C)
$$

(Distributivity of $\land$ over $\lor$)

#### 1.6 De Morgan's Laws

$$
\neg (A \lor B) \iff \neg A \land \neg B, \\ \neg (A \land B) \iff \neg A \lor \neg B
$$

#### 1.7 Absorption Laws

$$
A \lor (A \land B) \iff A, \\ A \land (A \lor B) \iff A
$$

#### 1.8 Zero Laws

$$
A \lor 1 \iff 1, \\ A \land 0 \iff 0
$$

#### 1.9 Identity Laws

$$
A \lor 0 \iff A, \\ A \land 1 \iff A
$$

#### 1.10 Law of the Excluded Middle

$$
A \lor \neg A \iff 1
$$

#### 1.11 Law of Non-Contradiction

$$
A \land \neg A \iff 0
$$

#### 1.12 Implication Identity

$$
A \rightarrow B \iff \neg A \lor B
$$

#### 1.13 Biconditional Identity

$$
A \leftrightarrow B \iff (A \rightarrow B) \land (B \rightarrow A)
$$

#### 1.14 Contrapositive

$$
A \rightarrow B \iff \neg B \rightarrow \neg A
$$

#### 1.15 Biconditional Negation

$$
A \leftrightarrow B \iff \neg A \leftrightarrow \neg B
$$

#### 1.16 Reductio ad Absurdum

$$
(A \rightarrow B) \land (A \rightarrow \neg B) \iff \neg A
$$

The 16 equivalence identities presented above contain 24 important equivalences. The variables $A$, $B$, and $C$ in the equivalence patterns can be replaced by any propositional formulas, resulting in infinitely many specific equivalences of the same type.

These are typical examples of equivalences derived from truth tables and form the foundational formulas of logical operations.


Here are some basic mathematical logic definitions, theorems, and other notes.

MY BLOG

MATHEMATICS

The Catalan Series

The Master Theorem For Time Complexity

Fast Inverse Square Root

How Do Computer Calculate the Log-Function

Use Differential Equation Method and Matrix Method to Find Fibonacci Sequence General Formula

The Basic Formal Logic (1)