Problem: cup_counterbalancing (2022 finals)

Courtney has an avant-garde kitchen counter made out of \(N\) rectangular sheets of metal of equal height, standing perpendicular to the ground. When viewed directly from the top, the counter looks like an (unfilled) polygon on a plane with \(N\) edges, the \(i\)th of which extends from point \((X_i, Y_i)\) to point \((X_{(i \text{ mod } N)+1},\, Y_{(i \text{ mod } N)+1})\). These line segments only intersect at their endpoints.

Courtney's friends often make fun of the counter, saying that the thin metal "surface" is useless for holding anything up, but Courtney wants to prove them wrong.

Courtney considers placing a cup with a circular rim of radius \(R\) onto the counter (with the rim facing down), centered at a point chosen uniformly randomly from the square \([0, L] \times [0, L]\). What is the chance that the cup stays balanced on the counter without falling? Your output will be accepted if it is within \(10^{-6}\) of the jury's answer.

*Note: The face-down cup is modeled as an unfilled ring—not as a disk.*


# Constraints

\(1 \le T \le 7\)
\(1 \le N \le 7\)
\(1 \le R, L \le 15\)
\(0 \le X_i, Y_i \le 15\)


# Input Format

Input begins with a single integer \(T\), the number of test cases. For each test case, there is first a line containing three space-separated integers \(N\), \(R\), and \(L\). Then, \(N\) lines follow, the \(i\)th of which contains two space-separated integers \(X_i\) and \(Y_i\),


# Output Format

For the \(i\)th test case, print a line containing `"Case #i: "` followed by a single real number between \(0\) and \(1\), the chance that the cup will balance if placed with its rim centered at a uniformly random point in the square \([0, L] \times [0, L]\).


# Sample Explanation

The counters in the first and second sample cases are depicted below, with the blue regions representing points at which centering the cup will lead to it staying balanced.

{{PHOTO_ID:556593865835592|WIDTH:650}}




Sample input:
2
3 1 4
0 0
2 2
4 0
5 1 4
0 0
1 2
0 3
4 3
3 1


Sample output:
Case #1: 0.1169663730642699
Case #2: 0.1353445414060363


Solution: Let's call a position where the circle doesn't fall *stable*. First, let's see how we can check whether a given point is stable. Consider all intersection points between that cup and surrounding positions (we can find these naively in linear time). If all of these points lie on a half-circle, then the ring will fall. Otherwise, there will be three points that create a triangle which contains the center of mass of the ring.

**Solution 1:**

As it turns out, any stable point in general form (we can ignore points that are edge cases because those points are ~\(0\%\) of the area) will have at least \(3\) intersections. In fact, if a point is stable it's true that there are a particular \(3\) intersection points that make it stable. This idea is informally supported with what we'll call the "Jenga Principle": if a ring is balancing on \(4\) or more points, one of them will be "loose" in the same way a Jenga block resting on three blocks below it will always be resting on at least one "loose" block which could be removed.

We can case out all of the possible ways that two or three line segments could support a ring, brute force those segments, and then compute the overlap of these shapes. This is tricky to do mathematically for the area on which a ring could balance on two segments (the shapes include more than just circles and line segments), so some method of curve approximation will be useful here.

**Solution 2:**

Actually, if we were to just pick random points \(10^6\) times, we'll get a reasonable approximation of the answer. To get even closer, there are other tricks we can try: first, we can better distribute our points to decrease variance, such as overlaying them evenly in a grid. Once we're there, there are several options to improve our sampling, and we can follow one or more to get a solution which is accurate enough and runs in time:

- Spend less time looking at spaces that are surrounded by points that are either all stable or all unstable: Because the bounds are small, these points are very unlikely to differ from their neighbors once the neighbors are close enough.
- Use a *k*-d tree or similar data structure to break up the search space, and stop searching if large areas are \(100\%\) similar after many samples.
- Multithread our solution and optimize the constant factor overhead.
- Analyze the variance of the solution ahead of time to determine how many samples are needed for a sufficiently accurate estimate.

This optimization can be a tedious but rather straightforward process, and the main challenge is not underestimating the complexity of the potential solutions, or spending too much time on implementation. Notably, the stable region isn't necessarily convex, isn't necessarily entirely located within the polygon, isn't necessarily connected, and can even be connected but with holes.


Code: import java.io.PrintWriter;
import java.text.DecimalFormat;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Random;
import java.util.Scanner;

public class CupCounterbalancing {
  static PrintWriter out;

