cp-algorithms
diff --git a/‎src/algebra/module-inverse.md
Lines changed: 26 additions & 29 deletions b/‎src/algebra/module-inverse.md
Lines changed: 26 additions & 29 deletions
diff --git a/‎test/test_modular_inverse.cpp
Lines changed: 28 additions & 0 deletions b/‎test/test_modular_inverse.cpp
Lines changed: 28 additions & 0 deletions
@@ -76,45 +76,42 @@ From these results, we can easily find the modular inverse using the [binary exp
 
 Even though this method is easier to understand than the method described in previous paragraph, in the case when $m$ is not a prime number, we need to calculate Euler phi function, which involves factorization of $m$, which might be very hard. If the prime factorization of $m$ is known, then the complexity of this method is $O(\log m)$.
 
-## Finding the modular inverse for every number modulo $m$ { #mod-inv-all-num data-toc-label="Finding the modular inverse for every number modulo m"}
+## Finding the modular inverse using Euclidean Division
 
-The problem is the following: 
-we want to compute the modular inverse for every number in the range $[1, m-1]$.
+Given that $m > i$ (or we can modulo to make it smaller in 1 step), according to [Euclidean Division](https://en.wikipedia.org/wiki/Euclidean_division)
 
-Applying the algorithms described in the previous sections, we can obtain a solution with complexity $O(m \log m)$. 
+$$m = k \cdot i + r$$
 
-Here we present a better algorithm with complexity $O(m)$.
-However for this specific algorithm we require that the modulus $m$ is prime.
+where $k = \left\lfloor \frac{m}{i} \right\rfloor$ and $r = m \bmod i$, then
 
-We denote by $\text{inv}[i]$ the modular inverse of $i$. Then for $i > 1$ the following equation is valid:
-
-$$\text{inv}[i] = - \left\lfloor \frac{m}{i} \right\rfloor \cdot \text{inv}[m \bmod i] \bmod m$$
+$$
+\begin{align*}
+& \implies & 0          & \equiv k \cdot i + r   & \mod m \\
+& \iff & r              & \equiv -k \cdot i      & \mod m \\
+& \iff & r \cdot i^{-1} & \equiv -k              & \mod m \\
+& \iff & i^{-1}         & \equiv -k \cdot r^{-1} & \mod m
+\end{align*}
+$$
 
-Thus the implementation is very simple:
+From there we can have the following recursive function (in C++) for computing the modular inverse for number $i$ with respect to module $m$
 
-```cpp
-inv[1] = 1;
-for(int i = 2; i < m; ++i)
-    inv[i] = m - (m/i) * inv[m%i] % m;
+```{.cpp file=modular_inverse_euclidean_division}
+int inv(int i) {
+  return i <= 1 ? i : m - (long long)(m/i) * inv(m % i) % m;
+}
 ```
 
-### Proof
-
-We have:
+The exact time complexity of the this recursion is not known. It's is somewhere between $O(\frac{\log m}{\log\log m})$ and $O(n^{\frac{1}{3} - \frac{2}{177} + \epsilon})$.
+See [On the length of Pierce expansions](https://arxiv.org/abs/2211.08374).
+In practice this implementation is fast, e.g. for the modulus $10^9 + 7$ it will always finish in less than 50 iterations.
 
-$$m \bmod i = m -  \left\lfloor \frac{m}{i} \right\rfloor \cdot i$$
+Applying this formula, we can also precompute the modular inverse for every number in the range $[1, m-1]$ in $O(m)$.
 
-Taking both sides modulo $m$ yields:
-
-$$m \bmod i \equiv - \left\lfloor \frac{m}{i} \right\rfloor \cdot i \mod m$$
-
-Multiply both sides by $i^{-1} \cdot (m \bmod i)^{-1}$ yields
-
-$$(m \bmod i) \cdot i^{-1} \cdot (m \bmod i)^{-1} \equiv -\left\lfloor \frac{m}{i} \right\rfloor \cdot i \cdot i^{-1} \cdot (m \bmod i)^{-1} \mod m,$$
-
-which simplifies to:
-
-$$i^{-1} \equiv -\left\lfloor \frac{m}{i} \right\rfloor \cdot (m \bmod i)^{-1} \mod m,$$
+```{.cpp file=modular_inverse_euclidean_division_all}
+inv[1] = 1;
+for(int i = 2; i < m; ++i)
+    inv[i] = m - (long long)(m/i) * inv[m%i] % m;
+```
 
 ## Finding the modular inverse for array of numbers modulo $m$
 
 
@@ -0,0 +1,28 @@
+#include <cassert>
+#include <initializer_list>
+#include <vector>
+using namespace std;
+
+namespace EuclideanDivision {
+int m;
+#include "modular_inverse_euclidean_division.h"
+} // namespace EuclideanDivision
+
+int main() {
+  for (int m : {307, 1'000'000'007}) {
+    for (int x : {1, 5, 100, 152, 299}) {
+      EuclideanDivision::m = m;
+      int x_inv = EuclideanDivision::inv(x);
+      assert((long long)x * x_inv % m == 1);
+    }
+  }
+
+  for (int m : {307, 12347}) {
+    vector<int> inv(m);
+#include "modular_inverse_euclidean_division_all.h"
+
+    for (int x = 1; x < m; x++) {
+      assert((long long)x * inv[x] % m == 1);
+    }
+  }
+}