<|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.

387. Lazy Judges
Time limit per test: 1.5 second(s)
Memory limit: 262144 kilobytes
input: standard
output: standard



It's no secret that in order to build a baseball stadium one needs to create a project first. The chief engineer in the town NN is preparing such a project now and he has met some difficulties. There is only one brigade of baseball judges in NN, and they are very lazy, so they have some restrictions for the project. As you know, a baseball field has a form of square with bases in its corners. There is one chief judge in the brigade and he wants the center of the baseball field to be located at a fixed point (0,0) (he will stay at this point and see all the course of the game). All the other n judges will also stay at some predefined points. Each of them will see only one segment with ends in the points with coordinates (xi1,yi1) and (xi2,yi2). It is also important that every base should belong to at least one of these segments. If there are several ways to build the baseball field, the chief engineer acts in the following way. He considers the set of all possible positions of the first base and chooses one of them randomly with a uniform distribution.

Your task is to calculate the expected area of the field.

Input
The first line of the input file contains an integer number n, 1 ≤ n ≤ 50. Each of the following n lines contains four integer numbers xi1, yi1, xi2, yi2 — coordinates of ends of the i-th segment. All coordinates do not exceed 100 by absolute value. It is guaranteed that there is at least one way to build the baseball field. All segments are non-degenerate. They may intersect, but it is guaranteed that they do not overlap.

Output
Output should contain one real number — the expected area of the field with relative or absolute error 10-9.

Example(s)
sample input
sample output
3
-3 -1 3 -1
-3 -1 0 2
3 -1 0 2
4.0000000000

<|response|>
## 1) Abridged problem statement (concise)

You are given **n (1…50)** line segments on the plane (they may intersect, but **no two overlap**).  
We want to place a **square centered at (0,0)**. Its **4 vertices** (bases) must each lie on **at least one** of the given segments.

Let the **first base** be one chosen vertex \(P\). Among all valid squares, the engineer considers all feasible positions of \(P\) and chooses **uniformly at random** from that feasible set:

- if it has **positive total length** (a union of subsegments), uniform by **arc length**,
- otherwise (only isolated points), uniform over those **points**.

Compute the **expected area** of the square.

---

## 2) Key observations

1. **Vertices are rotations around the origin**  
   If the first vertex is \(P=(x,y)\), the other vertices are
   \[
   R(P)=(-y,x),\quad R^2(P)=(-x,-y),\quad R^3(P)=(y,-x)
   \]
   where \(R\) is 90° counterclockwise rotation.

2. **Area depends only on \(|P|\)**  
   The diagonal endpoints are \(P\) and \(-P\), so diagonal length is \(2|P|\).  
   Square area:
   \[
   \text{area}=\frac{(2|P|)^2}{2}=2|P|^2=2(x^2+y^2)
   \]
   So we need \(\mathbb{E}[2|P|^2]\) over feasible \(P\).

3. **Feasible set is an intersection of rotated unions of segments**  
   Let \(U\) be the union of all segments. Validity requires:
   \[
   P\in U,\ R(P)\in U,\ R^2(P)\in U,\ R^3(P)\in U
   \]
   Equivalently:
   \[
   P \in V = U \cap R^{-1}(U) \cap R^{-2}(U) \cap R^{-3}(U)
   \]
   Thus feasible first-base positions are exactly \(V\), which is again a subset of segments (plus possibly isolated points).

4. **Parameterize per original segment and reduce to interval intersection on \(t\in[0,1]\)**  
   For each input segment \(S=[A,B]\), any point on it is \(P(t)=A+t(B-A)\).  
   For each rotation \(k=1,2,3\), compute the set of \(t\) for which \(P(t)\) lies on **some** rotated segment \(R^{-k}(S_j)\). This produces a union of \(t\)-intervals \(E_k\).  
   Then feasible points on this segment are exactly:
   \[
   t \in E_1 \cap E_2 \cap E_3
   \]

5. **Integration is easy because \(|P(t)|^2\) is quadratic in \(t\)**  
   If \(P(t)=A+tD\) with \(D=B-A\), then:
   \[
   |P(t)|^2 = (A+tD)\cdot(A+tD) = q_a t^2 + q_b t + q_c
   \]
   so we can integrate analytically over feasible \(t\)-intervals.

---

## 3) Full solution approach

### Step A — Geometry primitives
We need:
- rotate a point **clockwise** by 90° \(k\) times (this equals \(R^{-k}\)),
- intersection of two segments:
  - empty,
  - a single point,
  - or a collinear overlapping segment (possible with rotated copies).

