<|instruction|>
Solve the below problem. The solution should start with an abridged problem statement. Then key observations. Then full solution based on the observations. Then C++ and Python implementations with comments.

332. Largest Circle
Time limit per test: 0.25 second(s)
Memory limit: 65536 kilobytes
input: standard
output: standard



Yay! You've finally managed to buy a tiny plot of land just about 400 kilometers from Moscow. Only... only it was so expensive that you are no more able to afford building a house there. So you decide to go for a swimming pool. It should have the form of circle, and be as big as possible inside your plot. Having carefully measured the border of your plot, you know now that it is a convex N-gon.

What is the largest possible radius of a circular pool in it?

Input
The first line of input contains an integer N, . The next N lines contain two integers each, xi and yi, not exceeding 107 by absolute value — the coordinates of the vertices of the plot (a convex polygon) in the counter-clockwise direction. No three vertices lie on the same line.

Output
Output the sought radius. Your solution will be accepted if it is within 10-3 of the correct one.

Example(s)
sample input
sample output
4
0 0
1 0
1 1
0 1
0.5

<|response|>
## 1. Abridged problem statement

You are given a **convex polygon** with \(N\) vertices in **counter-clockwise** order (no three collinear). Find the **maximum radius** \(R\) of a circle that can be placed **entirely inside** the polygon. Output \(R\) with absolute/relative error up to \(10^{-3}\).

---

## 2. Key observations needed to solve the problem

1. **Largest inscribed circle = maximize minimum distance to sides**
   A circle of radius \(R\) fits in the polygon iff there exists a center point \(C\) such that for every polygon side \(i\),
   \[
   \text{dist}(C, \text{side}_i) \ge R.
   \]

2. **Distance-to-line constraint becomes a half-plane constraint**
   For a directed edge \(A \to B\) of a CCW convex polygon, the polygon interior lies on the **left** side of the directed line.
   Requiring distance \(\ge R\) from that side is equivalent to shifting the line **inward** by \(R\) along the inward unit normal, and requiring \(C\) to be on the left of the shifted line.

