
British economist John Maynard Keynes and Austrian economist Friedrich Hayek are two of the most renowned figures in the history of 20th-century economic thought. Their ideas not only had a profound impact on the past economic and social landscape, but they will continue to influence current and future economic and social development until a more convincing economic ideology takes hold. The debate between Keynes and Hayek is known as “the most classic and famous duel in economic history.” American journalist Nicholas Wapshott, who was a senior editor at *The Times* of London and *The New York Sun*, systematically examined this debate in his book *Keynes Hayek: The Clash that Defined Modern Economics*. As we revisit this classic academic feast, we also reflect on the serious economic question: “To what extent should the government intervene in the market?”

Today, people are already quite familiar with the starkly opposing economic views of Keynes and Hayek. Keynes believed that the root cause of an economic crisis is the lack of confidence among capitalists in investment prospects, leading to reduced investment, and thus advocated for effective government intervention and control over social economic activities. Hayek, on the other hand, believed that the root cause of economic crises was over-investment by capitalists and argued that government intervention would only damage the mechanisms of market development and functioning. Therefore, he firmly defended the idea of a free market order. Keynes argued that the government could expand public spending through fiscal deficits and national debt, increase state investment and consumption, and indirectly advocated for the "nationalization of currency." Hayek, by contrast, proposed that currency issuance should be handed over to private banks, which would provide sound and stable money based on their own interests, advocating for the "denationalization of currency."

Keynes supported inflationary policies to increase corporate profits, stimulate investment, expand employment, and achieve stable economic development. Hayek argued that inflation would widen the income gap, cause asset prices to soar, and eventually lead to greater unemployment and economic vicious cycles. He believed that business cycles would naturally return to full employment. Keynes argued for increased investment in public works, improving social welfare, and providing relief spending for those living below the poverty line. Hayek, however, believed that public goods do not directly satisfy any specific demand but provide certain conditions under which individuals and smaller groups can enjoy better opportunities to meet their respective needs.

The debate between Keynes and Hayek is characterized by a vast ideological chasm, and throughout their lives, they could never reconcile with each other’s viewpoints. Their debates were later carried on by their passionate disciples, such as John Kenneth Galbraith and Milton Friedman. In terms of personality, the two were also very different. Keynes was highly articulate and charismatic, often painting an optimistic vision of the world, which was embraced by an entire generation of politicians and economists on both sides of the Atlantic. In contrast, Hayek was meticulous and persistent, swimming against the tide, but found support among market advocates and libertarians. To this day, the two distinct schools of thought in economics continue their “battle.”

This economic battle has left government decision-makers oscillating between the two scholars' ideas, influencing the ups and downs of the world economy, affecting the lives and livelihoods of billions. The global economic crisis of 1929–1933 led to the widespread adoption of Keynesianism, represented by the "New Deal" of Franklin D. Roosevelt, and it seemed that Keynes had triumphed over Hayek. However, the oil crisis of the 1970s plunged the global economy into the dilemma of "stagflation," and the neoliberalism represented by the "Washington Consensus" became the dominant policy-making ideology, giving Hayek a comeback. By the end of the 20th century, the failure of neoliberal economic policies in Latin America, Russia, Southeast Asia, and Eastern Europe, coupled with the impact of the 2008 global financial crisis, suggested that Keynesianism might once again be on the rise.

In terms of policy practice, the struggle between Keynesianism and neoliberalism may continue, and it is still too early to conclude who will ultimately prevail. Interestingly, the Nobel Prize in Economics, which is seen as a "weather vane" for academic trends, has also been swayed by economic crises and shifts in government decision-making. Since the Nobel Prize in Economics was first awarded in 1968, Keynes, who had passed away by then, was naturally ineligible. Hayek, however, won the prize in 1974. During the 1970s and 1980s, when neoliberalism was on the rise, figures like Milton Friedman, Robert Lucas, Finn Kydland, and Edward Prescott also won the prize. Later, however, new Keynesians like Joseph Stiglitz and Peter Diamond, who carried on Keynes' legacy, were awarded the prize as well. In 2012, with the resurgence of neoliberalism, Alvin Roth and Lloyd Shapley took home the Nobel Prize in Economics.

