Add Chinese Remainder Theorem

Signed-off-by: Bruno Freitas Tissei <bft15@inf.ufpr.br>

algorithms/math/chinese_remainder_theorem.cpp

+/// Chinese Remainder Theorem
+///
+/// Description:
+///   Given $t$ linear congruences in the format: \[ x \equiv a_i \pmod{m_i} \]
+/// the Chinese Remainder Theorem (CRT) finds $x$ that satisfies every 
+/// given congruence or states that there are no solutions. \par
+///   This implementation does not require all $m_i$ to be coprime, it works
+/// with any value and returns the smallest possible $x$; infinitely many 
+/// solutions can be obtained by incrementing or decrementing 
+/// $LCM(m_1,m_2,...,m_t)$ from $x$.
+///
+/// Time: O(t log LCM(m_1,m_2,...,m_t))
+/// Space: O(t)
+///
+/// Caution:
+///   - It is very easy to get overflow, since the $LCM$ is computed from
+/// all $m_i$, BigInt or Python should be considered if inputs are too large
+
+ll norm(ll a, ll b) { return (a + b) % b; }
+
+pair<ll,ll> crt_single(ll a, ll n, ll b, ll m) {
+  ans_t e = ext_gcd(n, m);
+
+  if ((a - b) % e.d != 0)
+    return {-1,-1}; // No solution
+
+  ll ans = norm(a + e.x*(b-a) / e.d % (m/e.d)*n, n *m / e.d);
+  ll lcm = (n*m) / e.d;
+  return {norm(ans, lcm), lcm};
+}
+
+ll crt(vector<ll> a, vector<ll> m) {
+  ll res = a[0];
+  ll lcm = m[0];
+
+  int t = s.size();
+  for (int i = 1; i < t; ++i) {
+    auto ss = crt(res, lcm, a[i], m[i]);
+    if (ss.fi == -1)
+      return -1; // No solution
+
+    res = ss.fi;
+    lcm = ss.se;
+  }
+
+  return res;
+}

algorithms/math/extended_euclidean.cpp

@@ -3,17 +3,10 @@
 /// Time: O(log min(a,b))
 /// Space: O(1)
 
