Grand Finale: Circles
Contest: Codeforces Round 930 (Div. 1)
Difficulty rating: 3300
Tags: binary search, geometry

You are given $$$n$$$ circles on the plane. The $$$i$$$-th of these circles is given by a tuple of integers $$$(x_i, y_i, r_i)$$$, where $$$(x_i, y_i)$$$ are the coordinates of its center, and $$$r_i$$$ is the radius of the circle.

Please find a circle $$$C$$$ which meets the following conditions:

- $$$C$$$ is contained inside all $$$n$$$ circles given in the input.
- Among all circles $$$C$$$ that meet the first condition, the radius of the circle is maximum.

Let the largest suitable circle have the radius of $$$a$$$.

Your output $$$C$$$, described as $$$(x,y,r)$$$, will be accepted if it meets the following conditions:

- For each $$$i$$$, $$$\sqrt{(x_i-x)^2+(y_i-y)^2}+ r \le r_i+\max(a,1)\cdot 10^{-7}$$$.
- The absolute or relative error of $$$r$$$ does not exceed $$$10^{-7}$$$. Formally, your answer is accepted if and only if $$$\frac{\left|r - a\right|}{\max(1, a)} \le 10^{-7}$$$.

Input format:
The first line contains a single integer $$$n$$$ ($$$1 \le n \le 10^5$$$) — the number of circles.

The $$$i$$$-th of the following $$$n$$$ lines contains three integers $$$x_i$$$, $$$y_i$$$, $$$r_i$$$ ($$$-10^6 \le x_i,y_i \le 10^6$$$, $$$1 \le r_i \le 2 \cdot 10^6$$$).

It is guaranteed that there is a circle with a radius of at least $$$10^{-6}$$$ which is contained inside all $$$n$$$ circles.

Output format:
Output three real values, $$$x$$$, $$$y$$$, and $$$r$$$ — the coordinates of the center and the radius of the circle.

Constraints:
Time limit: 2.0 s
Memory limit: 512.0 MB

Note:
A two-dimensional plot depicting the first test case is given below. The output circle $$$C$$$ is dashed with blue lines.

A two-dimensional plot depicting the second test case is given below. The output circle $$$C$$$ is dashed with blue lines.

Example input:
4
1 1 3
-1 1 3
1 -1 2
-1 -1 2

Example output:
0.0000000000000000 -0.7637626158259733 0.9724747683480533

Example input:
4
41580 -23621 95642
-41580 -23621 95642
0 47821 95642
0 0 109750

Example output:
0.0000000000000000 0.0000000000000000 47821.0000000000000000