In fact, despite their dramatic disagreements, Keynes and Hayek had similar aims. They both studied the same subject: the cycles of boom and bust in commercial activity that emerged from the ruins of World War I. Both sought to establish a theoretical framework that combined real and monetary phenomena to explain the dynamics of the economy, especially the business cycle. Their economic theories both originated from Knut Wicksell's theory of monetary equilibrium and the concept of "natural interest rates." Furthermore, they shared a deep commitment to the basic principles of a market economy governed by freedom, democracy, and the rule of law. If we look purely at logical reasoning, both of their theoretical conclusions and policy recommendations could explain the causes of economic recessions and offer solutions. The historical choices made by policymakers ultimately explain the root cause of their disagreements—it was a matter of experience and practical judgment.

It’s worth mentioning that although their debate was often heated, with sharp criticisms exchanged, the mutual respect they maintained beyond academia, and their personal friendship, is even more commendable. Although Hayek had differing views on credit from Keynes, he praised Keynes' *A Treatise on Money* as "an outstanding description and classification of different forms of money." Similarly, Keynes was deeply moved by Hayek's great work *The Road to Serfdom*. Keynes not only provided housing for Hayek when he first arrived in London, but in their debate letters, they openly discussed their friendship. Hayek once said, “Although I still disagreed with Keynes' views and had fierce debates with him, we maintained the best of personal friendships. And as a man, I held him in the highest esteem in many ways.”

Their ideas not only had a profound impact on the past economic and social landscape, but they will continue to influence current and future economic and social development until a more convincing economic ideology takes hold.

The Historical Debate Between Hayek and Keynes


## The Largest Ear of Wheat

It is said that Socrates once led his students for a walk in a wheat field. His students asked him: how can one find the most suitable life partner?

Socrates looked at the wheat field and gave them a task: "Go into the wheat field and find the largest ear of wheat, but as you walk, you can only move forward and cannot turn back."

The students began walking while choosing, and the field was full of large ears of wheat. But which one was the largest? They walked with their heads down, checking one ear of wheat after another. They shook their heads at this one; they shook their heads at that one. They always believed that the largest ear was still ahead. Although they picked a few, they were not satisfied and threw them away. They thought there were plenty of chances and no need to decide too early.

The students continued walking with their heads down, carefully selecting as they went, for quite a long time. Suddenly, they heard Socrates' deep, resounding voice: "You have reached the end." At that moment, the students, empty-handed, realized what had happened and felt disheartened.

Socrates said to them: "There is certainly a largest ear of wheat in this field, but you may not encounter it; and even if you do, you may not make the correct judgment. Therefore, the largest ear is the one you just picked."

This story tells us that life is like walking in a wheat field, also searching for the largest ear of wheat. Some people, upon seeing a full and ripe ear, seize the opportunity and pick it without hesitation; others keep looking around and miss their chances again and again. Of course, the goal should be the largest ear, but holding the one in hand is the most tangible reality.

**So, the question arises: Is there a method that can help us find the largest ear of wheat as efficiently as possible?**

## Treasure Box Game

Let’s consider another example. Imagine there is a gambling game with the following rules: There are $m$ treasure boxes, each containing a random prize amount between $0$ and $x$ dollars, uniformly distributed. **However, you don’t know what the maximum amount $x$ in the boxes is.** Once you open a treasure box and see the prize, you have two choices: either take the prize and stop, or discard the box and open another one. However, once discarded, a box cannot be chosen again.

Now the question arises: how should you strategize to get the most money possible? Naturally, the more money, the better.

What would happen if we didn’t use any strategy and randomly picked a treasure box?

Let’s assume the amounts in the treasure boxes range from $1000, 999, 998, ..., 0$, with the top prize being $1000$. What is the probability of drawing the top prize?

Since each specific amount (including the top prize of $1000) has an equal chance of appearing, and the prize is uniformly distributed, the probability of drawing the top prize can be calculated as follows:

$$
P_{Jackpot} = \frac{1}{(\text{MaxMoney} - \text{MinMoney} + 1)}
$$

In this case, the maximum prize is $1000$ and the minimum is $0$, so:

$$
P_{Jackpot} = \frac{1}{(1000 - 0 + 1)} = \frac{1}{1001}
$$