3. **Feasibility is monotone in \(R\)**
   As \(R\) increases, the "allowed region" for the center (intersection of shifted half-planes) only shrinks.
   Therefore: if radius \(R\) is feasible, any \(R' < R\) is also feasible → use **binary search** on \(R\).

4. **Feasibility test = half-plane intersection (HPI)**
   For a given \(R\), compute the intersection of all shifted half-planes.
   Non-empty intersection ⇒ feasible.

---

## 3. Full solution approach

### Step A: Build half-planes for a candidate radius \(R\)
For each polygon edge \(A \to B\):

- Direction vector: \(d = B - A\)
- Left normal (points inward for CCW polygon): \(\text{perp}(d) = (-d_y, d_x)\)
- Unit inward normal: \(\hat{n} = \text{perp}(d) / \|d\|\)
- Shift the line inward by \(R\):
  \[
  A' = A + \hat{n}R,\quad B' = B + \hat{n}R
  \]
Now the valid centers must satisfy: "point is on the left of each directed line \(A' \to B'\)".

So for radius \(R\), we have \(N\) half-planes.

### Step B: Check if intersection is non-empty (Half-Plane Intersection)
Use a standard HPI algorithm (deque-based) that:
- sorts lines by direction angle,
- maintains a deque of candidate half-planes,
- removes planes that make the intersection empty.

If after processing, fewer than 3 constraints remain (or a contradiction is found), intersection is empty.

### Step C: Binary search the maximum radius
- low = 0
- high = \(10^7\) (safe bound from coordinate constraints; answer cannot exceed plot size scale)
- repeat ~100 iterations:
  - mid = (low+high)/2
  - if feasible(mid): low = mid else high = mid

Print `low` (maximum feasible radius found).

### Complexity
- Each feasibility check: \(O(N \log N)\) due to sorting, but we can sort the original edges once and keep the order (as in the reference), making checks close to \(O(N)\).
- Binary search: 100 checks.
- Works fast enough in typical constraints for convex polygons.

---

## 4. C++ Solution

```cpp
#include <bits/stdc++.h>
// #include <coding_library/geometry/point.hpp>
// #include <coding_library/geometry/polygon.hpp>

using namespace std;

template<typename T1, typename T2>
ostream& operator<<(ostream& out, const pair<T1, T2>& x) {
    return out << x.first << ' ' << x.second;
}

template<typename T1, typename T2>
istream& operator>>(istream& in, pair<T1, T2>& x) {
    return in >> x.first >> x.second;
}

template<typename T>
istream& operator>>(istream& in, vector<T>& a) {
    for(auto& x: a) {
        in >> x;
    }
    return in;
};

template<typename T>
ostream& operator<<(ostream& out, const vector<T>& a) {
    for(auto x: a) {
        out << x << ' ';
    }
    return out;
};

using coord_t = double;
const double INF = 1e7;

struct Point {
    static constexpr coord_t eps = 1e-9;

    coord_t x, y;
    Point(coord_t x = 0, coord_t y = 0) : x(x), y(y) {}

    Point operator+(const Point& p) const { return Point(x + p.x, y + p.y); }
    Point operator-(const Point& p) const { return Point(x - p.x, y - p.y); }
    Point operator*(coord_t c) const { return Point(x * c, y * c); }
    Point operator/(coord_t c) const { return Point(x / c, y / c); }

    coord_t operator*(const Point& p) const { return x * p.x + y * p.y; }
    coord_t operator^(const Point& p) const { return x * p.y - y * p.x; }

    bool operator==(const Point& p) const { return x == p.x && y == p.y; }
    bool operator!=(const Point& p) const { return x != p.x || y != p.y; }
    bool operator<(const Point& p) const {
        return x != p.x ? x < p.x : y < p.y;
    }
    bool operator>(const Point& p) const {
        return x != p.x ? x > p.x : y > p.y;
    }
    bool operator<=(const Point& p) const {
        return x != p.x ? x < p.x : y <= p.y;
    }
    bool operator>=(const Point& p) const {
        return x != p.x ? x > p.x : y >= p.y;
    }

    coord_t norm2() const { return x * x + y * y; }
    coord_t norm() const { return sqrt(norm2()); }
    coord_t angle() const { return atan2(y, x); }

    Point rotate(coord_t a) const {
        return Point(x * cos(a) - y * sin(a), x * sin(a) + y * cos(a));
    }

    Point perp() const { return Point(-y, x); }
    Point unit() const { return *this / norm(); }
    Point normal() const { return perp().unit(); }
    Point project(const Point& p) const {
        return *this * (*this * p) / norm2();
    }
    Point reflect(const Point& p) const {
        return *this * 2 * (*this * p) / norm2() - p;
    }

    friend ostream& operator<<(ostream& os, const Point& p) {
        return os << p.x << ' ' << p.y;
    }
    friend istream& operator>>(istream& is, Point& p) {
        return is >> p.x >> p.y;
    }

    friend int ccw(const Point& a, const Point& b, const Point& c) {
        coord_t v = (b - a) ^ (c - a);
        if(-eps <= v && v <= eps) {
            return 0;
        } else if(v > 0) {
            return 1;
        } else {
            return -1;
        }
    }

    friend bool point_on_segment(
        const Point& a, const Point& b, const Point& p
    ) {
        return ccw(a, b, p) == 0 && p.x >= min(a.x, b.x) - eps &&
               p.x <= max(a.x, b.x) + eps && p.y >= min(a.y, b.y) - eps &&
               p.y <= max(a.y, b.y) + eps;
    }

    friend bool point_in_triangle(
        const Point& a, const Point& b, const Point& c, const Point& p
    ) {
        int d1 = ccw(a, b, p);
        int d2 = ccw(b, c, p);
        int d3 = ccw(c, a, p);
        return (d1 >= 0 && d2 >= 0 && d3 >= 0) ||
               (d1 <= 0 && d2 <= 0 && d3 <= 0);
    }

    friend Point line_line_intersection(
        const Point& a1, const Point& b1, const Point& a2, const Point& b2
    ) {
        return a1 +
               (b1 - a1) * ((a2 - a1) ^ (b2 - a2)) / ((b1 - a1) ^ (b2 - a2));
    }

    friend bool collinear(const Point& a, const Point& b) {
        return abs(a ^ b) < eps;
    }

    friend Point circumcenter(const Point& a, const Point& b, const Point& c) {
        Point mid_ab = (a + b) / 2.0;
        Point mid_ac = (a + c) / 2.0;
        Point perp_ab = (b - a).perp();
        Point perp_ac = (c - a).perp();
        return line_line_intersection(
            mid_ab, mid_ab + perp_ab, mid_ac, mid_ac + perp_ac
        );
    }
};

using Line = pair<Point, Point>;

class HalfPlaneIntersection {
  private:
    vector<Line> lines;
    deque<int> dq;
    bool empty_intersection = false;

    Point dir(int i) const { return lines[i].second - lines[i].first; }

    bool outside(int i, const Point& pt) const {
        return ccw(lines[i].first, lines[i].second, pt) < 0;
    }

    Point inter(int i, int j) const {
        return line_line_intersection(
            lines[i].first, lines[i].second, lines[j].first, lines[j].second
        );
    }

    bool is_parallel(int i, int j) const {
        return abs(dir(i) ^ dir(j)) < Point::eps;
    }

    bool same_direction(int i, int j) const { return (dir(i) * dir(j)) > 0; }

  public:
    static vector<Line> sort_by_angle(const vector<Line>& lines) {
        vector<Line> sorted = lines;
        sort(sorted.begin(), sorted.end(), [](const Line& a, const Line& b) {
            return (a.second - a.first).angle() < (b.second - b.first).angle();
        });
        return sorted;
    }

    HalfPlaneIntersection(const vector<Line>& lines, bool is_sorted = false)
        : lines(is_sorted ? lines : sort_by_angle(lines)) {
        int n = this->lines.size();

        for(int i = 0; i < n; i++) {
            while(dq.size() > 1 &&
                  outside(i, inter(dq.back(), dq[dq.size() - 2]))) {
                dq.pop_back();
            }
            while(dq.size() > 1 && outside(i, inter(dq.front(), dq[1]))) {
                dq.pop_front();
            }

            if(!dq.empty() && is_parallel(i, dq.back())) {
                if(!same_direction(i, dq.back())) {
                    empty_intersection = true;
                    return;
                }
                if(outside(i, this->lines[dq.back()].first)) {
                    dq.pop_back();
                } else {
                    continue;
                }
            }

            dq.push_back(i);
        }

        while(dq.size() > 2 &&
              outside(dq.front(), inter(dq.back(), dq[dq.size() - 2]))) {
            dq.pop_back();
        }
        while(dq.size() > 2 && outside(dq.back(), inter(dq.front(), dq[1]))) {
            dq.pop_front();
        }

        if(dq.size() < 3) {
            empty_intersection = true;
        }
    }

    bool is_non_empty() const { return !empty_intersection; }

    vector<Point> get_polygon() const {
        if(empty_intersection) {
            return {};
        }
        vector<Point> result(dq.size());
        for(size_t i = 0; i + 1 < dq.size(); i++) {
            result[i] = inter(dq[i], dq[i + 1]);
        }
        result.back() = inter(dq.back(), dq.front());
        return result;
    }
};

int n;
vector<Point> pts;

void read() {
    cin >> n;
    pts.resize(n);
    cin >> pts;
}

void solve() {
    // The is a direct application of half-plane intersection. We can binary
    // search for the answer R, and then shrink each of the polygon sides
    // towards their perpendiculars. The check is then to make sure that there
    // is a point inside of each of these N half-planes. There is a standard
    // randomized algorithm (add planes in random order and slice the polygon),
    // but there is also a simple algorithm introduced by Zeyuan Zhu in his
    // 2006 CTSC thesis. More details can be found here:
    //     https://cp-algorithms.com/geometry/halfplane-intersection.html

    vector<Line> edges(n);
    for(int i = 0; i < n; i++) {
        edges[i] = {pts[i], pts[(i + 1) % n]};
    }

    edges = HalfPlaneIntersection::sort_by_angle(edges);

    auto check = [&](double R) {
        vector<Line> pushed(n);
        for(int i = 0; i < n; i++) {
            Point offset = (edges[i].second - edges[i].first).normal() * R;
            pushed[i] = {edges[i].first + offset, edges[i].second + offset};
        }
        return HalfPlaneIntersection(pushed, true).is_non_empty();
    };

    double low = 0, high = 1e7;
    for(int ops = 0; ops < 100; ops++) {
        double mid = (low + high) / 2.0;
        if(check(mid)) {
            low = mid;
        } else {
            high = mid;
        }
    }

    cout << fixed << setprecision(4) << low << "\n";
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    int T = 1;
    // cin >> T;
    for(int test = 1; test <= T; test++) {
        read();
        // cout << "Case #" << test << ": ";
        solve();
    }

    return 0;
}
```

---

## 5. Python Solution

```python
import sys
import math
from collections import deque

EPS = 1e-9

class Point:
    __slots__ = ("x", "y")
    def __init__(self, x=0.0, y=0.0):
        self.x = float(x)
        self.y = float(y)

    def __add__(self, other): return Point(self.x + other.x, self.y + other.y)
    def __sub__(self, other): return Point(self.x - other.x, self.y - other.y)
    def __mul__(self, k):     return Point(self.x * k, self.y * k)
    def __truediv__(self, k): return Point(self.x / k, self.y / k)

    # dot product
    def dot(self, other): return self.x * other.x + self.y * other.y
    # cross product (2D scalar)
    def cross(self, other): return self.x * other.y - self.y * other.x

    def norm(self): return math.hypot(self.x, self.y)

    def angle(self):
        # angle of direction vector for sorting half-planes
        return math.atan2(self.y, self.x)

    def perp(self):
        # left perpendicular vector
        return Point(-self.y, self.x)

    def unit(self):
        n = self.norm()
        return self / n

    def normal(self):
        # unit left normal
        return self.perp().unit()

def ccw(a: Point, b: Point, c: Point) -> int:
    # orientation of triangle (a,b,c)
    v = (b - a).cross(c - a)
    if -EPS <= v <= EPS:
        return 0
    return 1 if v > 0 else -1

def line_intersection(a1: Point, b1: Point, a2: Point, b2: Point) -> Point:
    # intersection of infinite lines (a1->b1) and (a2->b2)
    d1 = b1 - a1
    d2 = b2 - a2
    denom = d1.cross(d2)
    # For half-plane intersection usage, we assume not parallel when called.
    t = (a2 - a1).cross(d2) / denom
    return a1 + d1 * t

class HalfPlaneIntersection:
    # Each half-plane is "left of directed line (a->b)"
    def __init__(self, lines, already_sorted=False):
        # line = (Point a, Point b)
        if not already_sorted:
            lines = sorted(lines, key=lambda ln: (ln[1] - ln[0]).angle())
        self.lines = lines
        self.dq = deque()
        self.empty = False

        def direction(i):
            a, b = self.lines[i]
            return b - a

        def outside(i, p: Point) -> bool:
            a, b = self.lines[i]
            return ccw(a, b, p) < 0  # right side => outside

        def inter(i, j) -> Point:
            a1, b1 = self.lines[i]
            a2, b2 = self.lines[j]
            return line_intersection(a1, b1, a2, b2)

        def is_parallel(i, j) -> bool:
            return abs(direction(i).cross(direction(j))) < EPS

        def same_dir(i, j) -> bool:
            return direction(i).dot(direction(j)) > 0

        n = len(self.lines)

        for i in range(n):
            # Pop from back while last intersection violates new half-plane
            while len(self.dq) > 1 and outside(i, inter(self.dq[-1], self.dq[-2])):
                self.dq.pop()

            # Pop from front similarly
            while len(self.dq) > 1 and outside(i, inter(self.dq[0], self.dq[1])):
                self.dq.popleft()

            # Handle parallel lines
            if self.dq and is_parallel(i, self.dq[-1]):
                if not same_dir(i, self.dq[-1]):
                    self.empty = True
                    return
                # Keep the stricter half-plane
                a_old, _ = self.lines[self.dq[-1]]
                if outside(i, a_old):
                    self.dq.pop()
                else:
                    continue

            self.dq.append(i)

        # Final cleanup for cyclic validity
        while len(self.dq) > 2 and outside(self.dq[0], inter(self.dq[-1], self.dq[-2])):
            self.dq.pop()
        while len(self.dq) > 2 and outside(self.dq[-1], inter(self.dq[0], self.dq[1])):
            self.dq.popleft()

        if len(self.dq) < 3:
            self.empty = True

    def is_non_empty(self) -> bool:
        return not self.empty

def solve():
    data = sys.stdin.read().strip().split()
    if not data:
        return
    it = iter(data)
    n = int(next(it))
    pts = [Point(int(next(it)), int(next(it))) for _ in range(n)]

    # Build polygon edges (CCW)
    edges = [(pts[i], pts[(i + 1) % n]) for i in range(n)]

    # Sort edges by direction angle once (important for speed)
    edges = sorted(edges, key=lambda ln: (ln[1] - ln[0]).angle())

    def check(R: float) -> bool:
        pushed = []
        for a, b in edges:
            # unit inward normal (left normal for CCW polygon)
            offset = (b - a).normal() * R
            pushed.append((a + offset, b + offset))
        return HalfPlaneIntersection(pushed, already_sorted=True).is_non_empty()

    low, high = 0.0, 1e7
    for _ in range(100):
        mid = (low + high) / 2.0
        if check(mid):
            low = mid
        else:
            high = mid

    # Print with enough precision
    sys.stdout.write(f"{low:.4f}\n")

if __name__ == "__main__":
    solve()
```
