From 568886b27d1f4d8691e9eccaa000a06ee8c7d024 Mon Sep 17 00:00:00 2001
From: Ariel <ariel@arielherself.xyz>
Date: Mon, 11 Mar 2024 20:12:32 +0800
Subject: [PATCH] Update exgcd.cc

---
 number/exgcd.cc | 44 ++++++++++++++++++++++++++------------------
 1 file changed, 26 insertions(+), 18 deletions(-)

diff --git a/number/exgcd.cc b/number/exgcd.cc
index 4ca819e..d52476f 100644
--- a/number/exgcd.cc
+++ b/number/exgcd.cc
@@ -1,19 +1,23 @@
-
 namespace Exgcd {
+    template <typename T> T abs(T x) { return x < 0 ? -x : x; }
+
+    template <typename T>
     struct exgcd_solution_t {
-        ll x, y, gcd;
+        T x, y, gcd;
     };
 
+    template <typename T>
     struct diophantine_solution_t {
-        exgcd_solution_t x_min, y_min;
-        ll range;
+        exgcd_solution_t<T> x_min, y_min;
+        T range;
     };
 
     // solve `ax + by = gcd(a, b)`
-    optional<exgcd_solution_t> exgcd(ll a, ll b) {
+    template <typename T>
+    optional<exgcd_solution_t<T>> exgcd(T a, T b) {
         if (a < 0 || b < 0 || a == 0 && b == 0) return nullopt;
-        ll x, y, g;
-        function<void(ll, ll)> __exgcd = [&__exgcd, &x, &y, &g] (ll a, ll b) -> void {
+        T x, y, g;
+        function<void(T, T)> __exgcd = [&__exgcd, &x, &y, &g] (T a, T b) -> void {
             if (b == 0) {
                 g = a, x = 1, y = 0;
             } else {
@@ -26,7 +30,8 @@ namespace Exgcd {
         return {{ x, y, g }};
     };
 
-    optional<ll> inverse(ll a, ll b) {
+    template <typename T>
+    optional<T> inverse(T a, T b) {
         auto raw = exgcd(a, b);
         if (raw == nullopt || raw.value().gcd != 1) {
             return nullopt;
@@ -36,14 +41,15 @@ namespace Exgcd {
     }
 
     // solve { x = a_i (mod n_i) } if n_i's are coprime
-    optional<ll> crt(const vector<pll>& equations) {
-        ll prod = 1;
+    template <typename T>
+    optional<T> crt(const vector<pair<T, T>>& equations) {
+        T prod = 1;
         for (auto&& [a, n] : equations) {
             prod *= n;
         }
-        ll res = 0;
+        T res = 0;
         for (auto&& [a, n] : equations) {
-            ll m = prod / n;
+            T m = prod / n;
             auto m_rev = inverse(m, n);
             if (m_rev == nullopt) return nullopt;
             res = mod(res + a * mod(m * m_rev.value(), prod), prod);
@@ -51,22 +57,24 @@ namespace Exgcd {
         return res;
     }
 
-    // find minimal non-negative integral solutions of `ax + by = c`
-    optional<diophantine_solution_t> diophantine(ll a, ll b, ll c, bool force_positive = false) {
+    // find minimal non-negative integral solutions of `ax + by = c`. It's not guaranteed that the other variable is non-negative.
+    template <typename T>
+    optional<diophantine_solution_t<T>> diophantine(T a, T b, T c, bool force_positive = false) {
         if (a < 0 || b < 0 || a == 0 && b == 0) return nullopt;
         auto raw = exgcd(a, b).value();
         if (c % raw.gcd) {
             return nullopt;
         } else {
-            ll x = raw.x * c / raw.gcd, y = raw.y * c / raw.gcd;
-            ll kx = force_positive ? (x <= 0 ? (-x) * raw.gcd / b + 1 : 1 - (x + b / raw.gcd - 1) * raw.gcd / b) : (x <= 0 ? ((-x) + b / raw.gcd - 1) * raw.gcd / b : (- x * raw.gcd / b));
-            ll ky = force_positive ? (y <= 0 ? (- 1 - (-y) * raw.gcd / a) : (y + a / raw.gcd - 1) * raw.gcd / a - 1) : (y <= 0 ? (- ((-y) + a / raw.gcd - 1) * raw.gcd / a) : y * raw.gcd / a);
+            T x = raw.x * c / raw.gcd, y = raw.y * c / raw.gcd;
+            T kx = force_positive ? (x <= 0 ? (-x) * raw.gcd / b + 1 : 1 - (x + b / raw.gcd - 1) * raw.gcd / b) : (x <= 0 ? ((-x) + b / raw.gcd - 1) * raw.gcd / b : (- x * raw.gcd / b));
+            T ky = force_positive ? (y <= 0 ? (- 1 - (-y) * raw.gcd / a) : (y + a / raw.gcd - 1) * raw.gcd / a - 1) : (y <= 0 ? (- ((-y) + a / raw.gcd - 1) * raw.gcd / a) : y * raw.gcd / a);
             return {{ { x + b * kx / raw.gcd , y - a * kx / raw.gcd , raw.gcd }, { x + b * ky / raw.gcd , y - a * ky / raw.gcd, raw.gcd }, abs(kx - ky) + 1 }};
         }
     }
 
     // find the minimal non-negative integral solution of `ax = b (mod n)`
-    optional<ll> congruential(ll a, ll b, ll n) {
+    template <typename T>
+    optional<T> congruential(T a, T b, T n) {
         if (a == 0) {
             if (b != 0) return nullopt;
             return 0;