geo/range_closest_pair_query.hpp

• View this file on GitHub
• Last update: 2023-11-22 02:53:14+09:00
• Include: #include "geo/range_closest_pair_query.hpp"
• Link: https://codeforces.com/gym/104172/attachments/download/18933/Hong_Kong_Tutorial.pdf
• Link: https://qoj.ac/problem/5463

Depends on

• alg/monoid/min.hpp
• ds/hashmap.hpp
• ds/segtree/dual_segtree.hpp

Verified with

• test/mytest/range_closest_pair.test.cpp

Code

#include "ds/hashmap.hpp"
#include "ds/segtree/dual_segtree.hpp"
#include "alg/monoid/min.hpp"

// 点群 {p_i | i in [l, r)} に対する最近点対の計算を行うクエリ
// O(KNlogKN + QlogN)
// https://qoj.ac/problem/5463
// https://codeforces.com/gym/104172/attachments/download/18933/Hong_Kong_Tutorial.pdf
struct Range_Closest_Pair_Query {
  /*
  ・R を増やしながら、L ごとの答を管理する
  ・2^{k-1} <= ANS[L:R] < 2^{k} となる L ：レベル k （レベル 0：距離 0）
  ・レベル 0, 1, 2, ..., 29 のグリッドを用意する
  　・幅は 2^k
　  ・一辺 1.99 の正方形内で点対距離が 1 以上 → 8 個までありうる
  ・レベル 29, 28, ..., 0 の順に探索する：9 近傍
  　・答が見つかったらレベルを下げる。左向きに伝搬。
  　・レベルの減少は 30N 回までしか起きない
  */
  const int LOG = 30;
  vc<pair<int, int>> point;
  vc<pair<int, int>> query;
  void add_point(int x, int y) {
    assert(0 <= x && x < (1 << LOG));
    assert(0 <= y && y < (1 << LOG));
    point.eb(x, y);
  }
  void add_query(int L, int R) {
    assert(R - L >= 2);
    query.eb(L, R);
  }
  ll dist(int i, int j) {
    ll dx = point[i].fi - point[j].fi;
    ll dy = point[i].se - point[j].se;
    return dx * dx + dy * dy;
  }

  vc<ll> calc() {
    static HashMap<int, 20, true> MP;
    MP.reset();
    const int K = LOG;
    const int N = len(point), Q = len(query);
    using A9 = array<int, 9>;
    // それぞれのレベルのときのセル番号
    vv(int, IDX, K, N, -1);
    // 各セル番号に対する近傍
    vc<A9> nbd;
    FOR(k, 1, K) {
      MP.reset();
      auto to_64 = [&](int x, int y) -> u64 { return u64(x) << 30 | y; };
      int off = len(nbd);
      int p = off;
      FOR(i, N) {
        int x = point[i].fi >> (k);
        int y = point[i].se >> (k);
        u64 key = to_64(x, y);
        int idx = MP.index(key);
        if (!MP.used[idx]) {
          MP.used[idx] = 1, MP.IDS.eb(idx), MP.key[idx] = key,
          MP.val[idx] = p++;
        }
        IDX[k][i] = MP.val[idx];
      }
      nbd.resize(p);
      FOR(i, N) {
        int x = point[i].fi >> (k);
        int y = point[i].se >> (k);
        int me = MP[to_64(x, y)];
        int s = 0;
        FOR(dx, -1, 2) FOR(dy, -1, 2) {
          u64 key = to_64(x + dx, y + dy);
          nbd[me][s++] = (MP.count(key) ? MP[key] : -1);
        }
      }
    }

    vc<array<int, 8>> dat(len(nbd), {-1, -1, -1, -1, -1, -1, -1, -1});
    auto add = [&](int k, int i) -> void {
      int idx = IDX[k][i];
      for (auto&& j: dat[idx]) {
        if (j == -1) {
          j = i;
          return;
        }
      }
    };
    auto rm = [&](int k, int i) -> void {
      int idx = IDX[k][i];
      for (auto&& j: dat[idx]) {
        if (j == i) {
          j = -1;
          return;
        }
      }
    };

    auto solve_level = [&](int k, int i) -> vc<pair<int, ll>> {
      // レベル k の点群に対する答の計算
      vc<pair<int, ll>> res;
      int me = IDX[k][i];
      for (auto&& idx: nbd[me]) {
        if (idx == -1) continue;
        for (auto&& j: dat[idx]) {
          if (j == -1) continue;
          res.eb(j, dist(i, j));
        }
      }
      return res;
    };
    Dual_SegTree<Monoid_Min<ll>> seg(N);
    vc<int> LEVEL(N, -1);
    auto get_lv = [&](ll d) -> int {
      if (d == 0) return 0;
      return topbit(d) / 2 + 1;
    };

    vc<int> left(Q);
    vvc<int> query_at(N);
    FOR(qid, Q) {
      auto [L, R] = query[qid];
      left[qid] = L;
      query_at[--R].eb(qid);
    }

    vi ANS(Q);

    FOR(R, N) {
      // R 番目の点を用いた答の更新
      vc<pair<int, ll>> upd;
      FOR(k, 1, K) {
        auto res = solve_level(k, R);
        upd.insert(upd.end(), all(res));
      }

      for (auto [i, d]: upd) {
        int lv = get_lv(d);
        if (seg.get(i) < d) continue;
        // 答えの更新
        seg.apply(0, i + 1, d);
        // レベルの更新
        while (i >= 0 && LEVEL[i] > lv) {
          rm(LEVEL[i], i);
          LEVEL[i] = lv;
          if (lv) add(lv, i);
          --i;
        }
      }
      LEVEL[R] = K - 1;
      add(K - 1, R);
      for (auto&& qid: query_at[R]) { ANS[qid] = seg.get(left[qid]); }
    }
    return ANS;
  }
};