  static void solve(Scanner fs, PrintWriter out, int tt) throws InterruptedException {
    int n = fs.nextInt(), r = fs.nextInt(), max = fs.nextInt();
    Vec[] points = new Vec[n];
    for (int i = 0; i < n; i++) {
      points[i] = new Vec(fs.nextInt(), fs.nextInt());
    }
    double ans = solveNaive1e8(points, r, max);
    out.println("Case #" + tt + ": " + ans);
    out.flush();
  }

  static double solveNaive1e8(Vec[] points, int r, int max) {
    long time = System.currentTimeMillis();
    Vec test = new Vec(1, 0).rotate(random.nextDouble()).scale(max / Math.sqrt(1e8));
    double res = solveNaive(points, r, max, test);
    return res;
  }
  static double solveNaive1e8Multithreaded(
      Vec[] points, int r, int max, boolean print, int nThreads) throws InterruptedException {
    long time = System.currentTimeMillis();
    Thread[] threads = new Thread[nThreads];
    final double[] answers = new double[nThreads];
    for (int i = 0; i < threads.length; i++) {
      final int ii = i;
      threads[i] = new Thread(null, null, "T") {
        public void run() {
          answers[ii] = solveNaive1e8(points, r, max);
        }
      };
      threads[i].start();
    }
    for (Thread t : threads) {
      t.join();
    }
    double average = 0;
    for (double d : answers) {
      average += d;
    }
    return average / nThreads;
  }

  static Random random = new Random();

  static double solveNaive(Vec[] points, int r, int max, Vec test) {
    test.x = Math.abs(test.x);
    test.y = Math.abs(test.y);
    Vec start = new Vec(0, 0);
    int row = 0;
    Vec rotated = test.rot90();
    Seg[] segs = new Seg[points.length];
    for (int i = 0; i < points.length; i++) {
      segs[i] = new Seg(points[i], points[(i + 1) % points.length]);
    }
    int worksCount = 0;
    int allCount = 0;
    while (true) {
      Vec rowStart = start.add(rotated.scale(row));
      if (rowStart.y > max)
        break;
      Vec at = rowStart;
      while (at.x <= max && at.y <= max) {
        if (at.x >= 0 && at.y >= 0) {
          if (worksNaive(segs, at, r)) {
            worksCount++;
          }
          allCount++;
        }
        at = at.add(test);
      }
      row++;
    }
    row = -1;
    while (true) {
      Vec rowStart = start.add(rotated.scale(row));
      Vec at = rowStart;
      boolean hit = false;
      while (at.x <= max && at.y <= max) {
        if (at.x >= 0 && at.y >= 0) {
          hit = true;
          if (worksNaive(segs, at, r)) {
            worksCount++;
          }
          allCount++;
        }
        at = at.add(test);
      }
      row--;
      if (!hit) {
        break;
      }
    }
    return worksCount / (double) allCount;
  }

  static boolean worksNaive(Seg[] segs, Vec query, long r) {
    ArrayList<Vec> intersections = new ArrayList<>();
    Circle c = new Circle(query, r);
    for (Seg s : segs) {
      Vec[] ans = c.intersectSeg(s);
      if (ans != null) {
        for (Vec v : ans) {
          intersections.add(v);
        }
      }
    }
    double[] angles = new double[intersections.size()];
    for (int i = 0; i < angles.length; i++) angles[i] = intersections.get(i).sub(query).angle();
    Arrays.sort(angles);
    if (angles.length <= 2) {
      return false;
    }
    for (int i = 1; i < angles.length; i++) {
      if (angles[i] - angles[i - 1] >= Math.PI) {
        return false;
      }
    }
    if (angles[0] - (angles[angles.length - 1] - Math.PI * 2) > Math.PI) {
      return false;
    }
    return true;
  }

  static class Vec {
    static final double EPS = 1e-6;
    double x, y;
    public Vec(double x, double y) {
      this.x = x;
      this.y = y;
    }
    public Vec add(Vec o) { return new Vec(x + o.x, y + o.y); }
    public Vec sub(Vec o) { return new Vec(x - o.x, y - o.y); }
    public Vec scale(double s) { return new Vec(x * s, y * s); }
    public double dot(Vec o) { return x * o.x + y * o.y; }
    public double cross(Vec o) { return x * o.y - y * o.x; }
    public double mag2() { return dot(this); }
    public double mag() { return Math.sqrt(mag2()); }
    public Vec unit() { return scale(1 / mag()); }
    public Vec rot90() { return new Vec(-y, x); }
    public Vec rot270() { return new Vec(y, -x); }

