### pacha2880's blog

By pacha2880, history, 2 months ago,

What is the best way to solve this problem?, I have a solution in $O(n^4)$

Given are $N$ points $(x_i, y_i)$ in a two-dimensional plane. Find the minimum radius of a circle such that all the points are inside or on it.

 » 2 months ago, # | ← Rev. 3 →   +15 I suppose you can use inversion (geometry).The idea of inversion is that if you invert about a point $P$(using radius=1, but that only comes down to a constant difference), then a circle passing through that $P$ with radius $R$ becomes a straight line with distance $\frac{1}{2R}$ away from $P$, and the points inside the circle becomes the points on the opposite side of the line as $P$.Find the convex hull of the $N$ points. Without loss of generality, suppose that the $N$ points are their own convex hull. For each point $P$ on the convex hull, we carry out inversion on the remaining $N-1$ points. Now we need to find a line farthest from the inversion center $P$ such that the line separates $P$ from the inverted $N-1$ points. To do this, we find the convex hull of these points again, and process each edge of the resulting hull. We seek the largest distance after using each of the $N$ points as inversion centers.The time complexity should be $O(n^2 \log{n})$. Correct me if the algorithm is wrong; I'm not very familiar with this either.Edit: As pointed out below, there are a variety of solutions as this is a known problem.
 » 2 months ago, # | ← Rev. 3 →   +3 This is a well-known problem called minimum enclosing circle (it has also other names).I personally use Welzl's algorithm for it with expected time complexity of $O(n)$. Code// intersect two lines bool intersection(complex a, complex b, complex c, complex d, complex &inter) { double d1 = cp(a - b, d - c), d2 = cp(a - c, d - c), d3 = cp(a - b, a - c); if (fabs(d1) < EPS) return false; double t1 = d2 / d1, t2 = d3 / d1; inter = a + (b - a) * t1; return true; } // generate circle from 3 points pair, double> circle3p(point a, point b, point c) { complex A = complex(a.X, a.Y); complex B = complex(b.X, b.Y); complex C = complex(c.X, c.Y); complex ABmid = (A + B) / 2.0, BCmid = (B + C) / 2.0; complex ABnorm = complex((A - B).Y, -(A - B).X); complex BCnorm = complex((B - C).Y, -(B - C).X); complex center; bool valid = intersection(ABmid, ABmid + ABnorm, BCmid, BCmid + BCnorm, center); assert(valid); double r = length(A - center); return {center, r}; } // get circle from 1, 2, or 3 points pair, double> get_circle(vector cir) { double r; complex c; if (cir.size() == 1) { c = complex(cir[0].X, cir[0].Y); r = 0; } if (cir.size() == 2) { c = complex(cir[0].X + cir[1].X, cir[0].Y + cir[1].Y); c.real(c.X / 2), c.imag(c.Y / 2); r = length(cir[0] - cir[1]) / 2; } else { assert(cir.size() == 3); return circle3p(cir[0], cir[1], cir[2]); } return {c, r}; } // check if a points is inside a circle bool inside_circle(point p, vector cir) { if (cir.size() == 0) return false; if (cir.size() == 1) return same_vec(p, cir[0]); auto [c, r] = get_circle(cir); return cmp(length(complex(p.X, p.Y) - c), r) != 1;// lte (<=) } // welzl algorithm implementation // WARNING: don't forget to randomly shuffle points vector vector welzl(vector &points, int i = 0, vector cir = {}) { if (cir.size() == 3 || i == points.size()) return cir; auto new_cir = welzl(points, i + 1, cir); if (inside_circle(points[i], new_cir)) return new_cir; cir.push_back(points[i]); return welzl(points, i + 1, cir); }