ALGORITHM

COMPUTATION

CRYPTOGRAPHY

FRONTEND

GAME-THEORY

INVESTMENT

JS&TS

LINUX

MATHEMATICS

NODE

OTHERS

PROJECT-PRACTICE

REACT

THOUGHTS

古文整理

文史哲艺


## Algorithm Overview

Programmers are generally familiar with random numbers, and it is well-known that the random numbers generated by computers are often pseudo-random. What exactly is a pseudo-random sequence?

> **Pseudo-random numbers** are a sequence of numbers that appear random but are generated using a **deterministic** algorithm, meaning they are not truly random.

Since randomness is simulated through an algorithm, what kind of algorithm can achieve an effect close to true randomness, and how does it work?

One of the simpler methods is the linear congruential generator (LCG).

The LCG is based on the recursive formula:

$$
N_{j+1} = (A * N_{j} + B) \% M
$$

The most important elements in the LCG are three integers: the multiplier $A$, the increment $B$, and the modulus $M$.

These values, $A$, $B$, and $M$, are constants set by the generator. The period of the LCG is at most $M$, though in most cases, it will be less than $M$. To achieve the maximum period, the following conditions must be met:

- $B$ and $M$ are coprime;
- All prime factors of $M$ divide $A-1$;
- If $M \% 4 = 0$, then $A-1$ must also be divisible by 4;
- $A$, $B$, and $N_{0}$ must be smaller than $M$, and $A$, $B$ must be positive integers.

The advantage of the LCG is that it is **very easy to implement and generates numbers quickly**. However, its drawback is also clear: for a 32-bit number, the maximum period is limited to $2^{32}$, which is inadequate for applications that require high-quality random numbers, such as cryptography.

## Randomness Analysis

With this understanding, one might wonder, how does the LCG simulate a random process? True randomness implies that any number can appear, meaning **all outcomes are equally probable, following a uniform distribution**. So, how does the LCG achieve uniform distribution?

First, we need to understand the period of the LCG. From the definition above, it is clear that for a modulus $M$, the sequence generated by the LCG has a maximum period of $M$. Suppose the recursive formula is $N_{j+1} = (A * N_{j} + B) \mod M$. If the period is $T$, then $N_{k+T} = N_{k}$, and since $N_{k+1} = (A * N_{k} + B) \mod M$, this implies that each value of $N_{k}$ uniquely determines the value of $N_{k+1}$. Thus, the period $T$ must be less than or equal to $M$, because there are only $M$ possible distinct integers when taken modulo $M$. The $(M+1)$-th number will necessarily be the same as one of the previous numbers, and since the recursive relationship is one-to-one, the subsequent sequence will repeat, leading to a period $T \leq M$.

Next, it is important that the multiplier and modulus be coprime. This is easy to understand: consider $N_{j+1} = (A * N_{j} + B) \mod M$. If $A$ and $M$ are not coprime, suppose $A = a * d$ and $M = m * d$ where $d > 1$. Expanding the recursive formula yields $N_{j+1} = A * N_{j} + B + k * M$. Since $B$ is constant, it can be treated as a shift, and by subtracting $B$ from both sides, we get $N_{j+1} - B = (a * (N_{j} - B) + a * B + k * m) * d$. This shows that all values of $N_{j+1} - B$ contain the factor $d$. For example, if $d = 2$, the sequence of $M$ possible numbers will only yield even numbers plus the shift $B$, meaning the period $T \leq M / 2$. This disrupts the uniform distribution of the sequence.

A few simple examples can illustrate this. Suppose $A = 3$, $M = 5$, and $N_0 = 2$, the sequence would be $\{2, 1, 3, 4, 2, 1, \dots\}$, with a period of 5. If we set $A = 6$, $M = 10$, and $N_0 = 2$, the sequence would be $\{2, 2, 2, \dots\}$, with a period of 1.

Returning to the initial question, how does the LCG achieve uniform distribution? As seen from the discussion above, with carefully chosen parameters, it is possible to construct a sequence with a period of $M$, where each number appears exactly once in each cycle. This ensures that the frequency of appearance for each number in the sequence is the same, thus satisfying the uniform distribution requirement.

## Solving Linear Congruences

Having understood the relationship between random distribution and the linear congruential method, let's extend our knowledge to learn how to solve linear congruences. This is not just a theoretical mathematical problem; it has numerous practical applications. For instance:

> On a circular track, two frogs, $A$ and $B$, start from positions $x$ and $y$ respectively, with $A$ jumping $m$ units and $B$ jumping $n$ units each time. The total circumference of the track is $L$. If both frogs start at the same time, they are said to meet if they land on the same point. The question is, after how many jumps will the two frogs meet at the earliest?

Based on these conditions, we can formulate the following equation:

$$
(x + m*k) ≡ (y + n*k) \pmod{L}
$$

Where $≡$ denotes congruence modulo $L$.

Expanding this gives:

$$
x + m*k = y + n*k + L*k'
$$

Rearranging terms, we get:

$$
x - y = (n - m)*k + L*k'
$$

Thus,

$$
(x - y) ≡ (n - m)*k \pmod{L}
$$

Let $n - m = a$ and $x - y = b$, \
yielding the standard linear congruence:

$$
a*k ≡ b \pmod{L}
$$

> Definition: If $a$ and $b$ are integers, the expression $a*x ≡ b \pmod{M}$, where $x$ is an unknown integer, is called a linear congruence. Here, $≡$ and $(\pmod{M})$ indicate that both sides are taken modulo $M$.

To solve a linear congruence, we proceed with the following transformations:

1. The standard linear congruence $a*x ≡ b \pmod{M}$ is equivalent to the Diophantine equation $a*x + M*y = b$.
2. Let $d$ be the greatest common divisor (gcd) of $a$ and $M$, denoted $gcd(a,M) = d$. If the equation has a solution, $d$ must divide $b$. We rewrite the equation as $a_0*x + M_0*y = b_0$, where $a = a_0*d$, $M = M_0*d$, and $b = b_0*d$. Here, $a_0$ and $M_0$ are coprime, meaning $gcd(a_0, M_0) = 1$.
3. Let $x = x_0*b_0$ and $y = y_0*b_0$. The equation becomes $a_0*x_0 + M_0*y_0 = 1$, or equivalently $a_0*x_0 ≡ 1 \pmod{M_0}$.

Finally, we use modular inversion to compute $x_0$.

> **Theorem:** If $a$ and $M$ are coprime integers and $M > 1$, there exists a unique inverse of $a$ modulo $M$.

This theorem can be proven as follows:

Since $a$ and $M$ are coprime, there exist integers $s$ and $t$ such that $s*a + t*M = 1$. (If $a$ and $M$ have a common divisor $d > 1$, then $s*a + t*M = d*k$, which cannot equal 1). Taking the equation modulo $M$, we get $s*a ≡ 1 \pmod{M}$. Thus, $s$ is the modular inverse of $a$ modulo $M$.

Let $\bar{a}$ denote the modular inverse of $a$ modulo $M$, so $x_0 ≡ \bar{a} \pmod{M_0}$, meaning $x_0 ≡ \bar{a} + k*M_0$.

For example:

To solve $5*x ≡ 7 \pmod{9}$:

1. Since $gcd(5,9) = 1$, we know that $5$ has an inverse modulo $9$. Using the extended Euclidean algorithm, we find $2*5 + (-1)*9 = 1$.
2. Therefore, $2$ is the inverse of $5$ modulo $9$. Multiplying both sides of the original equation by $2$ gives $2*5*x ≡ 2*7 \pmod{9}$.
3. Simplifying, we get $x ≡ 14 \pmod{9} ≡ 5 \pmod{9}$.
4. Thus, the solution is $x = 5 + 9k$.

