Merge pull request #114 from NyaanNyaan/add_misc

add misc

data-structure/parallel-union-find.hpp

@@ -31,6 +31,27 @@ struct ParallelUnionFind {
       });
     }
   }
+  // [l1, r1) と [l2, r2) を unite する
+  // f(x, y) : x に y をマージ
+  template <typename F>
+  void unite(int l1, int r1, int l2, int r2, const F& f) {
+    assert(0 <= l1 and l1 <= r1 and r1 <= n);
+    assert(0 <= l2 and l2 <= r2 and r2 <= n);
+    assert(r1 - l1 == r2 - l2);
+    while (1) {
+      if (seg.same(l1, r1, l2, r2)) break;
+      int ok = 0, ng = r1 - l1;
+      while (ok + 1 < ng) {
+        int m = (ok + ng) / 2;
+        (seg.same(l1, l1 + m, l2, l2 + m) ? ok : ng) = m;
+      }
+      uf.unite(l1 + ok, l2 + ok, [&](int x, int y) {
+        for (int z : uf.enumerate(y)) seg.update(z, x);
+        f(x, y);
+      });
+    }
+  }
+
   void unite(int l, int r) { unite(l, l + 1, r, r + 1); }
   int find(int i) { return uf.find(i); }
   int same(int l, int r) { return uf.same(l, r); }

math/enumerate-convex.hpp

@@ -0,0 +1,82 @@
+#pragma once
+
+#include <functional>
+#include <utility>
+#include <vector>
+using namespace std;
+
+#include "stern-brocot-tree.hpp"
+
+// 下向き凸包の頂点列挙
+// (xl, yl) 始点, x in [xl, xr]
+// inside(x, y) : (x, y) が凸包内部か？
+// candicate(x, y, c, d) : (x, y) が凸包外部にあるとする。
+// 凸包内部の点 (x + sc, y + sd) が存在すればそのような s を返す
+// 存在しなければ任意の値 (-1 でもよい) を返す
+template <typename Int>
+vector<pair<Int, Int>> enumerate_convex(
+    Int xl, Int yl, Int xr, const function<bool(Int, Int)>& inside,
+    const function<Int(Int, Int, Int, Int)>& candicate) {
+  assert(xl <= xr);
+
+  // inside かつ x <= xr
+  auto f = [&](Int x, Int y) { return x <= xr && inside(x, y); };
+
+  // (a, b) から (c, d) 方向に進めるだけ進む
+  auto go = [&](Int a, Int b, Int c, Int d) -> Int {
+    assert(f(a, b));
+    Int r = 1, s = 0;
+    while (f(a + r * c, b + r * d)) r *= 2;
+    while ((r /= 2) != 0) {
+      if (f(a + r * c, b + r * d)) s += r, a += r * c, b += r * d;
+    }
+    return s;
+  };
+
+  // (a, b) が out, (a + c * k, b + d * k) が in とする
+  // out の間進めるだけ進む
+  auto go2 = [&](Int a, Int b, Int c, Int d, Int k) {
+    assert(!inside(a, b) and inside(a + c * k, b + d * k));
+    Int ok = 0, ng = k;
+    while (ok + 1 < ng) {
+      Int m = (ok + ng) / 2;
+      (inside(a + c * m, b + d * m) ? ng : ok) = m;
+    }
+    return ok;
+  };
+
+  vector<pair<Int, Int>> ps;
+  Int x = xl, y = yl;
+  assert(inside(x, y) and go(x, y, 0, -1) == 0);
+  ps.emplace_back(x, y);
+  while (x < xr) {
+    Int a, b;
+    if (f(x + 1, y)) {
+      a = 1, b = 0;
+    } else {
+      SternBrocotTreeNode<Int> sb;
+      while (true) {
+        assert(f(x + sb.lx, y + sb.ly));
+        assert(!f(x + sb.rx, y + sb.ry));
+        if (f(x + sb.lx + sb.rx, y + sb.ly + sb.ry)) {
+          Int s = go(x + sb.lx, y + sb.ly, sb.rx, sb.ry);
+          assert(s > 0);
+          sb.go_right(s);
+        } else {
+          Int s = candicate(x + sb.rx, y + sb.ry, sb.lx, sb.ly);
+          if (s <= 0 || !inside(x + sb.lx * s + sb.rx, y + sb.ly * s + sb.ry)) {
+            a = sb.lx, b = sb.ly;
+            break;
+          } else {
+            Int t = go2(x + sb.rx, y + sb.ry, sb.lx, sb.ly, s);
+            sb.go_left(t);
+          }
+        }
+      }
+    }
+    Int s = go(x, y, a, b);
+    x += a * s, y += b * s;
+    ps.emplace_back(x, y);
+  }
+  return ps;
+}

math/sum-of-floor.hpp → math/floor-sum.hpp

@@ -1,7 +1,5 @@
 #pragma once
 