#line 2 "ds/hashmap.hpp"

// u64 -> Val

template <typename Val, int LOG = 20, bool KEEP_IDS = false>
struct HashMap {
  static constexpr int N = (1 << LOG);
  u64* key;
  Val* val;
  vc<int> IDS;
  bitset<N> used;
  const int shift;
  const u64 r = 11995408973635179863ULL;
  HashMap() : key(new u64[N]), val(new Val[N]), shift(64 - LOG) {}
  u32 hash(u64 x) {
    static const u64 FIXED_RANDOM
        = std::chrono::steady_clock::now().time_since_epoch().count();
    return (u64(x + FIXED_RANDOM) * r) >> shift;
  }

  int index(const u64& k) {
    int i = 0;
    for (i = hash(k); used[i] && key[i] != k; (i += 1) &= (N - 1)) {}
    return i;
  }

  Val& operator[](const u64& k) {
    int i = index(k);
    if (!used[i]) {
      used[i] = 1, key[i] = k, val[i] = Val{};
      if constexpr (KEEP_IDS) IDS.eb(i);
    }
    return val[i];
  }

  Val get(const u64& k, Val default_value) {
    int i = index(k);
    if (!used[i]) return default_value;
    return val[i];
  }

  bool count(const u64& k) {
    int i = index(k);
    return used[i] && key[i] == k;
  }

  void reset() {
    static_assert(KEEP_IDS);
    for (auto&& i: IDS) used[i] = 0;
    IDS.clear();
  }

  // f(key, val)

  template <typename F>
  void enumerate_all(F f) {
    static_assert(KEEP_IDS);
    for (auto&& i: IDS) f(key[i], val[i]);
  }
};
#line 2 "ds/segtree/dual_segtree.hpp"

template <typename Monoid>
struct Dual_SegTree {
  using MA = Monoid;
  using A = typename MA::value_type;
  int n, log, size;
  vc<A> laz;

  Dual_SegTree() : Dual_SegTree(0) {}
  Dual_SegTree(int n) { build(n); }

  void build(int m) {
    n = m;
    log = 1;
    while ((1 << log) < n) ++log;
    size = 1 << log;
    laz.assign(size << 1, MA::unit());
  }

  A get(int p) {
    assert(0 <= p && p < n);
    p += size;
    for (int i = log; i >= 1; i--) push(p >> i);
    return laz[p];
  }

  vc<A> get_all() {
    FOR(i, size) push(i);
    return {laz.begin() + size, laz.begin() + size + n};
  }

  void apply(int l, int r, const A& a) {
    assert(0 <= l && l <= r && r <= n);
    if (l == r) return;
    l += size, r += size;
    if (!MA::commute) {
      for (int i = log; i >= 1; i--) {
        if (((l >> i) << i) != l) push(l >> i);
        if (((r >> i) << i) != r) push((r - 1) >> i);
      }
    }
    while (l < r) {
      if (l & 1) all_apply(l++, a);
      if (r & 1) all_apply(--r, a);
      l >>= 1, r >>= 1;
    }
  }

private:
  void push(int k) {
    if (laz[k] == MA::unit()) return;
    all_apply(2 * k, laz[k]), all_apply(2 * k + 1, laz[k]);
    laz[k] = MA::unit();
  }
  void all_apply(int k, A a) { laz[k] = MA::op(laz[k], a); }
};
#line 2 "alg/monoid/min.hpp"

template <typename E>
struct Monoid_Min {
  using X = E;
  using value_type = X;
  static constexpr X op(const X &x, const X &y) noexcept { return min(x, y); }
  static constexpr X unit() { return infty<E>; }
  static constexpr bool commute = true;
};
#line 4 "geo/range_closest_pair_query.hpp"

