Modular Multiplicative Inverse

Author: xirc

lib/cpalgo/algebra/README.md

@@ -109,4 +109,25 @@ We can compute the solution in `O(log(min(a,b)))`.
 - [Extended Euclidean Algorithm - Competitive Programming Algorithms](https://cp-algorithms.com/algebra/extended-euclid-algorithm.html)
 
 ### Challenges
-- [Extended Euclidean Algorithm - AOJ NTL1E](https://onlinejudge.u-aizu.ac.jp/courses/library/6/NTL/1/NTL_1_E)
+- [Extended Euclidean Algorithm - AOJ NTL1E](https://onlinejudge.u-aizu.ac.jp/courses/library/6/NTL/1/NTL_1_E)
+
+
+## Modular Multiplcative Inverse
+A *Modular Multiplicative Inverse* of an integer `a` is an integer `x` such that the product `ax` is congruent to 1 in modulus `m`.
+```math
+ax = 1 (mod m)
+x = a^{-1} (mod m)
+```
+
+- [Modular Multiplicative Inverse (using Fermat's little theorem and Binary Exponentation) | C++ code](modinv_binexp.hpp)
+- [Modular Multiplicative Inverse (using Extended Euclidean Algorithm) | C++ code](modinv_extgcd.hpp)
+- [All Modular Multiplicative Inverse | C++ code](modinvall.hpp)
+
+### References in English
+- [Modular multiplicative inverse - Wikipedia](https://en.wikipedia.org/wiki/Modular_multiplicative_inverse)
+- [Modular Inverse - Competitive Programming Algorithms](https://cp-algorithms.com/algebra/module-inverse.html)
+
+### Challenges
+- [D - 動的計画法](https://atcoder.jp/contests/abc024/tasks/abc024_d)
+- [E - Flatten](https://atcoder.jp/contests/abc152/tasks/abc152_e)
+- [E - Throne](https://atcoder.jp/contests/abc186/tasks/abc186_e)

lib/cpalgo/algebra/modinv_binexp.hpp

@@ -0,0 +1,36 @@
+#pragma once
+
+#include "binexp.hpp"
+
+
+// Modular Multiplicative Inverse
+//
+// (a ** -1) mod m
+//   where `m` is an prime number
+// undefined behavior if `m` is not an prime number
+//
+// Complexity:
+//   Time: O(log(m))
+//   Space: O(1)
+//
+// Algorithm:
+//   Fermat's little theorem and Binary Exponentation
+// { a ** (m - 1) = 1 } mod m
+// { a ** (m - 2) = a ** -1 } mod m
+//
+// NOTE:
+//   If `m` is not an prime but `a` and `m` is coprime,
+//   there is the Modular Multiplicative Inverse of `a`.
+//   Use an another solution using Extended Euclidean Algorithm.
+//
+//   We can solve it using Euler's theorem and Binary Exponentation, but it might be hard.
+//   Because we have to calculate Euler's Totient function that requires prime factorization of `m`.
+//   { a ** phi(m) = 1 } mod m
+//   { a ** (phi(m) - 1) = a ** -1 } mod m
+//     where phi(m) is Euler's Totient function.
+//
+int64_t modinv_binexp(int64_t a, int64_t m) {
+    return binexp<int64_t>(a, m - 2, 1, [&](auto const& lhs, auto const& rhs){
+        return lhs * rhs % m;
+    });
+}

lib/cpalgo/algebra/modinv_extgcd.hpp

@@ -0,0 +1,31 @@
+#pragma once
+
+#include "extgcd_iterative.hpp"
+
+
+// Modular Multiplicative Inverse
+//
+// returns (a ** -1) mod m
+//   where gcd(a,m) == 1
+// returns -1 if gcd(a,m) != 1
+//
+// Complexity:
+//   Time: O( log(min(a,b)) )
+//   Space: O(1)
+//
+// Algorithm:
+//   Extended Euclidean Algorithm
+// a * x + m * y = 1
+// (a * x + m * y = 1) mod m
+// (a * x = 1) mod m
+// (x = a ** -1) mod m
+//
+int64_t modinv(int64_t a, int64_t m) {
+    int64_t x, y;
+    auto g = extgcd(a, m, x, y);
+    if (g != 1) {
+        // There is no modular multiplicative inverse
+        return -1;
+    }
+    return (x % m + m) % m;
+}

lib/cpalgo/algebra/modinvall.hpp

@@ -0,0 +1,39 @@
+#pragma once
+
+#include <vector>
+
+
+// Modular Multiplicative Inverse for every number modulo `m`
+//   `m` must be an prime number
+// undefined behavior if `m` is not the prime number
+//
+// Complexity:
+//   Time: O(M)
+//   Space: O(M)
+//
+// Algorithm:
+//   m mod i = m - floor(m / i) * i
+//   m mod i = { m - floor(m / i) * i }   mod m
+//   m mod i = { - floor(m / i) * i }   mod m
+//   (m mod i) * (i ** -1) * ((m mod i) ** -1) =
+//       { - floor(m / i) * i * (i ** -1) * ((m mod i) ** 0) }   mod m
+//   i ** -1 =
+//     - floor(m / i) * { (m mod i) ** - 1 }   mod m
+//   i ** -1 =
+//     m - floor(m / i) * { (m mod i) ** -1 }   mod m
+//
+// NOTE:
+//  inverse[M % i] can be calculated using inverse[1..i-1]
+//  - i < M
+//  - M % i < i
+//  - M % i != 0 if M is the prime number
+//
+std::vector<int> modinvall(int M) {
+    std::vector<int> inverse(M, 0);
+    // inserve[0] is undefined
+    inverse[1] = 1;
+    for (int i = 2; i < M; ++i) {
+        inverse[i] = M - (M / i) * inverse[M % i] % M;
+    }
+    return inverse;
+}

lib/cpalgo/algebra/tests/modinv_binexp_test.cpp

@@ -0,0 +1,27 @@
+#include <bits/stdc++.h>
+#include "gtest/gtest.h"
+#include "cpalgo/algebra/modinv_binexp.hpp"
+
+TEST(ModinvBinexpTest, ComputeInverseInModulo7) {
+    EXPECT_EQ(1, modinv_binexp(1, 7));
+    EXPECT_EQ(4, modinv_binexp(2, 7));
+    EXPECT_EQ(5, modinv_binexp(3, 7));
+    EXPECT_EQ(2, modinv_binexp(4, 7));
+    EXPECT_EQ(3, modinv_binexp(5, 7));
+    EXPECT_EQ(6, modinv_binexp(6, 7));
+}
+
+TEST(ModinvBinexpTest, ComputeInverseInModulo13) {
+    EXPECT_EQ( 1, modinv_binexp( 1, 13));
+    EXPECT_EQ( 7, modinv_binexp( 2, 13));
+    EXPECT_EQ( 9, modinv_binexp( 3, 13));
+    EXPECT_EQ(10, modinv_binexp( 4, 13));
+    EXPECT_EQ( 8, modinv_binexp( 5, 13));
+    EXPECT_EQ(11, modinv_binexp( 6, 13));
+    EXPECT_EQ( 2, modinv_binexp( 7, 13));
+    EXPECT_EQ( 5, modinv_binexp( 8, 13));
+    EXPECT_EQ( 3, modinv_binexp( 9, 13));
+    EXPECT_EQ( 4, modinv_binexp(10, 13));
+    EXPECT_EQ( 6, modinv_binexp(11, 13));
+    EXPECT_EQ(12, modinv_binexp(12, 13));
+}

lib/cpalgo/algebra/tests/modinv_extgcd_test.cpp

@@ -0,0 +1,39 @@
+#include <bits/stdc++.h>
+#include "gtest/gtest.h"
+#include "cpalgo/algebra/modinv_extgcd.hpp"
+
+TEST(ModinvExtgcdTest, ComputeInverseInModulo7) {
+    EXPECT_EQ(1, modinv(1, 7));
+    EXPECT_EQ(4, modinv(2, 7));
+    EXPECT_EQ(5, modinv(3, 7));
+    EXPECT_EQ(2, modinv(4, 7));
+    EXPECT_EQ(3, modinv(5, 7));
+    EXPECT_EQ(6, modinv(6, 7));
+}
+
+TEST(ModinvExtgcdTest, ComputeInverseInModulo9) {
+    EXPECT_EQ( 1, modinv(1, 9));
+    EXPECT_EQ( 5, modinv(2, 9));
+    EXPECT_EQ(-1, modinv(3, 9));
+    EXPECT_EQ( 7, modinv(4, 9));
+    EXPECT_EQ( 2, modinv(5, 9));
+    EXPECT_EQ(-1, modinv(6, 9));
+    EXPECT_EQ( 4, modinv(7, 9));
+    EXPECT_EQ( 8, modinv(8, 9));
+}
+
+TEST(ModinvExtgcdTest, ComputeInverseInModulo13) {
+    EXPECT_EQ( 1, modinv( 1, 13));
+    EXPECT_EQ( 7, modinv( 2, 13));
+    EXPECT_EQ( 9, modinv( 3, 13));
+    EXPECT_EQ(10, modinv( 4, 13));
+    EXPECT_EQ( 8, modinv( 5, 13));
+    EXPECT_EQ(11, modinv( 6, 13));
+    EXPECT_EQ( 2, modinv( 7, 13));
+    EXPECT_EQ( 5, modinv( 8, 13));
+    EXPECT_EQ( 3, modinv( 9, 13));
+    EXPECT_EQ( 4, modinv(10, 13));
+    EXPECT_EQ( 6, modinv(11, 13));
+    EXPECT_EQ(12, modinv(12, 13));
+}
+

lib/cpalgo/algebra/tests/modinvall_test.cpp

@@ -0,0 +1,30 @@
+#include <bits/stdc++.h>
+#include "gtest/gtest.h"
+#include "cpalgo/algebra/modinvall.hpp"
+
+TEST(ModInvAllTest, ComputeAllInverseInModulo7) {
+    auto inv = modinvall(7);
+    EXPECT_EQ(1, inv[1]);
+    EXPECT_EQ(4, inv[2]);
+    EXPECT_EQ(5, inv[3]);
+    EXPECT_EQ(2, inv[4]);
+    EXPECT_EQ(3, inv[5]);
+    EXPECT_EQ(6, inv[6]);
+}
+
+TEST(ModinvExtgcdTest, ComputeAllInverseInModulo13) {
+    auto inv = modinvall(13);
+    EXPECT_EQ( 1,  inv[1]);
+    EXPECT_EQ( 7,  inv[2]);
+    EXPECT_EQ( 9,  inv[3]);
+    EXPECT_EQ(10,  inv[4]);
+    EXPECT_EQ( 8,  inv[5]);
+    EXPECT_EQ(11,  inv[6]);
+    EXPECT_EQ( 2,  inv[7]);
+    EXPECT_EQ( 5,  inv[8]);
+    EXPECT_EQ( 3,  inv[9]);
+    EXPECT_EQ( 4, inv[10]);
+    EXPECT_EQ( 6, inv[11]);
+    EXPECT_EQ(12, inv[12]);
+}
+