Problem: antisocial (2021 finals)

Cole is a programmer in the coal division of Metal (previously Facebulk), working across \(T\) office buildings with very pleasant coworkers and unlimited free ice cream. However, as a self-proclaimed antisocial code monkey, he still prefers to "mine" his own business when possible.

A given office building can be represented by an axis-aligned rectangle with opposite corners \((0, 0)\) and \((X_R, Y_R)\) on a 2D plane. Inside, there are \(N\) friendly coworkers, the \(i\)th of whom is permanently seated at position \((X_i, Y_i)\), distinct from all other coworkers.

Cole is currently at position \((X_A, Y_A)\), and would like to reach the unlimited free ice cream at a different position \((X_B, Y_B)\). Neither position is located at a coworker's seat.

Please help Cole determine the closest distance he'll need to get to another coworker in order to travel from \((X_A, Y_A)\) to \((X_B, Y_B)\) without leaving the building. Your answer will be accepted if it is within either \(10^{-6}\) or \(0.000001\%\) of the right answer.

# Constraints

\(1 \le T \le 800\)
\(1 \le N \le 800{,}000\)
\(2 \le X_R, Y_R \le 1{,}000{,}000{,}000\)
\(0 < X_i, X_A, X_B < X_R\)
\(0 < Y_i, Y_A, Y_B < Y_R\)

The sum of \(N\) across all buildings is at most \(3{,}000{,}000\).

# Input

Input begins with a single integer \(T\), the number of office buildings. For each building, there is first a line containing \(2\) space-separated integers, \(X_R\) and \(Y_R\). Then, there is a line containing \(4\) space-separated integers, \(X_A\), \(Y_A\), \(X_B\), and \(Y_B\). Then, there a line containing a single integer \(N\). Then, \(N\) lines follow, the \(i\)th of which consists of \(2\) space-separated integers, \(X_i\) and \(Y_i\).

# Output

For the \(i\)th building, print a line containing *"Case #i: "* followed by a single real number, the closest that Cole will need to get to another coworker in order to reach the free ice cream.

# Sample Explanation

In the first building, walking directly from \((3, 1)\) to \((17, 1)\) would take Cole within just 2 units of the coworker at \((10, 3)\). Instead, one optimal route might involve walking upward first, then across, and then back down, as illustrated below:

{{PHOTO_ID:441536274154579|WIDTH:450}}