In summary, we can solve a general linear congruence $a*x ≡ b \pmod{L}$ by transforming it into a standard congruence $a_0*x_0 ≡ 1 \pmod{M_0}$. After solving the standard congruence, we map the solution back. If the solution to the standard equation is $x_0 \in \{x = p + q*M_0\}$, then the solution to the original equation is $x \in \{x = \frac{b}{d}*p + q'*M_0\}$.

Through this, we have outlined the process of solving linear congruences and deepened our understanding of the linear congruential method.

## Insecurity of the Linear Congruential Method

The linear congruential method (LCG) is a pseudo-random number generation technique based on deterministic mathematical algorithms. Although it can meet the needs for general pseudo-random numbers, it is not suitable for high-security applications, such as cryptography or applications requiring strong randomness. Below are some key considerations regarding the insecurity of LCG:

1. **Determinism**: LCG is deterministic, meaning that the same seed and parameters will always produce the same sequence of random numbers. This poses a security risk, as an attacker could reproduce the sequence by guessing the seed and parameters, breaking the randomness.
2. **Periodicity**: The random number sequences generated by LCG are periodic. Once the sequence reaches its period, it starts repeating. This periodicity can cause problems in applications that require long-term randomness.
3. **Non-uniformity**: In some cases, LCG-generated sequences may not be uniformly distributed, especially when parameters are poorly chosen, leading to statistical biases.
4. **Security**: LCG is unsuitable for cryptographic or high-randomness security applications. For such use cases, cryptographically secure pseudo-random number generators (CSPRNGs) must be employed.

In conclusion, while LCG can be used for general pseudo-random number generation in simulations, games, or random testing, it should not be used in applications that require high security and strong randomness.

## Code Implementation

Below is the implementation of a Linear Congruential Generator (LCG) using TypeScript, along with a simple program to test the effectiveness of its randomness.

### LCG Implementation

```typescript
class LinearCongruentialGenerator {
  private seed: number;
  private readonly a: number;
  private readonly c: number;
  private readonly m: number;

  constructor(
    seed: number,
    a: number = 1664525,
    c: number = 1013904223,
    m: number = 2 ** 32
  ) {
    this.seed = seed;
    this.a = a;
    this.c = c;
    this.m = m;
  }

  public next(): number {
    this.seed = (this.a * this.seed + this.c) % this.m;
    return this.seed;
  }

  public nextFloat(): number {
    return this.next() / this.m;
  }
}

// Example usage
const lcg = new LinearCongruentialGenerator(12345);
console.log(lcg.next());
// Outputs a random integer
console.log(lcg.nextFloat());
// Outputs a random floating-point number between 0 and 1
```

### Code to Test Randomness

We can test the randomness of the generator by generating a set of random numbers and plotting their distribution. Below is a simple test program:

```typescript
function testRandomness(
  generator: LinearCongruentialGenerator,
  iterations: number = 10000
): void {
  const results: number[] = new Array(10).fill(0);

  for (let i = 0; i < iterations; i++) {
    const randomValue = Math.floor(generator.nextFloat() * 10);
    results[randomValue]++;
  }

  console.log("Randomness Test Results:");
  results.forEach((count, index) => {
    console.log(`${index}: ${"*".repeat(count / 100)}`);
  });
}

// Example usage
const lcgTest = new LinearCongruentialGenerator(12345);
testRandomness(lcgTest);
```

This test program generates a set of random numbers and counts the occurrences of each number. By observing the distribution of the output, we can make an initial assessment of the randomness of the generator.


Linear congruential method is a pseudo-random number generation method that generates a random number sequence based on a deterministic mathematical algorithm.

Linear Congruential for Pseudorandom Number


**Definition ($k$-th Order Minor):** In an $n$-th order determinant $D$, a $k$-th order minor $M$ is a $k \times k$ determinant formed by selecting any $k$ rows and $k$ columns from $D$ ($k \le n$), and taking the $k^2$ elements at the intersections of these rows and columns in their original order.

When $k < n$, after removing these $k$ rows and $k$ columns from $D$, the remaining $(n-k)$-th order determinant $M'$ formed by the leftover elements is called the cofactor of the $k$-th order minor $M$.

**Example:**

$$
D=\begin{vmatrix}  a_{11} & a_{12} & a_{13} & a_{14} & a_{15}\\  a_{21} & a_{22} & a_{23} & a_{24} & a_{25}\\ \vdots   & \vdots  & \vdots  & \vdots  & \vdots \\ a_{51} & a_{52} & a_{53} & a_{54} & a_{55} \end{vmatrix}, \\ M=\begin{vmatrix} a_{12}  & a_{13} & a_{15} \\ a_{22}  &  a_{23} & a_{25} \\ a_{42}  & a_{43} & a_{45} \end{vmatrix}, \\ M'=\begin{vmatrix} a_{31} & a_{34} \\ a_{51} & a_{54} \end{vmatrix}
$$

**Note:** In the above example, the algebraic cofactor of $M$ is

$$
(-1)^{(1+2+4)+(2+3+5)}M'=-M'
$$

**Definition (Algebraic Cofactor):** Let $M$ be a $k$-th order minor of $D$ with row and column indices $i_{1},i_{2},\cdots,i_{k}; \ j_{1},j_{2},\cdots,j_{k}$, respectively. The cofactor $M'$ of $M$ is multiplied by the sign $(-1)^{(i_{1}+i_{2}+\cdots+i_{k})+(j_{1}+j_{2}+\cdots+j_{k})}$ to form the algebraic cofactor of $M$.

**Lemma:** Every term in the product of any minor $M$ of determinant $D$ and its algebraic cofactor $A$ corresponds to a term in the expansion of determinant $D$, with the same sign.

**Proof:** First, consider the case where $M$ is located in the upper-left corner of determinant $D$:

$$
\begin{vmatrix}  a_{11} & a_{12} & \cdots & a_{1k} & |& a_{1,k+1} & \cdots & a_{1n}\\  \vdots

& \vdots & M & \cdots & | & \vdots &  &\vdots \\  a_{k1} & a_{k2} &\cdots  &  a_{kk}& | & a_{k,k+1} & \cdots & a_{kn}\\ -&- &- &- &-|-&- &- &- \\   a_{k+1,1}& a_{k+1,2} &\cdots  & a_{k+1,k} &  |& a_{k+1,k+1} & \cdots & a_{k+1,n}\\  \vdots & \vdots &  &  \vdots&  |& \vdots & M' & \vdots\\  a_{n1} & a_{n2} & \cdots & a_{nk} &  |& a_{n,k+1} & \cdots & a_{nn}\end{vmatrix}
$$

In this case, the algebraic cofactor of $M$ is:

$$
A=(-1)^{(1+2+\cdots +k)+(1+2+\cdots+k)}M'=M'
$$

Each term in $M$ can be written as $a_{1\alpha_{1}}a_{2\alpha_{2}}\cdots a_{k\alpha_{k}}$, where $\alpha_{1}\alpha_{2}\cdots\alpha_{k}$ is a permutation of $1,2,\cdots,k$. The sign of this term is $(-1)^{\tau(\alpha_1 \alpha_2\cdots\alpha_k)}$.

Each term in $M'$ can be written as $a_{k+1,\beta_{k+1}}a_{k+2,\beta_{k+2}}\cdots a_{n\beta_{n}}$, where $\beta_{k+1}\beta_{k+2}\cdots\beta_{n}$ is a permutation of $k+1,k+2,\cdots,n$, and its sign is $(-1)^{\tau((\beta_{k+1}-k)(\beta_{k+2}-k)\cdots(\beta_{n}-k))}$.

The product of these two terms is $a_{1\alpha_{1}}a_{2\alpha_{2}}\cdots a_{k\alpha_{k}}a_{k+1,\beta_{k+1}}\cdots a_{n\beta_{n}}$, and the combined sign is:

$$
(-1)^{\tau(\alpha_{1}\alpha_{2}\cdots\alpha_{k})+\tau((\beta_{k+1}-k)(\beta_{k+2}-k)\cdots(\beta_{n}-k))}= (-1)^{\tau(\alpha_{1}\alpha_{2}\cdots \alpha_{k}\beta_{k+1}\cdots \beta_{n}) }
$$

Thus, this product is a term in the expansion of determinant $D$ with the same sign.

Now, we will prove the general case.

Let minor $M$ be located in rows $i_1, i_2, \cdots, i_k$ and columns $j_1, j_2, \cdots, j_k$ of $D$, where $i_1 < i_2 < \cdots < i_k$, and $j_1 < j_2 < \cdots < j_k$.

We rearrange the rows and columns of determinant $D$ such that $M$ is moved to the top-left corner.

First, swap the $i_1$-th row with rows $i_1-1, i_1-2, \cdots, 2, 1$, which takes $i_1-1$ swaps to move the $i_1$-th row to the first row. Then, swap the $i_2$-th row $i_2-2$ times to move it to the second row, and so on. In total, this requires

$$
(i_1-1)+(i_2-2)+\cdots+(i_k-k)=(i_1+i_2+\cdots+i_k)-(1+2+\cdots+k)
$$

row swaps to move rows $i_1, i_2, \cdots, i_k$ to rows $1, 2, \cdots, k$.

Using similar column swaps, we move the columns of $M$ to columns $1, 2, \cdots, k$, requiring

$$
(j_1-1)+(j_2-2)+\cdots+(j_k-k)=(j_1+j_2+\cdots+j_k)-(1+2+\cdots+k)
$$

column swaps.

Let $D_1$ denote the new determinant after these transformations. We have:

$$
D_1 = (-1)^{(i_1+i_2+\cdots+i_k)-(1+2+\cdots+k)+(j_1+j_2+\cdots+j_k)-(1+2+\cdots+k)}D \\ = (-1)^{i_1+i_2+\cdots+i_k+j_1+j_2+\cdots+j_k}D
$$

Thus, the expansion of $D_1$ contains the same terms as that of $D$, but with each term having a sign difference of $(-1)^{i_1+i_2+\cdots+i_k+j_1+j_2+\cdots+j_k}$.

Now, $M$ is located in the top-left corner of $D_1$, and the cofactor and algebraic cofactor of $M$ in $D_1$ are both $M'$. Therefore, each term in the product $MM'$ corresponds to a term in $D_1$ with the same sign.

The algebraic cofactor $A$ of $M$ in $D$ gives $MA = (-1)^{i_1+i_2+\cdots+i_k+j_1+j_2+\cdots+j_k}MM'$, which is a term in $D_1$ with a sign difference of $(-1)^{i_1+i_2+\cdots+i_k+j_1+j_2+\cdots+j_k}$. Hence, $MA$ is a term in $D$ with the same sign.

**Theorem (Laplace's Theorem):** For any $k$ ($1 \le k \le n-1$) rows chosen from determinant $D$, the sum of the products of all $k$-th order minors formed from these rows and their respective algebraic cofactors equals the determinant $D$.

**Proof:** Let the minors formed from the chosen $k$ rows of $D$ be $M_1, M_2, \cdots, M_t$, with algebraic cofactors $A_1, A_2, \cdots, A_t$. The theorem claims that:

$$
D = M_1A_1 + M_2A_2 + \cdots + M_tA_t
$$

where

$$
t = \text{C}_{n}^{k} = \dfrac{n!}{k!(n-k)!}
$$

By the lemma, each term in $M_iA_i$ corresponds to a term in $D$ with the same sign, and there are no common terms between $M_iA_i$ and $M_jA_j$ ($i \ne j$).

The left-hand side of the equation contains $n!$ terms, and the right-hand side contains

$$\sum\limits_{i=1}^{t}k!(n-k)! = k!(n-k)!\dfrac{n!}{k!(n-k)!} = n!$$

terms. Hence, the theorem is proven.

**Example:** In the determinant $D=\begin{vmatrix} 1  & 2 & 1 &4 \\  0 & -1 & 2 & 1\\  1 & 0 & 1 & 3\\ 0  &  1& 3 & 1\end{vmatrix}$, selecting the first two rows yields six minors:

$$
M_{1}= \begin{vmatrix}   1& 2\\  0 &-1 \end{vmatrix}=-1, \quad M_{2}= \begin{vmatrix}  1 &1 \\   0& 2\end{vmatrix}=2, \quad M_{3}= \begin{vmatrix} 1  & 4\\  0 &1 \end{vmatrix}=1,
$$

$$
M_{4}= \begin{vmatrix}   2& 1\\   -1&2 \end{vmatrix}=5, \quad M_{5}= \begin{vmatrix}   2& 4\\  -1 &1 \end{vmatrix}=6, \quad M_{6}= \begin{vmatrix}   1& 4\\   2&1 \end{vmatrix}=-7.
$$

Their algebraic cofactors are:

$$
A_1=-8, \quad A_2=3, \quad A_3=-1, \quad A_4=1, \quad A_5=-3, \quad A_6=1 \Rightarrow D=-7.
$$

**Theorem:** The product of two $n$-order determinants

$$
D_{1}=\begin{vmatrix}  a_{11} & a_{12} & \cdots &a_{1n} \\ a_{21}  &a_{22}  &\cdots  & a_{2n}\\ \vdots  & \vdots &  &\vdots \\  a_{n1} & a_{n2} &\cdots  & a_{nn}\end{vmatrix}, \\ D_{2}=\begin{vmatrix}  b_{11} & b_{12} & \cdots &b_{1n} \\ b_{21}  &b_{22}  &\cdots  & b_{2n}\\ \vdots  & \vdots &  &\vdots \\  b_{n1} & b_{n2} &\cdots  & b_{nn}\end{vmatrix}
$$

is equal to a determinant

$$
C = \begin{vmatrix}  c_{11} &c_{12} & \cdots &c_{1n} \\ c_{21}  &c_{22}  &\cdots  & c_{2n}\\ \vdots  & \vdots &  &\vdots \\  c_{n1} & c_{n2} &\cdots  & c_{nn}\end{vmatrix}
$$

where each element $c_{ij}$ is the sum of the products of elements from the $i$-th row of $D_1$ and the corresponding elements from the $j$-th row of $D_2$,

i.e.

$$
c_{ij} = a_{i1}b_{1j} + a_{i2}b_{2j} + \cdots + a_{in}b_{nj}.
$$

**Proof:** Construct a $2n$-order determinant

$$
D = \begin{vmatrix}  a_{11} & a_{12} & \cdots &a_{1n} &0&0&0&0\\ a_{21}  &a_{22}  &\cdots  & a_{2n}&0&0&0&0\\ \vdots  & \vdots &  &\vdots &0&0&0&0\\  a_{n1} & a_{n2} &\cdots  & a_{nn} &0&0&0&0\\ -1&0&\cdots &0&b_{11}&b_{12}&\cdots&b_{1n}\\ 0&-1&\dotso&0&b_{21}&b_{22}&\cdots&b_{2n}\\ \vdots & \vdots & &\vdots&\vdots & \vdots &&\vdots \\ 0&0&\cdots&-1 &b_{n1}&b_{n2}&\cdots &b_{nn} \end{vmatrix}
$$

According to Laplace's theorem, expanding $D$ by the first $n$ rows, all of the minors formed by the first $n$ rows and other elements will be zero except for the minor at the top left, so

$$
D = \begin{vmatrix}  a_{11} & a_{12} & \cdots &a_{1n} \\ a_{21}  &a_{22}  &\cdots  & a_{2n}\\ \vdots  & \vdots &  &\vdots \\  a_{n1} & a_{n2} &\cdots  & a_{nn}\end{vmatrix}\cdot \begin{vmatrix}  b_{11} & b_{12} & \cdots &b_{1n} \\ b_{21}  &b_{22}  &\cdots  & b_{2n}\\ \vdots  & \vdots &  &\vdots \\  b_{n1} & b_{n2} &\cdots  & b_{nn}\end{vmatrix}=D_1D_2.
$$

Next, we show that $D = C$. Perform elementary row operations on $D$: add the $a_{11}$-multiple of the $(n+1)$-th row, the $a_{12}$-multiple of the $(n+2)$-th row, $\dots$, and the $a_{1n}$-multiple of the $(2n)$-th row to the first row to get

$$
D = \begin{vmatrix}  0&0&\cdots&0&c_{11}&c_{12}&\cdots&c_{1n}\\ a_{21}  &a_{22}  &\cdots  & a_{2n}&0&0&0&0\\ \vdots  & \vdots &  &\vdots &0&0&0&0\\  a_{n1} & a_{n2} &\cdots  & a_{nn} &0&0&0&0\\ -1&0&\cdots &0&b_{11}&b_{12}&\cdots&b_{1n}\\ 0&-1&\dotso&0&b_{21}&b_{22}&\cdots&b_{2n}\\ \vdots & \vdots & &\vdots&\vdots & \vdots &&\vdots \\ 0&0&\cdots&-1 &b_{n1}&b_{n2}&\cdots &b_{nn} \end{vmatrix}.
$$

Similarly, successively add the $a_{k1}$-multiple of the $(n+1)$-th row, the $a_{k2}$-multiple of the $(n+2)$-th row, $\dots$, and the $a_{kn}$-multiple of the $(2n)$-th row to the $k$-th row to obtain

$$
D = \begin{vmatrix}  0&0&\cdots&0& a_{11} & a_{12} & \cdots &a_{1n} \\ 0&0&\cdots&0&a_{21}  &a_{22}  &\cdots  & a_{2n}\\ \vdots &\vdots & &\vdots &\vdots  & \vdots &  &\vdots \\  0&0&\cdots&0&a_{n1} & a_{n2} &\cdots  & a_{nn} \\ -1&0&\cdots &0&b_{11}&b_{12}&\cdots&b_{1n}\\ 0&-1&\dotso&0&b_{21}&b_{22}&\cdots&b_{2n}\\ \vdots & \vdots & &\vdots&\vdots & \vdots &&\vdots \\ 0&0&\cdots&-1 &b_{n1}&b_{n2}&\cdots &b_{nn} \end{vmatrix}.
$$

In this determinant, the first $n$ rows contain only one non-zero $n$-order minor. Therefore, by Laplace's theorem,

$$
D = \begin{vmatrix}  c_{11} &c_{12} & \cdots &c_{1n} \\ c_{21}  &c_{22}  &\cdots  & c_{2n}\\ \vdots  & \vdots &  &\vdots \\  c_{n1} & c_{n2} &\cdots  & c_{nn}\end{vmatrix} ~\\ ~\\ \ast (-1)^{(1+2+\cdots+n)+(n+1+n+2+\cdots +2n)} ~\\ ~\\ \ast \begin{vmatrix} -1  & 0 & \cdots & 0\\  0 & -1 & \cdots & 0\\ \vdots  & \vdots &  &\vdots \\   0& 0 & \cdots &-1 \end{vmatrix} = C.
$$


Laplace's theorem in determinants and its applications

The Laplace Theorem for Determinant

This article mainly introduces the basic operations and examples of relational algebra such as union, difference, selection, projection and join.

The Basic Relation Calculation in Database


## Compound Interest Formula

Compound interest refers to the process of calculating interest not only on the principal amount but also on the interest that has already been added to it. The formula for compound interest is:

$$
A = P \left(1 + \frac{r}{n}\right)^{nt}
$$

where:

- $A$ is the final amount of the investment or loan (principal + interest).
- $P$ is the initial principal.
- $r$ is the annual interest rate (expressed as a decimal, e.g., 5% is 0.05).
- $n$ is the number of times interest is compounded per year (if compounded annually, $n = 1$; if compounded quarterly, $n = 4$).
- $t$ is the time the money is invested or borrowed for, in years.

**Example:** If you invest $1,000 at an annual interest rate of 5%, compounded annually, after 5 years, your investment will grow to:

$$
A = 1000 \times \left(1 + \frac{0.05}{1}\right)^{1 \times 5} \\ = 1000 \times (1.05)^5 \\ = 1276.28 \text{ dollars}
$$

## The Rule of 72

The Rule of 72 is a simple heuristic used to estimate the number of years required to double the investment under compound interest. The formula is:

$$
t \approx \frac{72}{r}
$$

where:

- $t$ is the time required to double the investment, in years.
- $r$ is the annual interest rate (as a percentage).

**Example:** If the annual interest rate is 6%, then according to the Rule of 72, the time required to double the investment is approximately:

$$
t \approx \frac{72}{6} = 12 \text{ years}
$$

The Rule of 72 is particularly convenient for quickly estimating the doubling time, especially when the interest rate is not very high (it is generally recommended to use the Rule of 72 for annual interest rates between 6% and 10% for minimal error).

## Proof of the Rule of 72

The mathematical derivation of the Rule of 72 is based on an approximation of the compound interest formula. The following is the mathematical proof of the Rule of 72:

### 1. Compound Interest Formula

We start with the compound interest formula:

$$
A = P \left(1 + \frac{r}{n}\right)^{nt}
$$

Here, we assume interest is compounded annually (i.e., $n = 1$), and we are interested in the time $t$ when the investment doubles, so the formula simplifies to:

$$
2P = P \left(1 + r\right)^t
$$

### 2. Solving for Time $t$

Dividing both sides by $P$, we get:

$$
2 = \left(1 + r\right)^t
$$

Next, we take the natural logarithm (ln) of both sides:

$$
\ln(2) = t \cdot \ln(1 + r)
$$

Solving for $t$:

$$
t = \frac{\ln(2)}{\ln(1 + r)}
$$

### 3. Approximation for Small $r$

For small $r$, we can use the approximation for the natural logarithm:

$$
\ln(1 + r) \approx r
$$

Therefore, the formula can be further simplified to:

$$
t \approx \frac{\ln(2)}{r}
$$

### 4. Calculation and Approximation

We know $\ln(2) \approx 0.693$, so the formula becomes:

$$
t \approx \frac{0.693}{r}
$$

For ease of estimation, 0.693 is often approximated as 0.72 in finance, leading to the Rule of 72:

$$
t \approx \frac{72}{r}
$$

### 5. Conclusion

The Rule of 72, $t \approx \frac{72}{r}$, is an approximate result derived from the above calculations and is suitable for situations with small annual interest rates (typically between 6% and 10%). While the Rule of 72 is not an exact formula, it provides a very practical tool for estimating doubling times in real-world applications.


Compound Interest Formula and The Rule of 72 are commonly used tools in finance for calculating investment growth and estimating the time it takes for an investment to double.

Compound Interest and 72 Formulation


## Problem Statement

The _weighing coin problem_ can be summarized into the following three variations:

1. Given $N$ coins, one of which is heavier (or lighter) than the others. Using a balance scale, what is the minimum number of weighings needed to guarantee finding the heavier (or lighter) coin?
2. Given $N$ coins, one of which has a different weight than the others. Using a balance scale, what is the minimum number of weighings needed to guarantee finding the coin with the different weight and determining whether it is heavier or lighter?
3. Given $N$ coins, one of which has a different weight than the others. Using a balance scale, what is the minimum number of weighings needed to guarantee finding the coin with the different weight without determining whether it is heavier or lighter?

For these three problems, the answers are different. The answer to the first problem is $⌈\log_{3}N⌉$ (where $⌈ x ⌉$ denotes rounding $x$ up to the nearest integer). The answer to the second problem is $⌈\log_{3}(2N+3)⌉$, and the answer to the third problem is $⌈\log_{3}(2N+1)⌉$.

In other words, using the balance scale $k$ times, we can distinguish between $3^k$, $\frac{3^k-3}{2}$, and $\frac{3^k-1}{2}$ coins respectively.

The first problem is relatively simple. Since we know whether the problem coin is heavier or lighter, by putting approximately $\frac{1}{3}$ (rounding up or down) of the coins on each side of the balance scale in each weighing, we can reduce the range of the problem coin to one-third of its previous range. We won't go into further detail here. We mainly discuss the other two problems.

## Strict Proof

Here we strictly prove: **By using the balance scale $k$ times, we can distinguish among at most $\frac{3^k-3}{2}$ coins and determine whether the problem coin is heavier or lighter.**

Let $L_i$ be the set of possibly lighter coins after the $i$th weighing, and $H_i$ be the set of possibly heavier coins.

If we know the problem coin and its weight after $k$ weighings, then:

$$
\text{card}(L_k) + \text{card}(H_k) \leq 1
$$

where $\text{card}(X)$ denotes the number of elements in set $X$.

Considering the $k$th weighing, there are only three possible outcomes: **left side heavier, right side heavier, or balanced**. So, before the $k$th weighing, the number of possibly problem coins can be at most $3$, which means:

$$
\text{card}(L_{k-1}) + \text{card}(H_{k-1}) \leq 3
$$

This conclusion can be recursively generalized:

$$
\text{card}(L_{k-i}) + \text{card}(H_{k-i}) \leq 3^{i}
$$

where $1\leq i\leq k-1$.

Therefore, after the first weighing, we have:

$$
\text{card}(L_{1}) + \text{card}(H_{1}) \leq 3^{k-1}
$$

After the first weighing, suppose we put $a$ coins on each side of the balance scale, and $b$ coins are not on the scale, making a total of $2a+b$ coins.

- If the scale is unbalanced, then either the heavier side or the lighter side contains the problem coin, so $a+a\leq 3^{k-1}$.
- If the scale is balanced, among the remaining $b$ coins, there is a lighter or heavier coin, so $b+b\leq 3^{k-1}$.

Thus, we have:

$$
a \leq \frac{3^{k-1}}{2}     \\
b \leq \frac{3^{k-1}}{2}
$$

Considering $a$ and $b$ are both integers, we have:

$$
a \leq \frac{3^{k-1}-1}{2} \\
b \leq \frac{3^{k-1}-1}{2}
$$

Therefore:

$$
N(k) = 2a+b \leq \frac{3^k-3}{2}
$$

$Q.E.D.$

Furthermore, we consider the third question: if we don't need to determine whether the problem coin is heavier or lighter.

In this case, $L_k$ and $H_k$ can each only contain one coin.

In this situation, it's still true that after $k$ weighings, if we haven't placed the problem coin on the balance scale (otherwise we would have already known if it's heavier or lighter), we can still conclude whether each coin is problematic. Therefore, the number of coins put on the scale cannot exceed the upper limit of the second question. So, the answer to the third question is obviously $\frac{3^k-1}{2}$.

## Solution Construction

We have proven the lower bound of this type of problem, and the other difficulty lies in constructing the solution.

Here is a general method to solve the problem. The process is consistent, but if you need a more practical construction method, you can refer to the end of this post.

Let's take the second question as an example. We use a method called "ternary numbering" to construct the solution.

**Taking $12$ coins as an example again:** We know from the proof that we only need to weigh $3$ times to find the problem coin. Therefore, our numbering will be $3$ digits. According to the following method:

- The first coin is numbered as $001$;
- The second coin's number is generated by moving each digit of the first coin's number $1$ place forward, i.e., $0\rightarrow1$, $1\rightarrow2$, $2\rightarrow0$. So the number of the second coin is $112$;
- Similarly, the third coin is $220$.

The new numbers of these three coins form a "rotation". They meet the following characteristic: **In the new numbers, when we encounter the first different digit from the left, it, and the preceding digit, together, form $01$, $12$, or $20$. We call them "normal" numbers.**

For the fourth coin, according to the new number of the first coin, we find the next ternary number. It needs to satisfy: in the new number from left to right, **the first non-0 digit is 1**. So, the fourth coin's number is $010$. Similarly, we continue this "rotation". The fifth coin's number is $121$, and the sixth coin's number is $202 \dots$

As an additional feature, the conversion between ternary and decimal is as follows:

When the original number is $3a-2$, the new number (normal order) in ternary is: $a+\frac{3^{⌈\log_{3}{2a}⌉}-1}{2}$;

The new number (in reverse order) and the cases when the original number is $3a-1$, $3a$, etc., follow the same pattern.

Here are the normal order and reverse order of the new numbers:

| Original Number | New Number (Normal Order) | Corresponding New Number (Reverse Order, interchange 0 and 2) |
| --------------- | ------------------------- | ------------------------------------------------------------- |
| 1               | 001                       | 221                                                           |
| 2               | 112                       | 110                                                           |
| 3               | 220                       | 002                                                           |
| 4               | 010                       | 212                                                           |
| 5               | 121                       | 101                                                           |
| 6               | 202                       | 020                                                           |
| 7               | 011                       | 211                                                           |
| 8               | 122                       | 100                                                           |
| 9               | 200                       | 022                                                           |
| 10              | 012                       | 210                                                           |
| 11              | 120                       | 102                                                           |
| 12              | 201                       | 021                                                           |

Ternary to decimal conversion:

When the original number is $3a-2$, the new number (normal order) in ternary is: $a+\frac{3^{⌈\log_{3}{2a}⌉}-1}{2}$;

The new number (reverse order) and the cases when the original number is $3a-1$, $3a$, etc., follow the same pattern.

The strategy for weighing is straightforward:

- In the $i$th weighing, we put the coins with the $i$th digit of their new number being $0$ or $2$ on both sides of the balance scale. If the result is left heavy, we record $0$; if it's right heavy, we record $2$; if it's balanced, we record $1$.
- Finally, we record the results of the three weighings, which gives us the number of the problem coin. If it is in normal order, then the problem coin is heavier than normal. If it is in reverse order, then the problem coin is lighter than normal.

This method is easily verifiable. We will give the solutions to the three sub-problems mentioned earlier.

## If the weight of the problem coin is not to be determined

We just need to number the extra last coin with $k$ digits of $1$. For example, for $13$ coins, the number of the last coin is $111$ and we can make $3$ weighings.

## If the amount is less than $\frac{3^k-3}{2}$

1. If the number of coins is a multiple of $3$, we number them as described above. The method of operation with the balance scale remains the same.
2. If the number of coins is one or two more than a multiple of $3$, we add $1$ to the original number of the last coin.

- For example, if there are $10$ coins, we leave out number $10$ and number the last coin as $11$.
- For example, if there are $11$ coins, we leave out number $11$ and number the last coin as $12$;

It is easy to prove that after the first operation with the balance scale, the number of coins with the first digit of the new number being $0$ or $2$ is the same; and since at least $\frac{1}{3}$ of the normal weight coins can be found after the first operation, if the operation encounters a situation where the number of coins with the $i$th digit being $0$ or $2$ is not the same, we just need to add the normal weight coins to make them the same.

## More Intuitive Construction

We have given the complete solution to the coin weighing problem, but in practical operations, it is very cumbersome to assign ternary numbers to coins without the aid of a computer. Is there a more intuitive method?

Yes, there is. We can dynamically construct the solution. Recalling our previous proof, we have the following conclusion:

$$
\text{card}(L_{k-i}) + \text{card}(H_{k-i}) \leq 3^{i}
$$

where $1\leq i\leq k-1$

The core of the problem is: **After each weighing, we need to reduce the solution space to one-third of its original size.**

The method is: Take out $\frac{2}{3}$ (rounding up) of the coins contained in $L_i$ and $H_i$, place them on both sides of the balance scale, and add normal weight coins to make the number of coins on both sides of the scale the same. This ensures that no matter which side of the balance scale is heavier or lighter, the solution space can be reduced to one-third of its original size.

To illustrate this method, let's take the example of $12$ coins again:

- **First step**:
  - Divide the coins into two groups: $1\ 2\ 3\ 4$ and $5\ 6\ 7\ 8$.
    - If the balance is even, then the problem coin is among $9\ 10\ 11\ 12$, go to step $(1a)$.
    - If the balance is uneven, then the problem coin is among $1\ 2\ 3\ 4$ or $5\ 6\ 7\ 8$, go to step $(1b)$.
- **Step $(1a)$**:
  - Put $9\ 10$ on the left side of the balance scale, and $1\ 11$ on the right side.
    - If the balance is even, then the problem coin is $12$, done.
    - If the balance is uneven, then the problem coin is among $9\ 10\ 11$.
- **Step $(1b)$**:
  - Divide the coins into two groups: $1\ 2\ 5$ and $3\ 4\ 6$.
    - If the balance is even, then the problem coin is $7$ or $8$, go to step $(2ba)$.
    - If the balance is uneven, then the problem coin is among $1\ 2$ or $5\ 6$, go to steps $(2bb)$ or $(2bc)$.
- **Step $(2aa)$**:
  - Put $12$ on the left side of the balance scale, and $1$ on the right side.
    - If the balance is even, then the problem coin is $11$, done.
    - If the balance is uneven, then the problem coin is $9$ or $10$.
- **Step $(2ab)$**:
  - Put $9$ on the left side of the balance scale, and $10$ on the right side.
    - If the balance is even, then the problem coin is $11$, done.
    - If the balance is uneven, then the problem coin is among $9\ 10$.
- **Step $(2ba)$**:
  - Put $7$ on the left side of the balance scale, and $8$ on the right side.
    - If the balance is even, then the problem coin is $7$ or $8$, done.
    - If the balance is uneven, then the problem coin is among $7$ or $8$.
- **Step $(2bb)$**:
  - Put $1$ on the left side of the balance scale, and $2$ on the right side.
    - If the balance is even, then the problem coin is $1$ or $2$, done.
    - If the balance is uneven, then the problem coin is among $1$ or $2$.
- **Step $(2bc)$**:
  - Put $3$ on the left side of the balance scale, and $4$ on the right side.
    - If the balance is even, then the problem coin is $3$ or $4$, done.
    - If the balance is uneven, then the problem coin is among $3$ or $4$.

This gives us a more intuitive dynamic solution.


The 'Coin Weight Problem' aims to solve the following mathematical problem: Among N coins, one coin has a different weight. with a balance scale, ensure finding the different coin with the fewest number of weighings.

The Problem of Weighing Coins with Balance.


Today, we mainly write about some problems related to combinatorial mathematics.

The Bernoulli numbers, denoted as $B_n$, are a sequence of rational numbers closely related to number theory. They appear in expressions for sums of powers of consecutive integers and have various applications in mathematics, including number theory, combinatorics, and calculus.

## Bernoulli Numbers

The Bernoulli numbers, denoted as $B_n$, form a rational number sequence closely related to number theory. The first few discovered Bernoulli numbers are:

$$B_0 = 1, \quad B_1 = -\frac{1}{2}, \quad B_2 = \frac{1}{6}, \quad B_3 = 0, \quad B_4 = -\frac{1}{30}, \quad \dots $$

## Sums of Equal Powers

The Bernoulli numbers are named after Jacob Bernoulli, who discovered remarkable relationships while investigating the formulas for sums of $m$th powers. Let's denote

$$S_{m}(n) = \sum_{k=0}^{n-1} k^m = 0^m + 1^m + \dots + (n-1)^m $$

Bernoulli observed the following sequence of formulas, outlining a pattern:

$$
\begin{align*}
S_0(n) &= n \\
S_1(n) &= \frac{1}{2}n^2 - \frac{1}{2}n \\
S_2(n) &= \frac{1}{3}n^3 - \frac{1}{2}n^2 + \frac{1}{6}n \\
S_3(n) &= \frac{1}{4}n^4 - \frac{1}{2}n^3 + \frac{1}{4}n^2 \\
S_4(n) &= \frac{1}{5}n^5 - \frac{1}{2}n^4 + \frac{1}{3}n^3 - \frac{1}{30}n
\end{align*}
$$

It can be observed that in $S_m(n)$, the coefficient of $n^{m+1}$ is always $\frac{1}{m+1}$, the coefficient of $n^m$ is always $-\frac{1}{2}$, the coefficient of $n^{m-1}$ is always $\frac{m}{12}$, the coefficient of $n^{m-3}$ is always $-\frac{m(m-1)(m-2)}{720}$, and the coefficient of $n^{m-4}$ is always zero.

Moreover, the coefficient of $n^{m-k}$ is always a constant times $m^{\underline{k}}$, where $m^{\underline{k}}$ denotes the falling factorial power, i.e., $\frac{m!}{(m-k)!}$.

## Recurrence Formula

$$
\begin{align*}
S_m{(n)} &= \frac{1}{m+1}(B_0n^{m+1}+\binom{m+1}{1}B_1 n^m+\dots+\binom{m+1}{m}B_m n) \\
&= \frac{1}{m+1}\sum_{k=0}^{m}\binom{m+1}{k}B_kn^{m+1-k}
\end{align*}
$$

The Bernoulli numbers are defined by the implicit recurrence relation:

$$
\begin{align*}
\sum_{j=0}^{m}\binom{m+1}{j}B_j &= 0, \quad (m>0) \\
B_0 &= 1
\end{align*}
$$

For example, $\binom{2}{0}B_0+\binom{2}{1}B_1=0$. The first few values are:

$$
\begin{array}{|c|c|c|c|c|c|c|c|c|c|c|}
\hline
n & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & \dots \\
\hline
B_n & 1 & -\frac{1}{2} & \frac{1}{6} & 0 & -\frac{1}{30} & 0 & \frac{1}{42} & 0 & -\frac{1}{30} & \dots \\
\hline
\end{array}
$$

### Proof

#### Proof by Induction

This proof method is adapted from Concrete Mathematics 6.5 BERNOULLI NUMBER.

Using binomial coefficient identities and induction:

$$
\begin{align*}
S_{m+1}(n)+n^{m+1} &= \sum_{k=0}^{n-1}(k+1)^{m+1} \\
&= \sum_{k=0}^{n-1}\sum_{j=0}^{m+1}\binom{m+1}{j}k^j \\
&= \sum_{j=0}^{m+1}\binom{m+1}{j}S_j(n)
\end{align*}
$$

Let $\hat{S}_{m}(n)=\frac{1}{m+1} \sum_{k=0}^{m} \binom{m+1}{k}B_kn^{m+1-k}$, we want to prove $S_m(n)=\hat{S}_m(n)$. Assume that for $j\in[0,m)$, $S_j(n)=\hat{S}_j(n)$.

Subtracting $S_{m+1}(n)$ from both sides:

$$
\begin{align*}
S_{m+1}(n)+n^{m+1} &=\sum_{j=0}^{m+1}\binom{m+1}{j}S_j(n) \\
n^{m+1} &=\sum_{j=0}^{m}\binom{m+1}{j}S_j(n) \\
&=\sum_{j=0}^{m-1}\binom{m+1}{j}\hat{S}_j(n)+\binom{m+1}{m}S_m(n)
\end{align*}
$$

By adding $\binom{m+1}{m}\hat{S}_m(n)-\binom{m+1}{m}\hat{S}_m(n)$ and simplifying:

$$
\begin{align*}
n^{m+1}&=\sum_{k=0}^{m}\frac{n^{k+1}}{k+1}\sum_{j=k}^{m}\binom{m+1}{j}\binom{j}{k}B_{j-k}+(m+1)\Delta
\end{align*}
$$

Let $\Delta = S_m(n)-\hat{S}_m(n)$, and expand $\hat{S}_j(n)$:

$$
\begin{align*}
n^{m+1}&=\sum_{k=0}^{m}\frac{n^{k+1}}{k+1}\sum_{j=k}^{m

}\binom{m+1}{j}\binom{j}{k}B_{j-k}+(m+1)\Delta\\
&=\sum_{k=0}^{m}\frac{n^{k+1}}{k+1}\binom{m+1}{k}\sum_{j=k}^{m}\binom{m-k+1}{j-k}B_{j-k}+(m+1)\Delta\\
\end{align*}
$$

Change the order of summation and apply binomial coefficient identity:

$$
\begin{align*}
n^{m+1}&=\sum_{k=0}^{m}\frac{n^{k+1}}{k+1}\binom{m+1}{k}\sum_{j=0}^{m-k}\binom{m-k+1}{j}B_{j}+(m+1)\Delta
\end{align*}
$$

Using the recursion relation:

$$
\begin{align*}
\sum_{j=0}^{m}\binom{m+1}{j}B_j &= 0, \quad (m>0) \\
B_0 &= 1 \\
\sum_{j=0}^{m}\binom{m+1}{j}B_j &= [m = 0]
\end{align*}
$$

Substituting:

$$
\begin{align*}
n^{m+1}&=\sum_{k=0}^{m}\frac{n^{k+1}}{k+1}\binom{m+1}{k}[m - k = 0]+(m+1)\Delta\\
&=\sum_{k=0}^{m}\frac{n^{k+1}}{k+1}\binom{m+1}{k}+(m+1)\Delta\\
&=\frac{n^{m+1}}{m+1}\binom{m+1}{m}+(m+1)\Delta\\
&=n^{m+1}+(m+1)\Delta
\end{align*}
$$

Thus, $\Delta=0$, and $S_m(n)=\hat{S}_m(n)$.

#### Proof Using Exponential Generating Functions

For the recurrence $\sum_{j=0}^{m}\binom{m+1}{j}B_j=[m=0]$,

adding $B_{m + 1}$ to both sides:

$$
\begin{align*}
\sum_{j=0}^{m+1}\binom{m+1}{j}B_j&=[m=0]+B_{m+1}\\
\sum_{j=0}^{m}\binom{m}{j}B_j&=[m=1]+B_{m}\\
\sum_{j=0}^{m}\dfrac{B_j}{j!}\cdot\dfrac{1}{(m-j)!}&=[m=1]+\dfrac{B_{m}}{m!}
\end{align*}
$$

Let $B(z) = \sum\limits_{i\ge 0}\dfrac{B_i}{i!}z^i$, notice the left side is in convolution form, thus:

$$
\begin{align*}
B(z)\mathrm{e}^z &= z+B(z)\\
B(z)&=\dfrac{z}{\mathrm{e}^z - 1}
\end{align*}
$$

Let $F_n(z) = \sum_{m\ge 0}\dfrac{S_m(n)}{m!}z^m$, then:

$$
\begin{align*}
F_n(z) &= \sum_{m\ge 0}\dfrac{S_m(n)}{m!}z^m\\
&= \sum_{m\ge 0}\sum_{i=0}^{n-1}\dfrac{i^mz^m}{m!}\\
\end{align*}
$$

Changing the order of summation:

$$
\begin{align*}
F_n(z) &=\sum_{i=0}^{n-1}\sum_{m\ge 0}\dfrac{i^mz^m}{m!}\\
       &=\sum_{i=0}^{n-1}\mathrm{e}^{iz}\\
       &=\dfrac{\mathrm{e}^{nz} - 1}{\mathrm{e}^z - 1}\\
       &=\dfrac{z}{\mathrm{e}^z - 1}\cdot\dfrac{\mathrm{e}^{nz} - 1}{z}
\end{align*}
$$

Substituting $B(z)=\dfrac{z}{\mathrm{e}^z - 1}$:

$$
\begin{align*}
F_n(z) &= B(z)\cdot\dfrac{\mathrm{e}^{nz} - 1}{z}\\
&= \left(\sum_{i\ge 0}\dfrac{B_i}{i!} \right)\left(\sum_{i\ge 1}\dfrac{n^i z^{i - 1}}{i!}\right)\\
&= \left(\sum_{i\ge 0}\dfrac{B_i}{i!} \right)\left(\sum_{i\ge 0}\dfrac{n^{i+1} z^{i}}{(i+1)!}\right)
\end{align*}
$$

As $F_n(z) = \sum_{m\ge 0}\dfrac{S_m(n)}{m!}z^m$, i.e., $S_m(n)=m![z^m]F_n(z)$:

$$
\begin{align*}
S_m(n)&=m![z^m]F_n(z)\\
             &= m!\sum_{i=0}^{m}\dfrac{B \times i}{i!}\cdot\dfrac{n^{m-i+1}}{(m-i+1)!}\\
             &=\dfrac{1}{m+1}\sum_{i=0}^{m}\binom{m+1}{i}B_in^{m-i+1}
\end{align*}
$$

Q.E.D


Introduction, derivation and proof of Bernoulli numbers

The Introduction to Bernoulli Number


## Concepts

The computer system's fault always be caused by instruments' and components' aging and inefficient. We trait their lifespan and present statistics information, conclude that every component in system have three stages: starting stage, normal working stage and aging stage.

At first starting stage, the component's working status may not stable, and has higher failure rate. At second stage it works normally and enters into stable, also has lower failure rate. But if it enters the third stage, system reliability will decrease due to aging factors. Its failure rate will re-increase again. This is also known as the _bathtub curve_. So the component should be kept in the second stage as much as possible.

The definition of reliability in computer systems refers to the probability of it operating normally from its initiation at time $t = 0$ to a specified time $t$. This probability can be represented by $R(t)$. The failure rate of a computer system refers to the proportion of components failing per unit time compared to the total number of components. It is denoted by $\lambda$. When $\lambda$ is constant, the relationship between reliability and failure rate can be expressed as

$$
R(t) = e^{-\lambda t}
$$

In computer systems, the time during which the system operates normally between two failure occurrences is termed the system's uptime. The average of all such uptime durations is referred to as the Mean Time Between Failures (MTBF). Its relationship with the failure rate is expressed as:

$$
MTBF = \frac{1}{\lambda}
$$

To express the maintainability of a computer, we commonly use the Mean Time to Repair (MTTR), which represents the efficiency of maintenance. It is the average duration from the occurrence of a failure to the restoration of operation. In computer systems, we also use "availability" to denote the system's operational efficiency, which refers to the probability of the computer system functioning normally at any given time. It can be expressed as:

$$
A = \frac{MTBF}{MTBF + MTTR}
$$

The reliability measurement of a computer typically relies on the aforementioned concepts.

## Reliability of series components

In a series system composed of $N$ components connected sequentially end-to-end, the entire system operates normally only when all components function properly. If the reliability of these components is $R_{1}, R_{2}, R_{3}, ..., R_{N}$, and their failure rates are $\lambda_{1}, \lambda_{2}, \lambda_{3}, ..., \lambda_{N}$, then:

The overall reliability of the system is:

$$
R = R_{1} \times R_{2} \times ... \times R_{N}
$$

The overall failure rate of the system is:

$$
\lambda = \lambda_{1} + \lambda_{2} + \lambda_{3} + ... + \lambda_{N}
$$

## Reliability of parallel components

In a parallel system composed of $N$ components connected in parallel, the entire system operates normally as long as at least one component functions properly. In a parallel system, only one component is essential, while the remaining $N-1$ components are redundant. The more redundant components there are, the higher the system reliability. Assuming the reliability of these components is $R_{1}, R_{2}, R_{3}, ..., R_{N}$:

The overall reliability of the system is calculated as:

$$
R = 1 - (1-R_{1})(1-R_{2})...(1-R_{N})
$$

Assuming the failure rate of all components in the system is $\lambda$, the overall failure rate of the system is:

$$
\mu = \frac{1}{\frac{1}{\lambda}\sum_{j=1}^{N}\frac{1}{j}}
$$

## N-modular redundancy system

An N-modular redundancy system consists of N (where $N = 2n + 1$) identical subsystems and one voting unit. When the system is operational, all subsystems simultaneously compute results, and the voting unit outputs the result that is in the majority, thus serving as the final outcome. This setup helps to mitigate failures caused by random factors, ensuring the reliability of the system. If at least $n+1$ components function properly, then the entire system operates normally.

When a module in the system fails, the remaining modules continue to provide service, ensuring the system's availability. The level of redundancy in an N-modular redundancy system is determined by the value of N, with a higher N value indicating greater redundancy and reliability.

Assuming the voter is entirely reliable and each subsystem has a reliability of $R_{0}$ , the reliability of an N-modular redundancy system is calculated as:

$$
R = \sum_{i=n+1}^{N}\binom{i}{N}R_{0}^{i}(1-R_{0})^{N-i}
$$

In the United States' lunar exploration satellites, N-modular redundancy systems are utilized to construct the computer systems onboard the satellite to mitigate failures in the space environment.


Outlines computer system reliability, including series, parallel components, and N-modular redundancy.

System Reliability Model


For the integration of piecewise continuous functions, the results of "rectangular approximation" and "trapezoidal approximation" are the same, with the difference lying only in the speed of approximation. The essence of this problem lies in "error estimation".

We only need to analyze the "error reduction process" associated with the "approximation process".

To better understand the error, I will specifically calculate and analyze the integrations using rectangular and trapezoidal approximations for continuous functions on closed intervals.

Given $f \in \mathcal{C}[a, b]$, let's start with the expression for rectangular approximation integration:

$$
\lim_{n\to\infty}\sum_{k=1}^n f\left(\frac{(n-k)a+kb}{n}\right)\cdot\frac{(a-b)}{n} = \\
\lim_{n\to\infty}\sum_{k=1}^{n-1} f\left(\frac{(n-k)a+kb}{n}\right)\cdot\frac{(a-b)}{n}+f(b)\cdot\frac{(a-b)}{n} \quad {(1)}
$$

Now, let's look at the expression for trapezoidal approximation:

$$
\lim_{n\to\infty}\sum_{k=1}^n\frac{f\left(\frac{(n-k+1)a+(k-1)b}{n}\right)+f\left(\frac{(n-k)a+kb}{n}\right)}{2}\cdot\frac{(a-b)}{n}\dots \quad {(2)}
$$

Let's rewrite the above expression, **pay attention to the change in the upper limit of summation**, observe:

$$
\begin{align*}
&\lim_{n\to\infty}\sum_{k=1}^n\frac{f\left(\frac{(n-k+1)a+(k-1)b}{n}\right)+f\left(\frac{(n-k)a+kb}{n}\right)}{2}\cdot\frac{(a-b)}{n} \\
&=\lim_{n\to\infty}\sum_{k=1}^n\frac{f\left(\frac{(n-k+1)a+(k-1)b}{n}\right)}{2}\cdot\frac{(a-b)}{n} \\
&\qquad\quad +\lim_{n\to\infty}\sum_{k=1}^n\frac{f\left(\frac{(n-k)a+kb}{n}\right)}{2}\cdot\frac{(a-b)}{n} \\
&=\lim_{n\to\infty}\sum_{k=0}^{n-1}\frac{f\left(\frac{(n-k)a+kb}{n}\right)}{2}\cdot\frac{(a-b)}{n} \\
&\qquad\quad +\lim_{n\to\infty}\sum_{k=1}^n\frac{f\left(\frac{(n-k)a+kb}{n}\right)}{2}\cdot\frac{(a-b)}{n} \\
&=\lim_{n\to\infty}\left[\left(\frac{f(a)}{2}\cdot\frac{(a-b)}{n}+\sum_{k=1}^{n-1}\frac{f\left(\frac{(n-k)a+kb}{n}\right)}{2}\cdot\frac{(a-b)}{n}\right) \right. \\
&\qquad\quad \left. +\left(\sum_{k=1}^{n-1}\frac{f\left(\frac{(n-k)a+kb}{n}\right)}{2}\cdot\frac{(a-b)}{n}+\frac{f(b)}{2}\cdot\frac{(a-b)}{n}\right)\right] \\
&=\lim_{n\to\infty}(\frac{f(a)}{2}\cdot\frac{(a-b)}{n}+\frac{f(b)}{2}\cdot\frac{(a-b)}{n})+ \\
& \lim_{n\to\infty}(2\cdot\sum_{k=1}^{n-1}\frac{f\left(\frac{(n-k)a+kb}{n}\right)}{2}\cdot\frac{(a-b)}{n}) \\
&=\lim_{n\to\infty}(\frac{f(a)+f(b)}{2}\cdot\frac{(a-b)}{n}+\sum_{k=1}^{n-1}f(\frac{(n-k)a+kb}{n})\cdot\frac{(a-b)}{n}) \\
&=\lim_{n\to\infty}\sum_{k=1}^{n-1}f\left(\frac{(n-k)a+kb}{n}\right)\cdot\frac{(a-b)}{n}+\frac{f(a)+f(b)}{2}\cdot\frac{(a-b)}{n} \dots \quad {(3)}
\end{align*}
$$

Comparing equation ${3}$ with equation ${1}$, did you notice? When integrating using two different approximation methods, the difference lies in:

$$
\frac{f(a)+f(b)}{2}\cdot\frac{(a-b)}{n}-f(b)\cdot\frac{(a-b)}{n} \\
=\frac{f(a)-f(b)}{2}\cdot\frac{(a-b)}{n}
$$

We denote this expression as $\Delta(n)$. Notice, it only involves one variable $n$. And it's easy to see:

- $I.$ If $f(a)=f(b)$, then $\Delta(n)\equiv0$.

- $II.$ If $f(a)\neq f(b)$, then $\Delta(n)=O(\frac{1}{n})$.

The first case might give you a feeling of reading an O. Henry novel. Surprisingly, only if $f(a)=f(b)$, no matter how the graph of $f$ fluctuates within the interval, rectangular and trapezoidal approximations are completely consistent, resulting in equal errors. As expected, when you write out the two expressions and observe each term of the summation, you'll find that the terms in both expressions are roughly the same.

The second case leads to a deeper conclusion. Let's denote the error of equation $(1)$ as $\delta_1(n)$ and the error of equation $(3)$ as $\delta_2(n)$, we have:

$$
\delta_1(n)-\delta_2(n)=\Delta(n)=O(\frac{1}{n})
$$

What does this mean? It means that $\delta_1(n)$ and $\delta_2(n)$ are at most one of them has an error term smaller than $o(\frac{1}{n})$! From your intuition, you might have guessed that the error of the trapezoidal approximation is smaller. Correct, the conclusion was proved.


Discusses the equivalence and error estimation in integrating piecewise continuous functions using rectangular and trapezoidal approximations.

Equivalence of Rectangular and Trapezoidal Approximations


<div class="mt-10">
  <div class="flex justify-center my-5">
    <div>
      <span>Do not stand by my grave, and weep.</span><br/>
      <span>I am not there, I do not sleep—</span><br/><br/>
      <span>I am the thousand winds that blow </span><br/>
      <span>I am the diamond glints in snow </span><br/>
      <span>I am the sunlight on ripened grain, </span><br/>
      <span>I am the gentle, autumn rain. </span><br/>
      <span>As you awake with morning's hush, </span><br/>
      <span>I am the swift, up-flinging rush </span><br/>
      <span>Of quiet birds in circling flight, </span><br/>
      <span>I am the day transcending night. </span><br/><br/>
      <span>Do not stand by my grave, and cry— </span><br/>
      <span>I am not there, I did not die.</span><br/>
    </div>
  </div>
  <div class="text-right">
  — Clare Harner, The Gypsy,<br/>
  December 1934
  </div>
</div>

---

<div class="mb-10">
  <div class="flex justify-center my-5">
    <div>
      <span>不要在我的墓前哭泣，</span><br/>
      <span>我不在那里，我没有离去。</span><br/><br/>
      <span>我是扬起的万里清风，</span><br/>
      <span>我是白雪晶莹的光芒，</span><br/>
      <span>我是阳光下成熟的麦穗，</span><br/>
      <span>我是温柔的秋雨纷纷扬扬，</span><br/>
      <span>我是夜幕中闪烁的星光，</span><br/>
      <span>当你在晨曦中醒来时，</span><br/>
      <span>我是盘旋的鸟飞舞直上。</span><br/><br/>
      <span>不要在我的墓前哭泣，</span><br/>
      <span>我不在那里，我没有离去。</span><br/>
    </div>
  </div>
  <div class="text-right">
  — 程虹教授 译
  </div>
</div>


Do Not Stand at My Grave and Weep

我没有离去


Node.js adopts an event-driven and non-blocking I/O model, allowing it to efficiently handle a large number of I/O requests and concurrent connections, which is crucial for a streaming media server. Meanwhile, FFmpeg, as a powerful video processing tool, provides rich video processing capabilities.

Combining Node.js's powerful extension capabilities, we can easily integrate FFmpeg to create a feature-complete video streaming media server, meeting various streaming media application requirements.

Now let's start from scratch and implement a simple video server: implementing the functionality to upload and play videos. In the process of learning, we will gradually understand how the video server works.

Before we begin, you need to:

- Understand the basic API of Node.js, such as the fs module, process module, etc.
- Understand the usage and common APIs of Express.js.
- Install and understand the basic usage of FFMPEG.

Now we want to implement a streaming media server. There are many streaming media protocols and formats on the market (RTP/RTCP/RTSP/RTMP/MMS/HLS, etc.), each with different applications. However, the vast majority of websites currently use HLS/M3U8 format streaming media, which is also supported by modern browsers.

Therefore, this article does not elaborate on these streaming media protocols but implements a simple HLS/M3U8 format streaming media server.

**The streaming media server implemented in this article is not technically difficult.**

## MP4 Video Slicing

Open some mainstream video websites, switch to the "Network" tab in Chrome Devtools, and then play the video. You can see that the browser is constantly downloading small files from a URL. These small files are video segments.

In general, web videos do not directly download the entire video file because a video file is usually several gigabytes, and downloading the entire file would be too large. Instead, the videos are split into many small segments, which are then continuously transferred to the browser. The browser then sequentially plays these segments.

For HLS/M3U8 type streaming media, the video playback address is the address for downloading the m3u8 file. The m3u8 file is a text-based list that stores the names, addresses, and order of video segments. Therefore, when playing this type of streaming media in the browser:

We can directly use ffmpeg to slice the video. Run the command in the terminal.

```bash
ffmpeg -i input_video.mp4 \
-c:v libx264 \
-c:a aac \
-hls_time 10 \
-hls_segment_type mpegts \
-hls_list_size 0 \
-f hls \
-max_muxing_queue_size 1024 \
output.m3u8
```

In this command:

- `-i input_video.mp4` specifies the input video file.
- `-c:v libx264 -c:a aac` specifies the video and audio codecs.
- `-hls_time 10` specifies the duration of each M3U8 segment, in seconds. Here, the duration of each segment is set to 10 seconds.
- `-hls_segment_type mpegts` specifies the type of M3U8 segment as MPEG-TS.
- `-hls_list_size 0` sets the maximum number of segments included in the M3U8 file. Here, it is set to 0, indicating no limit.
- `-f hls` specifies the output format as HLS.
- `-max_muxing_queue_size 1024` sets the maximum multiplexing queue size to ensure that the output does not exceed the specified size.
- The final output file is `output.m3u8`.

Of course, there are more advanced uses of the ffmpeg slicing command, such as limiting the size of all segments to not exceed `500kb`, etc.

![As shown in the image, executing the above command, we generated some video segments and an m3u8 video list.](https://image-assets.dreams.plus/20240323234019.webp)

We open the m3u8 list file, which shows the order and duration of each segment as follows.

```plaintext
#EXTM3U
#EXT-X-VERSION:3
#EXT-X-TARGETDURATION:11
#EXT-X-MEDIA-SEQUENCE:0
#EXTINF:11.386378,
output0.ts
#EXTINF:11.011000,
output1.ts
#EXTINF:9.050711,
output2.ts
#EXTINF:8.591911,
output3.ts
#EXTINF:6.506933,
output4.ts
#EXT-X-ENDLIST
```

## Frontend Playback

After the browser gets the m3u8 file, it sequentially traverses this list and then loads the segments listed in the list in order, i.e., output0.ts → output1.ts ... → output4.ts. Since it loads and plays the video while downloading, it immediately starts playing the first segment after loading the first segment and then simultaneously loads the second segment, and so on.

Next, let's create a new HTML file in the same directory and use `video.js` to play the m3u8 list we just generated.

```html
<html>
  <head>
    <title>Play the video</title>
    <link href="https://vjs.zencdn.net/8.10.0/video-js.css" rel="stylesheet" />
  </head>
  <body>
    <video class="video-js" controls preload="auto" width="640" height="264">
      <source src="./output.m3u8" type="application/x-mpegURL" />
    </video>
    <script src="https://vjs.zencdn.net/8.10.0/video.min.js"></script>
  </body>
</html>
```

![](https://image-assets.dreams.plus/20240323235457.webp)

As you can see, the browser successfully loads and plays the segments listed in the list.

So far, the process of slicing and playing streaming media is roughly understood. Therefore, it is not difficult for us to implement the basic functionality of a streaming media server.

## Design API

Before we start, we need to design the upload and play APIs.

- Video upload: `/upload`, POST request, accepts a file in the `multipart/form-data` form.
- Video playback: `/play/{videoId}/video.m3u8`, GET request, returns a text data of type m3u8. Where `videoId` is the unique identifier for each video.

## Implement Upload Functionality

Now we will implement the video upload functionality for the server. Video uploading involves three stages: receiving the uploaded file → video slicing → video persistent storage. After uploading, return a video playback address to the client for video playback.

Front-end video uploading is generally done through form files. After the server receives the file, it slices the file and stores it in the public area, then returns the generated m3u8 file address.

First, initialize the node project and install the following packages: `express` is the server framework for Node.js, `multer` is a form file processing middleware developed for Express, `nanoid` generates a unique identifier for the uploaded resource, `fluent-ffmpeg` is one of the ffmpeg extensions set for Node.js, `serve` allows you to set up a

server access environment in the local directory.

```bash
npm i multer express nanoid fluent-ffmpeg cors
npm i serve --save-dev
```

Add the following fields to `package.json`:

```json
{
  "main": "main.js",
  "type": "module"
}
```

Then create the `main.js` file and write the following code to create a server:

```js
import express from "express";
import fs from "fs";
import multer from "multer";
import { nanoid } from "nanoid";
import ffmpeg from "fluent-ffmpeg";
import path from "path";
import cors from "cors";

const server = express();
server.listen(3300);
console.log("Server started.");
server.use(cors()); // For easy debugging, allow cross-origin
```

According to the designed API, add a POST route to upload videos:

```js
const copeUpload = multer({
  dest: "uploads-temp/",
  filename: function (req, file, cb) {
    const uniqueSuffix = nanoid();
    cb(null, file.fieldname + "-" + uniqueSuffix);
  },
}).single("video");

server.post("/upload", copeUpload, function (req, res, next) {
  const tempFilePath = path.resolve(req.file.path);
  const videoId = nanoid();
  const storageDirectory = path.resolve("storage", videoId);
  fs.mkdirSync(storageDirectory);

  ffmpeg(tempFilePath)
    .videoCodec("libx264")
    .audioCodec("aac")
    .addOption("-hls_time", 10)
    .addOption("-hls_segment_type", "mpegts")
    .addOption("-hls_list_size", 0)
    .format("hls")
    .addOption("-max_muxing_queue_size", 1024)
    .output(`${storageDirectory}/video.m3u8`)
    .on("start", function () {
      console.log("Starting video slice process.");
    })
    .on("end", function () {
      fs.rmSync(tempFilePath); // Delete the uploaded temporary file
      console.log("OK");
    })
    .on("error", function (err) {
      fs.rmSync(tempFilePath); // Delete the uploaded temporary file
      console.error("Error:", err);
    })
    .run();

  res.json(`http://localhost:3300/play/${videoId}/video.m3u8`); // Return the playback address
});
```

## Implement Playback Functionality

Next is to set up the playback API for the server.

After uploading the video, the server will return a playback address to the client, which is the download address of the m3u8 file. So this part is very simple.

When the player requests the video address `play/abcde/video.m3u8` with the ID `abcde`, since the m3u8 specifies that the slice file and the m3u8 file are in the same directory, the address for requesting the slice is also

```plaintext
play/abcde/output01.ts
play/abcde/output02.ts
play/abcde/output03.ts
......
```

So what the server actually needs to do is to retrieve the files requested by `/play/${videoId}/*` from the `storage` directory and return them to the client.

```js
server.get("/play/:videoId/:filename", (req, res) => {
  const videoId = req.params["videoId"];
  const storageDirectory = path.resolve("storage", videoId);
  if (!existsSync(storageDirectory)) {
    res.status(404).send();
  }
  const filename = req.params["filename"];
  const filepath = path.join(storageDirectory, filename);
  if (!existsSync(filepath)) {
    res.status(404).send();
  }
  const data = fs.readFileSync(filepath);
  res.send(data);
});
```

## Functional Testing

In the above process, our video server is completed in less than 100 lines of code. Next, we will write code for the client to test the upload and playback functions of the server.

Requirements: After uploading the video, return the video playback address and play the specified video address.

Create an `index.html` in the same directory:

```html
<html>
  <head>
    <title>Video Server</title>
    <link href="https://vjs.zencdn.net/8.10.0/video-js.css" rel="stylesheet" />
  </head>
  <body>
    <div>
      <h2>Upload Video Files</h2>
      <form id="uploadForm" method="POST" enctype="multipart/form-data">
        <input type="file" name="video" accept="video/*" />
        <button type="submit">Upload</button>
      </form>
      <div id="response"></div>
      <script>
        const form = document.getElementById("uploadForm");
        form.addEventListener("submit", function (event) {
          event.preventDefault();
          const formData = new FormData(form);
          fetch("http://localhost:3300/upload", {
            method: "POST",
            body: formData,
          })
            .then((response) => response.text())
            .then((data) => {
              document.getElementById("response").innerText = "Upload successful, video address is: " + data;
            })
            .catch((error) => {
              console.error("Request error:", error);
            });
        });
      </script>
    </div>

    <hr />

    <div>
      <h2>Play Video</h2>
      <label for="m3u8-url">Enter the m3u8 address of the video:</label>
      <input type="text" id="m3u8-url" name="m3u8-url" placeholder="For example: https://example.com/video.m3u8" />
      <button onclick="playVideo()">Play</button>

      <video id="my-video" class="video-js" controls preload="auto" width="640" height="360" data-setup="{}">
        <source src="" type="application/x-mpegURL" />
      </video>
      <script src="https://vjs.zencdn.net/8.10.0/video.min.js"></script>
      <script>
        function playVideo() {
          var m3u8Url = document.getElementById("m3u8-url").value;
          var videoPlayer = videojs("my-video");
          videoPlayer.src({
            src: m3u8Url,
            type: "application/x-mpegURL",
          });
          videoPlayer.play();
        }
      </script>
    </div>
  </body>
</html>
```

Then start the server:

```bash
node main.js
```

The server is running at `http://localhost:3300`.

Then start the user interface:

```bash
npx serve
```

You can access it in the browser at `http://localhost:3000`.

Test results are as follows:

![Test result](https://image-assets.dreams.plus/vdsafvdag.gif)

## Conclusion

In this article, we implemented a simple streaming media server and successfully implemented two core functions: video upload and video playback.

Of course, this project is just a demo. It is far from enough for

a video server used in a production environment to have only these two functions. It also requires other more powerful functions such as video resource management, resource health monitoring, network fluctuation monitoring, encryption, and security, etc. I won't go into details about these. If you are interested, you can implement them yourself.

## Download Source Code

The source code repository for this project is [here on GitHub](https://github.com/PrinOrange/simple-node-media-server). It's released under the CC 1.0 open-source license.


Node.js's non-blocking I/O and ffmpeg's video processing allow for building an efficient video streaming server. Let's implement a video server.

MY BLOG

ALL POSTS