-int ext_gcd(int a, int b, int &x, int &y) {
-  if (a == 0) {
-    x = 0, y = 1;
-    return b;
-  }
+struct ans_t { ll x, y, d; };
 
-  int x1, y1;
-  int g = ext_gcd(b % a, a, x1, y1);
-
-  x = y1 - (b / a) * x1;
-  y = x1;
-
-  return g;
+ans_t ext_gcd(ll a, ll b) {
+  if (a == 0) return {0, 1, b};
+  ans_t e = ext_gcd(b % a, a);
+  return {e.y - (b/a)*e.x, e.x, e.d};
 }

algorithms/math/linear_diophantine_equation.cpp

@@ -24,15 +24,13 @@ struct Diophantine {
   { init(); }
 
   void init() {
-    int w0, z0;
-    d = ext_gcd(a, b, w0, z0);
+    ans_t e = ext_gcd(a, b); d = e.d;
     if (c % d == 0) {
-      x0 = w0 * (c / d);
-      y0 = z0 * (c / d);
+      x0 = e.x * (c / d);
+      y0 = e.y * (c / d);
       has_solution = true;
-    } else {
+    } else
       has_solution = false;
-    }
   }
 
   ii get(int t) {

algorithms/math/modular_multiplicative_inverse.cpp

@@ -11,7 +11,6 @@ int fermat(int a, int m) {
 // Extended Euclidean Algorithm: Used when m 
 // and a are coprime
 int extended_euclidean(int a, int m) {
-  int x, y;
-  int g = ext_gcd(a, m, x, y);
-  return (x % m + m) % m;
+  ans_t g = ext_gcd(a, m);
+  return (g.x % m + m) % m;
 }

contests/ICPC_LA18/F.cpp

+#include <bits/stdc++.h>
+
+#define EPS 1e-6
+#define MOD 1000000007
+#define inf 0x3f3f3f3f
+#define llinf 0x3f3f3f3f3f3f3f3f
+
+#define fi first
+#define se second
+#define pb push_back
+#define ende '\n'
+
+#define all(x) (x).begin(), (x).end()
+#define rall(x) (x).rbegin(), (x).rend()
+#define mset(x, y) memset(&x, (y), sizeof(x))
+
+using namespace std; 
+
+using ll = long long;
+using ii = pair<ll,ll>;
+
+vector<ll> karatsuba(const vector<ll> &a, 
+    const vector<ll> &b) 
+{
+  int n = a.size();
+  vector<ll> res(n + n);
+
+  if (n <= 32) {
+    for (int i = 0; i < n; i++)
+      for (int j = 0; j < n; j++)
+        res[i + j] += a[i] * b[j];
+
+    return res;
+  }
+
+  int k = n >> 1;
+  vector<ll> a1(a.begin(), a.begin() + k);
+  vector<ll> a2(a.begin() + k, a.end());
+  vector<ll> b1(b.begin(), b.begin() + k);
+  vector<ll> b2(b.begin() + k, b.end());
+
+  vector<ll> a1b1 = karatsuba(a1, b1);
+  vector<ll> a2b2 = karatsuba(a2, b2);
+
+  for (int i = 0; i < k; i++)
+    a2[i] += a1[i];
+  for (int i = 0; i < k; i++)
+    b2[i] += b1[i];
+
+  vector<ll> r = karatsuba(a2, b2);
+  for (int i = 0; i < a1b1.size(); i++) 
+    r[i] -= a1b1[i];
+  for (int i = 0; i < a2b2.size(); i++) 
+    r[i] -= a2b2[i];
+
+  for (int i = 0; i < r.size(); i++) 
+    res[i + k] += r[i];
+  for (int i = 0; i < a1b1.size(); i++) 
+    res[i] += a1b1[i];
+  for (int i = 0; i < a2b2.size(); i++) 
+    res[i + n] += a2b2[i];
+
+  return res;
+}
+
+const int base = 1000000000;
+const int base_d = 9;
+
+struct BigInt {
+  int sign;
+  vector<int> num;
+
+  BigInt() : sign(1) {}
+  BigInt(ll x) { *this = x; }
+
+  void operator=(const BigInt &x) {
+    sign = x.sign;
+    num = x.num;
+  }
+
+  void operator=(ll x) {
+    sign = 1;
+    if (x < 0) sign = -1, x = -x;
+    for (; x > 0; x /= base)
+      pb(x % base);
+  }
+
+  BigInt operator+(const BigInt &x) const {
+    if (sign != x.sign) return *this - (-x);
+
+    int carry = 0;
+    BigInt res = x;
+
+    for (int i = 0; i < max(size(), x.size()) || carry; ++i) {
+      if (i == (int) res.size())
+        res.push_back(0);
+
+      res[i] += carry + (i < size() ? num[i] : 0);
+      carry = res[i] >= base;
+      if (carry) res[i] -= base;
+    }
+
+    return res;
+  }
+
+  BigInt operator-(const BigInt &x) const {
+    if (sign != x.sign) return *this + (-x);
+    if (abs() < x.abs()) return -(x - *this);
+
+    int carry = 0;
+    BigInt res = *this;
+
+    for (int i = 0; i < x.size() || carry; ++i) {
+      res[i] -= carry + (i < x.size() ? x[i] : 0);
+      carry = res[i] < 0;
+      if (carry) res[i] += base;
+    }
+
+    res.trim();
+    return res;
+  }
+
+  void operator*=(int x) {
+    if (x < 0) sign = -sign, x = -x;
+
+    int carry = 0;
+    for (int i = 0; i < size() || carry; ++i) {
+      if (i == size()) pb(0);
+      ll cur = num[i] * (ll) x + carry;
+
+      carry  = (int) (cur / base);
+      num[i] = (int) (cur % base);
+    }
+
+    trim();
+  }
+
+  BigInt operator*(int x) const {
+    BigInt res = *this;
+    res *= x;
+    return res;
+  }
+
+  friend pair<BigInt, BigInt> divmod(const BigInt &a1,
+      const BigInt &b1)
+  {
+    int norm = base / (b1.back() + 1);
+    BigInt a = a1.abs() * norm;
+    BigInt b = b1.abs() * norm;
+    BigInt q, r;
+    q.resize(a.size());
+
+    for (int i = a.size() - 1; i >= 0; i--) {
+      r *= base;
+      r += a[i];
+
+      int s1 = r.size() <= b.size() ? 0 : r[b.size()];
+      int s2 = r.size() <= b.size() - 1 ? 0 : r[b.size() - 1];
+      int d = ((ll) base * s1 + s2) / b.back();
+
+      r -= b * d;
+      while (r < 0) r += b, --d;
+      q[i] = d;
+    }
+
+    q.sign = a1.sign * b1.sign;
+    r.sign = a1.sign;
+    q.trim(); r.trim();
+
+    return make_pair(q, r / norm);
+  }
+
+  BigInt operator/(const BigInt &x) const {
+    return divmod(*this, x).fi;
+  }
+
+  BigInt operator%(const BigInt &x) const {
+    return divmod(*this, x).se;
+  }
+
+  void operator/=(int x) {
+    if (x < 0) sign = -sign, x = -x;
+
+    for (int i = size() - 1, rem = 0; i >= 0; --i) {
+      ll cur = num[i] + rem * (ll) base;
+      num[i] = (int) (cur / x);
+      rem = (int) (cur % x);
+    }
+
+    trim();
+  }
+
+  BigInt operator/(int x) const {
+    BigInt res = *this;
+    res /= x;
+    return res;
+  }
+
+  int operator%(int x) const {
+    if (x < 0) x = -x;
+
+    int m = 0;
+    for (int i = size() - 1; i >= 0; --i)
+      m = (num[i] + m * (ll) base) % x;
+
+    return m * sign;
+  }
+
+  void operator+=(const BigInt &x) { *this = *this + x; }
+  void operator-=(const BigInt &x) { *this = *this - x; }
+  void operator*=(const BigInt &x) { *this = *this * x; }
+  void operator/=(const BigInt &x) { *this = *this / x; }
+
+  bool operator<(const BigInt &x) const {
+    if (sign != x.sign)
+      return sign < x.sign;
+
+    if (size() != x.size())
+      return size() * sign < x.size() * x.sign;
+
+    for (int i = size() - 1; i >= 0; i--)
+      if (num[i] != x[i])
+        return num[i] * sign < x[i] * sign;
+
+    return false;
+  }
+
+  bool operator>(const BigInt &x) const {
+    return x < *this;
+  }
+  bool operator<=(const BigInt &x) const {
+    return !(x < *this);
+  }
+  bool operator>=(const BigInt &x) const {
+    return !(*this < x);
+  }
+  bool operator==(const BigInt &x) const {
+    return !(*this < x) && !(x < *this);
+  }
+  bool operator!=(const BigInt &x) const {
+    return *this < x || x < *this;
+  }
+
+  void trim() {
+    while (!empty() && !back()) pop_back();
+    if (empty()) sign = 1;
+  }
+
+  bool is_zero() const {
+    return empty() || (size() == 1 && !num[0]);
+  }
+
+  BigInt operator-() const {
+    BigInt res = *this;
+    res.sign = -sign;
+    return res;
+  }
+
+  BigInt abs() const {
+    BigInt res = *this;
+    res.sign *= res.sign;
+    return res;
+  }
+
+  ll to_long() const {
+    ll res = 0;
+    for (int i = size() - 1; i >= 0; i--)
+      res = res * base + num[i];
+    return res * sign;
+  }
+
+  void read(const string &s) {
+    sign = 1;
+    num.clear();
+
+    int pos = 0;
+    while (pos < s.size() &&
+        (s[pos] == '-' || s[pos] == '+')) {
+      if (s[pos] == '-')
+        sign = -sign;
+      ++pos;
+    }
+
+    for (int i = s.size() - 1; i >= pos; i -= base_d) {
+      int x = 0;
+      for (int j = max(pos, i - base_d + 1); j <= i; j++)
+        x = x * 10 + s[j] - '0';
+      num.push_back(x);
+    }
+
+    trim();
+  }
+
+  friend istream& operator>>(istream &stream, BigInt &v) {
+    string s; stream >> s;
+    v.read(s);
+    return stream;
+  }
+
+  friend ostream& operator<<(ostream &stream, const BigInt &x) {
+    if (x.sign == -1)
+      stream << '-';
+
+    stream << (x.empty() ? 0 : x.back());
+    for (int i = x.size() - 2; i >= 0; --i)
+      stream << setw(base_d) << setfill('0') << x.num[i];
+
+    return stream;
+  }
+
+  static vector<int> convert_base(
+      const vector<int> &a,
+      int oldd, int newd) {
+    vector<ll> p(max(oldd, newd) + 1);
+    p[0] = 1;
+    for (int i = 1; i < p.size(); i++)
+      p[i] = p[i - 1] * 10;
+
+    ll cur = 0;
+    int curd = 0;
+    vector<int> res;
+
+    for (int i = 0; i < a.size(); i++) {
+      cur += a[i] * p[curd];
+      curd += oldd;
+
+      while (curd >= newd) {
+        res.pb(int(cur % p[newd]));
+        cur /= p[newd];
+        curd -= newd;
+      }
+    }
+
+    res.pb((int) cur);
+    while (!res.empty() && !res.back())
+      res.pop_back();
+    return res;
+  }
+
+  BigInt operator*(const BigInt &x) const {
+    vector<int> a6 = convert_base(this->num, base_d, 6);
+    vector<int> b6 = convert_base(x.num, base_d, 6);
+
+    vector<ll> a(all(a6));
+    vector<ll> b(all(b6));
+
+    while (a.size() < b.size()) a.pb(0);
+    while (b.size() < a.size()) b.pb(0);
+    while (a.size() & (a.size() - 1))
+      a.pb(0), b.pb(0);
+
+    vector<ll> c = karatsuba(a, b);
+
+    BigInt res;
+    int carry = 0;
+    res.sign = sign * x.sign;
+
+    for (int i = 0; i < c.size(); i++) {
+      ll cur = c[i] + carry;
+      res.pb((int) (cur % 1000000));
+      carry = (int) (cur / 1000000);
+    }
+
+    res.num = convert_base(res.num, 6, base_d);
+    res.trim();
+    return res;
+  }
+
+  // Handles vector operations.
+  int back() const { return num.back(); }
+  bool empty() const { return num.empty(); }
+  size_t size() const { return num.size(); }
+
+  void pop_back() { num.pop_back(); }
+  void resize(int x) { num.resize(x); }
+  void push_back(int x) { num.push_back(x); }
+
+  int &operator[](int i) { return num[i]; }
+  int operator[](int i) const { return num[i]; }
+};
+
+using pb = pair<BigInt,BigInt>;
+
+template <typename T> 
+struct ans_t { T x, y, d; };
+
+template <typename T>
+ans_t<T> ext_gcd(T a, T b) {
+  if (a == 0) return {0, 1, b};
+  ans_t<T> e = ext_gcd(b % a, a);
+  return {e.y - (b/a) * e.x, e.x, e.d};
+}
+
+template <typename T>
+T norm(T a, T b) {
+  return (a + b) % b;
+}
+
+template <typename T>
+pb crt(T a, T n, T b, T m) {
+  ans_t<T> e = ext_gcd(n, m);
+
+  if ((a - b) % e.d != 0)
+    return {-1,-1};
+
+  T ans = norm(a + e.x*(b-a) / e.d % (m/e.d)*n, n *m / e.d);
+  T lcm = (n*m) / e.d;
+  return {norm(ans, lcm), lcm};
+}
+
+int main() {
+  ios::sync_with_stdio(0);
+  cin.tie(0);
+
+  int b, z; cin >> b >> z;
+  vector<vector<int>> mat(b, vector<int>(z+1));
+  for (int i = 0; i < b; ++i)
+    for (int j = 0; j <= z; ++j)
+      cin >> mat[i][j];
+
+  vector<vector<int>> time(z+1, vector<int>(501));
+  for (int i = 0; i < b; ++i) {
+    int curr = mat[i][0];
+    time[curr][0] |= (1 << i);
+    for (int t = 1; t <= 500; ++t) {
+      curr = mat[i][curr]; 
+      time[curr][t] |= (1 << i);
+    }
+  }
+
+  vector<vector<int>> a(z + 1, vector<int>(b, -1));
+  vector<vector<int>> m(z + 1, vector<int>(b, -1));
+  for (int i = 1; i <= z; ++i) {
+    for (int t = 0; t <= 500; ++t) {
+      if (time[i][t] == (1 << b) - 1)
+        return cout << i << " " << t << ende, 0;
+
+      for (int j = 0; j < b; ++j) {
+        if (time[i][t] & (1 << j)) {
+          if (a[i][j] == -1)
+            a[i][j] = t;
+          else if (m[i][j] == -1)
+            m[i][j] = t - a[i][j];
+        }
+      }
+    }
+  }
+
+  int zoo = 0;
+  BigInt ans = llinf*2;
+  for (int i = 1; i <= z; ++i) {
+    bool poss = true;
+    for (int j = 0; j < b; ++j)
+      if (a[i][j] == -1 || m[i][j] == -1) {
+        poss = false;
+        break;
+      }
+
+    if (!poss) continue;
+
+    BigInt res = a[i][0];
+    BigInt lcm = m[i][0];
+    for (int j = 1; j < b; ++j) {
+      pb ss = crt(res, lcm, BigInt(a[i][j]), BigInt(m[i][j]));
+      if (ss.fi == -1) {
+        poss = false;
+        break;
+      }
+
+      res = ss.fi;
+      lcm = ss.se;
+    }
+
+    if (!poss) continue;
+    else {
+      if (ans > res)
+        zoo = i, ans = res;
+    }
+  }
+
+  if (ans == llinf*2) cout << "*" << ende;
+  else cout << zoo << " " << ans << ende;
+  return 0;
+}