    public Vec rotate(double theta) {
      double PI = Math.PI;
      double newX = x * Math.cos(theta) + y * Math.cos(PI / 2 + theta);
      double newY = x * Math.sin(theta) + y * Math.sin(PI / 2 + theta);
      return new Vec(newX, newY);
    }

    // Angle between 0 and 2PI
    public double angle() {
      return (Math.atan2(y, x) + 2 * Math.PI) % (2 * Math.PI);
    }

    public String toString() {
      DecimalFormat df = new DecimalFormat("#.##");
      return "(" + df.format(x) + ", " + df.format(y) + ")";
    }

    static boolean eq(double a, double b) {
      return Math.abs(a - b) < EPS;
    }

    static boolean leq(double a, double b) {
      return a - EPS < b;
    }
    static boolean geq(double a, double b) {
      return a + EPS > b;
    }

    public boolean equals(Object oo) {
      Vec o = (Vec) oo;
      return eq(x, o.x) && eq(y, o.y);
    }
  }

  static class Seg {
    Vec from, to, dir;

    public Seg(Vec from, Vec to) {
      this.from = from;
      this.to = to;
      dir = to.sub(from);
    }

    // Line-line intersection
    public Vec lineIntersect(Seg o) {
      double det = o.dir.x * dir.y - dir.x * o.dir.y;
      if (Vec.eq(det, 0))
        return null;
      double dist = (o.dir.x * (o.from.y - from.y) - o.dir.y * (o.from.x - from.x)) / det;
      return from.add(dir.scale(dist));
    }

    public boolean containsPoint(Vec o) {
      double distFromLine = dir.unit().cross(o.sub(from));
      if (!Vec.eq(distFromLine, 0)) {
        return false;
      }
      return Vec.eq(dir.mag(), from.sub(o).mag() + to.sub(o).mag());
    }

    // Seg-seg intersection
    public Vec segIntersection(Seg o) {
      Vec intersect = lineIntersect(o);
      if (intersect == null) {
        return null;
      }
      return containsPoint(intersect) && o.containsPoint(intersect) ? intersect : null;
    }

    // Returns 1 if above, 0 if on, -1 if below.
    public int side(Vec o) {
      Vec oDir = o.sub(from);
      double distFromLine = dir.unit().cross(oDir);
      return Vec.eq(distFromLine, 0) ? 0 : (int)Math.signum(distFromLine);
    }

    public boolean intersects(Seg o) {
      return side(o.from) != side(o.to) && o.side(from) != o.side(to);
    }

    public Vec getClosestTo(Vec o) {
      double percentThere = o.sub(from).dot(dir) / dir.mag2();
      return from.add(dir.scale(Math.max(0, Math.min(1, percentThere))));
    }

    public Vec projectToLine(Vec o) {
      return dir.scale(o.sub(from).dot(dir) / dir.mag2()).add(from);
    }

    // Returns the shortest segment from the x-axis to the point to the y axis.
    // toContain must be a point with positive x and y coordinates.
    public static Seg getShortestSegFromAxesContainingQ1Point(Vec toContain) {
      double slope = -Math.pow(toContain.y / toContain.x, 1.0 / 3);
      double b = toContain.y - toContain.x * slope;
      double xInt = -b / slope;
      return new Seg(new Vec(0, b), new Vec(xInt, 0));
    }

    public String toString() {
      return from + " -> " + to;
    }
  }

  static class Circle {
    Vec c;
    double r;

    public Circle(Vec c, double r) {
      this.c = c;
      this.r = r;
    }

    public boolean contains(Vec v) {
      return c.sub(v).mag2() - Vec.EPS * Vec.EPS <= r * r;
    }

    // When standing at this circle, returns right tangent, then left tangent.
    public Vec[] intersect(Circle o) {
      if (c.equals(o.c)) {
        return null;
      }
      Vec dir = o.c.sub(c);
      double d2 = dir.mag2(), d = Math.sqrt(d2);
      if (r + o.r < d || r + d < o.r || o.r + d < r) {
        return null;
      }
      if (Vec.eq(r + o.r, d) || Vec.eq(o.r + d, r)) {
        return new Vec[] {c.add(dir.scale(r / d))};
      }
      if (Vec.eq(r + d, o.r)) {
        return new Vec[] {c.sub(dir.scale(r / d))};
      }
      double d1 = (r * r + d2 - o.r * o.r) / (2 * d);
      double h = Math.sqrt(r * r - d1 * d1);
      Vec unitDir = dir.unit();
      Vec rInt = c.add(unitDir.scale(d1).add(unitDir.rot270().scale(h)));
      Vec lInt = c.add(unitDir.scale(d1).add(unitDir.rot90().scale(h)));
      return new Vec[] {rInt, lInt};
    }