-
-
 // sum_{0 <= i < N} (ai + b) // m
 template <typename T>
 T sum_of_floor(T n, T m, T a, T b) {

math/gaussian-integer.hpp

@@ -0,0 +1,57 @@
+#pragma once
+
+// x + yi
+template <typename T>
+struct Gaussian_Integer {
+  T x, y;
+  using G = Gaussian_Integer;
+
+  Gaussian_Integer(T _x = 0, T _y = 0) : x(_x), y(_y) {}
+  Gaussian_Integer(const pair<T, T>& p) : x(p.fi), y(p.se) {}
+
+  T norm() const { return x * x + y * y; }
+  G conj() const { return G{x, -y}; }
+
+  G operator+(const G& r) const { return G{x + r.x, y + r.y}; }
+  G operator-(const G& r) const { return G{x - r.x, y - r.y}; }
+  G operator*(const G& r) const {
+    return G{x * r.x - y * r.y, x * r.y + y * r.x};
+  }
+  G operator/(const G& r) const {
+    G g = G{*this} * r.conj();
+    T n = r.norm();
+    g.x += n / 2, g.y += n / 2;
+    return G{g.x / n - (g.x % n < 0), g.y / n - (g.y % n < 0)};
+  }
+  G operator%(const G& r) const { return G{*this} - G{*this} / r * r; }
+
+  G& operator+=(const G& r) { return *this = G{*this} + r; }
+  G& operator-=(const G& r) { return *this = G{*this} - r; }
+  G& operator*=(const G& r) { return *this = G{*this} * r; }
+  G& operator/=(const G& r) { return *this = G{*this} / r; }
+  G& operator%=(const G& r) { return *this = G{*this} % r; }
+  G operator-() const { return G{-x, -y}; }
+  G operator+() const { return G{*this}; }
+  bool operator==(const G& g) const { return x == g.x && y == g.y; }
+  bool operator!=(const G& g) const { return x != g.x || y != g.y; }
+
+  G pow(__int128_t e) const {
+    G res{1}, a{*this};
+    while (e) {
+      if (e & 1) res *= a;
+      a *= a, e >>= 1;
+    }
+    return res;
+  }
+
+  friend G gcd(G a, G b) {
+    while (b != G{0, 0}) {
+      trc(a, b, a / b, a % b);
+      swap(a %= b, b);
+    }
+    return a;
+  }
+  friend ostream& operator<<(ostream& os, const G& rhs) {
+    return os << rhs.x << " " << rhs.y;
+  }
+};

math/two-square.hpp

@@ -0,0 +1,91 @@
+#pragma once
+
+#include "../internal/internal-math.hpp"
+#include "../prime/fast-factorize.hpp"
+#include "gaussian-integer.hpp"
+
+// 解が存在しない場合 (-1, -1) を返す
+Gaussian_Integer<__int128_t> solve_norm_equation_prime(long long p) {
+  if (p % 4 == 3) return {-1, -1};
+  if (p == 2) return {1, 1};
+  long long x = 1;
+  while (true) {
+    x++;
+    long long z = internal::modpow<long long, __int128_t>(x, (p - 1) / 4, p);
+    if (__int128_t(z) * z % p == p - 1) {
+      x = z;
+      break;
+    }
+  }
+  long long y = 1, k = (__int128_t(x) * x + 1) / p;
+  while (k > 1) {
+    long long B = x % k, D = y % k;
+    if (B < 0) B += k;
+    if (D < 0) D += k;
+    if (B * 2 > k) B -= k;
+    if (D * 2 > k) D -= k;
+    long long nx = (__int128_t(x) * B + __int128_t(y) * D) / k;
+    long long ny = (__int128_t(x) * D - __int128_t(y) * B) / k;
+    x = nx, y = ny;
+    k = (__int128_t(x) * x + __int128_t(y) * y) / p;
+  }
+  return {x, y};
+}
+
+// p^e が long long に収まる
+vector<Gaussian_Integer<__int128_t>> solve_norm_equation_prime_power(
+    long long p, long long e) {
+  using G = Gaussian_Integer<__int128_t>;
+  if (p % 4 == 3) {
+    if (e % 2 == 1) return {};
+    long long x = 1;
+    for (int i = 0; i < e / 2; i++) x *= p;
+    return {G{x}};
+  }
+  if (p == 2) return {G{1, 1}.pow(e)};
+  G pi = solve_norm_equation_prime(p);
+  vector<G> pows(e + 1);
+  pows[0] = 1;
+  for (int i = 1; i <= e; i++) pows[i] = pows[i - 1] * pi;
+  vector<G> res(e + 1);
+  for (int i = 0; i <= e; i++) res[i] = pows[i] * (pows[e - i].conj());
+  return res;
+}
+
+// 0 <= arg < 90 の範囲の解のみ返す
+vector<Gaussian_Integer<__int128_t>> solve_norm_equation(long long N) {
+  using G = Gaussian_Integer<__int128_t>;
+  if (N < 0) return {};
+  if (N == 0) return {G{0}};
+  auto pes = factor_count(N);
+  for (auto& [p, e] : pes) {
+    if (p % 4 == 3 && e % 2 == 1) return {};
+  }
+  vector<G> res{G{1}};
+  for (auto& [p, e] : pes) {
+    vector<G> cur = solve_norm_equation_prime_power(p, e);
+    vector<G> nxt;
+    for (auto& g1 : res) {
+      for (auto& g2 : cur) nxt.push_back(g1 * g2);
+    }
+    res = nxt;
+  }
+
+  for (auto& g : res) {
+    while (g.x <= 0 || g.y < 0) g = G{-g.y, g.x};
+  }
+  return res;
+}
+
+// x,y 両方非負のみ, 辞書順で返す
+vector<pair<long long, long long>> two_square(long long N) {
+  if (N < 0) return {};
+  if (N == 0) return {{0, 0}};
+  vector<pair<long long, long long>> ans;
+  for (auto& g : solve_norm_equation(N)) {
+    ans.emplace_back(g.x, g.y);
+    if (g.y == 0) ans.emplace_back(g.y, g.x);
+  }
+  sort(begin(ans), end(ans));
+  return ans;
+}

misc/doubling.hpp

@@ -9,40 +9,40 @@ struct BinaryLifting {
   T I;
 
   BinaryLifting(int n, uint64_t lim, const T I_ = T())
-      : N(n), LOG(__lg(lim) + 2), I(I_) {
-    table.resize(n, vector<Data>(LOG, Data(-1, I)));
+      : N(n), LOG(__lg(max<uint64_t>(lim, 1)) + 2), I(I_) {
+    table.resize(LOG, vector<Data>(n, Data(-1, I)));
   }
 
-  void set_next(int k, int nxt, const T& t) { table[k][0] = Data(nxt, t); }
+  void set_next(int k, int nxt, const T& t) { table[0][k] = Data(nxt, t); }
 
   void build() {
     for (int k = 0; k + 1 < LOG; ++k)
       for (int i = 0; i < N; ++i) {
-        int pre = table[i][k].first;
+        int pre = table[k][i].first;
         if (pre == -1) {
-          table[i][k + 1] = table[i][k];
+          table[k + 1][i] = table[k][i];
         } else {
-          table[i][k + 1].first = table[pre][k].first;
-          table[i][k + 1].second = table[i][k].second + table[pre][k].second;
+          table[k + 1][i].first = table[k][pre].first;
+          table[k + 1][i].second = table[k][i].second + table[k][pre].second;
         }
       }
   }
 
-  // from i, move t times
+  // i から t 回進んだ先, (地点, モノイド和) を返す
   Data query(int i, uint64_t t) const {
     T d = I;
     for (int k = LOG - 1; k >= 0; k--) {
       if ((t >> k) & 1) {
-        d = d + table[i][k].second;
-        i = table[i][k].first;
+        d = d + table[k][i].second;
+        i = table[k][i].first;
       }
       if (i == -1) break;
     }
     return Data(i, d);
   }
 
   // query(i, pow(2, k))
-  inline Data query_pow(int i, int k) const { return table[i][k]; }
+  inline Data query_pow(int i, int k) const { return table[k][i]; }
 
   // assuming graph is DAG ( edge(u, v) <-> u < v )
   // find max j | j <= t, path from i to j exists
@@ -51,9 +51,9 @@ struct BinaryLifting {
     T d = I;
     uint64_t times = 0;
     for (int k = LOG - 1; k >= 0; k--) {
-      int nxt = table[thres][k].first;
+      int nxt = table[k][thres].first;
       if (nxt != -1 && nxt <= t) {
-        d = d + table[thres][k].second;
+        d = d + table[k][thres].second;
         thres = nxt;
         times += 1LL << k;
       }
@@ -68,9 +68,9 @@ struct BinaryLifting {
     T d = I;
     uint64_t times = 0;
     for (int k = LOG - 1; k >= 0; k--) {
-      int nxt = table[thres][k].first;
+      int nxt = table[k][thres].first;
       if (nxt != -1 && nxt >= t) {
-        d = d + table[thres][k].second;
+        d = d + table[k][thres].second;
         thres = nxt;
         times += 1LL << k;
       }

misc/mo-fast.hpp

@@ -11,6 +11,7 @@ struct Fast_Mo {
   int N, Q, width;
   vector<int> L, R, order;
   bool is_build;
+  int nl, nr;
 
   Fast_Mo(int _n, int _q) : N(_n), Q(_q), order(Q), is_build(false) {
     width = max<int>(1, 1.0 * N / max<double>(1.0, sqrt(Q / 2.0)));
@@ -28,7 +29,7 @@ struct Fast_Mo {
   void run(const AL &add_left, const AR &add_right, const DL &delete_left,
            const DR &delete_right, const REM &rem) {
     if (!is_build) build();
-    int nl = 0, nr = 0;
+    nl = nr = 0;
     for (auto idx : order) {
       while (nl > L[idx]) add_left(--nl);
       while (nr < R[idx]) add_right(nr++);
@@ -61,7 +62,7 @@ struct Fast_Mo {
   }
 
   int dist(int i, int j) { return abs(L[i] - L[j]) + abs(R[i] - R[j]); }
-  
+
   void climb(int iter = 3, int interval = 5) {
     vector<int> d(Q - 1);
     for (int i = 0; i < Q - 1; i++) d[i] = dist(order[i], order[i + 1]);