ALGORITHM

COMPUTATION

CRYPTOGRAPHY

FRONTEND

GAME-THEORY

INVESTMENT

JS&TS

LINUX

MATHEMATICS

NODE

OTHERS

PROJECT-PRACTICE

REACT

THOUGHTS

古文整理

文史哲艺


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


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

The Master Theorem For Time Complexity


## Problem Citation

A pile of stones, two people take turns to take at least 1 stone and at most 2 stones. Whoever gets the last stone loses. Is there a sure-win strategy of playing first and playing second? That's a famous **Bash Game Problem**

## Analysis

There is a pile of n items in total, and two players take turns picking up items from it. It is stipulated that each time at least one item must be taken, and at most m items can be taken. The one who takes the last item wins.

Let

$n=(m+1)q+r \quad( 0\leq r \leq m ) \\$

$\quad (\text{i})\quad$If $r=0$, the second player will win. The strategy is as follows:

If the first player takes away $k$ items, the second player takes away $m+1-k$ items. The result is $(m+1)(q-1)$ items left. By maintaining this method, the second player will win.

$\quad (\text{ii})\quad$If $r≠0$, the first player will win. The strategy is as follows:

The first player takes away $r$ items first. If the second player takes away $k$ items, the first player takes away $m+1-k$ items. The result is $(m+1)(q-1)$ items left. By maintaining this method, the first player will win.

In short, to ensure that the opponent is left with a multiple of $(m+1)$, one can eventually win.

## Extension

If it is stipulated that the one who takes the last item loses, let

$n-1=(m+1)q+r \quad( 0≤r≤m ) \\$

$\quad (\text{i})\quad$If $r=0$, the second player will win. The strategy is as follows: If the first player takes away $k$ items, the second player takes away $m+1-k$ items. The result is $(m+1)(q-1)+1$ item left. By maintaining this method, the first player will take the last item.

$\quad (\text{ii})\quad$If $r≠0$, the first player will win. The strategy is as follows: The first player takes away $r$ items first. If the second player takes away $k$ items, the first player takes away $m+1-k$ items. The result is $(m+1)(q-1)+1$ item left. By maintaining this method, the second player will take the last item.

## Mini Game

Two people take turns counting, each time reporting at least one, at most ten. Whoever reaches $100$ first wins.

Attached code:

```cpp
#include<cstdio>
#include<string.h>
#include<queue>
#include<algorithm>
#include<iostream>
using namespace std;
int main()
{
    int T;
    scanf("%d",&T);
    while(T--){
        int n,m;
        scanf("%d%d",&n,&m);
        if(n%(m+1)) printf("first\n");
        else printf("second\n");
    }
    return 0;
}
```

References:

- [1] [Summary of Bash Game](https://suool.blog.csdn.net/article/details/15777651?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-2.control&dist_request_id=ea4f2841-87a2-4c48-a15d-b529f81dffd9&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-2.control)
- [2] [Bash Game hdu 1846](https://www.cnblogs.com/pangbi/p/11509328.html)
- [3] [Bash Game - Baidu Baike](https://baike.baidu.com/item/巴什博弈/1819345?fr=aladdin)


A pile of stones, two people take turns to take at least 1 stone and at most 2 stones. Whoever gets the last stone loses.

The Bash Game


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


Magic Number - 0x5f3759df

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

Fast Inverse Square Root


余幼時即嗜學。家貧，無從致書以觀，每假借於藏書之家，手自筆錄，計日以還。天大寒，硯冰堅，手指不可屈伸，弗之怠。錄畢，走送之，不敢稍逾約。以是人多以書假余，余因得遍觀羣書。既加冠，益慕聖賢之道，又患無碩師名人與游，嘗趨百里外，從鄉之先達執經叩問。先達德隆望尊，門人弟子填其室，未嘗稍降辭色。余立侍左右，援疑質理，俯身傾耳以請；或遇其叱咄，色愈恭，禮愈至，不敢出一言以復；俟其欣悅，則又請焉。故余雖愚，卒獲有所聞。

當余之從師也，嘗負篋曳屣，行深山巨谷中。窮冬烈風，大雪深數尺，足膚皸裂而不知。至舍，四支僵勁不能動，媵人持湯沃灌，以衾擁覆，久而乃和。寓逆旅，主人日再食，無鮮肥滋味之享。同捨生皆被綺綉，戴珠纓寶飾之帽，腰白玉之環，左佩刀，右備容臭，燁然若神人。余則縕袍敝衣處其間，畧無慕豔意。以中有足樂者，不知口體之奉不若人也。蓋余之勤且艱若此。今雖耄老，未有所成，猶幸預君子之列，而承天子之寵光，綴公卿之後，日侍坐備顧問，四海亦謬稱其氏名，況才之過於余者乎？

今諸生學於太學，縣官日有廩稍之供，父母歲有裘葛之遺，無凍餒之患矣；坐大廈之下而誦《詩》《書》，無奔走之勞矣；有司業、博士為之師，未有問而不告、求而不得者也；凡所宜有之書，皆集於此，不必若余之手錄、假諸人而後見也。其業有不精、德有不成者，非天質之卑，則心不若余之專耳，豈他人之過哉。

東陽馬生君則，在太學已二年，流輩甚稱其賢。余朝京師，生以鄉人子謁余，譔長書以為贄，辭甚暢達；與之論辯，言和而色夷。自謂少時用心於學甚勞，是可謂善學者矣。其將歸見其親也，余故道為學之難以告之。謂余勉鄉人以學者，余之志也；詆我夸際遇之盛而驕鄉人者，豈知予者哉。


余朝京師，生以鄉人子謁余，譔長書以為贄，辭甚暢達；與之論辯，言和而色夷。自謂少時用心於學甚勞，是可謂善學者矣。其將歸見其親也，余故道為學之難以告之。謂余勉鄉人以學者，余之志也；詆我夸際遇之盛而驕鄉人者，豈知予者哉。

送東陽馬生序


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


Take the Glibc's Mindset for Example

How Do Computer Calculate the Log-Function


## Finding the Majority Element

Moore's Voting Algorithm is an algorithm used to find the element that appears more than half the time in an array. The basic idea of the algorithm is to find the element with a count exceeding half through the cancellation of different elements.

The algorithm steps are as follows:

1. Initialize two variables `candidate` and `count`, representing the candidate element and the counter. Initially, set `candidate` to the first element of the array and `count` to 1.
2. Traverse the array, starting from the second element:
   - If the current element is the same as `candidate`, increment the `count` by 1.
   - If the current element is different from `candidate`, decrement the `count` by 1.
     - If `count` becomes 0, set the current element as the new `candidate` and reset `count` to 1.
3. The final `candidate` is the element that may appear more than half the time.
4. Traverse the array, count the occurrences of the `candidate`, and ensure its count indeed exceeds half.

Below is an example code demonstrating how to use Moore's Voting Algorithm to find the element that appears more than half the time:

```c
#include <stdio.h>

int findMajorityElement(int arr[], int size) {
    int candidate = arr[0];
    int count = 1;

    for (int i = 1; i < size; i++) {
        if (arr[i] == candidate) {
            count++;
        } else {
            count--;
            if (count == 0) {
                candidate = arr[i];
                count = 1;
            }
        }
    }

    // Validate whether the candidate is the element appearing more than half the time
    int freq = 0;
    for (int i = 0; i < size; i++) {
        if (arr[i] == candidate) {
            freq++;
        }
    }
    if (freq > size / 2) {
        return candidate;
    } else {
        return -1;  // No element found appearing more than half the time
    }
}

int main() {
    int arr[] = {2, 4, 2, 2, 5, 2, 2, 6, 2, 2, 8};
    int size = sizeof(arr) / sizeof(arr[0]);

    int majorityElement = findMajorityElement(arr, size);

    if (majorityElement != -1) {
        printf("The element appearing more than half the time is: %d\n", majorityElement);
    } else {
        printf("No element appearing more than half the time\n");
    }

    return 0;
}
```

In the above example, for the array `{2, 4, 2, 2, 5, 2, 2, 6, 2, 2, 8}`, the element `2` appears more than half the time, so the algorithm outputs `2`.

## Extension: Finding Elements Appearing more than n-th of the Total

Moore's Voting Algorithm can be modified to find elements that appear more than n-th of the total count. The key modification lies in determining the required counting threshold.

Assuming the array length is `size`, to find elements that appear more than n-th of the total count, set the counting threshold to `size / n`. This means that when the count of an element reaches this threshold, it can be considered a valid element.

When finding elements that appear more than n-th of the total count, the following modifications can be made to Moore's Voting Algorithm:

1. Initialize two arrays `candidate` and `count`, both with a length of n-1. The `candidate` array is used to store candidate elements, and the `count` array is used to store corresponding counters. Initially, set all elements of the `candidate` array to -1 and all elements of the `count` array to 0.
2. Traverse the array for each element:
   - Check if the same candidate element exists in the `candidate` array for the current element.
     - If found, increment the corresponding counter by 1.
     - If not found, find the first position in the `candidate` array where the counter is 0. Set the current element as a new candidate in that position and set the corresponding counter to 1.
       - If no position with a counter of 0 is found, it means the `candidate` array is full. In this case, decrement all counters by 1 to cancel out one element.
3. Traverse the `candidate` array and validate whether each candidate element appears more than `size / n` times. Count the occurrences of each candidate in the original array, and if it exceeds `size / n`, output that candidate.
4. If no candidate element appears more than n-th of the total count, there is no element satisfying the condition.

The modified algorithm steps above involve traversing the array, using candidate elements and counters to cancel out, and ultimately finding elements that appear more than n-th of the total count.

Below is the modified example code:

```c
#include <stdio.h>

int findMajorityElement(int arr[], int size, int n) {
    int candidate[n - 1];
    int count[n - 1];

    for (int i = 0; i < n - 1; i++) {
        candidate[i] = -1;
        count[i] = 0;
    }

    for (int i = 0; i < size; i++) {
        int j;
        for (j = 0; j < n - 1; j++) {
            if (arr[i] == candidate[j]) {
                count[j]++;
                break;
            }
        }
        if (j == n - 1) {
            int k;
            for (k = 0; k < n - 1; k++) {
                if (count[k] == 0) {
                    candidate[k] = arr[i];
                    count[k] = 1;
                    break;
                }
            }
            if (k == n - 1) {
                for (k = 0; k < n - 1; k++) {
                    count[k]--;
                }
            }
        }
    }

    // Validate whether each candidate appears more than size/n times
    int freq = 0;
    for (int i = 0; i < n - 1; i++) {
        int count = 0;
        for (int j = 0; j < size; j++) {
            if (arr[j] == candidate[i]) {
                count++;
            }
        }
        if (count > size / n) {
            freq++;
            printf("The element appearing more than %d-th of the total is: %d\n", n, candidate[i]);
        }
    }

    if (freq == 0) {
        printf("No element appearing more than %d-th of the total\n", n);
    }

    return 0;
}

int main() {
    int arr[] = {2, 4, 2, 2, 5, 2, 2, 6, 2, 2, 8, 3, 3, 3, 3};


 int size = sizeof(arr) / sizeof(arr[0]);
    int n = 4;

    findMajorityElement(arr, size, n);

    return 0;
}
```

In the above example, for the array `{2, 4, 2, 2, 5, 2, 2, 6, 2, 2, 8, 3, 3, 3, 3}`, both elements `2` and `3` appear more than 4-th of the total count, so the algorithm outputs:

```plaintext
The element appearing more than 4-th of the total is: 2
The element appearing more than 4-th of the total is: 3
```


Moore's Voting Algorithm is an algorithm used to find the element that appears more than half the time in an array.

Moore's Voting Algorithm for Finding Majority Element


A*search algorithm is a **path finding** algorithm that finds **the** **single-pair** **shortest path** between the **start node(source)** and the **target node(destination)** of a **weighted** graph. The algorithm not only considers the actual cost from the start node to the current node(g) but also tries to estimate the cost will take from the current node to the target node using **heuristics** (h). Then it selects the node that has the **lowest f-value**(**f=g+h**) to be the next node to move until it hits the target node. **Dijkstra's algorithm** is a special case of A* algorithm where **heuristic is 0** for all nodes.

![A Star Search](https://miro.medium.com/v2/resize:fit:875/1*QvF-25wEHhQiWowGeJ8JjQ.png)

## How A\* Algorithm Works?

The formula for A\* algorithm considers two functions, **g(n) and h(n)**. Imagine that we are currently at node n, **g(n)** tells us the **actual cost** we took from the start node to the current node. **h(n)** is a **heuristic function** that **estimates** **the cost** we will take to get to the target node from the current node. Therefore, h(n) is an **educated guess** and it is crucial to determine the performance of A* algorithm. Based on the sum of the g(n) and h(n), **f(n)**, the algorithm is able to estimate the **total cost** from the start node to the end node.

**f(n) = g(n) + h(n)** , f(n): the estimate of the total cost from start node to target node through node n  
g(n): actual cost from start node to node n  
h(n): estimated cost from node n to target node

The algorithm selects the **minimum f-value node** as the next current node to explore and continuously update g(n) and h(n) if a better path is found for nodes it encountered. This process continues until the algorithm reaches the target node.

### A\* algorithm vs. Dijkstra's algorithm

**A\* Algorithm** only finds the shortest path between **the start node and the target node**, whereas **Dijkstra's algorithm** finds the shortest path from the **start node to all other nodes** because every node in a graph is a target node for it.

**A\* algorithm** runs **faster** than Dijkstra's algorithm because it uses a **heuristic** to direct **in the correct direction** towards the target node. However, **Dijkstra's algorithm expands evenly** in **all of the different available directions** because it has no knowledge regarding the target node before hand, it always processes the node that is the closest to the start node based on the cost or distance it already took.

![A* algorithm vs. Dijkstra's algorithm](https://miro.medium.com/v2/resize:fit:875/1*FzCo76pN8UJt2kIeOFlCzA.png)

## Graphical Explanation

> Heuristic is crucial to the performance of A\* algorithm. The heuristic function we use here just for the demo purpose.

![graph example](https://miro.medium.com/v2/resize:fit:875/1*5UTGfoZLPUA9dXJ_v_CsWg.png)

- **_Step 1_**  
   - **_a_**: Set the distance to the **start node** itself to **0** and the distance to **all other nodes** to **infinity**.  
   - **_b_**: Calculate the **f value** for the start node and set the **previous node** to none/nil/null.  
   - **_c_**: Initialize an **open list** and **close list** which are empty initially.

![step 1](https://miro.medium.com/v2/resize:fit:875/1*W6CaGse0EnoDi5RNW7Kb1g.png)

- **_Step 2_**: Place the **start node** into the **open list**

![step 2](https://miro.medium.com/v2/resize:fit:875/1*MJwjLSxWe_n1yW7rxkuu7w.png)

- **_Step 3  
   _\***-\* **_a_\***:\* Find the node with **the** **minimum f-value** in the **open list** and **removed** from the list. Denote this node as the **current node**.  
   - **b**: Check if the **current node** is the **target node** or not  
   - **c**: **Populate** all the current node's **neighboring nodes** and do following checks for each neighboring nodes

**1. Is the neighboring node in the close list?**

- Yes. Skip this node.
- No. Go to 2.**2. Calculate g value for the neighboring node and check if a better  
  path(lower g value) is found?** - No, Skip this node
- Yes. Go to 3.**3. Calculate g, h and f value for the neighboring node and set  
  the previous node to the current node\*\***4. Check if the neighbor node in the open list?\*\*
- Yes. Update g, h and f value of the neighboring node in the  
  open list
- No. Insert this neighboring node to the open list

- **d**: Place the **current node** to the **close list** because we have expanded this node.

![step 3](https://miro.medium.com/v2/resize:fit:875/1*ALCIZaJVV_iZ22aJDn8KwA.png)

- **_Step 4_**: Repeat _Step 3_ until reaches the target node

### demo iteration 1

![iteration 1–1](https://miro.medium.com/v2/resize:fit:875/1*3daYK311AZ5EhMZ7SlfMyA.png)

![iteration 1–2](https://miro.medium.com/v2/resize:fit:875/1*yxdTd3QM0kQiZ-Xa5Rdinw.png)

![iteration 1–3](https://miro.medium.com/v2/resize:fit:875/1*bs4TplpOIkxhzx_3TKv7-A.png)

### demo iteration 2

![iteration 2–1](https://miro.medium.com/v2/resize:fit:875/1*eX758M4Z1pHgXiVf2Re1pQ.png)

![iteration 2–2](https://miro.medium.com/v2/resize:fit:875/1*8qG-dpRBcGksbFNUql6cLg.png)

![iteration 2–3](https://miro.medium.com/v2/resize:fit:875/1*wGwIDY1BdPU4fLYd03E_tw.png)

### demo iteration 3

![iteration 3–1](https://miro.medium.com/v2/resize:fit:875/1*WEIx9GyNjStRCgMmd08ksg.png)

![iteration 3–2](https://miro.medium.com/v2/resize:fit:875/1*1o9QqFM6kjeVRD96Wzk_wg.png)

![iteration 3–3](https://miro.medium.com/v2/resize:fit:875/1*AEvC7BRbwoY7OBXwbIGk4A.png)

### demo iteration 4

![iteration 4–1](https://miro.medium.com/v2/resize:fit:875/1*nOouyfpxSOaOOP2nuL9TAA.png)

![iteration 4–2](https://miro.medium.com/v2/resize:fit:875/1*xzenBGkWstsSNEoNtGBHmQ.png)

![iteration 4–3](https://miro.medium.com/v2/resize:fit:875/1*N3nVtRAcjclt4DRBSLnkSQ.png)

### demo iteration 5

![iteration 5–1](https://miro.medium.com/v2/resize:fit:875/1*LfoBualez23-DK5xMWJ7zQ.png)

### final result

![](https://miro.medium.com/v2/resize:fit:875/1*L6s9xmzwxjWjOEe-tvTPKg.png)

At the end, we can **backtrack** the shortest path using the previous node. In this demo, we use A\* algorithm to search the shortest path **without** necessarily **search all the graph.**(We haven't examined node 4 yet, so you can see the node 4 is till on the open list). Unlike Dijkstra's algorithm, it has to expand all the nodes in the graph.

You can check [here](https://yuminlee2.medium.com/dijkstras-algorithm-6ad1e4bf9876) to know more about Dijkstra's algorithm.

## Code Implementation

### Data Structure

Use **Priority queue** to get the lowest f-value node for each run. In this demo, we implement code using **min heap** data structure .

### Complexity

Time: O(nlog(n)), Space: O(n)  
n: the total number of nodes in the graph

### Pseudocode

### initialize a min Heap contains the start node

minHeap = MinHeap([start])While minHeap is not empty:

### grab the minimum f-value node from minHeap and denote as current

current = minHeap.pop() # check if current is the target node  
 if current == target:  
 break

### populate all current node's neighbors

for neighbor in current.neighbors:  
 compute neighbor's g, h, f value

      if neighbor in minHeap:
         if neighbor.g < neighbor.g in minHeap:
             update minHeap
      else:
         insert neighbor into minHeap

### Python

- code

```python
from min_heap import MinHeap
class Node:
    def __init__(self, nodeId, hValue):
        self.id = nodeId # use node idx as id
        self.g = float("inf") # distance from the start node to the current node
        self.h = hValue # estimate distance from the current node to the target node
        self.f = float("inf") # total cost from the start node to the end node
        self.previousNode = None # use for backtrack

def AstartAlgorithm(graph, start, target, hValues):
    nodes = initializeNodes(graph, hValues)

    startNode = nodes[start]
    targetNode = nodes[target]

    # set distance to start node itself to 0
    startNode.g = 0
    startNode.f = startNode.g + startNode.h

   # init open list and close list
    openList = MinHeap([startNode])  # nodes to be expanded
    # closeList = set()  # nodes have expanded

    # repeat until the openList is empty
    while not openList.isEmpty():
        # remove the node with the lowest f-value
        currentNode = openList.remove()
        if currentNode == targetNode:
            break

        # populate all current nodes neighbors
        neighbors = graph[currentNode.id]
        for neighbor in neighbors:
            neighborIdx, distanceToNeighbor = neighbor
            neighborNode = nodes[neighborIdx]

            # check if neighbor in close list
            # if neighborNode in closeList:
            #     continue

            # check if find a better path
            newNeighborG = currentNode.g + distanceToNeighbor
            if newNeighborG >= neighborNode.g:
                continue

            # update neighbor's g, h, f and previousNode
            neighborNode.previousNode = currentNode
            neighborNode.g = newNeighborG
            neighborNode.f = neighborNode.g + neighborNode.h

            # check if the neighbor in the openList
            if not openList.containsNode(neighborNode):
                openList.insert(neighborNode)
            else:
                openList.update(neighborNode)

        # put current node to close list
        # closeList.add(currentNode)

    return backtrackPath(targetNode)

def initializeNodes(graph, hValues):
    nodes = []
    for i in range(len(graph)):
        nodes.append(Node(i, hValues[i]))
    return nodes

def backtrackPath(targetNode):
    if targetNode.previousNode is None:
        return []

    currentNode = targetNode
    path = []

    while currentNode is not None:
        path.append(currentNode)
        currentNode = currentNode.previousNode

    return path[::-1]


if __name__ == "__main__":
    hValues = [20, 16, 6, 10, 4, 0]
    graph = [
            # for vertex 0
            [
                [1, 2], [3, 6]
            ],
            # for vertex 1
            [
                [0, 2], [2, 5]
            ],
            # for vertex 2
            [
                [1, 5], [3, 7], [4, 6], [5, 9]
            ],
            # for vertex 3
            [
                [0, 6], [2, 7], [4, 10]
            ],
            # for vertex 4
            [
                [2, 6], [3, 10], [5, 6]
            ],
            # for vertex 5
            [
                [2, 9], [4, 6]
            ],
        ]

    # startNode = 1, targetNode = 5
    path = AstartAlgorithm(graph, 0, 5, hValues)
    print("The shortest path")
    print("order, nodeId, f-value")
    for idx, node in enumerate(path):
        print(idx+1, node.id, node.f)

"""
The shortest path
order, nodeId, f-value
(1, 0, 20)
(2, 1, 18)
(3, 2, 13)
(4, 5, 16)
"""
```

- min heap

```python
class MinHeap:
    def __init__(self, array):
        self.nodePositionInHeap = {node.id: idx for idx, node in enumerate(array) }
        self.heap = self.buildHeap(array)

    # Time: O(n) | Space: O(1)
    def buildHeap(self, array):
        lastParentNodeIdx = (len(array) - 2) // 2
        for currentIdx in range(lastParentNodeIdx, -1, -1):
            self.siftDown(currentIdx, len(array)- 1, array)
        return array

    # Time: O(logn) | Space: O(1)
    def remove(self):
        if self.isEmpty():
            return

        self.swap(0, len(self.heap)-1, self.heap)
        nodeToRemove = self.heap.pop()
        del self.nodePositionInHeap[nodeToRemove.id]
        self.siftDown(0, len(self.heap)-1, self.heap)
        return nodeToRemove

    # Time: O(logn) | Space: O(1)
    def insert(self, node):
        self.heap.append(node)
        self.nodePositionInHeap[node.id]= len(self.heap)-1
        self.siftUp(len(self.heap)-1, self.heap)

    # Time: O(logn) | Space: O(1)
    def siftDown(self, currentIdx, endIdx, heap):
        childOneIdx = currentIdx*2 + 1
        while childOneIdx <= endIdx:
            childTwoIdx = currentIdx*2 + 1 if currentIdx*2 + 1 <= endIdx else -1
            idxToSwap = childOneIdx
            if childTwoIdx != -1 and heap[childTwoIdx].f < heap[childOneIdx].f:
                idxToSwap = childTwoIdx

            if heap[idxToSwap].f < heap[currentIdx].f:
                self.swap(idxToSwap, currentIdx, heap)
                currentIdx = idxToSwap
                childOneIdx = currentIdx*2 + 1
            else:
                return

    # Time: O(logn) | Space: O(1)
    def siftUp(self, currentIdx, heap):
        parentIdx = (currentIdx - 1) // 2
        while currentIdx > 0 and heap[currentIdx].f < heap[parentIdx].f:
                self.swap(currentIdx, parentIdx, heap)
                currentIdx = parentIdx
                parentIdx = (currentIdx - 1) // 2

    def containsNode(self, node):
        return node.id in self.nodePositionInHeap

    def update(self, node):
        # update is occured when a slower f is found
        self.siftUp(self.nodePositionInHeap[node.id], self.heap)

    def isEmpty(self):
        return len(self.heap) == 0

    def swap(self, i, j, heap):
        self.nodePositionInHeap[heap[i].id] = j
        self.nodePositionInHeap[heap[j].id] = i
        heap[i], heap[j] = heap[j], heap[i]
```


A* search algorithm is a path finding algorithm that finds the single-pair shortest path between the start node(source) and the target node(destination) of a weighted graph.

A Star Algorithm


## Problem Introduction

Let's first consider a question: **Given an array of size 100 with elements arranged from 1 to 100 in order, how can we randomly select one number from it?**

The simplest way is to use the system method `Math.random() * 100`, which will give a random number between 0 and 99, and then use it as an index to select the corresponding element in the array.

Next, consider another question: **Given the same array, how can we randomly select 50 numbers?**

**The numbers must not repeat.** If you simply select 50 numbers randomly, the issue of duplicate numbers may arise.

One possible modification is to store each randomly generated number in an array, and for each subsequent random selection, check if the number is already in the array. If it is, continue to randomly select until 50 unique numbers are chosen.

This method works but has a problem. When the number of selections approaches the size of the array, the probability of repeated numbers increases significantly. For example, **if you need to select 99 numbers from an array of size 100**, following the above approach, as the number of selected numbers increases, the probability of selecting a duplicate number also increases, making the process less efficient.

In this case, we can consider another approach: **shuffle the array first, and then directly select the first 50 numbers.**

This approach works, but the key is to ensure that the "shuffle" is truly random. For example, a deck of cards has 54 cards, and there are $54!$ possible arrangements. A shuffle should be able to **generate any of these $54!$ arrangements with equal probability**.

## Shuffle Algorithm

The **Fisher-Yates shuffle** algorithm was proposed by Ronald Fisher and Frank Yates in 1938 and later adapted for computer programming by Richard Durstenfeld in 1964.

The algorithm is very simple: starting from the last element of the array, swap it with any one of the previous elements (including itself). Then, swap the second-to-last element with one of the remaining elements, and so on.

This process ensures that every element has an **equal probability** of appearing in any position.

## Program Implementation

```typescript
function knuthShuffle<T = any>(array: T[]): T[] {
  const result = array.slice();
  for (let i = result.length - 1; i > 0; i--) {
    const j = Math.floor(Math.random() * (i + 1));
    [result[i], result[j]] = [result[j], result[i]];
  }
  return result;
}

const arr = [1, 2, 3, 4, 5];
console.log(knuthShuffle(arr));
```

## Proof

The Knuth Shuffle algorithm (also known as the Fisher-Yates Shuffle) can generate all possible permutations with equal probability. To prove that the algorithm is effective in shuffling any sequence, the key is to ensure that **each permutation appears with equal probability**. Below is the mathematical proof.

For a sequence of $n$ elements, starting from the last element, randomly select one of the previous elements (including the current one) to swap. The steps are as follows:

1. For the $i$-th element (starting from $i = n$ down to $i = 2$), randomly select an index $j$ from the range [0, $i-1$] and swap element $a[i]$ with element $a[j]$.
2. Finally, all elements in the sequence are swapped once.

### Objective

We need to prove that this algorithm generates all possible permutations with equal probability. For a sequence of length $n$, there should be $n!$ different permutations, and each permutation should have a probability of $\frac{1}{n!}$.

### Analysis

1. **Total number of permutations**  
   For an array of length $n$, the total number of possible permutations is $n!$. That is, for $n$ elements, the number of unique permutations is:

   $$
   P(n) = n! = n \times (n-1) \times \cdots \times 2 \times 1
   $$

2. **Random selection process in the algorithm**  
   In the $i$-th step, the Knuth Shuffle algorithm randomly selects an index between 0 and $i$ to swap. This selection is uniform, so the probability of selecting any index is $\frac{1}{i+1}$.

   - Step 1: Choose one element out of $n$, with probability $\frac{1}{n}$.
   - Step 2: Choose one of the remaining $n-1$ elements, with probability $\frac{1}{n-1}$.
   - And so on, until no elements remain to swap.

3. **Probability of generating a specific permutation**  
   To prove that Knuth Shuffle generates all permutations with equal probability, we need to show that the probability of generating any specific permutation is $\frac{1}{n!}$.

   The probability of any element being swapped into a certain position is independent for each step:

   - The $n$-th element can be swapped with any of the $n$ elements, with probability $\frac{1}{n}$.
   - The $(n-1)$-th element can be swapped with any of the remaining $n-1$ elements, with probability $\frac{1}{n-1}$.
   - And so on.

   Therefore, the probability of generating any particular permutation is:

   $$
   \frac{1}{n} \times \frac{1}{n-1} \times \cdots \times \frac{1}{1} = \frac{1}{n!}
   $$

4. **Equal probability for all permutations**  
   Since the probability of generating each permutation is $\frac{1}{n!}$ and there are $n!$ distinct permutations, the Knuth Shuffle generates all possible permutations with equal probability.

In conclusion, the Knuth Shuffle algorithm ensures that every element is swapped with equal probability in each step, and the resulting permutations are uniformly distributed across all $n!$ possible permutations. Therefore, the Knuth Shuffle is an **effective** algorithm for shuffling any sequence in a uniformly random manner.


The Knuth Shuffle algorithm (also known as the Fisher-Yates Shuffle) generates all possible permutations with equal probability.

The Knuth Shuffle Algorithm


The Rabin-Karp algorithm, also known as the Karp-Rabin algorithm, was introduced by _Richard M. Karp_ and _Michael O. Rabin_ in 1987. It is designed to address the multiple pattern string matching problem.

Its implementation is somewhat unconventional. It begins by computing the hash values of two strings and then determines whether there is a match by comparing these hash values.

## Algorithm Analysis and Implementation

Choosing an appropriate hash function is crucial. Assuming the text string is $t[0, n)$, and the pattern string is $p[0, m)$, where $0<m<n$, let $Hash(t[i,j])$ represent the hash value of the substring $t[i, j]$.

When $Hash(t[0, m-1])!=Hash(p[0,m-1])$, it is natural to compare $Hash(t[1, m])$. In this process, if we recalculate the hash value for the substring $t[1, m]$, it would require $O(n)$ time complexity, which is not cost-effective. Observing that there are m-1 overlapping characters between the substrings $t[0, m-1]$ and $t[1, m]$, we can use a rolling hash function. This reduces the time complexity of recalculation to $O(1)$.

The rolling hash function used in the Rabin-Karp algorithm primarily leverages the concept of [Rabin fingerprint](https://en.wikipedia.org/wiki/Rabin_fingerprint). For example, the formula to calculate the hash value of the substring $t[0, m-1]$ is as follows:

$$
Hash(t[0, m-1])=t[0]*b^{m-1}+t[1]*b^{m-2}+...+t[m-1]*b^0
\\
\tag{t[0] represents the ASCII code of the character}
$$

Here, $b$ is a constant. In Rabin-Karp, it is generally set to 256 because the maximum value of a character does not exceed 255. The formula above has an issue - hash values could overflow. To address this, we take the modulus, and the value should be as large as possible and preferably a prime number. Here, we take 101.

The formula to calculate the hash value of the substring $t[1, m]$ is then:

$$
Hash(t[1,m])=(Hash(t[0,m-1])-t[0]*b^{m-1})*b+t[m]*b^0\\\tag{Please compare with the previous formula}
$$

The complete code is as follows:

```cpp
#include <iostream>
#include <string.h>

using namespace std;

#define BASE 256
#define MODULUS 101

void RabinKarp(char t[], char p[])
{
    int t_len = strlen(t);
    int p_len = strlen(p);

    // For rolling hash
    int h = 1;
    for (int i = 0; i < p_len - 1; i++)
        h = (h * BASE) % MODULUS;

    int t_hash = 0;
    int p_hash = 0;
    for (int i = 0; i < p_len; i++)
    {
        t_hash = (BASE * t_hash + t[i]) % MODULUS;
        p_hash = (BASE * p_hash + p[i]) % MODULUS;
    }

    int i = 0;
    while (i <= t_len - p_len)
    {
        // Considering the possibility of hash collisions, we use memcmp for additional verification
        if (t_hash == p_hash && memcmp(p, t + i, p_len) == 0)
            cout << p << " is found at index " << i << endl;

        // Rolling hash
        t_hash = (BASE * (t_hash - t[i] * h) + t[i + p_len]) % MODULUS;

        // Avoiding negative values
        if (t_hash < 0)
            t_hash = t_hash + MODULUS;

        i++;
    }
}

int main()
{
    char t[100] = "It is a test, but not just a test";
    char p[10] = "test";

    RabinKarp(t, p);

    return 0;
}
```

The output is as follows:

```text
test is found at index 8
test is found at index 29
```

## Complexity Analysis

Let's examine the space complexity first, which is easily determined: $S(n)=O(1)$.

Now, consider the time complexity. Let the length of the text string be n and the pattern string be m. Preprocessing requires $O(m)$, and during matching, in the best case where there are no hash collisions, $T_{best}(n)=O(n-m)$. In the worst case, where there is a collision every time, $T_{worst}(n)=O((n-m)*m)$. In practical scenarios, n is often much larger than m, so the final complexity table is:

| $S_{n}$        | $O(1)$  |
| -------------- | ------- |
| $T_{best}(n)$  | $O(n)$  |
| $T_{worst}(n)$ | $O(mn)$ |

## Application Analysis

The primary application of the Rabin-Karp algorithm is in plagiarism detection for articles, such as the detection system used by [Semantic Scholar](https://www.semanticscholar.org/).

However, from the complexity data above, the Rabin-Karp algorithm does not seem to have a significant advantage. Is it practical for detecting text plagiarism? Feedback from actual usage results indicates that the time complexity for plagiarism detection is only $O(n)$. I believe this is mainly due to the following two points:

1. In real-life articles, the text data does not often exhibit as many hash collisions as we might imagine.
2. The original content in a submitted article is likely to be much larger than the plagiarized content. In other words, successful matches do not occur as frequently as we might imagine.

## References

- [Rabin–Karp algorithm](https://en.wikipedia.org/wiki/Rabin–Karp_algorithm)
- [Searching for Patterns | Set 3 (Rabin-Karp Algorithm)](https://www.geeksforgeeks.org/searching-for-patterns-set-3-rabin-karp-algorithm/)
- [Computer Algorithms: Rabin-Karp String Searching](http://www.stoimen.com/blog/2012/04/02/computer-algorithms-rabin-karp-string-searching/)


It is designed to address the multiple pattern string matching problem.

Rabin-Karp Algorithm


## Problem Description

The Josephus problem is a classic mathematical problem that describes the following scenario:

There are $n$ people standing in a circle. Starting from a certain person, they count off in sequence, and every $m$-th person steps out of the circle. The counting and elimination process continues, looping back to the next person when the end of the circle is reached. This continues until only one person remains. The problem is to determine the initial position of the last person remaining in the circle.

For example, when $n=7$ and $m=3$, the order of elimination is: $3, 6, 2, 7, 5, 1$. The last person remaining is at position $4$.

## Solution

The Josephus problem can be solved using either recursion or mathematical formulas. The recursive solution assumes knowledge of the solution for $n-1$ people, denoted as $f(n-1, m)$, and then calculates the solution for $n$ people using the formula:

$
f(n, m) = (f(n-1, m) + m) \% n
$

Here, $\%$ denotes the modulo operation. The base case is $f(1, m) = 0$, representing the solution when there is only one person.

The mathematical formula solution is derived based on the recurrence relation, resulting in the same formula:

$
f(n, m) = (f(n-1, m) + m) \% n
$

## Programming Implementation

When it comes to solving the Josephus problem, dynamic programming is an effective approach. Here is an example C code implementing dynamic programming to solve the Josephus problem:

```c
#include <stdio.h>

int josephus(int n, int m) {
    int dp[n + 1];
    dp[1] = 0; // Initial position when there is only one person

    for (int i = 2; i <= n; i++) {
        dp[i] = (dp[i - 1] + m) % i;
    }

    return dp[n];
}

int main() {
    int n = 7; // Total number of people
    int m = 3; // Count off to m for elimination

    int lastPerson = josephus(n, m);
    printf("The last person remaining is: %d\n", lastPerson + 1); // Adding 1 due to 0-based indexing

    return 0;
}
```

In the above code, we use a dynamic programming array `dp` to store the next position after each person is eliminated. Through iteration, we can find the initial position of the last person remaining in the circle. In the `main` function, we set the total number of people `n` and the count to `m` for elimination. Then, we call the `josephus` function to calculate the position of the last person remaining and print the result.

Note that we assume people's positions start from 1, while array indices start from 0. Therefore, when printing the position of the last person remaining, we need to add 1.

Using dynamic programming provides an efficient solution to the Josephus problem, avoiding the repetitive calculations of recursion and improving code performance.


A description, analysis and solution for joseph circle problem.

MY BLOG

ALL POSTS