Therefore, the probability of randomly picking a treasure box containing the top prize of $1000$ is $0.0999\% ≈ 0.1\%$. The odds are very low indeed.

Since chasing the jackpot is unrealistic, let’s lower our expectations a bit. What if we aim for a prize of at least $800$? What is the probability of winning that?

We can use the formula for uniform distribution. For a prize uniformly distributed between $0$ and $1000$, the probability density function is constant. The formula for probability is:

$$
P(a \leq X \leq b) = \frac{b - a}{\text{range}}
$$

In this case, the range is from $0$ to $1000$, so the total range is $1000$. We want to calculate the probability for the interval from $800$ to $1000$, so:

$$
P(800 \leq X \leq 1000) = \frac{1000 - 800}{1000} = \frac{200}{1000} = 0.2
$$

Thus, the probability of randomly selecting a treasure box with a prize between $800$ and $1000$ is $20\%$.

While the probability has increased significantly, $20\%$ is still relatively low.

So, if you randomly pick a treasure box and decide to keep it, the chances of getting a high prize are still slim, which clearly isn’t the optimal strategy.

## Optimization

**In most cases, we don’t know the maximum prize inside the treasure boxes.** Just like in the story of Socrates' students picking the largest ear of wheat, they didn’t know how big the largest ear was.

Let’s change our approach. Even though we don’t know the maximum prize, we do know that if we sample a portion of the treasure boxes, the probability of finding the jackpot in the sample is equal to its probability in the entire set. By sampling some treasure boxes, we can estimate the value of the jackpot based on the maximum prize in the sample.

Let’s denote the jackpot in the entire set as $M$. Since the jackpot is uniformly distributed among the treasure boxes, we can take a sample from the boxes and record the maximum prize in that sample, denoted as $M'$. Clearly, this sample maximum $M'$ will be close to the overall jackpot $M$.

Sampling a portion of the boxes to estimate the jackpot $M$ is feasible, and the more samples we take, the closer $M'$ will be to $M$.

However, the sample size cannot be too small, or the difference between $M'$ and $M$ will be large, making the estimate ineffective. On the other hand, the sample size cannot be too large either, because there will be fewer treasure boxes left to choose from, reducing the chance of finding the jackpot. Consider an extreme example: if we sample $999$ treasure boxes out of $1000$, $M'$ will likely be almost equal to $M$, but with only one treasure box remaining, the chance of hitting the jackpot becomes negligible.

Thus, it is important to balance the sample size with the total number of treasure boxes. Is there an optimal sample size where $M'$ is close to $M$ and the probability of winning a high prize is maximized afterward?

Next, we will use mathematical methods to derive this “balance point.” Suppose the total number of treasure boxes is $m$, and the sample size is $r$. We will find the expression for the optimal $r$.

## Mathematical Derivation

Assume we have a total of $N$ treasure boxes, and we select the first $r$ boxes as a sample. In this scenario, we aim to find the optimal sample size $r$ that maximizes the probability of selecting the highest prize after discarding the sample.

### Probability Derivation

For a fixed $r$, suppose the highest prize appears in the $k$-th treasure box. To win the highest prize after discarding the sample, it is necessary that the second-highest prize appears in the first $k-1$ boxes, and specifically, this second-highest prize must also appear within the sample (the first $r$ boxes). The probability of this happening is $\frac{r}{k-1}$.

This can be visualized with a diagram, where the red arrow indicates the positions where the second-highest prize occurs:

![](https://raw.githubusercontent.com/PrinOrange/pic-host/master/myblog/20240922115506.webp)

Thus, we derive the following equation:

$$
P = \frac{1}{N} \sum_{k=r+1}^{N} \frac{r}{k-1} = \frac{r}{N} \sum_{k=r+1}^{N} \frac{1}{k-1}
$$

In this equation, $\frac{1}{N}$ represents the probability of the highest prize appearing in the $k$-th treasure box, and $\frac{r}{k-1}$ represents the probability that the second-highest prize is within the sample of the first $r$ boxes. $P$ is the probability of winning the highest prize after discarding the sample.

### Summation Approximation

The above equation contains a harmonic series, which does not have a simple closed form but can be approximated using the following formula:

$$
\sum_{k=1}^{N} \frac{1}{k} = \ln(N) + \gamma
$$

where $\ln(N)$ is the natural logarithm and $\gamma$ is Euler's constant. Here, we want to approximate the summation from $r+1$ to $N$. To do this, we rewrite the sum as:

$$
\sum_{k=r+1}^{N} \frac{1}{k-1} = \sum_{k=1}^{N-(r+1)+1} \frac{1}{k} = \sum_{k=1}^{N-r} \frac{1}{k}
$$

Thus, the sum starts from $k=1$ and follows the same harmonic series approximation:

$$
\sum_{k=r+1}^{N} \frac{1}{k-1} = \left( \ln(N) + \gamma \right) - \left( \ln(r) + \gamma \right) = \ln(N) - \ln(r)
$$

Finally, using logarithmic properties, we simplify this to:

$$
\ln(N) - \ln(r) = \ln\left(\frac{N}{r}\right)
$$

### Substituting into the Original Equation

Now, we substitute the simplified sum into the original probability equation:

$$
P = \frac{r}{N} \ln\left(\frac{N}{r}\right)
$$

Next, we introduce a new variable $x = \frac{r}{N}$ to represent the proportion of the total boxes chosen as the sample, yielding:

$$
P = x \ln\left(\frac{1}{x}\right)
$$

### Maximizing the Probability

To maximize $P$, we take the derivative of $P = x \ln\left(\frac{1}{x}\right)$ with respect to $x$:

$$
\frac{dP}{dx} = \frac{d}{dx} \left(x \ln\left(\frac{1}{x}\right)\right)
$$

Using the product rule:

$$
\frac{d}{dx}(uv) = u \frac{dv}{dx} + v \frac{du}{dx}
$$

where $u = x$ and $v = \ln\left(\frac{1}{x}\right)$, we calculate the derivatives:

$$
\frac{du}{dx} = 1
$$

and

$$
\frac{dv}{dx} = \frac{d}{dx} \ln\left(\frac{1}{x}\right) = \frac{d}{dx} \left(-\ln(x)\right) = -\frac{1}{x}
$$

Substituting these into the product rule:

$$
\frac{dP}{dx} = x \left(-\frac{1}{x}\right) + \ln\left(\frac{1}{x}\right) \cdot 1 = -1 + \ln\left(\frac{1}{x}\right)
$$

To find the maximum probability $P$, we set the derivative equal to zero:

$$
-1 + \ln\left(\frac{1}{x}\right) = 0
$$

Thus, we have:

$$
\ln\left(\frac{1}{x}\right) = 1
$$

Exponentiating both sides, we get:

$$
\frac{1}{x} = e
$$

Therefore, $x = \frac{1}{e}$, which means the probability is maximized when $\frac{r}{N} = \frac{1}{e}$.

### Maximum Probability

At this optimal point, the maximum probability is:

$$
P = P_{max} = \frac{1}{e} \approx 0.368
$$

This shows that the optimal strategy is to sample approximately $\frac{1}{e}$ (about 36.8%) of the total treasure boxes and then select the next box that is greater than all the ones in the sample. This strategy maximizes the chance of selecting the highest prize.

## Program Simulation

We can use TypeScript to simulate the process described above.

To verify the probability of winning the highest prize using this strategy, we can calculate the ratio of successful trials where the highest prize is obtained to the total number of trials. Below is a TypeScript program that simulates this process:

```typescript
function simulateGamblingGame(m: number, x: number, trials: number): number {
  let successCount = 0;

  for (let trial = 0; trial < trials; trial++) {
    const boxes = Array.from({ length: m }, () => Math.random() * x);
    const firstPart = boxes.slice(0, Math.floor(m * 0.37));

    // Remember the highest prize in the first 37% of boxes
    const p = Math.max(...firstPart);

    // Find the maximum prize in all boxes
    const maxPrize = Math.max(...boxes);

    // Check from the 37% onward if a prize higher than p is found
    for (const prize of boxes.slice(Math.floor(m * 0.37))) {
      if (prize > p) {
        if (prize === maxPrize) {
          successCount++;
        }
        break; // Stop once a higher prize is found
      }
    }
  }

  return successCount / trials; // Return the success probability
}

// Set parameters
const m = 100; // Total number of boxes
const x = 1000; // Maximum prize in each box
const trials = 100000; // Number of trials

const successProbability = simulateGamblingGame(m, x, trials);
console.log(`Probability of winning the highest prize: ${successProbability}`);
```

### Explanation

1. **Parameters**:
   - `m`: Total number of treasure boxes.
   - `x`: Maximum prize in each box.
   - `trials`: Number of simulation trials. More trials increase accuracy.
2. **Simulation Process**:

   - Randomly generate the prize for each box.
   - Calculate the maximum prize in the first 37% of the boxes (`p`).
   - Find the highest prize (`maxPrize`) in all boxes.
   - Starting from the 38% box, check if there is a prize higher than `p` and if it equals `maxPrize`.

3. **Result**:
   - Return the probability of winning the highest prize after many trials.

Running this simulation produces a result like this:

```plaintext
Probability of winning the highest prize: 0.367999
```

As shown, the probability of winning the highest prize increases to around 37%, which is a significant improvement from 0.1%.

---

To further verify, let's simulate the probability of winning the top $t\%$ prizes by dividing the prize ranges into $10\%$ increments (e.g., top 10%, top 20%, etc.).

### Modified TypeScript Code for Different Prize Ranges

```typescript
function simulateGamblingGameWithRanges(
  m: number,
  x: number,
  trials: number
): number[] {
  const successCounts = Array(10).fill(0);
  // Array to store the number of successes in each range

  for (let trial = 0; trial < trials; trial++) {
    const boxes = Array.from({ length: m }, () => Math.random() * x);

    const k = Math.floor(m * 0.37);
    const firstPart = boxes.slice(0, k);

    const p = Math.max(...firstPart);

    let chosenPrize = 0;

    // Find the first box greater than p after the first 37% boxes
    for (const prize of boxes.slice(k)) {
      if (prize > p) {
        chosenPrize = prize;
        break;
      }
    }

    // If no box greater than p is found, select the last box
    if (chosenPrize === 0) {
      chosenPrize = boxes[boxes.length - 1];
    }

    const sortedBoxes = [...boxes].sort((a, b) => b - a);

    // Determine the prize range that the chosen prize falls into
    for (let i = 1; i <= 9; i++) {
      const thresholdPrize = sortedBoxes[Math.floor((m * i * 10) / 100)];
      if (chosenPrize >= thresholdPrize) {
        successCounts[i - 1]++;
      }
    }

    // The top 100% prize is always chosen, so increment this count
    successCounts[9]++;
  }

  // Convert counts to probabilities
  const successProbabilities = successCounts.map((count) => count / trials);

  return successProbabilities;
}

// Set parameters
const m = 100; // Number of boxes
const x = 1000; // Maximum prize
const trials = 100000; // Number of trials

const successProbabilities = simulateGamblingGameWithRanges(m, x, trials);

// Output the probabilities for each prize range
for (let i = 1; i <= 10; i++) {
  console.log(
    `Probability of winning a prize in the top ${i * 10}%: ${(
      successProbabilities[i - 1] * 100
    ).toFixed(2)}%`
  );
}
```

### Results from the Program

```plaintext
Probability of winning a prize in the top 10%: 66.96%
Probability of winning a prize in the top 20%: 70.67%
Probability of winning a prize in the top 30%: 74.42%
Probability of winning a prize in the top 40%: 78.06%
Probability of winning a prize in the top 50%: 81.70%
Probability of winning a prize in the top 60%: 85.45%
Probability of winning a prize in the top 70%: 89.14%
Probability of winning a prize in the top 80%: 92.84%
Probability of winning a prize in the top 90%: 96.60%
Probability of winning a prize in the top 100%: 100.00%
```

We see that, using this strategy, we maximize the chances of getting a high-value prize.

## Conclusion

The simulation verifies that selecting 37% of the boxes first and then choosing the next box with a prize greater than the maximum from the first 37% is the optimal strategy. This maximizes the probability of winning the highest prize, known as the **37% rule**.

This strategy can also be applied in real-life scenarios, such as making optimal choices in uncertain situations.


This article proposes the classic 37% rule in game theory and its mathematical derivation process.

Socrates: The Story of Searching for the Largest Ear of Wheat

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

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

NOTE & BLOG

✨秋刀鱼之愿

LATEST