
## Wallis's Formula

Let's start with a definite integral:

$$
\int_{0}^{\frac{\pi}{2}} \sin^n x \, dx
$$

We can derive the following result:

$$
\int_{0}^{\frac{\pi}{2}} \sin^n x \, dx = -\int_{0}^{\frac{\pi}{2}} \sin^{n-1} x \, d(\cos x)
$$

Using integration by parts, where $\int u \, dv = uv - \int v \, du$, we get:

$$
-\int_{0}^{\frac{\pi}{2}} \sin^{n-1} x \, d(\cos x) \\ = -\left( \sin^{n-1} x \cos x \right) \Big|_{0}^{\frac{\pi}{2}} + \int_{0}^{\frac{\pi}{2}} \cos x \, d\left(\sin^{n-1} x \right)
$$

$$
= \int_{0}^{\frac{\pi}{2}} \cos^2 x \cdot (n-1) \sin^{n-2} x \, dx \\ = (n-1) \int_{0}^{\frac{\pi}{2}} \left( \sin^{n-2} x - \sin^n x \right) \, dx
$$

Thus:

$$
\int_{0}^{\frac{\pi}{2}} \sin^n x \, dx = (n-1) \int_{0}^{\frac{\pi}{2}} \left( \sin^{n-2} x - \sin^n x \right) \, dx
$$

We get:

$$
n \int_{0}^{\frac{\pi}{2}} \sin^n x \, dx = (n-1) \int_{0}^{\frac{\pi}{2}} \sin^{n-2} x \, dx \\ \implies \int_{0}^{\frac{\pi}{2}} \sin^n x \, dx = \frac{n-1}{n} \int_{0}^{\frac{\pi}{2}} \sin^{n-2} x \, dx
$$

Therefore:

**For even $n$, i.e., $n = 2k$:**

$$
\int_{0}^{\frac{\pi}{2}} \sin^{2k} x \, dx = \frac{2k-1}{2k} \int_{0}^{\frac{\pi}{2}} \sin^{2k-2} x \, dx \\ = \frac{2k-1}{2k} \cdot \frac{2k-3}{2k-2} \cdots \frac{1}{2} \int_{0}^{\frac{\pi}{2}} \sin^0 x \, dx \\= \frac{2k-1}{2k} \cdot \frac{2k-3}{2k-2} \cdots \frac{1}{2} \cdot \frac{\pi}{2}
$$

**For odd $n$, i.e., $n = 2k + 1$:**

$$
\int_{0}^{\frac{\pi}{2}} \sin^{2k+1} x \, dx \\ = \frac{2k}{2k+1} \int_{0}^{\frac{\pi}{2}} \sin^{2k-1} x \, dx \\ = \frac{2k}{2k+1} \cdot \frac{2k-2}{2k-1} \cdots \frac{2}{3} \int_{0}^{\frac{\pi}{2}} \sin^1 x \, dx \\ = \frac{2k}{2k+1} \cdot \frac{2k-2}{2k-1} \cdots \frac{2}{3}
$$

Here, we need the double factorial:

$$
m!! = m (m-2) (m-4) \cdots
$$

Thus:

$$
\int_{0}^{\frac{\pi}{2}} \sin^{2k} x \, dx = \frac{(2k-1)!!}{(2k)!!} \cdot \frac{\pi}{2}
$$

$$
\int_{0}^{\frac{\pi}{2}} \sin^{2k+1} x \, dx = \frac{(2k)!!}{(2k+1)!!}
$$

Since $\sin^{2k+1} x < \sin^{2k} x < \sin^{2k-1} x$ for $x \in (0, \frac{\pi}{2})$, we have:

$$
\int_{0}^{\frac{\pi}{2}} \sin^{2k+1} x \, dx < \int_{0}^{\frac{\pi}{2}} \sin^{2k} x \, dx < \int_{0}^{\frac{\pi}{2}} \sin^{2k-1} x \, dx
$$

Substituting the above formulas:

$$
\frac{(2k)!!}{(2k+1)!!} < \frac{(2k-1)!!}{(2k)!!} \cdot \frac{\pi}{2} < \frac{(2k-2)!!}{(2k-1)!!}
$$

We can examine the difference between the bounds for $\int_{0}^{\frac{\pi}{2}} \sin^{2k} x \, dx$:

$$
\lim_{k \rightarrow \infty} \left( \frac{(2k-2)!!}{(2k-1)!!} - \frac{(2k)!!}{(2k+1)!!} \right) \\ = \lim_{k \rightarrow \infty} \left( \frac{(2k-2)!!}{(2k-1)!!} - \frac{2k}{2k+1} \cdot \frac{(2k-2)!!}{(2k-1)!!} \right) \\ = \lim_{k \rightarrow \infty} \frac{1}{2k+1} \cdot \frac{(2k-2)!!}{(2k-1)!!} = 0
$$

Note: Here $\frac{2k-2}{2k-1} < 1$

For the inequality:

$$
\frac{(2k)!!}{(2k+1)!!} < \frac{(2k-1)!!}{(2k)!!} \cdot \frac{\pi}{2} < \frac{(2k-2)!!}{(2k-1)!!}
$$

Multiply all terms by $\frac{(2k)!!}{(2k-1)!!}$:

$$
\frac{(2k)!!^2}{(2k+1)!! (2k-1)!!} < \frac{\pi}{2} < \frac{(2k-2)!! (2k)!!}{(2k-1)!!^2}
$$

We can further examine the difference between the bounds for $\frac{\pi}{2}$:

$$
\lim_{k \rightarrow \infty} \left( \frac{(2k-2)!! (2k)!!}{(2k-1)!!^2} - \frac{(2k)!!^2}{(2k+1)!! (2k-1)!!} \right) \\= \lim_{k \rightarrow \infty} \left( \frac{1}{2k} - \frac{1}{2k+1} \right) \cdot \frac{(2k)!!^2}{(2k+1)!! (2k-1)!!} \\= \lim_{k \rightarrow \infty} \frac{1}{2k} \cdot \frac{(2k)!!^2}{(2k+1)!!^2} = 0
$$

Thus, we obtain **Wallis's Formula**:

$$
\lim_{k \rightarrow \infty} \frac{(2k)!!^2}{(2k+1)!! (2k-1)!!} = \frac{\pi}{2}
$$

**This formula introduces Wallis's result in preparation for deriving Stirling's formula.**

## Stirling's Formula

Stirling's formula reveals a relationship between $n!$, $n^n$, and $e^n$.

> As $n$ approaches infinity, $n!$ is approximately $\sqrt{2 \pi n} \left(\frac{n}{e}\right)^n$.

First, introduce a limit:

$$
\lim_{n \rightarrow \infty} \left( \frac{n!}{n^n} \right)^{\frac{1}{n}} = \frac{1}{e}
$$

Proof:

$$
\lim_{n \rightarrow \infty} \left( \frac{n!}{n^n} \right)^{\frac{1}{n}} = \lim_{n \rightarrow \infty} e^{\frac{1}{n} \ln \left( \frac{n!}{n^n} \right)}
$$

which is equivalent to proving:

$$
\lim_{n \rightarrow \infty} \frac{1}{n} \ln \left( \frac{n!}{n^n} \right) = -1
$$

We have:

$$
\lim_{n \rightarrow \infty} \frac{1}{n} \ln \left( \frac{n!}{n^n} \right) \\ = \lim_{n \rightarrow \infty} \frac{1}{n} \left( \ln \frac{1}{n} + \ln \frac{2}{n} + \cdots + \ln \frac{n}{n} \right) \\ = \int_{0}^{1} \ln x \, dx = (x \ln x - x) \Big|_{0}^{1} = -1
$$

QED.

Thus,

$$
\frac{n!}{n^n} \sim \frac{1}{e^n}
$$

We would like to verify:

$$
\lim_{n \rightarrow \infty} \frac{n!}{n^n} e^n = 1
$$

Let’s test this.

**Using the ratio test:**

$$
a_n = \frac{n!}{n^n} e^n, \quad a_{n+1} = \frac{(n+1)!}{(n+1)^{n+1}} e^{n+1}
$$

$$
\lim_{n \rightarrow \infty} \frac{a_{n+1}}{a_n} = \frac{e}{\left(1 + \frac{1}{n}\right)^n} > 1
$$

Although the limit of the expression is $e$, it is always approaching $e$ from below, so it is still greater than 1, indicating that $a_n$ is divergent.

Mathematicians have attempted various approaches and found:

$$
a_n = \frac{n!}{n^{n + \frac{1}{2}}} e^n
$$

To prove that $a_n$ converges, use the **ratio test** again:

$$
a_n = \frac{n!}{n^{n + \frac{1}{2}}} e^n, \quad a_{n+1} = \frac{(n+1)!}{(n+1)^{n + \frac{3}{2}}} e^{n+1}
$$

$$
\lim_{n \rightarrow \infty} \frac{a_{n+1}}{a_n} = \frac{e}{\left(1 + \frac{1}{n}\right)^{n + \frac{1}{2}}}
$$

Compare this with:

$$
\left(1 + \frac{1}{n}\right)^{n + \frac{1}{2}} = e^{(n + \frac{1}{2}) \ln \left(1 + \frac{1}{n}\right)}
$$

and $e$.

Thus, compare:

$$
(n + \frac{1}{2}) \ln \left(1 + \frac{1}{n}\right)
$$

with 1.

Specifically, compare $\ln \frac{n+1}{n}$ with $\frac{1}{n + \frac{1}{2}}$:

Since:

$$
\frac{1}{n + 1} < \ln \left(1 + \frac{1}{n}\right) < \frac{1}{n}
$$

**Note:** The right side follows from the inequality $\ln(1 + x) < x$ (for $x \neq 0$).

**The left side follows from:** $1 - \frac{n}{n+1} < \ln \left(\frac{n+1}{n}\right)$, where $\frac{n}{n+1}$ is considered as a whole. This is high school knowledge.

Because $\frac{1}{x}$ is a concave function:

$$
\int_{n}^{n+1} \frac{1}{x} \, dx = \ln \frac{n+1}{n} > \frac{1}{n + \frac{1}{2}}
$$

Therefore:

$$
\lim_{n \rightarrow \infty} \frac{a_{n+1}}{a_n} = \frac{e}{\left(1 + \frac{1}{n}\right)^{n + \frac{1}{2}}} < 1
$$

Thus, $a_n$ is monotonically decreasing and $a_n > 0$, so it must converge, i.e., it has a limit.

Let:

$$
\lim_{n \rightarrow \infty} \frac{n!}{n^{n + \frac{1}{2}}} e^n = a
$$

Then, we need to find this limit.

From the earlier Wallis formula:

$$
\lim_{k \rightarrow \infty} \frac{(2k)!!^2}{(2k+1)!! (2k-1)!!} = \frac{\pi}{2}
$$

$$
(2n)!! = 2n (2n-2) \cdots 2 = 2^n n!
$$

$$
(2n-1)!! = \frac{(2n)!}{(2n)!!} = \frac{(2n)!}{2^n n!}
$$

$$
(2n+1)!! = \frac{(2n+2)!}{(2n+2)!!} = \frac{(2n+2)!}{2^{n+1} (n+1)!}
$$

Thus, Wallis' formula becomes:

$$
\lim_{n \rightarrow \infty} \frac{16^n (n!)^4 (2n+2)}{(2n)!(2n+2)!} = \frac{\pi}{2}
$$

Using:

$$
n! \sim a e^{-n} n^{n + \frac{1}{2}}
$$

Substitute into:

$$
\lim_{n \rightarrow \infty} \frac{a^2 e^2}{4 \left(\frac{n+1}{n}\right)^{2n + \frac{3}{2}}} = \frac{\pi}{2}
$$

Solve for:

$$
a = \sqrt{2 \pi}
$$

Thus, as $n$ approaches infinity:

$$
n! \sim \sqrt{2 \pi} e^{-n} n^{n + \frac{1}{2}}
$$

or equivalently:

$$
n! \sim \sqrt{2 \pi n} \left(\frac{n}{e}\right)^n
$$


Stirling's formula finds a relationship of n!,nⁿ,eⁿ while n trends to infinity.

The Stirling Formula


## Content

The Box-Muller transform is a method used to generate normally distributed random variables from uniformly distributed random variables.

Specifically, it involves selecting two random variables, $U_1$ and $U_2$, which are uniformly distributed over the interval $[0,1]$. The random variables $X$ and $Y$ are constructed as follows:

$$
X = \cos(2\pi U_1) \sqrt{-2\ln U_2}
$$

$$
Y = \sin(2\pi U_1) \sqrt{-2\ln U_2}
$$

Then, $X$ and $Y$ follow a normal distribution with a mean of 0 and a variance of 1.

## Derivation

Assume that $X$ and $Y$ are independent, normally distributed random variables with a mean of 0 and a variance of 1. Let $p(X)$ and $p(Y)$ be their probability density functions. Thus, we have:

$$
p(X) = \frac{1}{\sqrt{2\pi}} e^{-\frac{X^2}{2}}, \\ p(Y) = \frac{1}{\sqrt{2\pi}} e^{-\frac{Y^2}{2}}
$$

Since $X$ and $Y$ are independent, their joint probability density function is:

$$
p(X,Y) = \frac{1}{2\pi} e^{-\frac{X^2 + Y^2}{2}}
$$

We now perform a coordinate transformation by defining:

$$
X = R \cos(\theta), \\ Y = R \sin(\theta)
$$

Thus, the joint distribution in polar coordinates becomes:

$$
\int_{0}^{\infty} \int_{0}^{2\pi} \frac{1}{2\pi} e^{-\frac{R^2}{2}} R \, d\theta \, dR = 1
$$

From this, we can derive the distribution functions for $R$ and $\theta$, denoted as $P_R$ and $P_\theta$:

$$
P_R(R \leq r) = \\ \int_{0}^{r} \int_{0}^{2\pi} \frac{1}{2\pi} e^{-\frac{R^2}{2}} R \, d\theta \, dR = 1 - e^{-\frac{r^2}{2}}
$$

$$
P_\theta(\theta \leq \phi) = \\ \int_{0}^{\phi} \int_{0}^{\infty} \frac{1}{2\pi} e^{-\frac{R^2}{2}} R \, dR \, d\theta = \frac{\phi}{2\pi}
$$

It is evident that $\theta$ follows a uniform distribution over $[0, 2\pi]$. Let

$$
F_R = 1 - e^{-\frac{r^2}{2}}
$$

The inverse function is given by:

$$
R = F_R^{-1}(z) \\ = \sqrt{-2 \ln(1 - z)}
$$

When $z$ is uniformly distributed over $[0,1]$, the distribution function of $R$ is $F_R(r)$. Therefore, we can select two random variables $U_1$ and $U_2$, uniformly distributed over $[0,1]$, such that:

$$
\theta = 2\pi U_1, \\ 1 - z = U_2, \\ R = \sqrt{-2 \ln U_2}
$$

Substituting these into the earlier expressions, we obtain:

$$
X = R \cos(\theta), \\ Y = R \sin(\theta)
$$

Thus, the original expressions for $X$ and $Y$ are recovered, where they follow a normal distribution with mean 0 and variance 1.

## Verification

To generate a random variable $X \sim N(0,1)$ following a standard normal distribution, we note that the inverse of the error function $\mathrm{erf}^{-1}$ is not elementary. Using a series expansion to approximate this inverse can introduce significant computational error.

Elementary functions, on the other hand, have well-established high-precision algorithms. Thus, we aim to construct a normally distributed random variable from a uniformly distributed one using elementary functions.

Let $U, V \in (0,1)$ be uniformly distributed random variables. Consider the following transformations:

$$
X = \sqrt{-2\ln U} \cos(2\pi V), \\ Y = \sqrt{-2\ln U} \sin(2\pi V)
$$

Next, we verify that $X$ and $Y$ are independent and follow a standard normal distribution.

Define:

$$
R = \sqrt{X^2 + Y^2}, \\ \Theta = \arctan\left(\frac{Y}{X}\right)
$$

Note that $\Theta = 2\pi V$ and $R = \sqrt{-2\ln U}$. We now show that $R$ and $\Theta$, under the polar coordinate transformation, are independent and that they yield standard normal random variables.

Consider the probability density function of one variable when the other

is fixed. The probability density function of $\Theta$ is:

$$
P_\Theta(\Theta) = \frac{1}{2\pi}
$$

which is straightforward. For $R$, the cumulative distribution function is:

$$
P(R \leq R_0) = \\ P\left(U \leq \exp\left(\frac{-R_0^2}{2}\right)\right)
$$

Thus, the probability density function of $R$ is:

$$
P_R(R) = R e^{-\frac{R^2}{2}}
$$

Hence, the joint probability density of $R$ and $\Theta$ is:

$$
P(R, \Theta) = \frac{1}{2\pi} R e^{-\frac{R^2}{2}}
$$

Since the polar coordinates yield:

$$
R \, dR \, d\Theta = dX \, dY
$$

we obtain the joint probability density of $X$ and $Y$:

$$
P(X,Y) = \frac{1}{2\pi} e^{-\frac{X^2 + Y^2}{2}} \\ = \left(\frac{1}{\sqrt{2\pi}} e^{-\frac{X^2}{2}}\right) \left(\frac{1}{\sqrt{2\pi}} e^{-\frac{Y^2}{2}}\right)
$$

This is precisely the probability density function of two independent, normally distributed random variables $X$ and $Y$, each following a standard normal distribution.

## Code Implementation

```typescript
function randInterval(): number {
  return Math.random();
}

function boxMuller(mu: number, sigma: number): [number, number] {
  const u = randInterval();
  const v = randInterval();
  const x = Math.cos(2 * Math.PI * u) * Math.sqrt(-2 * Math.log(v));
  const y = Math.sin(2 * Math.PI * u) * Math.sqrt(-2 * Math.log(v));
  return [x * sigma + mu, y * sigma + mu];
}

// Usage example
const [value1, value2] = boxMuller(0, 1);
console.log(value1, value2);
// Outputs two normally distributed random numbers
```

To verify that the `boxMuller` function generates numbers that follow a normal distribution, one can generate a large set of random numbers and plot a histogram to compare with a standard normal distribution. You can use libraries such as `plotly` to visualize the histogram and evaluate the output.

Here is a test program that generates a dataset using the `boxMuller` function and performs basic statistical tests, along with plotting the histogram. Assuming you are running in a Node.js environment, you can use `plotly-nodejs` or directly use `plotly.js` in the browser to generate the graph.

First, install the required dependencies:

```bash
npm i plotly.js-dist
```

Then, execute the following code:

```typescript
import * as fs from "fs";
import * as plotly from "plotly.js-dist";

// Generate random numbers following a normal distribution using Box-Muller
function randInterval(): number {
  return Math.random();
}

function boxMuller(mu: number, sigma: number): [number, number] {
  const u = randInterval();
  const v = randInterval();
  const x = Math.cos(2 * Math.PI * u) * Math.sqrt(-2 * Math.log(v));
  const y = Math.sin(2 * Math.PI * u) * Math.sqrt(-2 * Math.log(v));
  return [x * sigma + mu, y * sigma + mu];
}

// Generate n normally distributed random numbers
function generateNormalDistribution(
  mu: number,
  sigma: number,
  n: number
): number[] {
  const values: number[] = [];
  for (let i = 0; i < n / 2; i++) {
    const [value1, value2] = boxMuller(mu, sigma);
    values.push(value1, value2);
  }
  return values;
}

// Generate the data
const mu = 0;
const sigma = 1;
const sampleSize = 10000; // Sample size
const values = generateNormalDistribution(mu, sigma, sampleSize);

// Calculate the mean and variance of the generated random numbers
const mean = values.reduce((acc, val) => acc + val, 0) / values.length;
const variance =
  values.reduce((acc, val) => acc + (val - mean) ** 2, 0) / values.length;

console.log(`Mean: ${mean}`);
console.log(`Variance: ${variance}`);

// Plot the histogram
const trace = {
  x: values,
  type: "histogram",
  xbins: {
    size: 0.1, // Bin size for the histogram
  },
  marker: {
    color: "blue",
  },
  opacity: 0.7,
  name: "Box-Muller Distribution",
};

const data = [trace];

const layout = {
  title: "Generated Normal Distribution",
  xaxis: { title: "Value" },
  yaxis: { title: "Frequency" },
  bargap: 0.05,
};

const graphOptions = {
  filename: "box-muller-distribution",
  fileopt: "overwrite",
};
fs.writeFileSync("box-muller.html", plotly.plot(data, layout, graphOptions));

console.log("Box-Muller distribution data saved as box-muller.html.");
```

This program uses the `boxMuller` function to generate 10,000 random numbers that follow a normal distribution. It calculates the mean and variance of the generated dataset. The random numbers are visualized using a histogram, allowing for a direct comparison with a standard normal distribution. The histogram is saved as a file named `box-muller.html`, which can be opened in a browser for viewing.


The Box-Muller transformation is a method to construct a random variable that follows a Gaussian distribution by using a random variable that follows a uniform distribution.

Box-Muller Transformation


## 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


## 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


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

This article mainly introduces the content and connection between the four concepts of Taylor's formula, Taylor's theorem, Taylor series and Taylor's expansion.

Taylor Formula, Taylor Theorem, Taylor Series, Taylor Expansion


## Theorem

Let $n, p \in \mathbb{Z}_{>0}$ be (strictly) positive integers.

Then:

$$
\sum_{k = 1}^n k^p = \frac{1}{p + 1} \sum_{i = 0}^p (-1)^i \binom{p + 1}{i} B_i n^{p + 1 - i} \\
= \frac{n^{p + 1}}{p + 1} - \frac{B_1 \cdot n^p}{1!} + \frac{B_2 \cdot p \cdot n^{p - 1}}{2!} + \frac{B_4 \cdot p \cdot (p - 1) \cdot (p - 2) \cdot n^{p - 3}}{4!} + \cdots
$$

where:

- $B_i$ denotes the $i$th Bernoulli number.

## Proof

Let $x \ge 0$.

$$
\sum_{k = 0}^{n - 1} e^{k x} = \sum_{k = 0}^{n - 1} \sum_{p = 0}^\infty \frac{(k x)^p}{p!}
$$

$$
= \sum_{p = 0}^\infty \left(\sum_{k = 0}^{n - 1} k^p\right) \frac{x^p}{p!}
$$

We also have:

$$
\sum_{k = 0}^{n - 1} e^{k x} = \frac{1 - e^{n x}}{1 - e^x}
$$

Sum of Geometric Sequence

$$
= \frac{e^{n x} - 1}{x} \frac{x}{e^x - 1}
$$

multiplying numerator and denominator by $x$

$$
= \frac{1}{x} \left(\sum_{p = 0}^\infty \frac{(n x)^p}{p!} - 1\right) \sum_{p = 0}^\infty \frac{B_p x^p}{p!}
$$

Definition of _Bernoulli Numbers_ and _Power Series Expansion_ for Exponential Function

$$
= \frac{1}{x} \left(\sum_{p = 1}^\infty \frac{(n x)^p}{p!}\right) \sum_{p = 0}^\infty \frac{B_p x^p}{p!}
$$

Factorial of $0$ and Zeroth Power

$$
= \frac{1}{x} \left(\sum_{p = 0}^\infty \frac{(n x)^{p + 1}}{(p + 1)!}\right) \sum_{p = 0}^\infty \frac{B_p x^p}{p!}
$$

Translation of Index Variable of Summation

$$
= \sum_{p = 0}^\infty \frac{n^{p + 1} x^p}{(p + 1)!} \sum_{p = 0}^\infty \frac{B_p x^p}{p!}
$$

Power of Product

$$
= \sum_{p = 0}^\infty \sum_{i = 0}^p \frac{n^{p + 1 - i} x^{p - i}}{(p + 1 - i)! i!} B_i n^{p + 1 - i} x^p
$$

Definition of Cauchy Product

$$
= \sum_{p = 0}^\infty \left(\frac{1}{p + 1} \sum_{i = 0}^p \binom{p + 1}{i} B_i n^{p + 1 - i}\right) \frac{x^p}{p!}
$$

multiplying by $1$ and Product of Powers

$$
\sum_{k = 0}^{n - 1} k^p = \frac{1}{p + 1} \sum_{i = 0}^p \binom{p + 1}{i} B_i n^{p + 1 - i}
$$

Definition of Binomial Coefficient

$$
\implies \sum_{k = 1}^n k^p = \frac{1}{p + 1} \sum_{i = 0}^p (-1)^i \binom{p + 1}{i} B_i n^{p + 1 - i}
$$

Equating coefficients:

$$
\sum_{k = 0}^{n - 1} k^p = \frac{1}{p + 1} \sum_{i = 0}^p \binom{p + 1}{i} B_i n^{p + 1 - i}
$$

$$
\implies \sum_{k = 0}^{n - 1} k^p + n^p = \frac{1}{p + 1} \sum_{i = 0}^p \binom{p + 1}{i} B_i n^{p + 1 - i} + \left(\frac{1}{p + 1} \binom{p + 1}{1} n^p\right)
$$

adding $n^{p}$ to both sides and Binomial Coefficient with One

$$
\implies \sum_{k = 1}^n k^p = \frac{1}{p + 1} \sum_{i = 0}^p (-1)^i \binom{p + 1}{i} B_i n^{p + 1 - i}
$$

as $B_{1}=-\frac{1}{2}$ and Odd Bernoulli Numbers Vanish

**Q.E.D.**


Faulhaber's Formula, which provides a method to compute the sum of the p-th powers of the first n positive integers.

Proof of Faulhaber's Formula


## 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.

NOTE & BLOG

Posts of MATHEMATICS