// 点群 {p_i | i in [l, r)} に対する最近点対の計算を行うクエリ
// O(KNlogKN + QlogN)
// https://qoj.ac/problem/5463
// https://codeforces.com/gym/104172/attachments/download/18933/Hong_Kong_Tutorial.pdf
struct Range_Closest_Pair_Query {
  /*
  ・R を増やしながら、L ごとの答を管理する
  ・2^{k-1} <= ANS[L:R] < 2^{k} となる L ：レベル k （レベル 0：距離 0）
  ・レベル 0, 1, 2, ..., 29 のグリッドを用意する
  　・幅は 2^k
　  ・一辺 1.99 の正方形内で点対距離が 1 以上 → 8 個までありうる
  ・レベル 29, 28, ..., 0 の順に探索する：9 近傍
  　・答が見つかったらレベルを下げる。左向きに伝搬。
  　・レベルの減少は 30N 回までしか起きない
  */
  const int LOG = 30;
  vc<pair<int, int>> point;
  vc<pair<int, int>> query;
  void add_point(int x, int y) {
    assert(0 <= x && x < (1 << LOG));
    assert(0 <= y && y < (1 << LOG));
    point.eb(x, y);
  }
  void add_query(int L, int R) {
    assert(R - L >= 2);
    query.eb(L, R);
  }
  ll dist(int i, int j) {
    ll dx = point[i].fi - point[j].fi;
    ll dy = point[i].se - point[j].se;
    return dx * dx + dy * dy;
  }

  vc<ll> calc() {
    static HashMap<int, 20, true> MP;
    MP.reset();
    const int K = LOG;
    const int N = len(point), Q = len(query);
    using A9 = array<int, 9>;
    // それぞれのレベルのときのセル番号
    vv(int, IDX, K, N, -1);
    // 各セル番号に対する近傍
    vc<A9> nbd;
    FOR(k, 1, K) {
      MP.reset();
      auto to_64 = [&](int x, int y) -> u64 { return u64(x) << 30 | y; };
      int off = len(nbd);
      int p = off;
      FOR(i, N) {
        int x = point[i].fi >> (k);
        int y = point[i].se >> (k);
        u64 key = to_64(x, y);
        int idx = MP.index(key);
        if (!MP.used[idx]) {
          MP.used[idx] = 1, MP.IDS.eb(idx), MP.key[idx] = key,
          MP.val[idx] = p++;
        }
        IDX[k][i] = MP.val[idx];
      }
      nbd.resize(p);
      FOR(i, N) {
        int x = point[i].fi >> (k);
        int y = point[i].se >> (k);
        int me = MP[to_64(x, y)];
        int s = 0;
        FOR(dx, -1, 2) FOR(dy, -1, 2) {
          u64 key = to_64(x + dx, y + dy);
          nbd[me][s++] = (MP.count(key) ? MP[key] : -1);
        }
      }
    }

    vc<array<int, 8>> dat(len(nbd), {-1, -1, -1, -1, -1, -1, -1, -1});
    auto add = [&](int k, int i) -> void {
      int idx = IDX[k][i];
      for (auto&& j: dat[idx]) {
        if (j == -1) {
          j = i;
          return;
        }
      }
    };
    auto rm = [&](int k, int i) -> void {
      int idx = IDX[k][i];
      for (auto&& j: dat[idx]) {
        if (j == i) {
          j = -1;
          return;
        }
      }
    };

    auto solve_level = [&](int k, int i) -> vc<pair<int, ll>> {
      // レベル k の点群に対する答の計算
      vc<pair<int, ll>> res;
      int me = IDX[k][i];
      for (auto&& idx: nbd[me]) {
        if (idx == -1) continue;
        for (auto&& j: dat[idx]) {
          if (j == -1) continue;
          res.eb(j, dist(i, j));
        }
      }
      return res;
    };
    Dual_SegTree<Monoid_Min<ll>> seg(N);
    vc<int> LEVEL(N, -1);
    auto get_lv = [&](ll d) -> int {
      if (d == 0) return 0;
      return topbit(d) / 2 + 1;
    };

    vc<int> left(Q);
    vvc<int> query_at(N);
    FOR(qid, Q) {
      auto [L, R] = query[qid];
      left[qid] = L;
      query_at[--R].eb(qid);
    }

    vi ANS(Q);

    FOR(R, N) {
      // R 番目の点を用いた答の更新
      vc<pair<int, ll>> upd;
      FOR(k, 1, K) {
        auto res = solve_level(k, R);
        upd.insert(upd.end(), all(res));
      }

      for (auto [i, d]: upd) {
        int lv = get_lv(d);
        if (seg.get(i) < d) continue;
        // 答えの更新
        seg.apply(0, i + 1, d);
        // レベルの更新
        while (i >= 0 && LEVEL[i] > lv) {
          rm(LEVEL[i], i);
          LEVEL[i] = lv;
          if (lv) add(lv, i);
          --i;
        }
      }
      LEVEL[R] = K - 1;
      add(K - 1, R);
      for (auto&& qid: query_at[R]) { ANS[qid] = seg.get(left[qid]); }
    }
    return ANS;
  }
};

Back to top page