### Step B — Build feasibility intervals on a segment
For each original segment \(S_i=[A,B]\):
1. Parameterize: \(P(t)=A+t(B-A)\), \(t\in[0,1]\). Let \(D=B-A\), length \(L=|D|\), \(D^2=|D|^2\).
2. For each \(k=1,2,3\):
   - For each segment \(S_j=[C,D]\), rotate both endpoints **clockwise** \(k\) times to get \(S'_j = R^{-k}(S_j)\).
   - Compute intersection of \(S_i\) with \(S'_j\):
     - If intersection is a point \(P\): map to \(t\) by projection
       \[
       t = \frac{(P-A)\cdot (B-A)}{|B-A|^2}
       \]
       store \([t,t]\).
     - If intersection is an overlapping segment \([P,Q]\): map both to \(t_P,t_Q\), store \([t_P,t_Q]\).
   - Union all these pieces over all \(j\), and **merge** into disjoint sorted intervals \(E_k\).

Now \(t\) is feasible for a valid square (with first base on \(S_i\)) iff it belongs to all three sets \(E_1,E_2,E_3\).

### Step C — Sweep to integrate over feasible parts
To avoid complex interval-intersection code, do this:
1. Collect all interval endpoints from \(E_1,E_2,E_3\) plus \(0\) and \(1\).
2. Sort + deduplicate to get candidates \(t_0<t_1<\dots<t_m\).
3. For each consecutive interval \([t_r,t_{r+1}]\):
   - take midpoint \(mid\),
   - test membership of \(mid\) in \(E_1,E_2,E_3\); if true, the whole interval is feasible.
   - Add to accumulators:
     - feasible arc length: \(L\cdot (t_{r+1}-t_r)\)
     - integral of area over arc length:
       \[
       \int 2|P|^2 \, ds = 2L \int_{t_r}^{t_{r+1}} |P(t)|^2\,dt
       \]
       where \(|P(t)|^2\) is quadratic and has a closed-form antiderivative.

### Step D — Handle “only points” case
It can happen that the feasible set \(V\) has **zero total length** (only isolated points).  
In that case the distribution is uniform over those points.

We can reuse the candidate list of \(t\) endpoints: for each candidate \(t\), if it is in all \(E_k\), record the point \(P(t)\). After processing all segments:
- deduplicate points (with small tolerance),
- average \(2|P|^2\) over them.

### Step E — Final expectation
- If total feasible length \(> 0\):  
  \[
  \mathbb{E}[\text{area}] = \frac{\int_V 2|P|^2\,ds}{\int_V ds}
  \]
- Else: arithmetic mean over feasible points.

Complexity: \(O(n^2)\) segment intersections per segment and rotation → well within limits for \(n\le 50\).

---

## 4) C++ implementation (detailed comments)

```cpp
#include <bits/stdc++.h>
using namespace std;

using coord_t = double;
static const coord_t EPS = 1e-9;
static const coord_t EPS_MERGE = 1e-12;

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

    Point operator + (const Point& o) const { return Point(x + o.x, y + o.y); }
    Point operator - (const Point& o) const { return Point(x - o.x, y - o.y); }
    Point operator * (coord_t k)     const { return Point(x * k, y * k); }

    coord_t dot(const Point& o)   const { return x*o.x + y*o.y; }
    coord_t cross(const Point& o) const { return x*o.y - y*o.x; }

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

int n;
vector<pair<Point,Point>> segs;

// Rotate point clockwise by 90 degrees k times.
// k=1: (x,y) -> (y,-x). This equals R^{-1}.
Point rotate_cw(Point p, int k) {
    k %= 4;
    for(int i=0;i<k;i++){
        p = Point(p.y, -p.x);
    }
    return p;
}

// Segment-segment intersection between [a,b] and [c,d].
// Returns:
//   empty            => no intersection
//   {P}              => single intersection point
//   {P,Q}            => overlapping segment endpoints (collinear overlap)
vector<Point> segseg_intersection(const Point& a, const Point& b,
                                  const Point& c, const Point& d) {
    Point ab = b - a;
    Point cd = d - c;
    coord_t denom = ab.cross(cd);

    // Non-parallel lines: unique intersection point of infinite lines,
    // then check if it lies within both segments via parameters.
    if (fabs(denom) > EPS) {
        coord_t t = (c - a).cross(cd) / denom; // along [a,b]
        coord_t s = (c - a).cross(ab) / denom; // along [c,d]
        if (t >= -EPS && t <= 1 + EPS && s >= -EPS && s <= 1 + EPS) {
            return { a + ab * t };
        }
        return {};
    }

    // Parallel: if not collinear -> no intersection.
    if (fabs((c - a).cross(ab)) > EPS) return {};

    // Collinear: project c and d onto parameter t of [a,b].
    coord_t len2 = ab.norm2();
    coord_t tc = (c - a).dot(ab) / len2;
    coord_t td = (d - a).dot(ab) / len2;
    if (tc > td) swap(tc, td);

    // Overlap interval in t with [0,1]
    coord_t lo = max(coord_t(0), tc);
    coord_t hi = min(coord_t(1), td);
    if (hi < lo - EPS) return {};

    Point P = a + ab * lo;
    Point Q = a + ab * hi;
    if (hi - lo < EPS) return {P};
    return {P, Q};
}

// Merge intervals [lo,hi] (already in [0,1]) with a small tolerance.
static vector<pair<coord_t,coord_t>> merge_intervals(vector<pair<coord_t,coord_t>> v) {
    sort(v.begin(), v.end());
    vector<pair<coord_t,coord_t>> res;
    for (auto [lo,hi] : v) {
        if (res.empty() || lo > res.back().second + EPS_MERGE) {
            res.push_back({lo,hi});
        } else {
            res.back().second = max(res.back().second, hi);
        }
    }
    return res;
}

static bool in_intervals(coord_t t, const vector<pair<coord_t,coord_t>>& v) {
    for (auto [lo,hi] : v) {
        if (t >= lo - EPS && t <= hi + EPS) return true;
    }
    return false;
}

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

    cin >> n;
    segs.resize(n);
    for (int i=0;i<n;i++) {
        cin >> segs[i].first.x >> segs[i].first.y
            >> segs[i].second.x >> segs[i].second.y;
    }

    // Accumulators for the 1D (positive-length) feasible set:
    // total_length = ∫_V ds
    // total_integral = ∫_V 2|P|^2 ds  (so expected area = total_integral/total_length)
    coord_t total_length = 0.0;
    coord_t total_integral = 0.0;

    // For the 0D (points-only) feasible set
    vector<Point> discrete_points;

    // Process each original segment S_i = [A,B] as the locus of the first base.
    for (auto [A,B] : segs) {
        Point D = B - A;
        coord_t L = D.norm();
        if (L < 1e-12) continue; // non-degenerate guaranteed, but be safe
        coord_t len2 = D.norm2();

        // E[k] (k=0..2 corresponds to rotation k+1) stores merged t-intervals
        // where P(t) lies in R^{-(k+1)}(U).
        vector<vector<pair<coord_t,coord_t>>> E(3);

        // Build E_1, E_2, E_3
        for (int k=1;k<=3;k++) {
            vector<pair<coord_t,coord_t>> raw;

            for (auto [C,D2] : segs) {
                // Rotate S_j clockwise by k*90° to get R^{-k}(S_j)
                Point C2 = rotate_cw(C, k);
                Point D2r = rotate_cw(D2, k);

                // Intersect S_i with rotated S_j
                auto inter = segseg_intersection(A, B, C2, D2r);
                if (inter.empty()) continue;

                // Map intersection geometry back to parameter t on S_i:
                // t = ((P-A)·(B-A)) / |B-A|^2
                auto to_t = [&](const Point& P) -> coord_t {
                    coord_t t = (P - A).dot(B - A) / len2;
                    // clamp due to numerical noise
                    return max(coord_t(0), min(coord_t(1), t));
                };

                if (inter.size() == 1) {
                    coord_t t = to_t(inter[0]);
                    raw.push_back({t,t});
                } else {
                    coord_t t1 = to_t(inter[0]);
                    coord_t t2 = to_t(inter[1]);
                    if (t1 > t2) swap(t1,t2);
                    raw.push_back({t1,t2});
                }
            }

            E[k-1] = merge_intervals(raw);
        }

        // Candidate t-values where feasibility can change:
        // endpoints of all E_k intervals plus 0 and 1.
        vector<coord_t> ts = {0.0, 1.0};
        for (int k=0;k<3;k++) {
            for (auto [lo,hi] : E[k]) {
                ts.push_back(lo);
                ts.push_back(hi);
            }
        }
        sort(ts.begin(), ts.end());

        // Dedupe candidates with tolerance
        vector<coord_t> cand;
        for (coord_t t : ts) {
            t = max(coord_t(0), min(coord_t(1), t));
            if (cand.empty() || fabs(t - cand.back()) > EPS_MERGE)
                cand.push_back(t);
        }

        // Prepare |P(t)|^2 quadratic coefficients.
        // P(t) = A + t*(B-A) = A + t*D
        // |P(t)|^2 = (A+tD)·(A+tD) = qa t^2 + qb t + qc
        coord_t qa = (B-A).norm2();
        coord_t qb = 2.0 * (A.dot(B-A));
        coord_t qc = A.norm2();

        // Antiderivative of |P(t)|^2:
        // ∫ (qa t^2 + qb t + qc) dt = qa/3 t^3 + qb/2 t^2 + qc t
        auto F = [&](coord_t t) -> coord_t {
            return (qa/3.0)*t*t*t + (qb/2.0)*t*t + qc*t;
        };

        // Sweep consecutive candidate intervals and test midpoint feasibility.
        for (int i=0;i+1<(int)cand.size();i++) {
            coord_t tlo = cand[i], thi = cand[i+1];
            if (thi - tlo < EPS_MERGE) continue;

            coord_t mid = 0.5*(tlo + thi);
            bool ok = in_intervals(mid, E[0]) && in_intervals(mid, E[1]) && in_intervals(mid, E[2]);
            if (!ok) continue;

            // ds = L dt, area = 2|P|^2
            total_length   += (thi - tlo) * L;
            total_integral += 2.0 * L * (F(thi) - F(tlo));
        }

        // Also collect feasible candidate points themselves for the 0D case.
        for (coord_t t : cand) {
            bool ok = in_intervals(t, E[0]) && in_intervals(t, E[1]) && in_intervals(t, E[2]);
            if (ok) {
                Point P = A + (B-A)*t;
                discrete_points.push_back(P);
            }
        }
    }

    cout.setf(std::ios::fixed);
    cout << setprecision(10);

    // If V has positive total length -> expected value is integral/length.
    if (total_length > 1e-9) {
        cout << (total_integral / total_length) << "\n";
        return 0;
    }

    // Otherwise, V is points-only: dedupe points and average 2|P|^2.
    vector<Point> uniq;
    for (auto p : discrete_points) {
        bool found = false;
        for (auto q : uniq) {
            if ((p - q).norm() < 1e-6) { found = true; break; }
        }
        if (!found) uniq.push_back(p);
    }

    coord_t sum_area = 0.0;
    for (auto p : uniq) sum_area += 2.0 * p.norm2();
    cout << (sum_area / (coord_t)uniq.size()) << "\n";
    return 0;
}
```

---

## 5) Python implementation (detailed comments)

```python
import sys
import math

EPS = 1e-9
EPS_MERGE = 1e-12

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, o): return Point(self.x + o.x, self.y + o.y)
    def __sub__(self, o): return Point(self.x - o.x, self.y - o.y)
    def __mul__(self, k): return Point(self.x * k, self.y * k)

    def dot(self, o):   return self.x * o.x + self.y * o.y
    def cross(self, o): return self.x * o.y - self.y * o.x

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

def rotate_cw(p: Point, k: int) -> Point:
    """Rotate clockwise by 90 degrees k times. k=1: (x,y)->(y,-x) = R^{-1}."""
    k %= 4
    x, y = p.x, p.y
    for _ in range(k):
        x, y = y, -x
    return Point(x, y)

def segseg_intersection(a: Point, b: Point, c: Point, d: Point):
    """
    Intersection of segments [a,b] and [c,d]:
      []      -> empty
      [P]     -> single intersection point
      [P,Q]   -> collinear overlap segment endpoints
    """
    ab = b - a
    cd = d - c
    denom = ab.cross(cd)

    # Non-parallel: compute intersection parameters on both segments
    if abs(denom) > EPS:
        t = (c - a).cross(cd) / denom
        s = (c - a).cross(ab) / denom
        if -EPS <= t <= 1 + EPS and -EPS <= s <= 1 + EPS:
            return [a + ab * t]
        return []

    # Parallel but not collinear
    if abs((c - a).cross(ab)) > EPS:
        return []

    # Collinear: project endpoints onto ab parameter
    len2 = ab.norm2()
    tc = (c - a).dot(ab) / len2
    td = (d - a).dot(ab) / len2
    if tc > td:
        tc, td = td, tc

    lo = max(0.0, tc)
    hi = min(1.0, td)
    if hi < lo - EPS:
        return []

    P = a + ab * lo
    Q = a + ab * hi
    if hi - lo < EPS:
        return [P]
    return [P, Q]

def merge_intervals(intervals):
    """Merge [lo,hi] intervals with small tolerance."""
    intervals.sort()
    res = []
    for lo, hi in intervals:
        if not res or lo > res[-1][1] + EPS_MERGE:
            res.append([lo, hi])
        else:
            res[-1][1] = max(res[-1][1], hi)
    return res

def in_intervals(t, intervals):
    """Check membership in union of intervals."""
    for lo, hi in intervals:
        if t >= lo - EPS and t <= hi + EPS:
            return True
    return False

def solve(segs):
    total_length = 0.0          # ∫ ds over feasible set V
    total_integral = 0.0        # ∫ 2|P|^2 ds over V
    discrete_points = []        # for zero-length feasible set

    for A, B in segs:
        D = B - A
        L = D.norm()
        if L < 1e-12:
            continue
        len2 = D.norm2()

        # Build E_k for k=1..3: t such that P(t) in R^{-k}(U)
        E = []
        for k in (1, 2, 3):
            raw = []
            for C, D2 in segs:
                C2 = rotate_cw(C, k)
                D2r = rotate_cw(D2, k)
                inter = segseg_intersection(A, B, C2, D2r)
                if not inter:
                    continue

                def to_t(P):
                    t = (P - A).dot(B - A) / len2
                    return max(0.0, min(1.0, t))

                t1 = to_t(inter[0])
                if len(inter) == 1:
                    raw.append((t1, t1))
                else:
                    t2 = to_t(inter[1])
                    if t1 > t2:
                        t1, t2 = t2, t1
                    raw.append((t1, t2))

            E.append(merge_intervals(raw))

        # Candidate t where feasibility can change: endpoints + {0,1}
        ts = [0.0, 1.0]
        for k in range(3):
            for lo, hi in E[k]:
                ts.append(lo)
                ts.append(hi)
        ts.sort()

        cand = []
        for t in ts:
            t = max(0.0, min(1.0, t))
            if not cand or abs(t - cand[-1]) > EPS_MERGE:
                cand.append(t)

        # |P(t)|^2 quadratic coefficients:
        # P(t)=A+tD, |P(t)|^2 = qa t^2 + qb t + qc
        qa = (B - A).norm2()
        qb = 2.0 * A.dot(B - A)
        qc = A.norm2()

        # Antiderivative of |P(t)|^2
        def F(t):
            return (qa/3.0)*t*t*t + (qb/2.0)*t*t + qc*t

        # Sweep intervals; use midpoint test for membership in all E_k
        for i in range(len(cand) - 1):
            tlo, thi = cand[i], cand[i+1]
            if thi - tlo < EPS_MERGE:
                continue
            mid = 0.5 * (tlo + thi)

            if (in_intervals(mid, E[0]) and
                in_intervals(mid, E[1]) and
                in_intervals(mid, E[2])):

                total_length += (thi - tlo) * L
                total_integral += 2.0 * L * (F(thi) - F(tlo))

        # Collect candidate endpoints as feasible points for 0D case
        for t in cand:
            if (in_intervals(t, E[0]) and
                in_intervals(t, E[1]) and
                in_intervals(t, E[2])):
                discrete_points.append(A + (B - A) * t)

    # Positive-length feasible set: expectation by length measure
    if total_length > 1e-9:
        return total_integral / total_length

    # Points-only feasible set: dedupe and average
    uniq = []
    for p in discrete_points:
        ok = True
        for q in uniq:
            if (p - q).norm() < 1e-6:
                ok = False
                break
        if ok:
            uniq.append(p)

    # guaranteed at least one solution
    s = 0.0
    for p in uniq:
        s += 2.0 * p.norm2()
    return s / len(uniq)

def main():
    data = sys.stdin.read().strip().split()
    if not data:
        return
    it = iter(data)
    n = int(next(it))
    segs = []
    for _ in range(n):
        x1 = int(next(it)); y1 = int(next(it))
        x2 = int(next(it)); y2 = int(next(it))
        segs.append((Point(x1, y1), Point(x2, y2)))

    ans = solve(segs)
    print(f"{ans:.10f}")

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

