Copyright © 2014 - 2022 by https://github.com/cp-algorithms +repo_name: cp-algorithms/cp-algorithms +edit_uri: edit/main/src/ +copyright: Text is available under the Creative Commons Attribution Share Alike 4.0 International License
Copyright © 2014 - 2025 by cp-algorithms contributors extra_javascript: - javascript/config.js - - https://polyfill.io/v3/polyfill.min.js?features=es6 - - https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js + - https://cdnjs.cloudflare.com/polyfill/v3/polyfill.min.js?features=es6 + - https://unpkg.com/mathjax@3/es5/tex-mml-chtml.js extra_css: - stylesheets/extra.css @@ -47,17 +49,20 @@ markdown_extensions: alternate_style: true - attr_list - pymdownx.emoji: - emoji_index: !!python/name:materialx.emoji.twemoji - emoji_generator: !!python/name:materialx.emoji.to_svg + emoji_index: !!python/name:material.extensions.emoji.twemoji + emoji_generator: !!python/name:material.extensions.emoji.to_svg - meta + - toc: + permalink: true plugins: + - toggle-sidebar: + toggle_button: all - mkdocs-simple-hooks: hooks: on_env: "hooks:on_env" - search - - tags: - tags_file: tags.md + - tags - literate-nav: nav_file: navigation.md - git-revision-date-localized: @@ -69,9 +74,11 @@ plugins: docs_path: src/ token: !ENV MKDOCS_GIT_COMMITTERS_APIKEY enabled: !ENV [MKDOCS_ENABLE_GIT_COMMITTERS, False] + branch: !ENV [MKDOCS_GIT_COMMITERS_BRANCH, main] - macros + - rss extra: analytics: provider: google - property: UA-85220282-1 + property: G-7FLS2HCYHH diff --git a/scripts/install-mkdocs.sh b/scripts/install-mkdocs.sh index d7ccb735f..9c4e73c3f 100755 --- a/scripts/install-mkdocs.sh +++ b/scripts/install-mkdocs.sh @@ -2,9 +2,11 @@ pip install \ "mkdocs-material>=9.0.2" \ + mkdocs-toggle-sidebar-plugin \ mkdocs-macros-plugin \ mkdocs-literate-nav \ mkdocs-git-authors-plugin \ mkdocs-git-revision-date-localized-plugin \ mkdocs-simple-hooks \ + mkdocs-rss-plugin \ plugins/mkdocs-git-committers-plugin-2 diff --git a/src/algebra/balanced-ternary.md b/src/algebra/balanced-ternary.md index 6e9cdfa73..12e4d4ddc 100644 --- a/src/algebra/balanced-ternary.md +++ b/src/algebra/balanced-ternary.md @@ -6,7 +6,7 @@ e_maxx_link: balanced_ternary # Balanced Ternary - + This is a non-standard but still positional **numeral system**. Its feature is that digits can have one of the values `-1`, `0` and `1`. Nevertheless, its base is still `3` (because there are three possible values). Since it is not convenient to write `-1` as a digit, diff --git a/src/algebra/big-integer.md b/src/algebra/big-integer.md index 38a4ef77d..7f8dd1816 100644 --- a/src/algebra/big-integer.md +++ b/src/algebra/big-integer.md @@ -166,7 +166,7 @@ This method is often used for calculations modulo non-prime number M; in this ca The idea is to choose a set of prime numbers (typically they are small enough to fit into standard integer data type) and to store an integer as a vector of remainders from division of the integer by each of those primes. -Chinese remainder theorem states that this representation is sufficient to uniquely restore any number from 0 to product of these primes minus one. [Garner algorithm](chinese-remainder-theorem.md) allows to restore the number from such representation to normal integer. +Chinese remainder theorem states that this representation is sufficient to uniquely restore any number from 0 to product of these primes minus one. [Garner algorithm](garners-algorithm.md) allows to restore the number from such representation to normal integer. This method allows to save memory compared to the classical approach (though the savings are not as dramatic as in factorization representation). Besides, it allows to perform fast addition, subtraction and multiplication in time proportional to the number of prime numbers used as modulos (see [Chinese remainder theorem](chinese-remainder-theorem.md) article for implementation). diff --git a/src/algebra/binary-exp.md b/src/algebra/binary-exp.md index 05dbacc76..52dfd27a5 100644 --- a/src/algebra/binary-exp.md +++ b/src/algebra/binary-exp.md @@ -96,7 +96,7 @@ Compute $x^n \bmod m$. This is a very common operation. For instance it is used in computing the [modular multiplicative inverse](module-inverse.md). **Solution:** -Since we know that the module operator doesn't interfere with multiplications ($a \cdot b \equiv (a \bmod m) \cdot (b \bmod m) \pmod m$), we can directly use the same code, and just replace every multiplication with a modular multiplication: +Since we know that the modulo operator doesn't interfere with multiplications ($a \cdot b \equiv (a \bmod m) \cdot (b \bmod m) \pmod m$), we can directly use the same code, and just replace every multiplication with a modular multiplication: ```cpp long long binpow(long long a, long long b, long long m) { @@ -144,13 +144,13 @@ vector
+ 
+
-There is still one big open question.
-We don't know $p$ yet, so how can we argue about the sequence $\{x_i \bmod p\}$?
+Yet, there is still an open question.
+How can we exploit the properties of the sequence $\{x_i \bmod p\}$ to our advantage without even knowing the number $p$ itself?
It's actually quite easy.
There is a cycle in the sequence $\{x_i \bmod p\}_{i \le j}$ if and only if there are two indices $s, t \le j$ such that $x_s \equiv x_t \bmod p$.
This equation can be rewritten as $x_s - x_t \equiv 0 \bmod p$ which is the same as $p ~|~ \gcd(x_s - x_t, n)$.
Therefore, if we find two indices $s$ and $t$ with $g = \gcd(x_s - x_t, n) > 1$, we have found a cycle and also a factor $g$ of $n$.
-Notice that it is possible that $g = n$.
-In this case we haven't found a proper factor, and we have to repeat the algorithm with different parameter (different starting value $x_0$, different constant $c$ in the polynomial function $f$).
+It is possible that $g = n$.
+In this case we haven't found a proper factor, so we must repeat the algorithm with a different parameter (different starting value $x_0$, different constant $c$ in the polynomial function $f$).
To find the cycle, we can use any common cycle detection algorithm.
### Floyd's cycle-finding algorithm
-This algorithm finds a cycle by using two pointers.
-These pointers move over the sequence at different speeds.
-In each iteration the first pointer advances to the next element, but the second pointer advances two elements.
-It's not hard to see, that if there exists a cycle, the second pointer will make at least one full cycle and then meet the first pointer during the next few cycle loops.
+This algorithm finds a cycle by using two pointers moving over the sequence at differing speeds.
+During each iteration, the first pointer will advance one element over, while the second pointer advances to every other element.
+Using this idea it is easy to observe that if there is a cycle, at some point the second pointer will come around to meet the first one during the loops.
If the cycle length is $\lambda$ and the $\mu$ is the first index at which the cycle starts, then the algorithm will run in $O(\lambda + \mu)$ time.
-This algorithm is also known as **tortoise and the hare algorithm**, based on the tale in which a tortoise (here a slow pointer) and a hare (here a faster pointer) make a race.
+This algorithm is also known as the [Tortoise and Hare algorithm](../others/tortoise_and_hare.md), based on the tale in which a tortoise (the slow pointer) and a hare (the faster pointer) have a race.
-It is actually possible to determine the parameter $\lambda$ and $\mu$ using this algorithm (also in $O(\lambda + \mu)$ time and $O(1)$ space), but here is just the simplified version for finding the cycle at all.
-The algorithm and returns true as soon as it detects a cycle.
-If the sequence doesn't have a cycle, then the function will never stop.
-However this cannot happen during Pollard's rho algorithm.
+It is actually possible to determine the parameter $\lambda$ and $\mu$ using this algorithm (also in $O(\lambda + \mu)$ time and $O(1)$ space).
+When a cycle is detected, the algorithm will return 'True'.
+If the sequence doesn't have a cycle, then the function will loop endlessly.
+However, using Pollard's Rho Algorithm, this can be prevented.
```text
function floyd(f, x0):
@@ -290,8 +288,8 @@ function floyd(f, x0):
### Implementation
-First here is a implementation using the **Floyd's cycle-finding algorithm**.
-The algorithm runs (usually) in $O(\sqrt[4]{n} \log(n))$ time.
+First, here is an implementation using the **Floyd's cycle-finding algorithm**.
+The algorithm generally runs in $O(\sqrt[4]{n} \log(n))$ time.
```{.cpp file=pollard_rho}
long long mult(long long a, long long b, long long mod) {
@@ -353,16 +351,15 @@ long long mult(long long a, long long b, long long mod) {
Alternatively you can also implement the [Montgomery multiplication](montgomery_multiplication.md).
-As already noticed above: if $n$ is composite and the algorithm returns $n$ as factor, you have to repeat the procedure with different parameter $x_0$ and $c$.
+As stated previously, if $n$ is composite and the algorithm returns $n$ as factor, you have to repeat the procedure with different parameters $x_0$ and $c$.
E.g. the choice $x_0 = c = 1$ will not factor $25 = 5 \cdot 5$.
-The algorithm will just return $25$.
-However the choice $x_0 = 1$, $c = 2$ will factor it.
+The algorithm will return $25$.
+However, the choice $x_0 = 1$, $c = 2$ will factor it.
### Brent's algorithm
-Brent uses a similar algorithm as Floyd.
-It also uses two pointer.
-But instead of advancing the pointers by one and two respectably, we advance them in powers of two.
+Brent implements a similar method to Floyd, using two pointers.
+The difference being that instead of advancing the pointers by one and two places respectively, they are advanced by powers of two.
As soon as $2^i$ is greater than $\lambda$ and $\mu$, we will find the cycle.
```text
@@ -380,12 +377,12 @@ function floyd(f, x0):
return true
```
-Brent's algorithm also runs in linear time, but is usually faster than Floyd's algorithm, since it uses less evaluations of the function $f$.
+Brent's algorithm also runs in linear time, but is generally faster than Floyd's, since it uses less evaluations of the function $f$.
### Implementation
-The straightforward implementation using Brent's algorithms can be speeded up by noticing, that we can omit the terms $x_l - x_k$ if $k < \frac{3 \cdot l}{2}$.
-Also, instead of performing the $\gcd$ computation at every step, we multiply the terms and do it every few steps and backtrack if we overshoot.
+The straightforward implementation of Brent's algorithm can be sped up by omitting the terms $x_l - x_k$ if $k < \frac{3 \cdot l}{2}$.
+In addition, instead of performing the $\gcd$ computation at every step, we multiply the terms and only actually check $\gcd$ every few steps and backtrack if overshot.
```{.cpp file=pollard_rho_brent}
long long brent(long long n, long long x0=2, long long c=1) {
@@ -422,7 +419,7 @@ long long brent(long long n, long long x0=2, long long c=1) {
}
```
-The combination of a trial division for small prime numbers together with Brent's version of Pollard's rho algorithm will make a very powerful factorization algorithm.
+The combination of a trial division for small prime numbers together with Brent's version of Pollard's rho algorithm makes a very powerful factorization algorithm.
## Practice Problems
diff --git a/src/algebra/fft.md b/src/algebra/fft.md
index c77d81a1d..38ed620a5 100644
--- a/src/algebra/fft.md
+++ b/src/algebra/fft.md
@@ -15,8 +15,7 @@ Some researchers attribute the discovery of the FFT to Runge and König in 1924.
But actually Gauss developed such a method already in 1805, but never published it.
Notice, that the FFT algorithm presented here runs in $O(n \log n)$ time, but it doesn't work for multiplying arbitrary big polynomials with arbitrary large coefficients or for multiplying arbitrary big integers.
-It can easily handle polynomials of size $10^5$ with small coefficients, or multiplying two numbers of size $10^6$, but at some point the range and the precision of the used floating point numbers will not no longer be enough to give accurate results.
-That is usually enough for solving competitive programming problems, but there are also more complex variations that can perform arbitrary large polynomial/integer multiplications.
+It can easily handle polynomials of size $10^5$ with small coefficients, or multiplying two numbers of size $10^6$, which is usually enough for solving competitive programming problems. Beyond the scale of multiplying numbers with $10^6$ bits, the range and precision of the floating point numbers used during the computation will not be enough to give accurate final results, though there are more complex variations that can perform arbitrary large polynomial/integer multiplications.
E.g. in 1971 Schönhage and Strasser developed a variation for multiplying arbitrary large numbers that applies the FFT recursively in rings structures running in $O(n \log n \log \log n)$.
And recently (in 2019) Harvey and van der Hoeven published an algorithm that runs in true $O(n \log n)$.
@@ -98,7 +97,7 @@ It is easy to see that
$$A(x) = A_0(x^2) + x A_1(x^2).$$
-The polynomials $A_0$ and $A_1$ are only half as much coefficients as the polynomial $A$.
+The polynomials $A_0$ and $A_1$ have only half as many coefficients as the polynomial $A$.
If we can compute the $\text{DFT}(A)$ in linear time using $\text{DFT}(A_0)$ and $\text{DFT}(A_1)$, then we get the recurrence $T_{\text{DFT}}(n) = 2 T_{\text{DFT}}\left(\frac{n}{2}\right) + O(n)$ for the time complexity, which results in $T_{\text{DFT}}(n) = O(n \log n)$ by the **master theorem**.
Let's learn how we can accomplish that.
diff --git a/src/algebra/fibonacci-numbers.md b/src/algebra/fibonacci-numbers.md
index 3653d4bd0..6d18015f0 100644
--- a/src/algebra/fibonacci-numbers.md
+++ b/src/algebra/fibonacci-numbers.md
@@ -22,6 +22,8 @@ Fibonacci numbers possess a lot of interesting properties. Here are a few of the
$$F_{n-1} F_{n+1} - F_n^2 = (-1)^n$$
+>This can be proved by induction. A one-line proof by Knuth comes from taking the determinant of the 2x2 matrix form below.
+
* The "addition" rule:
$$F_{n+k} = F_k F_{n+1} + F_{k-1} F_n$$
@@ -50,7 +52,7 @@ such that $k_1 \ge k_2 + 2,\ k_2 \ge k_3 + 2,\ \ldots,\ k_r \ge 2$ (i.e.: the
It follows that any number can be uniquely encoded in the Fibonacci coding.
And we can describe this representation with binary codes $d_0 d_1 d_2 \dots d_s 1$, where $d_i$ is $1$ if $F_{i+2}$ is used in the representation.
-The code will be appended by a $1$ do indicate the end of the code word.
+The code will be appended by a $1$ to indicate the end of the code word.
Notice that this is the only occurrence where two consecutive 1-bits appear.
$$\begin{eqnarray}
@@ -74,9 +76,8 @@ The encoding of an integer $n$ can be done with a simple greedy algorithm:
To decode a code word, first remove the final $1$. Then, if the $i$-th bit is set (indexing from 0 from the leftmost to the rightmost bit), sum $F_{i+2}$ to the number.
-## Formulas for the $n^{\text{th}}$ Fibonacci number { data-toc-label="Formulas for the -th Fibonacci number" }
-The $n$-th Fibonacci number can be easily found in $O(n)$ by computing the numbers one by one up to $n$. However, there are also faster ways, as we will see.
+## Formulas for the $n^{\text{th}}$ Fibonacci number { data-toc-label="Formulas for the -th Fibonacci number" }
### Closed-form expression
@@ -94,30 +95,159 @@ where the square brackets denote rounding to the nearest integer.
As these two formulas would require very high accuracy when working with fractional numbers, they are of little use in practical calculations.
-### Matrix form
+### Fibonacci in linear time
-It is easy to prove the following relation:
+The $n$-th Fibonacci number can be easily found in $O(n)$ by computing the numbers one by one up to $n$. However, there are also faster ways, as we will see.
-$$\begin{pmatrix}F_{n-1} & F_{n} \cr\end{pmatrix} = \begin{pmatrix}F_{n-2} & F_{n-1} \cr\end{pmatrix} \cdot \begin{pmatrix}0 & 1 \cr 1 & 1 \cr\end{pmatrix}$$
+We can start from an iterative approach, to take advantage of the use of the formula $F_n = F_{n-1} + F_{n-2}$, therefore, we will simply precalculate those values in an array. Taking into account the base cases for $F_0$ and $F_1$.
+
+```{.cpp file=fibonacci_linear}
+int fib(int n) {
+ int a = 0;
+ int b = 1;
+ for (int i = 0; i < n; i++) {
+ int tmp = a + b;
+ a = b;
+ b = tmp;
+ }
+ return a;
+}
+```
-Denoting $P \equiv \begin{pmatrix}0 & 1 \cr 1 & 1 \cr\end{pmatrix}$, we have:
+In this way, we obtain a linear solution, $O(n)$ time, saving all the values prior to $n$ in the sequence.
-$$\begin{pmatrix}F_n & F_{n+1} \cr\end{pmatrix} = \begin{pmatrix}F_0 & F_1 \cr\end{pmatrix} \cdot P^n$$
+### Matrix form
-Thus, in order to find $F_n$, we must raise the matrix $P$ to $n$. This can be done in $O(\log n)$ (see [Binary exponentiation](binary-exp.md)).
+To go from $(F_n, F_{n-1})$ to $(F_{n+1}, F_n)$, we can express the linear recurrence as a 2x2 matrix multiplication:
+
+$$
+\begin{pmatrix}
+1 & 1 \\
+1 & 0
+\end{pmatrix}
+\begin{pmatrix}
+F_n \\
+F_{n-1}
+\end{pmatrix}
+=
+\begin{pmatrix}
+F_n + F_{n-1} \\
+F_{n}
+\end{pmatrix}
+=
+\begin{pmatrix}
+F_{n+1} \\
+F_{n}
+\end{pmatrix}
+$$
+
+This lets us treat iterating the recurrence as repeated matrix multiplication, which has nice properties. In particular,
+
+$$
+\begin{pmatrix}
+1 & 1 \\
+1 & 0
+\end{pmatrix}^n
+\begin{pmatrix}
+F_1 \\
+F_0
+\end{pmatrix}
+=
+\begin{pmatrix}
+F_{n+1} \\
+F_{n}
+\end{pmatrix}
+$$
+
+where $F_1 = 1, F_0 = 0$.
+In fact, since
+
+$$
+\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}
+= \begin{pmatrix} F_2 & F_1 \\ F_1 & F_0 \end{pmatrix}
+$$
+
+we can use the matrix directly:
+
+$$
+\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n
+= \begin{pmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{pmatrix}
+$$
+
+Thus, in order to find $F_n$ in $O(\log n)$ time, we must raise the matrix to n. (See [Binary exponentiation](binary-exp.md))
+
+```{.cpp file=fibonacci_matrix}
+struct matrix {
+ long long mat[2][2];
+ matrix friend operator *(const matrix &a, const matrix &b){
+ matrix c;
+ for (int i = 0; i < 2; i++) {
+ for (int j = 0; j < 2; j++) {
+ c.mat[i][j] = 0;
+ for (int k = 0; k < 2; k++) {
+ c.mat[i][j] += a.mat[i][k] * b.mat[k][j];
+ }
+ }
+ }
+ return c;
+ }
+};
+
+matrix matpow(matrix base, long long n) {
+ matrix ans{ {
+ {1, 0},
+ {0, 1}
+ } };
+ while (n) {
+ if(n&1)
+ ans = ans*base;
+ base = base*base;
+ n >>= 1;
+ }
+ return ans;
+}
-### Fast Doubling Method
+long long fib(int n) {
+ matrix base{ {
+ {1, 1},
+ {1, 0}
+ } };
+ return matpow(base, n).mat[0][1];
+}
+```
-Using above method we can find these equations:
+### Fast Doubling Method
-$$ \begin{array}{rll}
- F_{2k} &= F_k \left( 2F_{k+1} - F_{k} \right). \\
- F_{2k+1} &= F_{k+1}^2 + F_{k}^2.
-\end{array}$$
+By expanding the above matrix expression for $n = 2\cdot k$
+
+$$
+\begin{pmatrix}
+F_{2k+1} & F_{2k}\\
+F_{2k} & F_{2k-1}
+\end{pmatrix}
+=
+\begin{pmatrix}
+1 & 1\\
+1 & 0
+\end{pmatrix}^{2k}
+=
+\begin{pmatrix}
+F_{k+1} & F_{k}\\
+F_{k} & F_{k-1}
+\end{pmatrix}
+^2
+$$
+
+we can find these simpler equations:
+
+$$ \begin{align}
+F_{2k+1} &= F_{k+1}^2 + F_{k}^2 \\
+F_{2k} &= F_k(F_{k+1}+F_{k-1}) = F_k (2F_{k+1} - F_{k})\\
+\end{align}.$$
Thus using above two equations Fibonacci numbers can be calculated easily by the following code:
-```cpp
+```{.cpp file=fibonacci_doubling}
pair
+
+ 
+
The idea behind is this:
A number is prime, if none of the smaller prime numbers divides it.
-Since we iterate over the prime numbers in order, we already marked all numbers, who are divisible by at least one of the prime numbers, as divisible.
+Since we iterate over the prime numbers in order, we already marked all numbers, which are divisible by at least one of the prime numbers, as divisible.
Hence if we reach a cell and it is not marked, then it isn't divisible by any smaller prime number and therefore has to be prime.
## Implementation
@@ -53,7 +55,7 @@ Using such implementation the algorithm consumes $O(n)$ of the memory (obviously
It's simple to prove a running time of $O(n \log n)$ without knowing anything about the distribution of primes - ignoring the `is_prime` check, the inner loop runs (at most) $n/i$ times for $i = 2, 3, 4, \dots$, leading the total number of operations in the inner loop to be a harmonic sum like $n(1/2 + 1/3 + 1/4 + \cdots)$, which is bounded by $O(n \log n)$.
Let's prove that algorithm's running time is $O(n \log \log n)$.
-The algorithm will perform $\frac{n}{p}$ operations for every prime $p \le n$ the inner loop.
+The algorithm will perform $\frac{n}{p}$ operations for every prime $p \le n$ in the inner loop.
Hence, we need to evaluate the next expression:
$$\sum_{\substack{p \le n, \\\ p \text{ prime}}} \frac n p = n \cdot \sum_{\substack{p \le n, \\\ p \text{ prime}}} \frac 1 p.$$
@@ -61,7 +63,7 @@ $$\sum_{\substack{p \le n, \\\ p \text{ prime}}} \frac n p = n \cdot \sum_{\subs
Let's recall two known facts.
- The number of prime numbers less than or equal to $n$ is approximately $\frac n {\ln n}$.
- - The $k$-th prime number approximately equals $k \ln k$ (that follows immediately from the previous fact).
+ - The $k$-th prime number approximately equals $k \ln k$ (this follows from the previous fact).
Thus we can write down the sum in the following way:
@@ -115,7 +117,7 @@ Such optimization doesn't affect the complexity (indeed, by repeating the proof
Since all even numbers (except $2$) are composite, we can stop checking even numbers at all. Instead, we need to operate with odd numbers only.
-First, it will allow us to half the needed memory. Second, it will reduce the number of operations performing by algorithm approximately in half.
+First, it will allow us to halve the needed memory. Second, it will reduce the number of operations performed by algorithm approximately in half.
### Memory consumption and speed of operations
@@ -139,13 +141,13 @@ Another drawback from `bitset` is that you need to know the size at compile time
### Segmented Sieve
-It follows from the optimization "sieving till root" that there is no need to keep the whole array `is_prime[1...n]` at all time.
+It follows from the optimization "sieving till root" that there is no need to keep the whole array `is_prime[1...n]` at all times.
For sieving it is enough to just keep the prime numbers until the root of $n$, i.e. `prime[1... sqrt(n)]`, split the complete range into blocks, and sieve each block separately.
Let $s$ be a constant which determines the size of the block, then we have $\lceil {\frac n s} \rceil$ blocks altogether, and the block $k$ ($k = 0 ... \lfloor {\frac n s} \rfloor$) contains the numbers in a segment $[ks; ks + s - 1]$.
We can work on blocks by turns, i.e. for every block $k$ we will go through all the prime numbers (from $1$ to $\sqrt n$) and perform sieving using them.
It is worth noting, that we have to modify the strategy a little bit when handling the first numbers: first, all the prime numbers from $[1; \sqrt n]$ shouldn't remove themselves; and second, the numbers $0$ and $1$ should be marked as non-prime numbers.
-While working on the last block it should not be forgotten that the last needed number $n$ is not necessary located in the end of the block.
+While working on the last block it should not be forgotten that the last needed number $n$ is not necessarily located at the end of the block.
As discussed previously, the typical implementation of the Sieve of Eratosthenes is limited by the speed how fast you can load data into the CPU caches.
By splitting the range of potential prime numbers $[1; n]$ into smaller blocks, we never have to keep multiple blocks in memory at the same time, and all operations are much more cache-friendlier.
@@ -254,6 +256,8 @@ However, this algorithm also has its own weaknesses.
## Practice Problems
+* [Leetcode - Four Divisors](https://leetcode.com/problems/four-divisors/)
+* [Leetcode - Count Primes](https://leetcode.com/problems/count-primes/)
* [SPOJ - Printing Some Primes](http://www.spoj.com/problems/TDPRIMES/)
* [SPOJ - A Conjecture of Paul Erdos](http://www.spoj.com/problems/HS08PAUL/)
* [SPOJ - Primal Fear](http://www.spoj.com/problems/VECTAR8/)
@@ -269,4 +273,4 @@ However, this algorithm also has its own weaknesses.
* [SPOJ - Prime Generator](http://www.spoj.com/problems/PRIME1/)
* [SPOJ - Printing some primes (hard)](http://www.spoj.com/problems/PRIMES2/)
* [Codeforces - Nodbach Problem](https://codeforces.com/problemset/problem/17/A)
-* [Codefoces - Colliders](https://codeforces.com/problemset/problem/154/B)
+* [Codeforces - Colliders](https://codeforces.com/problemset/problem/154/B)
diff --git a/src/combinatorics/binomial-coefficients.md b/src/combinatorics/binomial-coefficients.md
index 382e16b06..9f065c1a7 100644
--- a/src/combinatorics/binomial-coefficients.md
+++ b/src/combinatorics/binomial-coefficients.md
@@ -85,7 +85,9 @@ int C(int n, int k) {
### Improved implementation
-Note that in the above implementation numerator and denominator have the same number of factors ($k$), each of which is greater than or equal to 1. Therefore, we can replace our fraction with a product $k$ fractions, each of which is real-valued. However, on each step after multiplying current answer by each of the next fractions the answer will still be integer (this follows from the property of factoring in). C++ implementation:
+Note that in the above implementation numerator and denominator have the same number of factors ($k$), each of which is greater than or equal to 1. Therefore, we can replace our fraction with a product $k$ fractions, each of which is real-valued. However, on each step after multiplying current answer by each of the next fractions the answer will still be integer (this follows from the property of factoring in).
+
+C++ implementation:
```cpp
int C(int n, int k) {
@@ -101,6 +103,7 @@ Here we carefully cast the floating point number to an integer, taking into acco
### Pascal's Triangle
By using the recurrence relation we can construct a table of binomial coefficients (Pascal's triangle) and take the result from it. The advantage of this method is that intermediate results never exceed the answer and calculating each new table element requires only one addition. The flaw is slow execution for large $n$ and $k$ if you just need a single value and not the whole table (because in order to calculate $\binom n k$ you will need to build a table of all $\binom i j, 1 \le i \le n, 1 \le j \le n$, or at least to $1 \le j \le \min (i, 2k)$). The time complexity can be considered to be $\mathcal{O}(n^2)$.
+
C++ implementation:
```cpp
@@ -217,7 +220,6 @@ When $m$ is not square-free, a [generalization of Lucas's theorem for prime powe
* [LightOj - Necklaces](http://www.lightoj.com/volume_showproblem.php?problem=1419)
* [HACKEREARTH: Binomial Coefficient](https://www.hackerearth.com/problem/algorithm/binomial-coefficient-1/description/)
* [SPOJ - Ada and Teams](http://www.spoj.com/problems/ADATEAMS/)
-* [DevSkill - Drive In Grid](https://devskill.com/CodingProblems/ViewProblem/61)
* [SPOJ - Greedy Walking](http://www.spoj.com/problems/UCV2013E/)
* [UVa 13214 - The Robot's Grid](https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=5137)
* [SPOJ - Good Predictions](http://www.spoj.com/problems/GOODB/)
@@ -225,11 +227,11 @@ When $m$ is not square-free, a [generalization of Lucas's theorem for prime powe
* [SPOJ - Topper Rama Rao](http://www.spoj.com/problems/HLP_RAMS/)
* [UVa 13184 - Counting Edges and Graphs](https://uva.onlinejudge.org/index.php?option=onlinejudge&page=show_problem&problem=5095)
* [Codeforces - Anton and School 2](http://codeforces.com/contest/785/problem/D)
-* [DevSkill - Parandthesis](https://devskill.com/CodingProblems/ViewProblem/255)
* [Codeforces - Bacterial Melee](http://codeforces.com/contest/760/problem/F)
* [Codeforces - Points, Lines and Ready-made Titles](http://codeforces.com/contest/872/problem/E)
* [SPOJ - The Ultimate Riddle](https://www.spoj.com/problems/DCEPC13D/)
* [CodeChef - Long Sandwich](https://www.codechef.com/MAY17/problems/SANDWICH/)
+* [Codeforces - Placing Jinas](https://codeforces.com/problemset/problem/1696/E)
## References
* [Blog fishi.devtail.io](https://fishi.devtail.io/weblog/2015/06/25/computing-large-binomial-coefficients-modulo-prime-non-prime/)
diff --git a/src/combinatorics/bracket_sequences.md b/src/combinatorics/bracket_sequences.md
index c737e1aa5..e24a18cb0 100644
--- a/src/combinatorics/bracket_sequences.md
+++ b/src/combinatorics/bracket_sequences.md
@@ -29,10 +29,10 @@ For this case there exists a very simple algorithm.
Let $\text{depth}$ be the current number of open brackets.
Initially $\text{depth} = 0$.
We iterate over all character of the string, if the current bracket character is an opening bracket, then we increment $\text{depth}$, otherwise we decrement it.
-If at any time the variable $\text{depth}$ gets negative, or at the end it is different from $0$, than the string is not a balances sequence.
+If at any time the variable $\text{depth}$ gets negative, or at the end it is different from $0$, then the string is not a balanced sequence.
Otherwise it is.
-If there are several bracket types involved, then the algorithm needs to be changes.
+If there are several bracket types involved, then the algorithm needs to be changed.
Instead of a counter $\text{depth}$ we create a stack, in which we will store all opening brackets that we meet.
If the current bracket character is an opening one, we put it onto the stack.
If it is a closing one, then we check if the stack is non-empty, and if the top element of the stack is of the same type as the current closing bracket.
@@ -49,7 +49,7 @@ The number of balanced bracket sequences of length $2n$ ($n$ pairs of brackets)
$$\frac{1}{n+1} \binom{2n}{n}$$
-If we allow $k$ types of brackets, then each pair be of any of the $k$ types (independently of the others), thus the number of balanced bracket sequences is:
+If we allow $k$ types of brackets, then each pair can be of any of the $k$ types (independently of the others), thus the number of balanced bracket sequences is:
$$\frac{1}{n+1} \binom{2n}{n} k^n$$
@@ -82,7 +82,7 @@ When we meet an opening brackets, we will decrement $\text{depth}$, and when we
If we are at some point meet an opening bracket, and the balance after processing this symbol is positive, then we have found the rightmost position that we can change.
We change the symbol, compute the number of opening and closing brackets that we have to add to the right side, and arrange them in the lexicographically minimal way.
-If we find do suitable position, then this sequence is already the maximal possible one, and there is no answer.
+If we don't find a suitable position, then this sequence is already the maximal possible one, and there is no answer.
```{.cpp file=next_balanced_brackets_sequence}
bool next_balanced_sequence(string & s) {
@@ -117,7 +117,7 @@ To generate then, we can start with the lexicographically smallest sequence $((\
However, if the length of the sequence is not very long (e.g. $n$ smaller than $12$), then we can also generate all permutations conveniently with the C++ STL function `next_permutation`, and check each one for balanceness.
-Also they can be generate using the ideas we used for counting all sequences with dynamic programming.
+Also they can be generated using the ideas we used for counting all sequences with dynamic programming.
We will discuss the ideas in the next two sections.
## Sequence index
@@ -142,19 +142,19 @@ Thus we can compute this array in $O(n^2)$.
Now let us generate the index for a given sequence.
First let there be only one type of brackets.
-We will us the counter $\text{depth}$ which tells us how nested we currently are, and iterate over the characters of the sequence.
+We will use the counter $\text{depth}$ which tells us how nested we currently are, and iterate over the characters of the sequence.
If the current character $s[i]$ is equal to $($, then we increment $\text{depth}$.
If the current character $s[i]$ is equal to $)$, then we must add $d[2n-i-1][\text{depth}+1]$ to the answer, taking all possible endings starting with a $($ into account (which are lexicographically smaller sequences), and then decrement $\text{depth}$.
New let there be $k$ different bracket types.
-Thus, when we look at the current character $s[i]$ before recomputing $\text{depth}$, we have to go through all bracket types that are smaller than the current character, and try to put this bracket into the current position (obtaining a new balance $\text{ndepth} = \text{depth} \pm 1$), and add the number of ways to finish the sequence (length $2n-i-1$, balance $ndepth$) to the answer:
+Thus, when we look at the current character $s[i]$ before recomputing $\text{depth}$, we have to go through all bracket types that are smaller than the current character, and try to place this bracket into the current position (obtaining a new balance $\text{ndepth} = \text{depth} \pm 1$), and add the number of ways to finish the sequence (length $2n-i-1$, balance $ndepth$) to the answer:
$$d[2n - i - 1][\text{ndepth}] \cdot k^{\frac{2n - i - 1 - ndepth}{2}}$$
This formula can be derived as follows:
First we "forget" that there are multiple bracket types, and just take the answer $d[2n - i - 1][\text{ndepth}]$.
-Now we consider how the answer will change is we have $k$ types of brackets.
+Now we consider how the answer will change if we have $k$ types of brackets.
We have $2n - i - 1$ undefined positions, of which $\text{ndepth}$ are already predetermined because of the opening brackets.
But all the other brackets ($(2n - i - 1 - \text{ndepth})/2$ pairs) can be of any type, therefore we multiply the number by such a power of $k$.
@@ -169,10 +169,9 @@ First, we start with only one bracket type.
We will iterate over the characters in the string we want to generate.
As in the previous problem we store a counter $\text{depth}$, the current nesting depth.
-In each position we have to decide if we use an opening of a closing bracket.
-To have to put an opening bracket character, it $d[2n - i - 1][\text{depth}+1] \ge k$.
-We increment the counter $\text{depth}$, and move on to the next character.
-Otherwise we decrement $k$ by $d[2n - i - 1][\text{depth}+1]$, put a closing bracket and move on.
+At each position, we have to decide whether to place an opening or a closing bracket. To place an opening bracket, $d[2n - i - 1][\text{depth}+1] \ge k$ must be true.
+If so, we increment the counter $\text{depth}$, and move on to the next character.
+Otherwise, we decrement $k$ by $d[2n - i - 1][\text{depth}+1]$, place a closing bracket, and move on.
```{.cpp file=kth_balances_bracket}
string kth_balanced(int n, int k) {
diff --git a/src/combinatorics/burnside.md b/src/combinatorics/burnside.md
index fa799399c..894b1e87d 100644
--- a/src/combinatorics/burnside.md
+++ b/src/combinatorics/burnside.md
@@ -266,3 +266,8 @@ int solve(int n, int m) {
return sum / s.size();
}
```
+## Practice Problems
+* [CSES - Counting Necklaces](https://cses.fi/problemset/task/2209)
+* [CSES - Counting Grids](https://cses.fi/problemset/task/2210)
+* [Codeforces - Buildings](https://codeforces.com/gym/101873/problem/B)
+* [CS Academy - Cube Coloring](https://csacademy.com/contest/beta-round-8/task/cube-coloring/)
diff --git a/src/combinatorics/catalan-numbers.md b/src/combinatorics/catalan-numbers.md
index 8ede5e92c..94d559a78 100644
--- a/src/combinatorics/catalan-numbers.md
+++ b/src/combinatorics/catalan-numbers.md
@@ -9,7 +9,7 @@ Catalan numbers is a number sequence, which is found useful in a number of combi
This sequence was named after the Belgian mathematician [Catalan](https://en.wikipedia.org/wiki/Eug%C3%A8ne_Charles_Catalan), who lived in the 19th century. (In fact it was known before to Euler, who lived a century before Catalan).
-The first few numbers Catalan numbers, $C_n$ (starting from zero):
+The first few Catalan numbers $C_n$ (starting from zero):
$1, 1, 2, 5, 14, 42, 132, 429, 1430, \ldots$
@@ -47,7 +47,7 @@ You can also think it in this manner. By definition, $C_n$ denotes number of cor
$( ) ( ( ) )$ can be divided into $( )$ and $( ( ) )$, but cannot be divided into $( ) ($ and $( ) )$. Again summing over all admissible $k's$, we get the recurrence relation on $C_n$.
-#### C++ implementation Show/Hide
+#### C++ implementation
```cpp
const int MOD = ....
diff --git a/src/combinatorics/inclusion-exclusion.md b/src/combinatorics/inclusion-exclusion.md
index cb96ba96d..ed6767216 100644
--- a/src/combinatorics/inclusion-exclusion.md
+++ b/src/combinatorics/inclusion-exclusion.md
@@ -153,7 +153,7 @@ we want to break a sequence of $20$ units into $6$ groups, which is the same as
$$N_0 = \binom{25}{5}$$
-We will now calculate the number of "bad" solutions with the inclusion-exclusion principle. The "bad" solutions will be those in which one or more $x_i$ are greater than $9$.
+We will now calculate the number of "bad" solutions with the inclusion-exclusion principle. The "bad" solutions will be those in which one or more $x_i$ are greater than or equal to $9$.
Denote by $A_k ~ (k = 1,2\ldots 6)$ the set of solutions where $x_k \ge 9$, and all other $x_i \ge 0 ~ (i \ne k)$ (they may be $\ge 9$ or not). To calculate the size of $A_k$, note that we have essentially the same combinatorial problem that was solved in the two paragraphs above, but now $9$ of the units are excluded from the slots and definitely belong to the first group. Thus:
@@ -169,6 +169,14 @@ Combining all this into the formula of inclusions-exceptions and given that we s
$$\binom{25}{5} - \left(\binom{6}{1} \cdot \binom{16}{5} - \binom{6}{2} \cdot \binom{7}{5}\right) $$
+This easily generalizes to $d$ numbers that sum up to $s$ with the restriction $0 \le x_i \le b$:
+
+$$\sum_{i=0}^d (-1)^i \binom{d}{i} \binom{s+d-1-(b+1)i}{d-1}$$
+
+As above, we treat binomial coefficients with negative upper index as zero.
+
+Note this problem could also be solved with dynamic programming or generating functions. The inclusion-exclusion answer is computed in $O(d)$ time (assuming math operations like binomial coefficient are constant time), while a simple DP approach would take $O(ds)$ time.
+
### The number of relative primes in a given interval
Task: given two numbers $n$ and $r$, count the number of integers in the interval $[1;r]$ that are relatively prime to n (their greatest common divisor is $1$).
@@ -441,7 +449,6 @@ A list of tasks that can be solved using the principle of inclusions-exclusions:
* [TopCoder SRM 390 "SetOfPatterns" [difficulty: medium]](http://www.topcoder.com/stat?c=problem_statement&pm=8307)
* [TopCoder SRM 176 "Deranged" [difficulty: medium]](http://community.topcoder.com/stat?c=problem_statement&pm=2013)
* [TopCoder SRM 457 "TheHexagonsDivOne" [difficulty: medium]](http://community.topcoder.com/stat?c=problem_statement&pm=10702&rd=14144&rm=303184&cr=22697599)
-* [Test>>>thebest "HarmonicTriples" (in Russian) [difficulty: medium]](http://esci.ru/ttb/statement-62.htm)
* [SPOJ #4191 MSKYCODE "Sky Code" [difficulty: medium]](http://www.spoj.com/problems/MSKYCODE/)
* [SPOJ #4168 SQFREE "Square-free integers" [difficulty: medium]](http://www.spoj.com/problems/SQFREE/)
* [CodeChef "Count Relations" [difficulty: medium]](http://www.codechef.com/JAN11/problems/COUNTREL/)
@@ -451,3 +458,4 @@ A list of tasks that can be solved using the principle of inclusions-exclusions:
* [SPOJ - EASY MATH [difficulty: medium]](http://www.spoj.com/problems/EASYMATH/)
* [SPOJ - MOMOS - FEASTOFPIGS [difficulty: easy]](https://www.spoj.com/problems/MOMOS/)
* [Atcoder - Grid 2 [difficulty: easy]](https://atcoder.jp/contests/dp/tasks/dp_y/)
+* [Codeforces - Count GCD](https://codeforces.com/contest/1750/problem/D)
diff --git a/src/combinatorics/stars_and_bars.md b/src/combinatorics/stars_and_bars.md
index 0f1561852..a8e5c395f 100644
--- a/src/combinatorics/stars_and_bars.md
+++ b/src/combinatorics/stars_and_bars.md
@@ -39,6 +39,17 @@ E.g. the solution $1 + 3 + 0 = 4$ for $n = 4$, $k = 3$ can be represented using
It is easy to see, that this is exactly the stars and bars theorem.
Therefore the solution is $\binom{n + k - 1}{n}$.
+## Number of positive integer sums
+
+A second theorem provides a nice interpretation for positive integers. Consider solutions to
+
+$$x_1 + x_2 + \dots + x_k = n$$
+
+with $x_i \ge 1$.
+
+We can consider $n$ stars, but this time we can put at most _one bar_ between stars, since two bars between stars would represent $x_i=0$, i.e. an empty box.
+There are $n-1$ gaps between stars to place $k-1$ bars, so the solution is $\binom{n-1}{k-1}$.
+
## Number of lower-bound integer sums
This can easily be extended to integer sums with different lower bounds.
@@ -61,3 +72,11 @@ So we have reduced the problem to the simpler case with $x_i' \ge 0$ and again c
With some help of the [Inclusion-Exclusion Principle](./inclusion-exclusion.md), you can also restrict the integers with upper bounds.
See the [Number of upper-bound integer sums](./inclusion-exclusion.md#number-of-upper-bound-integer-sums) section in the corresponding article.
+
+## Practice Problems
+
+* [Codeforces - Array](https://codeforces.com/contest/57/problem/C)
+* [Codeforces - Kyoya and Coloured Balls](https://codeforces.com/problemset/problem/553/A)
+* [Codeforces - Colorful Bricks](https://codeforces.com/contest/1081/problem/C)
+* [Codeforces - Two Arrays](https://codeforces.com/problemset/problem/1288/C)
+* [Codeforces - One-Dimensional Puzzle](https://codeforces.com/contest/1931/problem/G)
diff --git a/src/data_structures/disjoint_set_union.md b/src/data_structures/disjoint_set_union.md
index d25f76900..2964a87cf 100644
--- a/src/data_structures/disjoint_set_union.md
+++ b/src/data_structures/disjoint_set_union.md
@@ -98,7 +98,7 @@ The trick is to make the paths for all those nodes shorter, by setting the paren
You can see the operation in the following image.
On the left there is a tree, and on the right side there is the compressed tree after calling `find_set(7)`, which shortens the paths for the visited nodes 7, 5, 3 and 2.
-
+
The new implementation of `find_set` is as follows:
@@ -375,13 +375,13 @@ If we add an edge $(a, b)$ that connects two connected components into one, then
Let's derive a formula, which computes the parity issued to the leader of the set that will get attached to another set.
Let $x$ be the parity of the path length from vertex $a$ up to its leader $A$, and $y$ as the parity of the path length from vertex $b$ up to its leader $B$, and $t$ the desired parity that we have to assign to $B$ after the merge.
-The path contains the of the three parts:
+The path consists of the three parts:
from $B$ to $b$, from $b$ to $a$, which is connected by one edge and therefore has parity $1$, and from $a$ to $A$.
Therefore we receive the formula ($\oplus$ denotes the XOR operation):
$$t = x \oplus y \oplus 1$$
-Thus regardless of how many joins we perform, the parity of the edges is carried from on leader to another.
+Thus regardless of how many joins we perform, the parity of the edges is carried from one leader to another.
We give the implementation of the DSU that supports parity. As in the previous section we use a pair to store the ancestor and the parity. In addition for each set we store in the array `bipartite[]` whether it is still bipartite or not.
diff --git a/src/data_structures/fenwick.md b/src/data_structures/fenwick.md
index c67044d87..439885b83 100644
--- a/src/data_structures/fenwick.md
+++ b/src/data_structures/fenwick.md
@@ -6,44 +6,48 @@ e_maxx_link: fenwick_tree
# Fenwick Tree
-Let, $f$ be some group operation (binary associative function over a set with identity element and inverse elements) and $A$ be an array of integers of length $N$.
+Let $f$ be some group operation (a binary associative function over a set with an identity element and inverse elements) and $A$ be an array of integers of length $N$.
+Denote $f$'s infix notation as $*$; that is, $f(x,y) = x*y$ for arbitrary integers $x,y$.
+(Since this is associative, we will omit parentheses for order of application of $f$ when using infix notation.)
-Fenwick tree is a data structure which:
+The Fenwick tree is a data structure which:
-* calculates the value of function $f$ in the given range $[l, r]$ (i.e. $f(A_l, A_{l+1}, \dots, A_r)$) in $O(\log N)$ time;
-* updates the value of an element of $A$ in $O(\log N)$ time;
-* requires $O(N)$ memory, or in other words, exactly the same memory required for $A$;
-* is easy to use and code, especially, in the case of multidimensional arrays.
+* calculates the value of function $f$ in the given range $[l, r]$ (i.e. $A_l * A_{l+1} * \dots * A_r$) in $O(\log N)$ time
+* updates the value of an element of $A$ in $O(\log N)$ time
+* requires $O(N)$ memory (the same amount required for $A$)
+* is easy to use and code, especially in the case of multidimensional arrays
-The most common application of Fenwick tree is _calculating the sum of a range_ (i.e. using addition over the set of integers $\mathbb{Z}$: $f(A_1, A_2, \dots, A_k) = A_1 + A_2 + \dots + A_k$).
+The most common application of a Fenwick tree is _calculating the sum of a range_.
+For example, using addition over the set of integers as the group operation, i.e. $f(x,y) = x + y$: the binary operation, $*$, is $+$ in this case, so $A_l * A_{l+1} * \dots * A_r = A_l + A_{l+1} + \dots + A_{r}$.
-Fenwick tree is also called **Binary Indexed Tree**, or just **BIT** abbreviated.
-
-Fenwick tree was first described in a paper titled "A new data structure for cumulative frequency tables" (Peter M. Fenwick, 1994).
+The Fenwick tree is also called a **Binary Indexed Tree** (BIT).
+It was first described in a paper titled "A new data structure for cumulative frequency tables" (Peter M. Fenwick, 1994).
## Description
### Overview
-For the sake of simplicity, we will assume that function $f$ is just a *sum function*.
+For the sake of simplicity, we will assume that function $f$ is defined as $f(x,y) = x + y$ over the integers.
-Given an array of integers $A[0 \dots N-1]$.
-A Fenwick tree is just an array $T[0 \dots N-1]$, where each of its elements is equal to the sum of elements of $A$ in some range $[g(i), i]$:
+Suppose we are given an array of integers, $A[0 \dots N-1]$.
+(Note that we are using zero-based indexing.)
+A Fenwick tree is just an array, $T[0 \dots N-1]$, where each element is equal to the sum of elements of $A$ in some range, $[g(i), i]$:
-$$T_i = \sum_{j = g(i)}^{i}{A_j},$$
+$$T_i = \sum_{j = g(i)}^{i}{A_j}$$
where $g$ is some function that satisfies $0 \le g(i) \le i$.
-We will define the function in the next few paragraphs.
+We will define $g$ in the next few paragraphs.
-The data structure is called tree, because there is a nice representation of the data structure as tree, although we don't need to model an actual tree with nodes and edges.
-We will only need to maintain the array $T$ to handle all queries.
+The data structure is called a tree because there is a nice representation of it in the form of a tree, although we don't need to model an actual tree with nodes and edges.
+We only need to maintain the array $T$ to handle all queries.
**Note:** The Fenwick tree presented here uses zero-based indexing.
-Many people will actually use a version of the Fenwick tree that uses one-based indexing.
-Therefore you will also find an alternative implementation using one-based indexing in the implementation section.
+Many people use a version of the Fenwick tree that uses one-based indexing.
+As such, you will also find an alternative implementation which uses one-based indexing in the implementation section.
Both versions are equivalent in terms of time and memory complexity.
-Now we can write some pseudo-code for the two operations mentioned above - get the sum of elements of $A$ in the range $[0, r]$ and update (increase) some element $A_i$:
+Now we can write some pseudo-code for the two operations mentioned above.
+Below, we get the sum of elements of $A$ in the range $[0, r]$ and update (increase) some element $A_i$:
```python
def sum(int r):
@@ -60,20 +64,21 @@ def increase(int i, int delta):
The function `sum` works as follows:
-1. first, it adds the sum of the range $[g(r), r]$ (i.e. $T[r]$) to the `result`
-2. then, it "jumps" to the range $[g(g(r)-1), g(r)-1]$, and adds this range's sum to the `result`
-3. and so on, until it "jumps" from $[0, g(g( \dots g(r)-1 \dots -1)-1)]$ to $[g(-1), -1]$; that is where the `sum` function stops jumping.
+1. First, it adds the sum of the range $[g(r), r]$ (i.e. $T[r]$) to the `result`.
+2. Then, it "jumps" to the range $[g(g(r)-1), g(r)-1]$ and adds this range's sum to the `result`.
+3. This continues until it "jumps" from $[0, g(g( \dots g(r)-1 \dots -1)-1)]$ to $[g(-1), -1]$; this is where the `sum` function stops jumping.
-The function `increase` works with the same analogy, but "jumps" in the direction of increasing indices:
+The function `increase` works with the same analogy, but it "jumps" in the direction of increasing indices:
-1. sums of the ranges $[g(j), j]$ that satisfy the condition $g(j) \le i \le j$ are increased by `delta` , that is `t[j] += delta`. Therefore we updated all elements in $T$ that correspond to ranges in which $A_i$ lies.
+1. The sum for each range of the form $[g(j), j]$ which satisfies the condition $g(j) \le i \le j$ is increased by `delta`; that is, `t[j] += delta`.
+Therefore, it updates all elements in $T$ that correspond to ranges in which $A_i$ lies.
-It is obvious that the complexity of both `sum` and `increase` depend on the function $g$.
-There are lots of ways to choose the function $g$, as long as $0 \le g(i) \le i$ for all $i$.
-For instance the function $g(i) = i$ works, which results just in $T = A$, and therefore summation queries are slow.
-We can also take the function $g(i) = 0$.
-This will correspond to prefix sum arrays, which means that finding the sum of the range $[0, i]$ will only take constant time, but updates are slow.
-The clever part of the Fenwick algorithm is, that there it uses a special definition of the function $g$ that can handle both operations in $O(\log N)$ time.
+The complexity of both `sum` and `increase` depend on the function $g$.
+There are many ways to choose the function $g$ such that $0 \le g(i) \le i$ for all $i$.
+For instance, the function $g(i) = i$ works, which yields $T = A$ (in which case, the summation queries are slow).
+We could also take the function $g(i) = 0$.
+This would correspond to prefix sum arrays (in which case, finding the sum of the range $[0, i]$ will only take constant time; however, updates are slow).
+The clever part of the algorithm for Fenwick trees is how it uses a special definition of the function $g$ which can handle both operations in $O(\log N)$ time.
### Definition of $g(i)$ { data-toc-label='Definition of ' }
@@ -116,14 +121,16 @@ h(31) = 63 &= 0111111_2 \\\\
Unsurprisingly, there also exists a simple way to perform $h$ using bitwise operations:
-$$h(j) = j ~\|~ (j+1),$$
+$$h(j) = j ~|~ (j+1),$$
-where $\|$ is the bitwise OR operator.
+where $|$ is the bitwise OR operator.
The following image shows a possible interpretation of the Fenwick tree as tree.
The nodes of the tree show the ranges they cover.
-
+
+ 
+
## Implementation
@@ -148,7 +155,7 @@ struct FenwickTree {
bit.assign(n, 0);
}
- FenwickTree(vector
+
+ 
+
The dot product holds some notable properties:
@@ -181,7 +183,9 @@ To see the next important property we should take a look at the set of points $\
You can see that this set of points is exactly the set of points for which the projection onto $\mathbf a$ is the point $C \cdot \dfrac{\mathbf a}{|\mathbf a|}$ and they form a hyperplane orthogonal to $\mathbf a$.
You can see the vector $\mathbf a$ alongside with several such vectors having same dot product with it in 2D on the picture below:
-
+
+ 
+
In 2D these vectors will form a line, in 3D they will form a plane.
Note that this result allows us to define a line in 2D as $\mathbf r\cdot \mathbf n=C$ or $(\mathbf r - \mathbf r_0)\cdot \mathbf n=0$ where $\mathbf n$ is vector orthogonal to the line and $\mathbf r_0$ is any vector already present on the line and $C = \mathbf r_0\cdot \mathbf n$.
@@ -192,14 +196,18 @@ In the same manner a plane can be defined in 3D.
### Definition
Assume you have three vectors $\mathbf a$, $\mathbf b$ and $\mathbf c$ in 3D space joined in a parallelepiped as in the picture below:
-
+
+ 
+
How would you calculate its volume?
From school we know that we should multiply the area of the base with the height, which is projection of $\mathbf a$ onto direction orthogonal to base.
That means that if we define $\mathbf b \times \mathbf c$ as the vector which is orthogonal to both $\mathbf b$ and $\mathbf c$ and which length is equal to the area of the parallelogram formed by $\mathbf b$ and $\mathbf c$ then $|\mathbf a\cdot (\mathbf b\times\mathbf c)|$ will be equal to the volume of the parallelepiped.
For integrity we will say that $\mathbf b\times \mathbf c$ will be always directed in such way that the rotation from the vector $\mathbf b$ to the vector $\mathbf c$ from the point of $\mathbf b\times \mathbf c$ is always counter-clockwise (see the picture below).
-
+
+ 
+
This defines the cross (or vector) product $\mathbf b\times \mathbf c$ of the vectors $\mathbf b$ and $\mathbf c$ and the triple product $\mathbf a\cdot(\mathbf b\times \mathbf c)$ of the vectors $\mathbf a$, $\mathbf b$ and $\mathbf c$.
diff --git a/src/geometry/chebyshev-transformation.png b/src/geometry/chebyshev-transformation.png
new file mode 100644
index 000000000..bcad9616f
Binary files /dev/null and b/src/geometry/chebyshev-transformation.png differ
diff --git a/src/geometry/convex-hull.md b/src/geometry/convex-hull.md
index 1117e1e64..f74061d37 100644
--- a/src/geometry/convex-hull.md
+++ b/src/geometry/convex-hull.md
@@ -23,7 +23,7 @@ with the same Y coordinate, the one with the smaller X coordinate is considered.
step takes $\mathcal{O}(N)$ time.
Next, all the other points are sorted by polar angle in clockwise order.
-If the polar angle between two points is the same, the nearest point is chosen instead.
+If the polar angle between two or more points is the same, the tie should be broken by distance from $P_0$, in increasing order.
Then we iterate through each point one by one, and make sure that the current
point and the two before it make a clockwise turn, otherwise the previous
@@ -47,6 +47,9 @@ of the algorithm, otherwise you wouldn't get the smallest convex hull.
```{.cpp file=graham_scan}
struct pt {
double x, y;
+ bool operator == (pt const& t) const {
+ return x == t.x && y == t.y;
+ }
};
int orientation(pt a, pt b, pt c) {
@@ -86,6 +89,9 @@ void convex_hull(vector
+
+ 
+
One has to keep points on the convex hull and normal vectors of the hull's edges.
When you have a $(x;1)$ query you'll have to find the normal vector closest to it in terms of angles between them, then the optimum linear function will correspond to one of its endpoints.
@@ -90,7 +92,9 @@ Assume we're in some vertex corresponding to half-segment $[l,r)$ and the functi
Here is the illustration of what is going on in the vertex when we add new function:
-
+
+ 
+
Let's go to implementation now. Once again we will use complex numbers to keep linear functions.
diff --git a/src/geometry/delaunay.md b/src/geometry/delaunay.md
index c5f5c22d1..68d67ce8f 100644
--- a/src/geometry/delaunay.md
+++ b/src/geometry/delaunay.md
@@ -30,7 +30,9 @@ Because of the duality, we only need a fast algorithm to compute only one of $V$
## Quad-edge data structure
During the algorithm $D$ will be stored inside the quad-edge data structure. This structure is described in the picture:
-
+
+ 
+
In the algorithm we will use the following functions on edges:
diff --git a/src/geometry/intersecting_segments.md b/src/geometry/intersecting_segments.md
index 77416975b..ac26c8fd5 100644
--- a/src/geometry/intersecting_segments.md
+++ b/src/geometry/intersecting_segments.md
@@ -16,18 +16,24 @@ The naive solution algorithm is to iterate over all pairs of segments in $O(n^2)
Let's draw a vertical line $x = -\infty$ mentally and start moving this line to the right.
In the course of its movement, this line will meet with segments, and at each time a segment intersect with our line it intersects in exactly one point (we will assume that there are no vertical segments).
-
+
+ 
+
Thus, for each segment, at some point in time, its point will appear on the sweep line, then with the movement of the line, this point will move, and finally, at some point, the segment will disappear from the line.
We are interested in the **relative order of the segments** along the vertical.
Namely, we will store a list of segments crossing the sweep line at a given time, where the segments will be sorted by their $y$-coordinate on the sweep line.
-
+
+ 
+
This order is interesting because intersecting segments will have the same $y$-coordinate at least at one time:
-
+
+ 
+
We formulate key statements:
@@ -162,7 +168,7 @@ pair
+
+ 
+
- $k < 1$ and $b < 1$.
@@ -51,7 +53,9 @@ We have two cases:
which returns us back to the case $k>1$.
You can see new reference point $O'$ and axes $X'$ and $Y'$ in the picture below:
-
+
+ 
+
As you see, in new reference system linear function will have coefficient $\tfrac 1 k$ and its zero will be in the point $\lfloor k\cdot n + b \rfloor-(k\cdot n+b)$ which makes formula above correct.
## Complexity analysis
diff --git a/src/geometry/manhattan-distance.md b/src/geometry/manhattan-distance.md
new file mode 100644
index 000000000..87514868a
--- /dev/null
+++ b/src/geometry/manhattan-distance.md
@@ -0,0 +1,189 @@
+---
+tags:
+ - Original
+---
+
+# Manhattan Distance
+
+## Definition
+For points $p$ and $q$ on a plane, we can define the distance between them as the sum of the differences between their $x$ and $y$ coordinates:
+
+$$d(p,q) = |x_p - x_q| + |y_p - y_q|$$
+
+Defined this way, the distance corresponds to the so-called [Manhattan (taxicab) geometry](https://en.wikipedia.org/wiki/Taxicab_geometry), in which the points are considered intersections in a well designed city, like Manhattan, where you can only move on the streets horizontally or vertically, as shown in the image below:
+
+
+
+ 
+
+
+This images show some of the smallest paths from one black point to the other, all of them with length $12$.
+
+There are some interesting tricks and algorithms that can be done with this distance, and we will show some of them here.
+
+## Farthest pair of points in Manhattan distance
+
+Given $n$ points $P$, we want to find the pair of points $p,q$ that are farther apart, that is, maximize $|x_p - x_q| + |y_p - y_q|$.
+
+Let's think first in one dimension, so $y=0$. The main observation is that we can bruteforce if $|x_p - x_q|$ is equal to $x_p - x_q$ or $-x_p + x_q$, because if we "miss the sign" of the absolute value, we will get only a smaller value, so it can't affect the answer. More formally, it holds that:
+
+$$|x_p - x_q| = \max(x_p - x_q, -x_p + x_q)$$
+
+So, for example, we can try to have $p$ such that $x_p$ has the plus sign, and then $q$ must have the negative sign. This way we want to find:
+
+$$\max\limits_{p, q \in P}(x_p + (-x_q)) = \max\limits_{p \in P}(x_p) + \max\limits_{q \in P}( - x_q ).$$
+
+Notice that we can extend this idea further for 2 (or more!) dimensions. For $d$ dimensions, we must bruteforce $2^d$ possible values of the signs. For example, if we are in $2$ dimensions and bruteforce that $p$ has both the plus signs we want to find:
+
+$$\max\limits_{p, q \in P} [(x_p + (-x_q)) + (y_p + (-y_q))] = \max\limits_{p \in P}(x_p + y_p) + \max\limits_{q \in P}(-x_q - y_q).$$
+
+As we made $p$ and $q$ independent, it is now easy to find the $p$ and $q$ that maximize the expression.
+
+The code below generalizes this to $d$ dimensions and runs in $O(n \cdot 2^d \cdot d)$.
+
+```cpp
+long long ans = 0;
+for (int msk = 0; msk < (1 << d); msk++) {
+ long long mx = LLONG_MIN, mn = LLONG_MAX;
+ for (int i = 0; i < n; i++) {
+ long long cur = 0;
+ for (int j = 0; j < d; j++) {
+ if (msk & (1 << j)) cur += p[i][j];
+ else cur -= p[i][j];
+ }
+ mx = max(mx, cur);
+ mn = min(mn, cur);
+ }
+ ans = max(ans, mx - mn);
+}
+```
+
+## Rotating the points and Chebyshev distance
+
+It's well known that, for all $m, n \in \mathbb{R}$,
+
+$$|m| + |n| = \text{max}(|m + n|, |m - n|).$$
+
+To prove this, we just need to analyze the signs of $m$ and $n$. And it's left as an exercise.
+
+We may apply this equation to the Manhattan distance formula to find out that
+
+$$d((x_1, y_1), (x_2, y_2)) = |x_1 - x_2| + |y_1 - y_2| = \text{max}(|(x_1 + y_1) - (x_2 + y_2)|, |(x_1 - y_1) - (x_2 - y_2)|).$$
+
+The last expression in the previous equation is the [Chebyshev distance](https://en.wikipedia.org/wiki/Chebyshev_distance) of the points $(x_1 + y_1, x_1 - y_1)$ and $(x_2 + y_2, x_2 - y_2)$. This means that, after applying the transformation
+
+$$\alpha : (x, y) \to (x + y, x - y),$$
+
+the Manhattan distance between the points $p$ and $q$ turns into the Chebyshev distance between $\alpha(p)$ and $\alpha(q)$.
+
+Also, we may realize that $\alpha$ is a [spiral similarity](https://en.wikipedia.org/wiki/Spiral_similarity) (rotation of the plane followed by a dilation about a center $O$) with center $(0, 0)$, rotation angle of $45^{\circ}$ in clockwise direction and dilation by $\sqrt{2}$.
+
+Here's an image to help visualizing the transformation:
+
+
+
+ 
+
+
+## Manhattan Minimum Spanning Tree
+
+The Manhattan MST problem consists of, given some points in the plane, find the edges that connect all the points and have a minimum total sum of weights. The weight of an edge that connects two points is their Manhattan distance. For simplicity, we assume that all points have different locations.
+Here we show a way of finding the MST in $O(n \log{n})$ by finding for each point its nearest neighbor in each octant, as represented by the image below. This will give us $O(n)$ candidate edges, which, as we show below, will guarantee that they contain the MST. The final step is then using some standard MST, for example, [Kruskal algorithm using disjoint set union](https://cp-algorithms.com/graph/mst_kruskal_with_dsu.html).
+
+
+
+ 
+ *The 8 octants relative to a point S*
+
+
+The algorithm shown here was first presented in a paper from [H. Zhou, N. Shenoy, and W. Nichollos (2002)](https://ieeexplore.ieee.org/document/913303). There is also another know algorithm that uses a Divide and conquer approach by [J. Stolfi](https://www.academia.edu/15667173/On_computing_all_north_east_nearest_neighbors_in_the_L1_metric), which is also very interesting and only differ in the way they find the nearest neighbor in each octant. They both have the same complexity, but the one presented here is easier to implement and has a lower constant factor.
+
+First, let's understand why it is enough to consider only the nearest neighbor in each octant. The idea is to show that for a point $s$ and any two other points $p$ and $q$ in the same octant, $d(p, q) < \max(d(s, p), d(s, q))$. This is important, because it shows that if there was a MST where $s$ is connected to both $p$ and $q$, we could erase one of these edges and add the edge $(p,q)$, which would decrease the total cost. To prove this, we assume without loss of generality that $p$ and $q$ are in the octanct $R_1$, which is defined by: $x_s \leq x$ and $x_s - y_s > x - y$, and then do some casework. The image below give some intuition on why this is true.
+
+
+ *The 8 octants relative to a point S*
+
+ 
+ *Intuitively, the limitation of the octant makes it impossible that $p$ and $q$ are both closer to $s$ than to each other*
+
+
+
+Therefore, the main question is how to find the nearest neighbor in each octant for every single of the $n$ points.
+
+## Nearest Neighbor in each Octant in O(n log n)
+
+For simplicity we focus on the NNE octant ($R_1$ in the image above). All other directions can be found with the same algorithm by rotating the input.
+
+We will use a sweep-line approach. We process the points from south-west to north-east, that is, by non-decreasing $x + y$. We also keep a set of points which don't have their nearest neighbor yet, which we call "active set". We add the images below to help visualize the algorithm.
+
+
+ *Intuitively, the limitation of the octant makes it impossible that $p$ and $q$ are both closer to $s$ than to each other*
+
+ 
+ *In black with an arrow you can see the direction of the line-sweep. All the points below this lines are in the active set, and the points above are still not processed. In green we see the points which are in the octant of the processed point. In red the points that are not in the searched octant.*
+
+
+
+ *In black with an arrow you can see the direction of the line-sweep. All the points below this lines are in the active set, and the points above are still not processed. In green we see the points which are in the octant of the processed point. In red the points that are not in the searched octant.*
+
+ 
+ *In this image we see the active set after processing the point $p$. Note that the $2$ green points of the previous image had $p$ in its north-north-east octant and are not in the active set anymore, because they already found their nearest neighbor.*
+
+
+When we add a new point point $p$, for every point $s$ that has it in its octant we can safely assign $p$ as the nearest neighbor. This is true because their distance is $d(p,s) = |x_p - x_s| + |y_p - y_s| = (x_p + y_p) - (x_s + y_s)$, because $p$ is in the north-north-east octant. As all the next points will not have a smaller value of $x + y$ because of the sorting step, $p$ is guaranteed to have the smaller distance. We can then remove all such points from the active set, and finally add $p$ to the active set.
+
+The next question is how to efficiently find which points $s$ have $p$ in the north-north-east octant. That is, which points $s$ satisfy:
+
+- $x_s \leq x_p$
+- $x_p - y_p < x_s - y_s$
+
+Because no points in the active set are in the $R_1$ region of another, we also have that for two points $q_1$ and $q_2$ in the active set, $x_{q_1} \neq x_{q_2}$ and their ordering implies $x_{q_1} < x_{q_2} \implies x_{q_1} - y_{q_1} \leq x_{q_2} - y_{q_2}$.
+
+You can try to visualize this on the images above by thinking of the ordering of $x - y$ as a "sweep-line" that goes from the north-west to the south-east, so perpendicular to the one that is drawn.
+
+This means that if we keep the active set ordered by $x$ the candidates $s$ are consecutively placed. We can then find the largest $x_s \leq x_p$ and process the points in decreasing order of $x$ until the second condition $x_p - y_p < x_s - y_s$ breaks (we can actually allow that $x_p - y_p = x_s - y_s$ and that deals with the case of points with equal coordinates). Notice that because we remove from the set right after processing, this will have an amortized complexity of $O(n \log(n))$.
+ Now that we have the nearest point in the north-east direction, we rotate the points and repeat. It is possible to show that actually we also find this way the nearest point in the south-west direction, so we can repeat only 4 times, instead of 8.
+
+In summary we:
+
+- Sort the points by $x + y$ in non-decreasing order;
+- For every point, we iterate over the active set starting with the point with the largest $x$ such that $x \leq x_p$, and we break the loop if $x_p - y_p \geq x_s - y_s$. For every valid point $s$ we add the edge $(s,p, d(s,p))$ in our list;
+- We add the point $p$ to the active set;
+- Rotate the points and repeat until we iterate over all the octants.
+- Apply Kruskal algorithm in the list of edges to get the MST.
+
+Below you can find a implementation, based on the one from [KACTL](https://github.com/kth-competitive-programming/kactl/blob/main/content/geometry/ManhattanMST.h).
+
+```{.cpp file=manhattan_mst}
+struct point {
+ long long x, y;
+};
+
+// Returns a list of edges in the format (weight, u, v).
+// Passing this list to Kruskal algorithm will give the Manhattan MST.
+vector
+ *In this image we see the active set after processing the point $p$. Note that the $2$ green points of the previous image had $p$ in its north-north-east octant and are not in the active set anymore, because they already found their nearest neighbor.*
+
+ 
+
## Distance between two polygons
One of the most common applications of Minkowski sum is computing the distance between two convex polygons (or simply checking whether they intersect).
@@ -97,9 +99,9 @@ vector
+
+ 
+
+
+## Implementation
+The following implementation returns a vector of vertices for each face, outer face goes first.
+Inner faces are returned in counter-clockwise orders and the outer face is returned in clockwise order.
+
+For simplicity we find the next edge by doing binary search by angle.
+```{.cpp file=planar}
+struct Point {
+ int64_t x, y;
+
+ Point(int64_t x_, int64_t y_): x(x_), y(y_) {}
+
+ Point operator - (const Point & p) const {
+ return Point(x - p.x, y - p.y);
+ }
+
+ int64_t cross (const Point & p) const {
+ return x * p.y - y * p.x;
+ }
+
+ int64_t cross (const Point & p, const Point & q) const {
+ return (p - *this).cross(q - *this);
+ }
+
+ int half () const {
+ return int(y < 0 || (y == 0 && x < 0));
+ }
+};
+
+std::vector
+
+ 
+
+
+## Implementation
+```{.cpp file=planar_implicit}
+using dbl = long double;
+
+const dbl eps = 1e-9;
+
+struct Point {
+ dbl x, y;
+
+ Point(){}
+ Point(dbl x_, dbl y_): x(x_), y(y_) {}
+
+ Point operator * (dbl d) const {
+ return Point(x * d, y * d);
+ }
+
+ Point operator + (const Point & p) const {
+ return Point(x + p.x, y + p.y);
+ }
+
+ Point operator - (const Point & p) const {
+ return Point(x - p.x, y - p.y);
+ }
+
+ dbl cross (const Point & p) const {
+ return x * p.y - y * p.x;
+ }
+
+ dbl cross (const Point & p, const Point & q) const {
+ return (p - *this).cross(q - *this);
+ }
+
+ dbl dot (const Point & p) const {
+ return x * p.x + y * p.y;
+ }
+
+ dbl dot (const Point & p, const Point & q) const {
+ return (p - *this).dot(q - *this);
+ }
+
+ bool operator < (const Point & p) const {
+ if (fabs(x - p.x) < eps) {
+ if (fabs(y - p.y) < eps) {
+ return false;
+ } else {
+ return y < p.y;
+ }
+ } else {
+ return x < p.x;
+ }
+ }
+
+ bool operator == (const Point & p) const {
+ return fabs(x - p.x) < eps && fabs(y - p.y) < eps;
+ }
+
+ bool operator >= (const Point & p) const {
+ return !(*this < p);
+ }
+};
+
+struct Line{
+ Point p[2];
+
+ Line(Point l, Point r){p[0] = l; p[1] = r;}
+ Point& operator [](const int & i){return p[i];}
+ const Point& operator[](const int & i)const{return p[i];}
+ Line(const Line & l){
+ p[0] = l.p[0]; p[1] = l.p[1];
+ }
+ Point getOrth()const{
+ return Point(p[1].y - p[0].y, p[0].x - p[1].x);
+ }
+ bool hasPointLine(const Point & t)const{
+ return std::fabs(p[0].cross(p[1], t)) < eps;
+ }
+ bool hasPointSeg(const Point & t)const{
+ return hasPointLine(t) && t.dot(p[0], p[1]) < eps;
+ }
+};
+
+std::vector
+This problem may arise when you need to locate some points in a Voronoi diagram or in some simple polygon. @@ -15,7 +15,9 @@ This problem may arise when you need to locate some points in a Voronoi diagram Firstly, for each query point $p\ (x_0, y_0)$ we want to find such an edge that if the point belongs to any edge, the point lies on the edge we found, otherwise this edge must intersect the line $x = x_0$ at some unique point $(x_0, y)$ where $y < y_0$ and this $y$ is maximum among all such edges. The following image shows both cases. -
+ 
+
We will solve this problem offline using the sweep line algorithm. Let's iterate over x-coordinates of query points and edges' endpoints in increasing order and keep a set of edges $s$. For each x-coordinate we will add some events beforehand.
@@ -27,7 +29,9 @@ Finally, for each query point we will add one _get_ event for its x-coordinate.
For each x-coordinate we will sort the events by their types in order (_vertical_, _get_, _remove_, _add_).
The following image shows all events in sorted order for each x-coordinate.
-
+
+ 
+
We will keep two sets during the sweep-line process.
A set $t$ for all non-vertical edges, and one set $vert$ especially for the vertical ones.
diff --git a/src/geometry/tangents-to-two-circles.md b/src/geometry/tangents-to-two-circles.md
index 4f4dfa4e9..546150a62 100644
--- a/src/geometry/tangents-to-two-circles.md
+++ b/src/geometry/tangents-to-two-circles.md
@@ -14,7 +14,9 @@ The described algorithm will also work in the case when one (or both) circles de
## The number of common tangents
The number of common tangents to two circles can be **0,1,2,3,4** and **infinite**.
Look at the images for different cases.
-
+
+ 
+
Here, we won't be considering **degenerate** cases, i.e *when the circles coincide (in this case they have infinitely many common tangents), or one circle lies inside the other (in this case they have no common tangents, or if the circles are tangent, there is one common tangent).*
diff --git a/src/geometry/vertical_decomposition.md b/src/geometry/vertical_decomposition.md
index 934db8044..96853994f 100644
--- a/src/geometry/vertical_decomposition.md
+++ b/src/geometry/vertical_decomposition.md
@@ -39,7 +39,9 @@ For simplicity we will show how to do this for an upper segment, the algorithm f
Here is a graphic representation of the three cases.
-
+
+ 
+
Finally we should remark on processing all the additions of $1$ or $-1$ on all stripes in $[x_1, x_2]$. For each addition of $w$ on $[x_1, x_2]$ we can create events $(x_1, w),\ (x_2, -w)$
and process all these events with a sweep line.
diff --git a/src/graph/01_bfs.md b/src/graph/01_bfs.md
index 15022c141..5d7aee6f4 100644
--- a/src/graph/01_bfs.md
+++ b/src/graph/01_bfs.md
@@ -7,7 +7,7 @@ tags:
It is well-known, that you can find the shortest paths between a single source and all other vertices in $O(|E|)$ using [Breadth First Search](breadth-first-search.md) in an **unweighted graph**, i.e. the distance is the minimal number of edges that you need to traverse from the source to another vertex.
We can interpret such a graph also as a weighted graph, where every edge has the weight $1$.
-If not all edges in graph have the same weight, that we need a more general algorithm, like [Dijkstra](dijkstra.md) which runs in $O(|V|^2 + |E|)$ or $O(|E| \log |V|)$ time.
+If not all edges in graph have the same weight, then we need a more general algorithm, like [Dijkstra](dijkstra.md) which runs in $O(|V|^2 + |E|)$ or $O(|E| \log |V|)$ time.
However if the weights are more constrained, we can often do better.
In this article we demonstrate how we can use BFS to solve the SSSP (single-source shortest path) problem in $O(|E|)$, if the weight of each edge is either $0$ or $1$.
diff --git a/src/graph/2SAT.md b/src/graph/2SAT.md
index f22da73f6..fbd1a6a35 100644
--- a/src/graph/2SAT.md
+++ b/src/graph/2SAT.md
@@ -14,7 +14,7 @@ Find an assignment of $a, b, c$ such that the following formula is true:
$$(a \lor \lnot b) \land (\lnot a \lor b) \land (\lnot a \lor \lnot b) \land (a \lor \lnot c)$$
-SAT is NP-complete, there is no known efficient solution known for it.
+SAT is NP-complete, there is no known efficient solution for it.
However 2SAT can be solved efficiently in $O(n + m)$ where $n$ is the number of variables and $m$ is the number of clauses.
## Algorithm:
@@ -39,7 +39,9 @@ b \Rightarrow a & \lnot b \Rightarrow \lnot a & b \Rightarrow \lnot a & c \Right
You can see the implication graph in the following image:
-
+
+ 
+
It is worth paying attention to the property of the implication graph:
if there is an edge $a \Rightarrow b$, then there also is an edge $\lnot b \Rightarrow \lnot a$.
@@ -59,7 +61,9 @@ The following image shows all strongly connected components for the example.
As we can check easily, neither of the four components contain a vertex $x$ and its negation $\lnot x$, therefore the example has a solution.
We will learn in the next paragraphs how to compute a valid assignment, but just for demonstration purposes the solution $a = \text{false}$, $b = \text{false}$, $c = \text{false}$ is given.
-
+
+ 
+
Now we construct the algorithm for finding the solution of the 2-SAT problem on the assumption that the solution exists.
@@ -102,64 +106,81 @@ Below is the implementation of the solution of the 2-SAT problem for the already
In the graph the vertices with indices $2k$ and $2k+1$ are the two vertices corresponding to variable $k$ with $2k+1$ corresponding to the negated variable.
```{.cpp file=2sat}
-int n;
-vector
+
+ 
+
The value of the flow of a network is the sum of all the flows that get produced in the source $s$, or equivalently to the sum of all the flows that are consumed by the sink $t$.
A **maximal flow** is a flow with the maximal possible value.
@@ -53,7 +55,9 @@ In the visualization with water pipes, the problem can be formulated in the foll
how much water can we push through the pipes from the source to the sink?
The following image shows the maximal flow in the flow network.
-
+
+ 
+
## Ford-Fulkerson method
@@ -66,7 +70,7 @@ We can create a **residual network** from all these edges, which is just a netwo
The Ford-Fulkerson method works as follows.
First, we set the flow of each edge to zero.
Then we look for an **augmenting path** from $s$ to $t$.
-An augmenting path is a simple path in the residual graph, i.e. along the edges whose residual capacity is positive.
+An augmenting path is a simple path in the residual graph where residual capacity is positive for all the edges along that path.
If such a path is found, then we can increase the flow along these edges.
We keep on searching for augmenting paths and increasing the flow.
Once an augmenting path doesn't exist anymore, the flow is maximal.
@@ -79,19 +83,30 @@ we update $f((u, v)) ~\text{+=}~ C$ and $f((v, u)) ~\text{-=}~ C$ for every edge
Here is an example to demonstrate the method.
We use the same flow network as above.
Initially we start with a flow of 0.
-
+
+
+
We can find the path $s - A - B - t$ with the residual capacities 7, 5, and 8.
Their minimum is 5, therefore we can increase the flow along this path by 5.
This gives a flow of 5 for the network.
-
+
+ 
+ 
+
Again we look for an augmenting path, this time we find $s - D - A - C - t$ with the residual capacities 4, 3, 3, and 5.
Therefore we can increase the flow by 3 and we get a flow of 8 for the network.
-
+
+
+ 
+ 
+
This time we find the path $s - D - C - B - t$ with the residual capacities 1, 2, 3, and 3, and hence, we increase the flow by 1.
-
+
+
+ 
+ 
+
This time we find the augmenting path $s - A - D - C - t$ with the residual capacities 2, 3, 1, and 2.
We can increase the flow by 1.
@@ -101,7 +116,10 @@ In the original flow network, we are not allowed to send any flow from $A$ to $D
But because we already have a flow of 3 from $D$ to $A$, this is possible.
The intuition of it is the following:
Instead of sending a flow of 3 from $D$ to $A$, we only send 2 and compensate this by sending an additional flow of 1 from $s$ to $A$, which allows us to send an additional flow of 1 along the path $D - C - t$.
-
+
+
+ 
+
+
Now, it is impossible to find an augmenting path between $s$ and $t$, therefore this flow of $10$ is the maximal possible.
We have found the maximal flow.
@@ -184,7 +202,7 @@ int maxflow(int s, int t) {
## Integral flow theorem ## { #integral-theorem}
-The theorem simply says, that if every capacity in the network is an integer, then the flow in each edge will be an integer in the maximal flow.
+The theorem says, that if every capacity in the network is an integer, then the size of the maximum flow is an integer, and there is a maximum flow such that the flow in each edge is an integer as well. In particular, Ford-Fulkerson method finds such a flow.
## Max-flow min-cut theorem
@@ -200,7 +218,9 @@ It says that the capacity of the maximum flow has to be equal to the capacity of
In the following image, you can see the minimum cut of the flow network we used earlier.
It shows that the capacity of the cut $\{s, A, D\}$ and $\{B, C, t\}$ is $5 + 3 + 2 = 10$, which is equal to the maximum flow that we found.
Other cuts will have a bigger capacity, like the capacity between $\{s, A\}$ and $\{B, C, D, t\}$ is $4 + 3 + 5 = 12$.
-
+
+
+ 
+
A minimum cut can be found after performing a maximum flow computation using the Ford-Fulkerson method.
One possible minimum cut is the following:
@@ -213,3 +233,4 @@ This partition can be easily found using [DFS](depth-first-search.md) starting a
- [CSES - Download Speed](https://cses.fi/problemset/task/1694)
- [CSES - Police Chase](https://cses.fi/problemset/task/1695)
- [CSES - School Dance](https://cses.fi/problemset/task/1696)
+- [CSES - Distinct Routes](https://cses.fi/problemset/task/1711)
diff --git a/src/graph/euler_path.md b/src/graph/euler_path.md
index 1846e8c87..a54906253 100644
--- a/src/graph/euler_path.md
+++ b/src/graph/euler_path.md
@@ -65,14 +65,14 @@ Reformulate the problem. Let the numbers written on the bottoms be the vertices
The program below searches for and outputs a Eulerian loop or path in a graph, or outputs $-1$ if it does not exist.
First, the program checks the degree of vertices: if there are no vertices with an odd degree, then the graph has an Euler cycle, if there are $2$ vertices with an odd degree, then in the graph there is only an Euler path (but no Euler cycle), if there are more than $2$ such vertices, then in the graph there is no Euler cycle or Euler path.
-To find the Euler path (not a cycle), let's do this: if $V1$ and $V2$ are two vertices of odd degree,then just add an edge $(V1, V2)$, in the resulting graph we find the Euler cycle (it will obviously exist), and then remove the "fictitious" edge $(V1, V2)$ from the answer.
+To find the Euler path (not a cycle), let's do this: if $V1$ and $V2$ are two vertices of odd degree, then just add an edge $(V1, V2)$, in the resulting graph we find the Euler cycle (it will obviously exist), and then remove the "fictitious" edge $(V1, V2)$ from the answer.
We will look for the Euler cycle exactly as described above (non-recursive version), and at the same time at the end of this algorithm we will check whether the graph was connected or not (if the graph was not connected, then at the end of the algorithm some edges will remain in the graph, and in this case we need to print $-1$).
Finally, the program takes into account that there can be isolated vertices in the graph.
Notice that we use an adjacency matrix in this problem.
Also this implementation handles finding the next with brute-force, which requires to iterate over the complete row in the matrix over and over.
A better way would be to store the graph as an adjacency list, and remove edges in $O(1)$ and mark the reversed edges in separate list.
-This way we can archive a $O(N)$ algorithm.
+This way we can achieve an $O(N)$ algorithm.
```cpp
int main() {
diff --git a/src/graph/finding-cycle.md b/src/graph/finding-cycle.md
index d629d2feb..af0ebda9f 100644
--- a/src/graph/finding-cycle.md
+++ b/src/graph/finding-cycle.md
@@ -129,5 +129,6 @@ void find_cycle() {
```
### Practice problems:
+- [AtCoder : Reachability in Functional Graph](https://atcoder.jp/contests/abc357/tasks/abc357_e)
- [CSES : Round Trip](https://cses.fi/problemset/task/1669)
- [CSES : Round Trip II](https://cses.fi/problemset/task/1678/)
diff --git a/src/graph/finding-negative-cycle-in-graph.md b/src/graph/finding-negative-cycle-in-graph.md
index 615c75c08..a9b87c7f9 100644
--- a/src/graph/finding-negative-cycle-in-graph.md
+++ b/src/graph/finding-negative-cycle-in-graph.md
@@ -19,6 +19,9 @@ Bellman-Ford algorithm allows you to check whether there exists a cycle of negat
The details of the algorithm are described in the article on the [Bellman-Ford](bellman_ford.md) algorithm.
Here we'll describe only its application to this problem.
+The standard implementation of Bellman-Ford looks for a negative cycle reachable from some starting vertex $v$ ; however, the algorithm can be modified to just looking for any negative cycle in the graph.
+For this we need to put all the distance $d[i]$ to zero and not infinity — as if we are looking for the shortest path from all vertices simultaneously; the validity of the detection of a negative cycle is not affected.
+
Do $N$ iterations of Bellman-Ford algorithm. If there were no changes on the last iteration, there is no cycle of negative weight in the graph. Otherwise take a vertex the distance to which has changed, and go from it via its ancestors until a cycle is found. This cycle will be the desired cycle of negative weight.
### Implementation
@@ -27,33 +30,33 @@ Do $N$ iterations of Bellman-Ford algorithm. If there were no changes on the las
struct Edge {
int a, b, cost;
};
-
-int n, m;
+
+int n;
vector
+
+ 
+
## Sample problems
@@ -183,3 +185,8 @@ int query(int a, int b) {
## Practice problems
- [SPOJ - QTREE - Query on a tree](https://www.spoj.com/problems/QTREE/)
+- [CSES - Path Queries II](https://cses.fi/problemset/task/2134)
+- [Codeforces - Subway Lines](https://codeforces.com/gym/101908/problem/L)
+- [Codeforces - Tree Queries](https://codeforces.com/contest/1254/problem/D)
+- [Codeforces - Tree or not Tree](https://codeforces.com/contest/117/problem/E)
+- [Codeforces - The Tree](https://codeforces.com/contest/1017/problem/G)
diff --git a/src/graph/hungarian-algorithm.md b/src/graph/hungarian-algorithm.md
new file mode 100644
index 000000000..1bc57c1f6
--- /dev/null
+++ b/src/graph/hungarian-algorithm.md
@@ -0,0 +1,313 @@
+---
+tags:
+ - Translated
+e_maxx_link: assignment_hungary
+---
+
+# Hungarian algorithm for solving the assignment problem
+
+## Statement of the assignment problem
+
+There are several standard formulations of the assignment problem (all of which are essentially equivalent). Here are some of them:
+
+- There are $n$ jobs and $n$ workers. Each worker specifies the amount of money they expect for a particular job. Each worker can be assigned to only one job. The objective is to assign jobs to workers in a way that minimizes the total cost.
+
+- Given an $n \times n$ matrix $A$, the task is to select one number from each row such that exactly one number is chosen from each column, and the sum of the selected numbers is minimized.
+
+- Given an $n \times n$ matrix $A$, the task is to find a permutation $p$ of length $n$ such that the value $\sum A[i]\left[p[i]\right]$ is minimized.
+
+- Consider a complete bipartite graph with $n$ vertices per part, where each edge is assigned a weight. The objective is to find a perfect matching with the minimum total weight.
+
+It is important to note that all the above scenarios are "**square**" problems, meaning both dimensions are always equal to $n$. In practice, similar "**rectangular**" formulations are often encountered, where $n$ is not equal to $m$, and the task is to select $\min(n,m)$ elements. However, it can be observed that a "rectangular" problem can always be transformed into a "square" problem by adding rows or columns with zero or infinite values, respectively.
+
+We also note that by analogy with the search for a **minimum** solution, one can also pose the problem of finding a **maximum** solution. However, these two problems are equivalent to each other: it is enough to multiply all the weights by $-1$.
+
+## Hungarian algorithm
+
+### Historical reference
+
+The algorithm was developed and published by Harold **Kuhn** in 1955. Kuhn himself gave it the name "Hungarian" because it was based on the earlier work by Hungarian mathematicians Dénes Kőnig and Jenő Egerváry.
++In 1957, James **Munkres** showed that this algorithm runs in (strictly) polynomial time, independently from the cost.
+Therefore, in literature, this algorithm is known not only as the "Hungarian", but also as the "Kuhn-Mankres algorithm" or "Mankres algorithm".
+However, it was recently discovered in 2006 that the same algorithm was invented **a century before Kuhn** by the German mathematician Carl Gustav **Jacobi**. His work, _About the research of the order of a system of arbitrary ordinary differential equations_, which was published posthumously in 1890, contained, among other findings, a polynomial algorithm for solving the assignment problem. Unfortunately, since the publication was in Latin, it went unnoticed among mathematicians. + +It is also worth noting that Kuhn's original algorithm had an asymptotic complexity of $\mathcal{O}(n^4)$, and only later Jack **Edmonds** and Richard **Karp** (and independently **Tomizawa**) showed how to improve it to an asymptotic complexity of $\mathcal{O}(n^3)$. + +### The $\mathcal{O}(n^4)$ algorithm + +To avoid ambiguity, we note right away that we are mainly concerned with the assignment problem in a matrix formulation (i.e., given a matrix $A$, you need to select $n$ cells from it that are in different rows and columns). We index arrays starting with $1$, i.e., for example, a matrix $A$ has indices $A[1 \dots n][1 \dots n]$. + +We will also assume that all numbers in matrix A are **non-negative** (if this is not the case, you can always make the matrix non-negative by adding some constant to all numbers). + +Let's call a **potential** two arbitrary arrays of numbers $u[1 \ldots n]$ and $v[1 \ldots n]$, such that the following condition is satisfied: + +$$u[i]+v[j]\leq A[i][j],\quad i=1\dots n,\ j=1\dots n$$ + +(As you can see, $u[i]$ corresponds to the $i$-th row, and $v[j]$ corresponds to the $j$-th column of the matrix). + +Let's call **the value $f$ of the potential** the sum of its elements: + +$$f=\sum_{i=1}^{n} u[i] + \sum_{j=1}^{n} v[j].$$ + +On one hand, it is easy to see that the cost of the desired solution $sol$ **is not less than** the value of any potential. + +!!! info "" + + **Lemma.** $sol\geq f.$ + +??? info "Proof" + + The desired solution of the problem consists of $n$ cells of the matrix $A$, so $u[i]+v[j]\leq A[i][j]$ for each of them. Since all the elements in $sol$ are in different rows and columns, summing these inequalities over all the selected $A[i][j]$, you get $f$ on the left side of the inequality, and $sol$ on the right side. + +On the other hand, it turns out that there is always a solution and a potential that turns this inequality into **equality**. The Hungarian algorithm described below will be a constructive proof of this fact. For now, let's just pay attention to the fact that if any solution has a cost equal to any potential, then this solution is **optimal**. + +Let's fix some potential. Let's call an edge $(i,j)$ **rigid** if $u[i]+v[j]=A[i][j].$ + +Recall an alternative formulation of the assignment problem, using a bipartite graph. Denote with $H$ a bipartite graph composed only of rigid edges. The Hungarian algorithm will maintain, for the current potential, **the maximum-number-of-edges matching** $M$ of the graph $H$. As soon as $M$ contains $n$ edges, then the solution to the problem will be just $M$ (after all, it will be a solution whose cost coincides with the value of a potential). + +Let's proceed directly to **the description of the algorithm**. + +**Step 1.** At the beginning, the potential is assumed to be zero ($u[i]=v[i]=0$ for all $i$), and the matching $M$ is assumed to be empty. + +**Step 2.** Further, at each step of the algorithm, we try, without changing the potential, to increase the cardinality of the current matching $M$ by one (recall that the matching is searched in the graph of rigid edges $H$). To do this, the usual [Kuhn Algorithm for finding the maximum matching in bipartite graphs](kuhn_maximum_bipartite_matching.md) is used. Let us recall the algorithm here. +All edges of the matching $M$ are oriented in the direction from the right part to the left one, and all other edges of the graph $H$ are oriented in the opposite direction. + +Recall (from the terminology of searching for matchings) that a vertex is called saturated if an edge of the current matching is adjacent to it. A vertex that is not adjacent to any edge of the current matching is called unsaturated. A path of odd length, in which the first edge does not belong to the matching, and for all subsequent edges there is an alternating belonging to the matching (belongs/does not belong) - is called an augmenting path. +From all unsaturated vertices in the left part, a [depth-first](depth-first-search.md) or [breadth-first](breadth-first-search.md) traversal is started. If, as a result of the search, it was possible to reach an unsaturated vertex of the right part, we have found an augmenting path from the left part to the right one. If we include odd edges of the path and remove the even ones in the matching (i.e. include the first edge in the matching, exclude the second, include the third, etc.), then we will increase the matching cardinality by one. + +If there was no augmenting path, then the current matching $M$ is maximal in the graph $H$. + +**Step 3.** If at the current step, it is not possible to increase the cardinality of the current matching, then a recalculation of the potential is performed in such a way that, at the next steps, there will be more opportunities to increase the matching. + +Denote by $Z_1$ the set of vertices of the left part that were visited during the last traversal of Kuhn's algorithm, and through $Z_2$ the set of visited vertices of the right part. + +Let's calculate the value $\Delta$: + +$$\Delta = \min_{i\in Z_1,\ j\notin Z_2} A[i][j]-u[i]-v[j].$$ + +!!! info "" + + **Lemma.** $\Delta > 0.$ + +??? info "Proof" + + Suppose $\Delta=0$. Then there exists a rigid edge $(i,j)$ with $i\in Z_1$ and $j\notin Z_2$. It follows that the edge $(i,j)$ must be oriented from the right part to the left one, i.e. $(i,j)$ must be included in the matching $M$. However, this is impossible, because we could not get to the saturated vertex $i$ except by going along the edge from j to i. So $\Delta > 0$. + +Now let's **recalculate the potential** in this way: + +- for all vertices $i\in Z_1$, do $u[i] \gets u[i]+\Delta$, + +- for all vertices $j\in Z_2$, do $v[j] \gets v[j]-\Delta$. + +!!! info "" + + **Lemma.** The resulting potential is still a correct potential. + +??? info "Proof" + + We will show that, after recalculation, $u[i]+v[j]\leq A[i][j]$ for all $i,j$. For all the elements of $A$ with $i\in Z_1$ and $j\in Z_2$, the sum $u[i]+v[j]$ does not change, so the inequality remains true. For all the elements with $i\notin Z_1$ and $j\in Z_2$, the sum $u[i]+v[j]$ decreases by $\Delta$, so the inequality is still true. For the other elements whose $i\in Z_1$ and $j\notin Z_2$, the sum increases, but the inequality is still preserved, since the value $\Delta$ is, by definition, the maximum increase that does not change the inequality. + +!!! info "" + + **Lemma.** The old matching $M$ of rigid edges is valid, i.e. all edges of the matching will remain rigid. + +??? info "Proof" + + For some rigid edge $(i,j)$ to stop being rigid as a result of a change in potential, it is necessary that equality $u[i] + v[j] = A[i][j]$ turns into inequality $u[i] + v[j] < A[i][j]$. However, this can happen only when $i \notin Z_1$ and $j \in Z_2$. But $i \notin Z_1$ implies that the edge $(i,j)$ could not be a matching edge. + +!!! info "" + + **Lemma.** After each recalculation of the potential, the number of vertices reachable by the traversal, i.e. $|Z_1|+|Z_2|$, strictly increases. + +??? info "Proof" + + First, note that any vertex that was reachable before recalculation, is still reachable. Indeed, if some vertex is reachable, then there is some path from reachable vertices to it, starting from the unsaturated vertex of the left part; since for edges of the form $(i,j),\ i\in Z_1,\ j\in Z_2$ the sum $u[i]+v[j]$ does not change, this entire path will be preserved after changing the potential. + Secondly, we show that after a recalculation, at least one new vertex will be reachable. This follows from the definition of $\Delta$: the edge $(i,j)$ which $\Delta$ refers to will become rigid, so vertex $j$ will be reachable from vertex $i$. + +Due to the last lemma, **no more than $n$ potential recalculations can occur** before an augmenting path is found and the matching cardinality of $M$ is increased. +Thus, sooner or later, a potential that corresponds to a perfect matching $M^*$ will be found, and $M^*$ will be the answer to the problem. +If we talk about the complexity of the algorithm, then it is $\mathcal{O}(n^4)$: in total there should be at most $n$ increases in matching, before each of which there are no more than $n$ potential recalculations, each of which is performed in time $\mathcal{O}(n^2)$. + +We will not give the implementation for the $\mathcal{O}(n^4)$ algorithm here, since it will turn out to be no shorter than the implementation for the $\mathcal{O}(n^3)$ one, described below. + +### The $\mathcal{O}(n^3)$ algorithm + +Now let's learn how to implement the same algorithm in $\mathcal{O}(n^3)$ (for rectangular problems $n \times m$, $\mathcal{O}(n^2m)$). + +The key idea is to **consider matrix rows one by one**, and not all at once. Thus, the algorithm described above will take the following form: + +1. Consider the next row of the matrix $A$. + +2. While there is no increasing path starting in this row, recalculate the potential. + +3. As soon as an augmenting path is found, propagate the matching along it (thus including the last edge in the matching), and restart from step 1 (to consider the next line). + +To achieve the required complexity, it is necessary to implement steps 2-3, which are performed for each row of the matrix, in time $\mathcal{O}(n^2)$ (for rectangular problems in $\mathcal{O}(nm)$). + +To do this, recall two facts proved above: + +- With a change in the potential, the vertices that were reachable by Kuhn's traversal will remain reachable. + +- In total, only $\mathcal{O}(n)$ recalculations of the potential could occur before an augmenting path was found. + +From this follow these **key ideas** that allow us to achieve the required complexity: + +- To check for the presence of an augmenting path, there is no need to start the Kuhn traversal again after each potential recalculation. Instead, you can make the Kuhn traversal in an **iterative form**: after each recalculation of the potential, look at the added rigid edges and, if their left ends were reachable, mark their right ends reachable as well and continue the traversal from them. + +- Developing this idea further, we can present the algorithm as follows: at each step of the loop, the potential is recalculated. Subsequently, a column that has become reachable is identified (which will always exist as new reachable vertices emerge after every potential recalculation). If the column is unsaturated, an augmenting chain is discovered. Conversely, if the column is saturated, the matching row also becomes reachable. + +- To quickly recalculate the potential (faster than the $\mathcal{O}(n^2)$ naive version), you need to maintain auxiliary minima for each of the columns: + +
$minv[j]=\min_{i\in Z_1} A[i][j]-u[i]-v[j].$
+ + It's easy to see that the desired value $\Delta$ is expressed in terms of them as follows: + +
$\Delta=\min_{j\notin Z_2} minv[j].$
+ + Thus, finding $\Delta$ can now be done in $\mathcal{O}(n)$. + + It is necessary to update the array $minv$ when new visited rows appear. This can be done in $\mathcal{O}(n)$ for the added row (which adds up over all rows to $\mathcal{O}(n^2)$). It is also necessary to update the array $minv$ when recalculating the potential, which is also done in time $\mathcal{O}(n)$ ($minv$ changes only for columns that have not yet been reached: namely, it decreases by $\Delta$). + +Thus, the algorithm takes the following form: in the outer loop, we consider matrix rows one by one. Each row is processed in time $\mathcal{O}(n^2)$, since only $\mathcal{O}(n)$ potential recalculations could occur (each in time $\mathcal{O}(n)$), and the array $minv$ is maintained in time $\mathcal{O}(n^2)$; Kuhn's algorithm will work in time $\mathcal{O}(n^2)$ (since it is presented in the form of $\mathcal{O}(n)$ iterations, each of which visits a new column). + +The resulting complexity is $\mathcal{O}(n^3)$ or, if the problem is rectangular, $\mathcal{O}(n^2m)$. + +## Implementation of the Hungarian algorithm + +The implementation below was developed by **Andrey Lopatin** several years ago. It is distinguished by amazing conciseness: the entire algorithm consists of **30 lines of code**. + +The implementation finds a solution for the rectangular matrix $A[1\dots n][1\dots m]$, where $n\leq m$. The matrix is 1-based for convenience and code brevity: this implementation introduces a dummy zero row and zero column, which allows us to write many cycles in a general form, without additional checks. + +Arrays $u[0 \ldots n]$ and $v[0 \ldots m]$ store potential. Initially, they are set to zero, which is consistent with a matrix of zero rows (Note that it is unimportant for this implementation whether or not the matrix $A$ contains negative numbers). + +The array $p[0 \ldots m]$ contains a matching: for each column $j = 1 \ldots m$, it stores the number $p[j]$ of the selected row (or $0$ if nothing has been selected yet). For the convenience of implementation, $p[0]$ is assumed to be equal to the number of the current row. + +The array $minv[1 \ldots m]$ contains, for each column $j$, the auxiliary minima necessary for a quick recalculation of the potential, as described above. + +The array $way[1 \ldots m]$ contains information about where these minimums are reached so that we can later reconstruct the augmenting path. Note that, to reconstruct the path, it is sufficient to store only column values, since the row numbers can be taken from the matching (i.e., from the array $p$). Thus, $way[j]$, for each column $j$, contains the number of the previous column in the path (or $0$ if there is none). + +The algorithm itself is an outer **loop through the rows of the matrix**, inside which the $i$-th row of the matrix is considered. The first _do-while_ loop runs until a free column $j0$ is found. Each iteration of the loop marks visited a new column with the number $j0$ (calculated at the last iteration; and initially equal to zero - i.e. we start from a dummy column), as well as a new row $i0$ - adjacent to it in the matching (i.e. $p[j0]$; and initially when $j0=0$ the $i$-th row is taken). Due to the appearance of a new visited row $i0$, you need to recalculate the array $minv$ and $\Delta$ accordingly. If $\Delta$ is updated, then the column $j1$ becomes the minimum that has been reached (note that with such an implementation $\Delta$ could turn out to be equal to zero, which means that the potential cannot be changed at the current step: there is already a new reachable column). After that, the potential and the $minv$ array are recalculated. At the end of the "do-while" loop, we found an augmenting path ending in a column $j0$ that can be "unrolled" using the ancestor array $way$. + +The constant INF is "infinity", i.e. some number, obviously greater than all possible numbers in the input matrix $A$. + +```{.cpp file=hungarian} +vector
+ To solve it, we simply build an assignment problem, putting the number "infinity" in place of the missing edges. After that, we solve the problem with the Hungarian algorithm, and remove edges of infinite weight from the answer (they could enter the answer if the problem does not have a solution in the form of a perfect matching). + +- Given a bipartite graph, it is required to find in it **the maximum matching with the maximum weight**.
+ The solution is again obvious, all weights must be multiplied by minus one. + +- The task of **detecting moving objects in images**: two images were taken, as a result of which two sets of coordinates were obtained. It is required to correlate the objects in the first and second images, i.e. determine for each point of the second image, which point of the first image it corresponded to. In this case, it is required to minimize the sum of distances between the compared points (i.e., we are looking for a solution in which the objects have taken the shortest path in total).
+ To solve, we simply build and solve an assignment problem, where the weights of the edges are the Euclidean distances between points. + +- The task of **detecting moving objects by locators**: there are two locators that can't determine the position of an object in space, but only its direction. Both locators (located at different points) received information in the form of $n$ such directions. It is required to determine the position of objects, i.e. determine the expected positions of objects and their corresponding pairs of directions in such a way that the sum of distances from objects to direction rays is minimized.
+ Solution: again, we simply build and solve the assignment problem, where the vertices of the left part are the $n$ directions from the first locator, the vertices of the right part are the $n$ directions from the second locator, and the weights of the edges are the distances between the corresponding rays. + +- Covering a **directed acyclic graph with paths**: given a directed acyclic graph, it is required to find the smallest number of paths (if equal, with the smallest total weight) so that each vertex of the graph lies in exactly one path.
+ The solution is to build the corresponding bipartite graph from the given graph and find the maximum matching of the minimum weight in it. See separate article for more details. + +- **Tree coloring book**. Given a tree in which each vertex, except for leaves, has exactly $k-1$ children. It is required to choose for each vertex one of the $k$ colors available so that no two adjacent vertices have the same color. In addition, for each vertex and each color, the cost of painting this vertex with this color is known, and it is required to minimize the total cost.
+ To solve this problem, we use dynamic programming. Namely, let's learn how to calculate the value $d[v][c]$, where $v$ is the vertex number, $c$ is the color number, and the value $d[v][c]$ itself is the minimum cost needed to color all the vertices in the subtree rooted at $v$, and the vertex $v$ itself with color $c$. To calculate such a value $d[v][c]$, it is necessary to distribute the remaining $k-1$ colors among the children of the vertex $v$, and for this, it is necessary to build and solve the assignment problem (in which the vertices of the left part are colors, the vertices of the right part are children, and the weights of the edges are the corresponding values of $d$).
+ Thus, each value $d[v][c]$ is calculated using the solution of the assignment problem, which ultimately gives the asymptotic $\mathcal{O}(nk^4)$. + +- If, in the assignment problem, the weights are not on the edges, but on the vertices, and only **on the vertices of the same part**, then it's not necessary to use the Hungarian algorithm: just sort the vertices by weight and run the usual [Kuhn algorithm](kuhn_maximum_bipartite_matching.md) (for more details, see a [separate article](http://e-maxx.ru/algo/vertex_weighted_matching)). + +- Consider the following **special case**. Let each vertex of the left part be assigned some number $\alpha[i]$, and each vertex of the right part $\beta[j]$. Let the weight of any edge $(i,j)$ be equal to $\alpha[i]\cdot \beta[j]$ (the numbers $\alpha[i]$ and $\beta[j]$ are known). Solve the assignment problem.
+ To solve it without the Hungarian algorithm, we first consider the case when both parts have two vertices. In this case, as you can easily see, it is better to connect the vertices in the reverse order: connect the vertex with the smaller $\alpha[i]$ to the vertex with the larger $\beta[j]$. This rule can be easily generalized to an arbitrary number of vertices: you need to sort the vertices of the first part in increasing order of $\alpha[i]$ values, the second part in decreasing order of $\beta[j]$ values, and connect the vertices in pairs in that order. Thus, we obtain a solution with complexity of $\mathcal{O}(n\log n)$. + +- **The Problem of Potentials**. Given a matrix $A[1 \ldots n][1 \ldots m]$, it is required to find two arrays $u[1 \ldots n]$ and $v[1 \ldots m]$ such that, for any $i$ and $j$, $u[i] + v[j] \leq a[i][j]$ and the sum of elements of arrays $u$ and $v$ is maximum.
+ Knowing the Hungarian algorithm, the solution to this problem will not be difficult: the Hungarian algorithm just finds such a potential $u, v$ that satisfies the condition of the problem. On the other hand, without knowledge of the Hungarian algorithm, it seems almost impossible to solve such a problem. + + !!! info "Remark" + + This task is also called the **dual problem** of the assignment problem: minimizing the total cost of the assignment is equivalent to maximizing the sum of the potentials. + +## Literature + +- [Ravindra Ahuja, Thomas Magnanti, James Orlin. Network Flows [1993]](https://books.google.it/books/about/Network_Flows.html?id=rFuLngEACAAJ&redir_esc=y) + +- [Harold Kuhn. The Hungarian Method for the Assignment Problem [1955]](https://link.springer.com/chapter/10.1007/978-3-540-68279-0_2) + +- [James Munkres. Algorithms for Assignment and Transportation Problems [1957]](https://www.jstor.org/stable/2098689) + +## Practice Problems + +- [UVA - Crime Wave - The Sequel](http://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=1687) + +- [UVA - Warehouse](http://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=1829) + +- [SGU - Beloved Sons](http://acm.sgu.ru/problem.php?contest=0&problem=210) + +- [UVA - The Great Wall Game](http://livearchive.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=1277) + +- [UVA - Jogging Trails](http://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=1237) diff --git a/src/graph/kuhn_maximum_bipartite_matching.md b/src/graph/kuhn_maximum_bipartite_matching.md index 64647cb5f..25b835ca1 100644 --- a/src/graph/kuhn_maximum_bipartite_matching.md +++ b/src/graph/kuhn_maximum_bipartite_matching.md @@ -15,10 +15,14 @@ that no selected edge shares a vertex with any other selected edge. ### Required Definitions * A **matching** $M$ is a set of pairwise non-adjacent edges of a graph (in other words, no more than one edge from the set should be incident to any vertex of the graph $M$). -The **cardinality** of a matching is the number of edges in it. The maximum (or largest) matching is a matching whose cardinality is maximum among all possible matchings -in a given graph. All those vertices that have an adjacent edge from the matching (i.e., which have degree exactly one in the subgraph formed by $M$) are called **saturated** +The **cardinality** of a matching is the number of edges in it. +All those vertices that have an adjacent edge from the matching (i.e., which have degree exactly one in the subgraph formed by $M$) are called **saturated** by this matching. +* A **maximal matching** is a matching $M$ of a graph $G$ that is not a subset of any other matching. + +* A **maximum matching** (also known as maximum-cardinality matching) is a matching that contains the largest possible number of edges. Every maximum matching is a maximal matching. + * A **path** of length $k$ here means a *simple* path (i.e. not containing repeated vertices or edges) containing $k$ edges, unless specified otherwise. * An **alternating path** (in a bipartite graph, with respect to some matching) is a path in which the edges alternately belong / do not belong to the matching. diff --git a/src/graph/lca.md b/src/graph/lca.md index db36c18a0..f577d4b2f 100644 --- a/src/graph/lca.md +++ b/src/graph/lca.md @@ -30,7 +30,9 @@ So the $\text{LCA}(v_1, v_2)$ can be uniquely determined by finding the vertex w Let's illustrate this idea. Consider the following graph and the Euler tour with the corresponding heights: -
+ 
+
$$\begin{array}{|l|c|c|c|c|c|c|c|c|c|c|c|c|c|}
\hline
diff --git a/src/graph/lca_binary_lifting.md b/src/graph/lca_binary_lifting.md
index e8f51dc94..e65d57d7f 100644
--- a/src/graph/lca_binary_lifting.md
+++ b/src/graph/lca_binary_lifting.md
@@ -95,4 +95,6 @@ void preprocess(int root) {
```
## Practice Problems
+* [LeetCode - Kth Ancestor of a Tree Node](https://leetcode.com/problems/kth-ancestor-of-a-tree-node)
* [Codechef - Longest Good Segment](https://www.codechef.com/problems/LGSEG)
+* [HackerEarth - Optimal Connectivity](https://www.hackerearth.com/practice/algorithms/graphs/graph-representation/practice-problems/algorithm/optimal-connectivity-c6ae79ca/)
diff --git a/src/graph/lca_farachcoltonbender.md b/src/graph/lca_farachcoltonbender.md
index 9c9b9e6ac..506509359 100644
--- a/src/graph/lca_farachcoltonbender.md
+++ b/src/graph/lca_farachcoltonbender.md
@@ -22,7 +22,9 @@ The LCA of two nodes $u$ and $v$ is the node between the occurrences of $u$ and
In the following picture you can see a possible Euler-Tour of a graph and in the list below you can see the visited nodes and their heights.
-
+
+
+
$$\begin{array}{|l|c|c|c|c|c|c|c|c|c|c|c|c|c|}
\hline
@@ -118,7 +120,7 @@ int min_by_h(int i, int j) {
}
void precompute_lca(int root) {
- // get euler tour & indices of first occurences
+ // get euler tour & indices of first occurrences
first_visit.assign(n, -1);
height.assign(n, 0);
euler_tour.reserve(2 * n);
diff --git a/src/graph/min_cost_flow.md b/src/graph/min_cost_flow.md
index cdc2a6a22..b13716d40 100644
--- a/src/graph/min_cost_flow.md
+++ b/src/graph/min_cost_flow.md
@@ -17,7 +17,7 @@ Sometimes the task is given a little differently:
you want to find the maximum flow, and among all maximal flows we want to find the one with the least cost.
This is called the **minimum-cost maximum-flow problem**.
-Both these problems can be solved effectively with the algorithm of sucessive shortest paths.
+Both these problems can be solved effectively with the algorithm of successive shortest paths.
## Algorithm
@@ -42,7 +42,7 @@ the residual network contains only unsaturated edges (i.e. edges in which $F_{i
Now we can talk about the **algorithms** to compute the minimum-cost flow.
At each iteration of the algorithm we find the shortest path in the residual graph from $s$ to $t$.
-In contrary to Edmonds-Karp we look for the shortest path in terms of the cost of the path, instead of the number of edges.
+In contrast to Edmonds-Karp, we look for the shortest path in terms of the cost of the path instead of the number of edges.
If there doesn't exists a path anymore, then the algorithm terminates, and the stream $F$ is the desired one.
If a path was found, we increase the flow along it as much as possible (i.e. we find the minimal residual capacity $R$ of the path, and increase the flow by it, and reduce the back edges by the same amount).
If at some point the flow reaches the value $K$, then we stop the algorithm (note that in the last iteration of the algorithm it is necessary to increase the flow by only such an amount so that the final flow value doesn't surpass $K$).
@@ -160,3 +160,9 @@ int min_cost_flow(int N, vector
+
+ 
+ 
+
It is easy to see that any spanning tree will necessarily contain $n-1$ edges.
diff --git a/src/graph/push-relabel-faster.md b/src/graph/push-relabel-faster.md
index 0c000d207..24b9e0ce4 100644
--- a/src/graph/push-relabel-faster.md
+++ b/src/graph/push-relabel-faster.md
@@ -91,9 +91,6 @@ int max_flow(int s, int t)
}
}
- int max_flow = 0;
- for (int i = 0; i < n; i++)
- max_flow += flow[i][t];
- return max_flow;
+ return excess[t];
}
```
diff --git a/src/graph/push-relabel.md b/src/graph/push-relabel.md
index c1e75b83c..569a1198a 100644
--- a/src/graph/push-relabel.md
+++ b/src/graph/push-relabel.md
@@ -101,8 +101,7 @@ vector
+
+
+ 
+
The range minimum query `[l, r]` is equivalent to the lowest common ancestor query `[l', r']`, where `l'` is the node corresponding to the element `A[l]` and `r'` the node corresponding to the element `A[r]`.
Indeed the node corresponding to the smallest element in the range has to be an ancestor of all nodes in the range, therefor also from `l'` and `r'`.
@@ -33,7 +35,9 @@ And is also has to be the lowest ancestor, because otherwise `l'` and `r'` would
In the following image you can see the LCA queries for the RMQ queries `[1, 3]` and `[5, 9]`.
In the first query the LCA of the nodes `A[1]` and `A[3]` is the node corresponding to `A[2]` which has the value 2, and in the second query the LCA of `A[5]` and `A[9]` is the node corresponding to `A[8]` which has the value 3.
-
+
+ 
+
Such a tree can be built in $O(N)$ time and the Farach-Colton and Benders algorithm can preprocess the tree in $O(N)$ and find the LCA in $O(1)$.
diff --git a/src/graph/search-for-connected-components.md b/src/graph/search-for-connected-components.md
index f18779055..cfc0f27ac 100644
--- a/src/graph/search-for-connected-components.md
+++ b/src/graph/search-for-connected-components.md
@@ -20,43 +20,88 @@ Given an undirected graph $G$ with $n$ nodes and $m$ edges. We are required to f
``` cpp
int n;
-vector
++
+ 
+ 
+
*In the image left is the MST and right is the second best MST.* - +
In the given graph suppose we root the MST at the blue vertex on the top, and then run our algorithm by start picking the edges not in MST.
diff --git a/src/graph/strong-orientation.md b/src/graph/strong-orientation.md
index 50397b287..69a7b37d9 100644
--- a/src/graph/strong-orientation.md
+++ b/src/graph/strong-orientation.md
@@ -78,7 +78,7 @@ void find_bridges(int v) {
bridge_cnt++;
}
} else {
- low[v] = min(low[v], low[nv]);
+ low[v] = min(low[v], tin[nv]);
}
}
}
diff --git a/src/graph/strongly-connected-components-tikzpicture/.gitignore b/src/graph/strongly-connected-components-tikzpicture/.gitignore
new file mode 100644
index 000000000..8164e021b
--- /dev/null
+++ b/src/graph/strongly-connected-components-tikzpicture/.gitignore
@@ -0,0 +1,3 @@
+*.log
+*.aux
+*.pdf
diff --git a/src/graph/strongly-connected-components-tikzpicture/cond_graph.svg b/src/graph/strongly-connected-components-tikzpicture/cond_graph.svg
new file mode 100644
index 000000000..507c1a071
--- /dev/null
+++ b/src/graph/strongly-connected-components-tikzpicture/cond_graph.svg
@@ -0,0 +1,138 @@
+
+
diff --git a/src/graph/strongly-connected-components-tikzpicture/cond_graph.tex b/src/graph/strongly-connected-components-tikzpicture/cond_graph.tex
new file mode 100644
index 000000000..3a62e5248
--- /dev/null
+++ b/src/graph/strongly-connected-components-tikzpicture/cond_graph.tex
@@ -0,0 +1,22 @@
+\documentclass[tikz]{standalone}
+\usepackage{xcolor}
+\usetikzlibrary{arrows,positioning,quotes,backgrounds,arrows.meta,bending,positioning,shapes,shapes.geometric}
+\begin{document}
+\begin{tikzpicture}
+ [scale=2.5,very thick,every node/.style={draw,inner sep=0.18cm,outer sep=1.5mm}, every path/.style={->}]
+
+ %\draw[help lines] (-1,-2) grid (7,2);
+ \node[rectangle,green,fill=green!20] (0) at (0,0) {\color{black}\textbf{\{0,7\}}};
+ \node[rectangle,red,fill=red!20] (1) at (1,0.5) {\color{black}\textbf{\{1,2,3,5,6\}}};
+ \node[rectangle,blue,fill=blue!20] (4) at (2,0) {\color{black}\textbf{\{4,9\}}};
+ \node[rectangle,yellow,fill=yellow!20] (8) at (1,-0.3) {\color{black}\textbf{\{8\}}};
+
+ \path (0) edge (1);
+ \path (0) edge (8);
+ \path (1) edge (4);
+ \path (8) edge (1);
+ \path (8) edge (4);
+
+\end{tikzpicture}
+
+\end{document}
diff --git a/src/graph/strongly-connected-components-tikzpicture/graph.svg b/src/graph/strongly-connected-components-tikzpicture/graph.svg
new file mode 100644
index 000000000..6dc60d3de
--- /dev/null
+++ b/src/graph/strongly-connected-components-tikzpicture/graph.svg
@@ -0,0 +1,166 @@
+
+
diff --git a/src/graph/strongly-connected-components-tikzpicture/graph.tex b/src/graph/strongly-connected-components-tikzpicture/graph.tex
new file mode 100644
index 000000000..c33dccec6
--- /dev/null
+++ b/src/graph/strongly-connected-components-tikzpicture/graph.tex
@@ -0,0 +1,50 @@
+\documentclass[tikz]{standalone}
+\usepackage{xcolor}
+\usetikzlibrary{arrows,positioning,quotes,backgrounds,arrows.meta,bending,positioning,shapes,shapes.geometric}
+\begin{document}
+\pgfdeclarelayer{background}
+\pgfsetlayers{background,main}
+\begin{tikzpicture}
+ [scale=1.3,very thick,every circle node/.style={draw,inner sep=0.18cm,outer sep=1.1mm,fill=brown!30}, every path/.style={->}]
+
+ %\draw[help lines] (-1,-2) grid (7,2);
+ \node[circle] (0) at (0,0) {\textbf0};
+ \node[circle] (1) at (1,1) {\textbf1};
+ \node[circle] (2) at (3,1) {\textbf2};
+ \node[circle] (3) at (5,1) {\textbf3};
+ \node[circle] (4) at (6,0) {\textbf4};
+ \node[circle] (5) at (4,0) {\textbf5};
+ \node[circle] (6) at (2,0) {\textbf6};
+ \node[circle] (7) at (1,-1) {\textbf7};
+ \node[circle] (8) at (3,-1) {\textbf8};
+ \node[circle] (9) at (5,-1) {\textbf9};
+
+ \path (0) edge (1);
+ \path (0) edge[bend left=20] (7);
+ \path (1) edge[loop below] (1);
+ \path (1) edge[bend left=20] (2);
+ \path (2) edge[bend left=20] (1);
+ \path (2) edge (5);
+ \path (3) edge (2);
+ \path (3) edge (4);
+ \path (4) edge[bend left=20] (9);
+ \path (5) edge (3);
+ \path (5) edge (6);
+ \path (5) edge (9);
+ \path (6) edge (2);
+ \path (7) edge[bend left=20] (0);
+ \path (7) edge (6);
+ \path (7) edge (8);
+ \path (8) edge (6);
+ \path (8) edge (9);
+ \path (9) edge[bend left=20] (4);
+
+ \begin{pgfonlayer}{background}
+ \draw[thick, red, fill=red!20!white] plot [smooth cycle,tension=0.65] coordinates{(0.65,1.3)(4.7,1.52)(5.4,0.7)(4.25,-.4)(3,-.33)(1.5,-.4)};
+ \draw[thick, green, fill=green!20!white] plot [smooth cycle,tension=0.7] coordinates{(0.2,0.4)(1.5,-0.9)(.9,-1.5)(-.45,-.1)};
+ \draw[thick, blue, fill=blue!20!white] plot [smooth cycle,tension=0.7] coordinates{(5.8,0.5)(4.45,-0.8)(5.1,-1.5)(6.6,-.1)};
+ \draw[thick, yellow, fill=yellow!20!white] plot [smooth cycle,tension=0.9] coordinates{(3,-0.5)(3.6,-1)(3,-1.5)(2.4,-1)};
+ \end{pgfonlayer}
+\end{tikzpicture}
+
+\end{document}
diff --git a/src/graph/strongly-connected-components-tikzpicture/info.txt b/src/graph/strongly-connected-components-tikzpicture/info.txt
new file mode 100644
index 000000000..8165f78ac
--- /dev/null
+++ b/src/graph/strongly-connected-components-tikzpicture/info.txt
@@ -0,0 +1,5 @@
+These are the pictures (graphs) used in the article about Strongly Connected Components and the Condensation graph.
+
+Compile the .tex file with pdflatex, and then use a command line tool like pdf2svg to convert the compiled pdf into an svg file. We want svg because it is scalable, i.e. still looks good when the user zooms in far.
+
+This svg can than be directly included in the markdown.
diff --git a/src/graph/strongly-connected-components.md b/src/graph/strongly-connected-components.md
index a09b36464..5fd7a525b 100644
--- a/src/graph/strongly-connected-components.md
+++ b/src/graph/strongly-connected-components.md
@@ -4,151 +4,143 @@ tags:
e_maxx_link: strong_connected_components
---
-# Finding strongly connected components / Building condensation graph
+# Strongly connected components and the condensation graph
## Definitions
-You are given a directed graph $G$ with vertices $V$ and edges $E$. It is possible that there are loops and multiple edges. Let's denote $n$ as number of vertices and $m$ as number of edges in $G$.
+Let $G=(V,E)$ be a directed graph with vertices $V$ and edges $E \subseteq V \times V$. We denote with $n=|V|$ the number of vertices and with $m=|E|$ the number of edges in $G$. It is easy to extend all definitions in this article to multigraphs, but we will not focus on that.
-**Strongly connected component** is a maximal subset of vertices $C$ such that any two vertices of this subset are reachable from each other, i.e. for any $u, v \in C$:
+A subset of vertices $C \subseteq V$ is called a **strongly connected component** if the following conditions hold:
-$$u \mapsto v, v \mapsto u$$
+- for all $u,v\in C$, if $u \neq v$ there exists a path from $u$ to $v$ and a path from $v$ to $u$, and
+- $C$ is maximal, in the sense that no vertex can be added without violating the above condition.
-where $\mapsto$ means reachability, i.e. existence of the path from first vertex to the second.
+We denote with $\text{SCC}(G)$ the set of strongly connected components of $G$. These strongly connected components do not intersect with each other, and cover all vertices in the graph. Thus, the set $\text{SCC}(G)$ is a partition of $V$.
-It is obvious, that strongly connected components do not intersect each other, i.e. this is a partition of all graph vertices. Thus we can give a definition of condensation graph $G^{SCC}$ as a graph containing every strongly connected component as one vertex. Each vertex of the condensation graph corresponds to the strongly connected component of graph $G$. There is an oriented edge between two vertices $C_i$ and $C_j$ of the condensation graph if and only if there are two vertices $u \in C_i, v \in C_j$ such that there is an edge in initial graph, i.e. $(u, v) \in E$.
+Consider this graph $G_\text{example}$, in which the strongly connected components are highlighted:
-The most important property of the condensation graph is that it is **acyclic**. Indeed, suppose that there is an edge between $C$ and $C'$, let's prove that there is no edge from $C'$ to $C$. Suppose that $C' \mapsto C$. Then there are two vertices $u' \in C$ and $v' \in C'$ such that $v' \mapsto u'$. But since $u$ and $u'$ are in the same strongly connected component then there is a path between them; the same for $v$ and $v'$. As a result, if we join these paths we have that $v \mapsto u$ and at the same time $u \mapsto v$. Therefore $u$ and $v$ should be at the same strongly connected component, so this is contradiction. This completes the proof.
+
+
+ *In the image left is the MST and right is the second best MST.* - +
+ 
+ 
+
Topological order can be **non-unique** (for example, if there exist three vertices $a$, $b$, $c$ for which there exist paths from $a$ to $b$ and from $a$ to $c$ but not paths from $b$ to $c$ or from $c$ to $b$).
The example graph also has multiple topological orders, a second topological order is the following:
-
+
+
+ 
+
A Topological order may **not exist** at all.
It only exists, if the directed graph contains no cycles.
-Otherwise because there is a contradiction: if there is a cycle containing the vertices $a$ and $b$, then $a$ needs to have a smaller index than $b$ (since you can reach $b$ from $a$) and also a bigger one (as you can reach $a$ from $b$).
+Otherwise, there is a contradiction: if there is a cycle containing the vertices $a$ and $b$, then $a$ needs to have a smaller index than $b$ (since you can reach $b$ from $a$) and also a bigger one (as you can reach $a$ from $b$).
The algorithm described in this article also shows by construction, that every acyclic directed graph contains at least one topological order.
-A common problem in which topological sorting occurs is the following. There are $n$ variables with unknown values. For some variables we know that one of them is less than the other. You have to check whether these constraints are contradictory, and if not, output the variables in ascending order (if several answers are possible, output any of them). It is easy to notice that this is exactly the problem of finding topological order of a graph with $n$ vertices.
+A common problem in which topological sorting occurs is the following. There are $n$ variables with unknown values. For some variables, we know that one of them is less than the other. You have to check whether these constraints are contradictory, and if not, output the variables in ascending order (if several answers are possible, output any of them). It is easy to notice that this is exactly the problem of finding the topological order of a graph with $n$ vertices.
## The Algorithm
-To solve this problem we will use [depth-first search](depth-first-search.md).
+To solve this problem, we will use [depth-first search](depth-first-search.md).
Let's assume that the graph is acyclic. What does the depth-first search do?
@@ -65,24 +65,26 @@ vector
+
+ 
+
++A visual representation of simulated annealing, searching for the maxima of this function with multiple local maxima. +
+
+ 
+
## Proofs
@@ -135,22 +137,76 @@ void build(int a = 0, int b = 1, int c = 1, int d = 0, int level = 1) {
The search algorithm was already described in the proof that all fractions appear in the tree, but we will repeat it here. The algorithm is a binary search algorithm. Initially we stand at the root of the tree and we compare our target with the current fraction. If they are the same we are done and stop the process. If our target is smaller we move to the left child, otherwise we move to the right child.
-Here is an implementation that returns the path to a given fraction $\frac{x}{y}$ as a sequence of `'L'` and `'R'` characters, meaning traversal to the left and right child respectively. This sequence of characters uniquely defines all positive fractions and is called the Stern-Brocot number system.
+### Naive search
+
+Here is an implementation that returns the path to a given fraction $\frac{p}{q}$ as a sequence of `'L'` and `'R'` characters, meaning traversal to the left and right child respectively. This sequence of characters uniquely defines all positive fractions and is called the Stern-Brocot number system.
```cpp
-string find(int x, int y, int a = 0, int b = 1, int c = 1, int d = 0) {
- int m = a + c, n = b + d;
- if (x == m && y == n)
- return "";
- if (x*n < y*m)
- return 'L' + find(x, y, a, b, m, n);
- else
- return 'R' + find(x, y, m, n, c, d);
+string find(int p, int q) {
+ int pL = 0, qL = 1;
+ int pR = 1, qR = 0;
+ int pM = 1, qM = 1;
+ string res;
+ while(pM != p || qM != q) {
+ if(p * qM < pM * q) {
+ res += 'L';
+ tie(pR, qR) = {pM, qM};
+ } else {
+ res += 'R';
+ tie(pL, qL) = {pM, qM};
+ }
+ tie(pM, qM) = pair{pL + pR, qL + qR};
+ }
+ return res;
}
```
Irrational numbers in the Stern-Brocot number system corresponds to infinite sequences of characters. Along the endless path towards the irrational number the algorithm will find reduced fractions with gradually increasing denominators that provides increasingly better approximations of the irrational number. So by taking a prefix of the infinite sequence approximations with any desired precision can be achieved. This application is important in watch-making, which explains why the tree was discovered in that domain.
+Note that for a fraction $\frac{p}{q}$, the length of the resulting sequence could be as large as $O(p+q)$, for example when the fraction is of form $\frac{p}{1}$. This means that the algorithm above **should not be used, unless this is an acceptable complexity**!
+
+### Logarithmic search
+
+Fortunately, it is possible to enhance the algorithm above to guarantee $O(\log (p+q))$ complexity. For this we should note that if the current boundary fractions are $\frac{p_L}{q_L}$ and $\frac{p_R}{q_R}$, then by doing $a$ steps to the right we move to the fraction $\frac{p_L + a p_R}{q_L + a q_R}$, and by doing $a$ steps to the left, we move to the fraction $\frac{a p_L + p_R}{a q_L + q_R}$.
+
+Therefore, instead of doing steps of `L` or `R` one by one, we can do $k$ steps in the same direction at once, after which we would switch to going into other direction, and so on. In this way, we can find the path to the fraction $\frac{p}{q}$ as its run-length encoding.
+
+As the directions alternate this way, we will always know which one to take. So, for convenience we may represent a path to a fraction $\frac{p}{q}$ as a sequence of fractions
+
+$$
+\frac{p_0}{q_0}, \frac{p_1}{q_1}, \frac{p_2}{q_2}, \dots, \frac{p_n}{q_n}, \frac{p_{n+1}}{q_{n+1}} = \frac{p}{q}
+$$
+
+such that $\frac{p_{k-1}}{q_{k-1}}$ and $\frac{p_k}{q_k}$ are the boundaries of the search interval on the $k$-th step, starting with $\frac{p_0}{q_0} = \frac{0}{1}$ and $\frac{p_1}{q_1} = \frac{1}{0}$. Then, after the $k$-th step we move to a fraction
+
+$$
+\frac{p_{k+1}}{q_{k+1}} = \frac{p_{k-1} + a_k p_k}{q_{k-1} + a_k q_k},
+$$
+
+where $a_k$ is a positive integer number. If you're familiar with [continued fractions](../algebra/continued-fractions.md), you would recognize that the sequence $\frac{p_i}{q_i}$ is the sequence of the convergent fractions of $\frac{p}{q}$ and the sequence $[a_1; a_2, \dots, a_{n}, 1]$ represents the continued fraction of $\frac{p}{q}$.
+
+This allows to find the run-length encoding of the path to $\frac{p}{q}$ in the manner which follows the algorithm for computing continued fraction representation of the fraction $\frac{p}{q}$:
+
+```cpp
+auto find(int p, int q) {
+ bool right = true;
+ vector
+
+ 
+
+
+Here we need to find out the point **C**, i.e the starting point of the cycle.
+
+## Proposed algorithm
+The algorithm is called **Floyd’s Cycle Algorithm or Tortoise And Hare algorithm**.
+In order to figure out the starting point of the cycle, we need to figure out of the the cycle even exists or not.
+So, it involved two steps:
+1. Figure out the presence of the cycle.
+2. Find out the starting point of the cycle.
+
+### Step 1: Presence of the cycle
+1. Take two pointers $slow$ and $fast$.
+2. Both of them will point to head of the linked list initially.
+3. $slow$ will move one step at a time.
+4. $fast$ will move two steps at a time. (twice as speed as $slow$ pointer).
+5. Check if at any point they point to the same node before any one(or both) reach null.
+6. If they point to any same node at any point of their journey, it would indicate that the cycle indeed exists in the linked list.
+7. If we get null, it would indicate that the linked list has no cycle.
+
+
+
+ 
+
+
+Now, that we have figured out that there is a cycle present in the linked list, for the next step we need to find out the starting point of cycle, i.e., **C**.
+### Step 2: Starting point of the cycle
+1. Reset the $slow$ pointer to the **head** of the linked list.
+2. Move both pointers one step at a time.
+3. The point they will meet at will be the starting point of the cycle.
+
+```java
+// Presence of cycle
+public boolean hasCycle(ListNode head) {
+ ListNode slow = head;
+ ListNode fast = head;
+
+ while(fast != null && fast.next != null){
+ slow = slow.next;
+ fast = fast.next.next;
+ if(slow==fast){
+ return true;
+ }
+ }
+
+ return false;
+}
+```
+
+```java
+// Assuming there is a cycle present and slow and fast are point to their meeting point
+slow = head;
+while(slow!=fast){
+ slow = slow.next;
+ fast = fast.next;
+}
+
+return slow; // the starting point of the cycle.
+```
+
+## Why does it work
+
+### Step 1: Presence of the cycle
+Since the pointer $fast$ is moving with twice as speed as $slow$, we can say that at any point of time, $fast$ would have covered twice as much distance as $slow$.
+We can also deduce that the difference between the distance covered by both of these pointers is increasing by $1$.
+```
+slow: 0 --> 1 --> 2 --> 3 --> 4 (distance covered)
+fast: 0 --> 2 --> 4 --> 6 --> 8 (distance covered)
+diff: 0 --> 1 --> 2 --> 3 --> 4 (difference between distance covered by both pointers)
+```
+Let $L$ denote the length of the cycle, and $a$ represent the number of steps required for the slow pointer to reach the entry of cycle. There exists a positive integer $k$ ($k > 0$) such that $k \cdot L \geq a$.
+When the slow pointer has moved $k \cdot L$ steps, and the fast pointer has covered $2 \cdot k \cdot L$ steps, both pointers find themselves within the cycle. At this point, there is a separation of $k \cdot L$ between them. Given that the cycle's length remains $L$, this signifies that they meet at the same point within the cycle, resulting in their encounter.
+
+### Step 2: Starting point of the cycle
+
+Lets try to calculate the distance covered by both of the pointers till they point they met within the cycle.
+
+
+
+ 
+
+
+$slowDist = a + xL + b$ , $x\ge0$
+
+$fastDist = a + yL + b$ , $y\ge0$
+
+- $slowDist$ is the total distance covered by slow pointer.
+- $fastDist$ is the total distance covered by fast pointer.
+- $a$ is the number of steps both pointers need to take to enter the cycle.
+- $b$ is the distance between **C** and **G**, i.e., distance between the starting point of cycle and meeting point of both pointers.
+- $x$ is the number of times the slow pointer has looped inside the cycle, starting from and ending at **C**.
+- $y$ is the number of times the fast pointer has looped inside the cycle, starting from and ending at **C**.
+
+$fastDist = 2 \cdot (slowDist)$
+
+$a + yL + b = 2(a + xL + b)$
+
+Resolving the formula we get:
+
+$a=(y-2x)L-b$
+
+where $y-2x$ is an integer
+
+This basically means that $a$ steps is same as doing some number of full loops in cycle and go $b$ steps backwards.
+Since the fast pointer already is $b$ steps ahead of the entry of cycle, if fast pointer moves another $a$ steps it will end up at the entry of the cycle.
+And since we let the slow pointer start at the start of the linked list, after $a$ steps it will also end up at the cycle entry. So, if they both move $a$ step they both will meet the entry of cycle.
+
+# Problems:
+- [Linked List Cycle (EASY)](https://leetcode.com/problems/linked-list-cycle/)
+- [Happy Number (Easy)](https://leetcode.com/problems/happy-number/)
+- [Find the Duplicate Number (Medium)](https://leetcode.com/problems/find-the-duplicate-number/)
+
diff --git a/src/others/tortoise_hare_algo.png b/src/others/tortoise_hare_algo.png
new file mode 100644
index 000000000..24d2f6882
Binary files /dev/null and b/src/others/tortoise_hare_algo.png differ
diff --git a/src/others/tortoise_hare_cycle_found.png b/src/others/tortoise_hare_cycle_found.png
new file mode 100644
index 000000000..f2f3d84aa
Binary files /dev/null and b/src/others/tortoise_hare_cycle_found.png differ
diff --git a/src/others/tortoise_hare_proof.png b/src/others/tortoise_hare_proof.png
new file mode 100644
index 000000000..bd62b57a2
Binary files /dev/null and b/src/others/tortoise_hare_proof.png differ
diff --git a/src/overrides/partials/content.html b/src/overrides/partials/content.html
index e0b148fb8..babde6339 100644
--- a/src/overrides/partials/content.html
+++ b/src/overrides/partials/content.html
@@ -1,9 +1,9 @@
{% if page.edit_url %}
-
+
{% include ".icons/material/pencil.svg" %}
diff --git a/src/overrides/partials/header.html b/src/overrides/partials/header.html
index c288bfafb..967330642 100644
--- a/src/overrides/partials/header.html
+++ b/src/overrides/partials/header.html
@@ -87,6 +87,13 @@
+
{% include "partials/source.html" %}
+
+
+
{% endif %}
{% if "navigation.tabs.sticky" in features %}
diff --git a/src/preview.md b/src/preview.md
index 775d4fe87..e0ffe3866 100644
--- a/src/preview.md
+++ b/src/preview.md
@@ -29,7 +29,7 @@
# Article Preview
-Information for contributors
+Information for contributors