The closest such a route gets to the coworker is \(7\) units (for example at coordinates \((10, 10)\).

In the second building, Cole's route must begin and end at a distance of \(\sqrt{4^2 + 2^2} \approx 4.47\) units from the coworker at \((10, 3)\). It's possible for Cole to get around them without getting any closer than that, for example with the following route:

{{PHOTO_ID:4816243665107323|WIDTH:450}}
  
In the third building, one optimal route for Cole might look as follows:

{{PHOTO_ID:440615224330860|WIDTH:280}}


Sample input:
7
20 10
3 1 17 1
1
10 3
20 10
6 1 14 1
1
10 3
5 10
1 1 1 9
2
1 5
4 5
5 5
2 3 3 2
2
2 2
3 3
5 5
2 3 3 2
14
2 2
3 3
1 1
1 2
1 3
1 4
2 4
3 4
4 4
4 3
4 2
4 1
3 1
2 1
100 100
83 92 34 18
5
47 41
58 90
86 53
25 76
65 31
400 700
363 34 96 512
20
247 522
182 289
165 653
305 392
265 186
346 513
189 605
379 328
357 628
255 58
5 319
270 583
195 74
374 651
121 169
152 162
392 645
88 42
183 562
393 640


Sample output:
Case #1: 7.00000000
Case #2: 4.47213595
Case #3: 1.50000000
Case #4: 1.00000000
Case #5: 0.70710678
Case #6: 20.67002661
Case #7: 80.21533519


Solution: We can imagine "danger circles" around all coworkers increase at a constant size until there is no path to the finish. This idea is very similar to the construction of a Voronoi diagram, which can be done in \(O(N \log(N))\) time using Fortune's algorithm. If we remove the restriction that you are not allowed to leave the building's borders, the optimal path maximizing the distance to the nearest person will be a sequence of edges on the Voronoi diagram, as well as a segment from the start position to a Voronoi node and a segment from the end position to a Voronoi node. For the first and last segment, we simply want to walk to a Voronoi node that isn't *closer* to our nearest neighbor than we are, of the Voronoi nodes tangential to this neighbor's corresponding node in the Delaunay triangulation. Valid start and end Voronoi nodes can be found by iterating through all neighboring Voronoi nodes to the nearest person, and simply taking the one which is most in the direction away from the nearest coworker. Then, we can handle the path of Voronoi edges by sorting the edges and using a disjoint set (similarly to Kruskal's minimum spanning tree algorithm, though in this case as a maximum spanning tree), stopping once the start and end Voronoi nodes are unioned.

Note that there are also 2 cases in which the path might not remain on the Voronoi diagram's edges:

  1. At the start (and similarly at the end), the path must move from the the starting point \(a = (X_A, Y_A)\) to some graph node. If it's in the Voronoi cell of some closest point \(x\), then it should move to some node on the boundary of that cell without getting closer to \(x\) (any such node will do). A valid choice of node is the one which is closest angle-wise (relative to \(x\)) to \(a\).
  
  2. The path can also travel along the borders of the building's rectangle — there should effectively be a graph node at each rectangle corner and at each intersection of an (infinite) Voronoi edge with the rectangle's border, with these nodes connected around the border by edges (weighted like the other edges). A good way of constructing this situation is to begin by mirroring all \(N\) points across the rectangle's \(4\) sides, and then computing the Voronoi diagram for all \(5N\) resulting points. Provided that we then include Voronoi edges lying inside or on the border of the rectangle (excluding ones lying outside it), the set of edges will then be exactly what we need. 
  This is because the introduction of the mirrored points doesn't affect the portion of the Voronoi diagram strictly inside the rectangle (as no mirrored point can be the closest point to a point inside the rectangle; the corresponding original point must be closer), while at each point along the rectangle's boundary, whatever original point was closest to it, the corresponding mirrored point must now be equidistant (resulting in a Voronoi node/edge at the border). 
  For a visualization of this concept with a 3-point example, check out [this tool](http://alexbeutel.com/webgl/voronoi.html) with the following input data: `{"sites":[220,240,270,250,290,220,180,240,220,160,380,240,220,360,270,150,270,350,130,250,330,250,290,180,290,380,110,220,310,220],"queries":[]}`.
  
The time complexity of the above algorithm is \(O(N \log(N))\).

It's worth noting that it's also possible to solve this problem for many (e.g. \(10^5\)) queries of start/end points, though that extension was not included in the contest as its implementation would be rather long.

Computation of which Voronoi cell and graph node each starting/ending point corresponds to must now be done more efficiently, in less than \(O(N)\) time per such point. It's possible to augment Fortune's algorithm to determine which Voronoi cell each such point lies within. Alternatively, you can do a line sweep from left to right, finding for each query point which Voronoi edge is above it (form which you can determine which Voronoi cell it is in). For each cell, if the angles of its surrounding Voronoi vertices (relative to the cell's central point) are precomputed and sorted, it's possible to binary search over them for each starting/ending point to find its enclosing triangle. These changes take a total of \(O((N+K) \log(N))\) time.

To solve the second part of the problem, one approach involves still generating a maximum spanning tree using Kruskal's algorithm. The answer for a query is then based on the minimum edge weight on the path connecting its starting/ending points' corresponding nodes in this tree. This can be computed for each query using binary lifting (or other methods). Thus, the time complexity of this portion also comes out to \(O((N+K) \log(N))\) time.

Another approach for the second part of the problem is to use a parallel binary search. For each query, we can store its minimum and maximum possible answer, and iterate through all edges until these ranges are sufficiently small. At each iteration, the range for each query's possible answer will be halved. An iteration consists of maintaining a disjoint set of all components, processing queries and edges in decreasing order of weight. When we get to a query, we simply check if its two nodes are in the same component, and adjust either the high or low value of that queries binary search accordingly (this will change its mid for where it is sorted in the next iteration). At a Voronoi edge, we simply union the two components.



Code: #include <algorithm>
#include <cassert>
#include <cmath>
#include <functional>
#include <iostream>
#include <queue>
#include <utility>
#include <vector>
using namespace std;

const int LIM = 6000006;
const long double EPS = 1e-9;

using LD = long double;
using pii = pair<int, int>;

// From https://github.com/zigui-ps/VoronoiDiagram
struct point {
  LD x, y;

  point(): x(0), y(0) {}
  point(LD _x, LD _y): x(_x), y(_y) {}

  LD operator/(point p) const { return x * p.y - y * p.x; }
  point operator+(point p) const { return point(x + p.x, y + p.y); }
  point operator-(point p) const { return point(x - p.x, y - p.y); }
  friend point operator*(LD b, point a) { return point(b * a.x, b * a.y); }

  LD size() const { return hypot(x, y); }
  LD sz2() const { return x * x + y * y; }
  point r90() const { return point(-y, x); }
  LD dot(point p) { return x * p.x + y * p.y; }
};

int dcmp(LD x) {
  return x < -EPS ? -1 : x > EPS ? 1 : 0;
}

point line_intersect(point a, point b, point u, point v) {
  return u + (((a - u) / b) / (v / b)) * v;
}

point get_circumcenter(point p0, point p1, point p2) {
  return line_intersect(
    0.5 * (p0 + p1), (p0 - p1).r90(), 0.5 * (p1 + p2), (p1 - p2).r90()
  );
}

point project_point_segment(point a, point b, point c) {
  LD r = (b - a).dot(b - a);
  if (fabs(r) < EPS) {
    return a;
  }
  r = (c - a).dot(b - a) / r;
  if (r < 0) {
    return a;
  }
  if (r > 1) {
    return b;
  }
  return a + r * (b - a);
}

bool point_in_polygon_incl(const vector<point> &p, point q) {
  bool c = 0;
  for (int i = 0; i < p.size(); i++) {
    int j = (i + 1) % p.size();
    if ((project_point_segment(p[i], p[j], q) - q).sz2() < EPS) {
      return true;
    }
    if (
      (p[i].y <= q.y && q.y < p[j].y || p[j].y <= q.y && q.y < p[i].y) &&
      q.x < p[i].x + (p[j].x - p[i].x) * (q.y - p[i].y) / (p[j].y - p[i].y)
    )  {
      c = !c;
    }
  }
  return c;
}

// https://www.youtube.com/watch?v=h_vvP4ah6Ck
LD parabola_intersect(point left, point right, LD sweepline) {
  auto f2 = [](point left, point right, LD sweepline) {
    int sign = left.x < right.x ? 1 : -1;
    point m = 0.5 * (left + right);
    point v = line_intersect(
      m, (right - left).r90(), point(0, sweepline), point(1, 0)
    );
    point w = line_intersect(m, (left - v).r90(), v, left - v);
    LD l1 = (v - w).size();
    LD l2 = sqrt(pow(sweepline - m.y, 2) - (m - w).sz2());
    LD l3 = (left - v).size();
    return v.x + (m.x - v.x) * l3 / (l1 + sign * l2);
  };
  if (fabs(left.y - right.y) < fabs(left.x - right.x) * EPS) {
    return f2(left, right, sweepline);
  }
  int sign = left.y < right.y ? -1 : 1;
  point v = line_intersect(
    left, right - left, point(0, sweepline), point(1, 0)
  );
  LD d1 = (0.5 * (left + right) - v).sz2(), d2 = (0.5 * (left - right)).sz2();
  return v.x + sign * sqrt(max(0.0L, d1 - d2));
}

struct Beachline {
  struct node {
    point pt;
    int idx;
    int end;
    node *link[2], *par, *prv, *nxt;

    node() {}
    node(point _pt, int idx):
      pt(_pt), idx(idx), end(0), link{0, 0}, par(0), prv(0), nxt(0) {}
  } *root;

  LD sweepline;

  Beachline() : sweepline(-1e20), root(NULL) { }

  inline int dir(node* x) {
    return x->par->link[0] != x;
  }

  //     p        n          p            n
  //    / \      / \        / \          / \
  //   n   d => a   p   or a   n   =>   p   d
  //  / \          / \        / \      / \
  // a   b        b   d      c   d    a   c
  void rotate(node* n) {
    node* p = n->par;
    int d = dir(n);
    p->link[d] = n->link[!d];
    if (n->link[!d]) {
      n->link[!d]->par = p;
    }
    n->par = p->par;
    if (p->par) {
      p->par->link[dir(p)] = n;
    }
    n->link[!d] = p;
    p->par = n;
  }

  void splay(node* x, node* f = NULL) {
    while (x->par != f) {
      if (x->par->par == f) {
        // no-op
      } else if (dir(x) == dir(x->par)) {
        rotate(x->par);
      } else {
        rotate(x);
      }
      rotate(x);
    }
    if (f == NULL) {
      root = x;
    }
  }

  void insert(node* n, node* p, int d) {
    splay(p);
    node *c = p->link[d];
    n->link[d] = c;
    if (c) {
      c->par = n;
    }
    p->link[d] = n;
    n->par = p;
    node *prv = !d ? p->prv : p, *nxt = !d ? p : p->nxt;
    n->prv = prv;
    if (prv) {
      prv->nxt = n;
    }
    n->nxt = nxt;
    if (nxt) {
      nxt->prv = n;
    }
  }

  void erase(node* n) {
    node *prv = n->prv, *nxt = n->nxt;
    if (!prv && !nxt) {
      if (n == root) {
        root = NULL;
      }
      return;
    }
    n->prv = NULL;
    if (prv) {
      prv->nxt = nxt;
    }
    n->nxt = NULL;
    if (nxt) {
      nxt->prv = prv;
    }
    splay(n);
    if (!nxt) {
      root->par = NULL;
      n->link[0] = NULL;
      root = prv;
    } else {
      splay(nxt, n);
      node* c = n->link[0];
      nxt->link[0] = c;
      c->par = nxt;
      n->link[0] = NULL;
      n->link[1] = NULL;
      nxt->par = NULL;
      root = nxt;
    }
  }

  bool get_event(node *cur, LD &next_sweep) {
    if (!cur->prv || !cur->nxt) {
      return false;
    }
    point u = (cur->pt - cur->prv->pt).r90();
    point v = (cur->nxt->pt - cur->pt).r90();
    if (dcmp(u / v) != 1) {
      return false;
    }
    point p = get_circumcenter(cur->pt, cur->prv->pt, cur->nxt->pt);
    next_sweep = p.y + (p - cur->pt).size();
    return true;
  }

  node* find_beachline(LD x) {
    node* cur = root;
    while (cur) {
      LD left = cur->prv
                  ? parabola_intersect(cur->prv->pt, cur->pt, sweepline)
                  : -1e30;
      LD right = cur->nxt
                  ? parabola_intersect(cur->pt, cur->nxt->pt, sweepline)
                  : 1e30;
      if (left <= x && x <= right) {
        splay(cur);
        return cur;
      }
      cur = cur->link[x > right];
    }
    return NULL;
  }
};

using BeachNode = Beachline::node;

static BeachNode arr[LIM];
static int sz;
static BeachNode* new_node(point point, int idx) {
  arr[sz] = BeachNode(point, idx);
  return arr + (sz++);
}

struct event {
  int type, idx, prv, nxt;
  BeachNode *cur;
  LD sweep;

  event(LD sweep, int idx) : type(0), sweep(sweep), idx(idx) {}
  event(LD sweep, BeachNode* cur) :
    type(1), sweep(sweep), prv(cur->prv->idx), cur(cur), nxt(cur->nxt->idx) {}

  bool operator>(const event& l) const { return sweep > l.sweep; }
};

void voronoi_diagram(
  vector<point> &input,
  vector<point> &vertex,
  vector<pii> &edge,
  vector<pii> &area
) {
  Beachline beachline = Beachline();
  priority_queue<event, vector<event>, greater<event>> events;

  auto add_edge = [&](int u, int v, int a, int b, BeachNode* c1, BeachNode* c2) {
    if (c1) {
      c1->end = edge.size() * 2;
    }
    if (c2) {
      c2->end = edge.size() * 2 + 1;
    }
    edge.emplace_back(u, v);
    area.emplace_back(a, b);
  };
  auto write_edge = [&](int idx, int v) {
    if (idx % 2 == 0) {
      edge[idx / 2].first = v;
    } else {
      edge[idx / 2].second = v;
    }
  };
  auto add_event = [&](BeachNode* cur) {
    LD nxt;
    if (beachline.get_event(cur, nxt)) {
      events.emplace(nxt, cur);
    }
  };

  int n = input.size(), cnt = 0;
  /*arr = new BeachNode[n * 4]; */sz = 0;
  sort(input.begin(), input.end(), [](const point &l, const point &r) {
    return l.y != r.y ? l.y < r.y : l.x < r.x;
  });

  BeachNode *tmp = beachline.root = new_node(input[0], 0), *t2;
  for (int i = 1; i < n; i++) {
    if (dcmp(input[i].y - input[0].y) == 0) {
      add_edge(-1, -1, i - 1, i, 0, tmp);
      beachline.insert(t2 = new_node(input[i], i), tmp, 1);
      tmp = t2;
    } else {
      events.emplace(input[i].y, i);
    }
  }

  while (events.size()) {
    event q = events.top();
    events.pop();
    BeachNode *prv, *cur, *nxt, *site;
    int v = vertex.size(), idx = q.idx;
    beachline.sweepline = q.sweep;
    if (q.type == 0) {
      point pt = input[idx];
      cur = beachline.find_beachline(pt.x);
      beachline.insert(site = new_node(pt, idx), cur, 0);
      beachline.insert(prv = new_node(cur->pt, cur->idx), site, 0);
      add_edge(-1, -1, cur->idx, idx, site, prv);
      add_event(prv);
      add_event(cur);
    } else {
      cur = q.cur;
      prv = cur->prv;
      nxt = cur->nxt;
      if (!prv || !nxt || prv->idx != q.prv || nxt->idx != q.nxt) {
        continue;
      }
      vertex.push_back(get_circumcenter(prv->pt, nxt->pt, cur->pt));
      write_edge(prv->end, v);
      write_edge(cur->end, v);
      add_edge(v, -1, prv->idx, nxt->idx, 0, prv);
      beachline.erase(cur);
      add_event(prv);
      add_event(nxt);
    }
  }
  // delete arr;
}

int XR, YR;
int N, M;
vector<pair<int, LD>> adj[LIM];
point KP[2];
int KN[2];
LD ans;

void process_voronoi_diagrams(vector<point> P) {
  int i, j, k;
  vector<point> vertex;
  vector<pii> edge, area;
  voronoi_diagram(P, vertex, edge, area);
  M = vertex.size();
  for (int i = 0; i < area.size(); i++) {
    int e1 = edge[i].first, e2 = edge[i].second;
    if (e1 < 0 || e2 < 0) {
      continue; // Infinite edge.
    }
    point v1 = vertex[e1], v2 = vertex[e2];
    if (
      min(v1.x, v2.x) < -EPS || max(v1.x, v2.x) > XR + EPS ||
      min(v1.y, v2.y) < -EPS || max(v1.y, v2.y) > YR + EPS
    ) {
      continue; // Edge outside rectangle.
    }
    // Add edge to graph with appropriate weight.
    int pi1 = area[i].first, pi2 = area[i].second;
    point p1 = P[pi1], p2 = P[pi2];
    LD d = (p1 - project_point_segment(v1, v2, p1)).size();
    adj[e1].emplace_back(e2, d);
    adj[e2].emplace_back(e1, d);
    // Consider any key points in triangles adjacent to this edge.
    point mid = 0.5 * (p1 + p2);
    for (int j : {0, 1}) {
      for (int k : {0, 1}) {
        point kp = KP[j], p = k ? p2 : p1;
        ans = min(ans, (kp - p).size());  // Must be this close to a point.
        if (point_in_polygon_incl({p, v1, mid}, kp)) {
          KN[j] = e1;
        }
        if (point_in_polygon_incl({p, mid, v2}, kp)) {
          KN[j] = e2;
        }
      }
    }
  }
}

LD solve() {
  vector<point> P;
  for (int i = 0; i < LIM; i++) {
    adj[i].clear();
  }
  // Input.
  cin >> XR >> YR;
  cin >> KP[0].x >> KP[0].y >> KP[1].x >> KP[1].y;
  cin >> N;
  for (int i = 0; i < N; i++) {
    int a, b;
    cin >> a >> b;
    P.emplace_back(a, b);
    // Additionally mirror each point across each side of the rectangle.
    P.emplace_back(-a, b);
    P.emplace_back(a, -b);
    P.emplace_back(2 * XR - a, b);
    P.emplace_back(a, 2 * YR - b);
  }
  // Compute Voronoi diagram and corresponding graph.
  ans = 1e18;
  KN[0] = KN[1] = -1;
  process_voronoi_diagrams(P);
  assert(KN[0] >= 0 && KN[1] >= 0);
  // MST.
  vector<bool> visit(M);
  priority_queue<pair<LD, int>> Q;
  Q.emplace(ans, KN[0]);
  while (!Q.empty()) {
    ans = min(ans, Q.top().first);
    int i = Q.top().second;
    Q.pop();
    if (i == KN[1]) {
      break;
    }
    if (visit[i]) {
      continue;
    }
    visit[i] = true;
    for (auto c : adj[i]) {
      if (!visit[c.first]) {
        Q.emplace(c.second, c.first);
      }
    }
  }
  return ans;
}

int main() {
  cout.precision(10);
  int T;
  cin >> T;
  for (int t = 1; t <= T; t++) {
    cout << "Case #" << t << ": " << fixed << solve() << endl;
  }
  return 0;
}