Editorial:
For a given center coordinate $$$(x,y)$$$, we can model the objective function $$$r=f(x,y)$$$ to maximize as follows.
$$$$$$f(x,y)=\min_{1 \le i \le n} {\left({r_i-\sqrt{(x_i-x)^2+(y_i-y)^2}}\right)}$$$$$$
Formally, this can be modeled as follows:
$$$$$$\begin{split} \max & \ r \\ \text{s.t.} & \ \sqrt{(x_i-x)^2+(y_i-y)^2} \le r_i-r \ \ \forall \ 1 \le i \le n \end{split}$$$$$$
This is a second order cone program. Generally, this is not as simple to solve as linear or quadratic programs. There are multiple solvers out there, but all of them are either commercial or simply unfit for use on Codeforces. And even if we use general approaches to second order cone programming, their time complexity is at least $$$\mathcal{O}(n^{3.5})$$$, so they cannot work. Almost all general approaches to convex programming do not help very much for this task also, because of the following reasons.
• Gradient Descent cannot work in this task. The objective function is not smooth.
• Coordinate Descent also cannot work in this task for the same reason, and it is not hard to find a countercase.
• Subgradient methods can solve this task, but they take $$$\mathcal{O}(1/\epsilon^2)$$$ iterations to converge. This is a very big dependency on the precision, so they take too much time.
• We do not expect anyone to write Adaptive Coordinate Descent during contest time, but we still tested the PRAXIS algorithm and checked that it does not pass.
• Using nested ternary searches does not help very much. The time complexity will be $$$\mathcal{O}(n \log^2 (1/\epsilon))$$$, but the constants are large due to multiple function calls and floating point operations.
• Heuristics such as the Nelder-Mead method are not proven to converge. In many cases we found it not converging successfully.
In the editorial, we will explain a solution of expected $$$\mathcal{O}(n)$$$ time complexity, with constants depending on the number of dimensions (which is $$$2$$$ in our case).
First, we observe that any answer can be described as a "basis" of at most $$$3$$$ selected circles. For $$$k=0,1,2,3$$$, the basis is as follows.
• Case 1: $$$k=0$$$
This is the case of no circles, which can be regarded as an "identity element" of bases. Some very huge circle that will enclose any given circle, for example the circle with $$$10^{18}$$$ as radius and $$$(0,0)$$$ as center, can represent this case.
• Case 2: $$$k=1$$$
This is the case of a single circle. The circle that represents this case is the one selected circle.
• Case 3: $$$k=2$$$
This is the case of two circles, for which the basis is the largest circle in the intersection of two circles. We can find this by choosing the center as the midpoint of the two points $$$A$$$ and $$$B$$$, where $$$A$$$ and $$$B$$$ are the points where the line segment between circle centers and the circle itself intersect.
• Case 4: $$$k=3$$$
This is the case of three circles, and the most tricky of all four cases. It can be found as one solution to the following system of quadratic equations:
$$$$$$(x_A-x)^2+(y_A-y)^2 = (r_A-r)^2 \\ (x_B-x)^2+(y_B-y)^2 = (r_B-r)^2 \\ (x_C-x)^2+(y_C-y)^2 = (r_C-r)^2$$$$$$
This can be solved algebraically. If you are used to solving geometry problems in Olympiads, you may know the "Problem of Apollonius". This is one special case of the problem, which is not very often mentioned because often the three circles do not intersect with each other (thus most of the time one solution for $$$r$$$ must be negative). Or if one does not know the problem, they can find that the locus of the circle center that meets with two circle edges (as in Case 3) is a conic section, precisely one side of a pair of hyperbola. By finding the intersections of two hyperbola, one can find the center of this circle. Thus, this circle can be found in $$$\mathcal{O}(1)$$$ time, with constants based on number of dimensions.
By the system of quadratic equations stated above, we see that the solution of the system must be unique (one is inside the intersection, one is outside the union, and we want the former). Thus, for some $$$n>3$$$, the bases can be represented with some $$$k \le 3$$$ circles.
Do note that, already, we have an $$$\mathcal{O}(n^3)$$$ solution based on this. We enumerate all possible bases, and find the smallest one out of them (since that one will be the one that satisfies the conditions).
Now, for some basis circle $$$(x,y,r)$$$, we can check whether some given circle $$$(x_i,y_i,r_i)$$$ is "violated" by this basis (i.e. this basis is not inside this circle) in $$$\mathcal{O}(1)$$$ by comparing the distance. We take an incremental approach. We add circles to the basis one by one, up to $$$3$$$ circles. We choose to change the basis if only the next added circle is violated by the current basis. Of course, still this takes $$$\mathcal{O}(n^3)$$$ time, and sorting the circles based on any argument will not help anyways.
However, adding the circles in random order will make it only take expected $$$\mathcal{O}(n)$$$ time. This can be analyzed by the probability that a new basis will be added. The precise analysis for the time complexity is left as a practice for the reader.
If you have seen the minimum enclosing circle problem, the Welzl's algorithm to solve that precise problem will seem very similar to this. In fact, the suggested algorithm itself is not very different from Welzl's algorithm. If you want, though, you may take an iterative approach instead of modifying Welzl's algorithm. Using three nested loops will do.
There is only two issues left. The first is the issue of numerical instability. Sometimes, due to numerical instability, the algorithm will return a NaN value. It is not that hard to solve, we can simply shuffle again and rerun the algorithm. The given time limit will be sufficient for doing this unless the constant is too big. The second is more tricky, and it is about hacks.
If you use a fixed seed for randomization (or time-based seeds to some extent), someone may hack you by adding an adversarial test case which makes the time complexity explode back to $$$\mathcal{O}(n^3)$$$. This is very bad. How do we solve this? To intervent this, we will make the algorithm halt when the number of iterations exceed a certain limit. Formally, let us choose a constant $$$c$$$, and the algorithm will halt and rerun when the number of iterations exceeds $$$c \cdot n$$$. Then, if the probability of the algorithm successfully running in this number of iterations is $$$p$$$, the expected number of reruns needed is $$$\mathcal{O}(1/p)$$$, and hacking is much harder. The algorithm will successfully find the solution in $$$\mathcal{O}((c/p) \cdot n)$$$ time. Of course, the probability $$$p$$$ depends on $$$c$$$, and you must tune the value as needed. Empirically $$$c=30$$$ runs very well, almost always taking no longer than $$$0.4$$$ seconds under proper optimizations, and $$$1.5$$$ seconds without.