    public double intersectionArea(Circle o) {
      double d = o.c.sub(c).mag();
      if (r + o.r < d) {
        return 0;
      }
      double minR = Math.min(r, o.r), maxR = Math.max(r, o.r), pi = Math.PI;
      if (Vec.leq(d + minR, maxR)) {
        return pi * minR * minR;
      }
      double theta1 = 2 * Math.acos((r * r + d * d - o.r * o.r) / (2 * r * d));
      double theta2 = 2 * Math.acos((o.r * o.r + d * d - r * r) / (2 * o.r * d));
      double part1Area = theta1 / 2 * r * r;
      double part2Area = theta2 / 2 * o.r * o.r;
      double tri1 = r * r * Math.sin(theta1) / 2;
      double tri2 = o.r * o.r * Math.sin(theta2) / 2;
      return part1Area + part2Area - tri1 - tri2;
    }

    // Returns right tangent, then left tangent from perspective of the point
    public Vec[] getTangentPoints(Vec p) {
      if (contains(p)) {
        return null;
      }
      double d2 = c.sub(p).mag2();
      return new Circle(p, Math.sqrt(d2 - r * r)).intersect(this);
    }

    // Line going from my left to his right, then my right to his left lines go
    // from me to him
    public Seg[] internalTangentLines(Circle o) {
      Vec[] tangentPoints = new Circle(c, r + o.r).getTangentPoints(o.c);
      Vec offset1 = tangentPoints[0].sub(o.c).rot90().unit().scale(o.r);
      Vec offset2 = tangentPoints[1].sub(o.c).rot270().unit().scale(o.r);
      return new Seg[] {
          new Seg(tangentPoints[0].add(offset1), o.c.add(offset1)),
          new Seg(tangentPoints[1].add(offset2), o.c.add(offset2))
      };
    }

    // Right external tangent, then left external tangent, from my perspective
    // lines go from me to him
    public Seg[] externalTangentLines(Circle o) {
      if (o.r > r) {
        Seg[] oAnswer = o.externalTangentLines(this);
        return new Seg[] {
            new Seg(oAnswer[1].to, oAnswer[1].from),
            new Seg(oAnswer[0].to, oAnswer[0].from)
        };
      }
      Vec[] tangentPoints = new Circle(c, r - o.r).getTangentPoints(o.c);
      Vec offset1 = tangentPoints[0].sub(o.c).rot270().unit().scale(o.r);
      Vec offset2 = tangentPoints[1].sub(o.c).rot90().unit().scale(o.r);
      return new Seg[] {
          new Seg(tangentPoints[1].add(offset2), o.c.add(offset2)),
          new Seg(tangentPoints[0].add(offset1), o.c.add(offset1))
      };
    }

    // Line (not line segment)-circle intersection in the order of line.dir
    public Vec[] intersectLine(Seg line) {
      Vec closest = line.projectToLine(c);
      double d2 = closest.sub(c).mag2();
      if (d2 > r * r) {
        return null;
      }
      double l = Math.sqrt(r * r - d2);
      if (Vec.eq(l, 0)) {
        return new Vec[] {closest};
      }
      Vec lVec = line.dir.unit().scale(l);
      return new Vec[] {closest.sub(lVec), closest.add(lVec)};
    }

    // Line segment-circle intersection
    public Vec[] intersectSeg(Seg seg) {
      Vec[] lineIntersections = intersectLine(seg);
      if (lineIntersections == null)
        return null;
      ArrayList<Vec> contained = new ArrayList<>();
      for (Vec v : lineIntersections) {
        if (seg.containsPoint(v)) {
          contained.add(v);
        }
      }
      if (contained.isEmpty()) {
        return null;
      }
      return contained.toArray(new Vec[contained.size()]);
    }

    public String toString() {
      DecimalFormat df = new DecimalFormat();
      return "center: " + c + ", r: " + df.format(r);
    }
  }

  public static void main(String[] args) throws InterruptedException {
    Scanner fs = new Scanner(System.in);
    out = new PrintWriter(System.out);
    int T = fs.nextInt();
    for (int t = 1; t <= T; t++) {
      solve(fs, out, t);
    }
    fs.close();
    out.close();
  }
}