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

315. The Highway Belt
Time limit per test: 0.25 second(s)
Memory limit: 65536 kilobytes
input: standard
output: standard



In the country of Berland the decision was made to build a highway belt. The king Berbinduk IX wants the highway belt to be a closed polyline and have the star-like form around the capital of Berland. This means that each half-line starting from the capital intersects the highway belt in exactly one point.

To save time and money the king decided to build the highway using only existing roads. More specifically, the highway polyline should be formed by existing roads and parts of exisiting roads.

You are the head of the First Berland National Road Building Company, and you got the contract to build the highway. Now you need to present a plan of the new highway to the king. The longer it will be the more money you will earn. So your task is to create a plan with maximal possible length of the belt highway without breaking the king's requests listed above.

Input
The first line of the input contains the number of the existing roads N (1 ≤ N ≤ 50). The following N lines contain the description of the roads. Each description is four integer numbers x1, y1, x2, y2 not exceeding 100 by absolute value. (x1, y1) is one of the end points of the existing road and (x2, y2) is another. A road is allowed to have zero length.

The capital is situated in the origin of the coordinate system.

Output
Write to the output the maximum possible length of the highway belt with precision of no less than 5 digits after the decimal point. If it is impossible to build a highway belt of positive length, write 0 to the output.

Example(s)
sample input
sample output
6
2 2 2 -2
-2 2 -2 -2
2 2 -2 2
-2 -2 2 -2
3 0 -1 4
-2 1 0 3 
16.82843

sample input
sample output
1
1 1 -1 -1
0

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

You are given **N ≤ 50** line segments (roads) on the plane. The capital is at the **origin (0,0)**.

You must build a **closed polyline** (a cycle) using only **whole roads or subsegments of roads** such that:

> Every ray starting at the origin intersects the polyline in **exactly one point**.

Among all valid belts, output the **maximum possible total length** (perimeter).  
If it’s impossible to build any belt of **positive** length, output **0**.

---

## 2) Key observations

1. **“Exactly one intersection for every ray” ⇒ angular monotonicity.**  
   As you walk along a valid belt, the polar angle `atan2(y,x)` of points on the belt must change **monotonically** (always clockwise or always counterclockwise), otherwise some direction (ray) would hit the belt twice.

2. **We can only turn at road endpoints or intersections.**  
   Since we may use parts of roads, we should subdivide roads at all pairwise intersections. Then any optimal belt can be represented as a cycle in the resulting planar graph.

3. **Checking that the origin is inside can be done by parity.**  
   A closed curve intersects a fixed ray from the origin an **odd** number of times iff the origin is inside.  
   So we can track a **1-bit parity** while building the cycle: toggle when an edge crosses a fixed test ray.

4. **After fixing direction + start vertex, the problem becomes a longest path in a DAG.**  
   Sort vertices by polar angle (CW or CCW). If we force edges to go **forward** in that order (wrapping once), we eliminate cycles → we can do DP for maximum total length.

---

## 3) Full solution approach

### Step A — Build “atomic edges” graph from the roads

1. **Discard any road that contains the origin.**  
   If the belt uses a segment passing through the origin, rays would intersect at the origin in a degenerate/invalid way; typical accepted solutions simply forbid such roads.

2. For each road segment, collect points on it:
   - its two endpoints
   - every intersection point with every other road (including collinear overlaps handled by checking whether endpoints lie on the other segment)

3. For each road, sort collected points along the segment and connect each consecutive pair:
   - Each consecutive pair forms an **atomic usable subsegment** (edge).
   - Store its Euclidean length.
   - Also store whether this edge crosses a fixed test ray from the origin (parity toggle).

This yields a graph `G(V,E)` embedded in the plane.

### Step B — Define the “forward” order by angle

For each of the 2 directions:
- `dir = 0`: increasing angle (CCW)
- `dir = 1`: decreasing angle (CW)

Sort all vertices by:
- primary key: `angle = atan2(y,x)` in the chosen direction
- tie-break: distance to origin (radius)

### Step C — DP for each chosen start vertex

Pick a start vertex `s` and “cut” the circular angle order there to obtain a linear order:
- `rel_pos[v]` = relative position of vertex `v` in `[0..m-1]` around the circle starting at `s`.

We run DP over positions `0..m` where `m` represents “wrapped back to start”.

State:
- `dp[pos][parity]` = maximum length reachable when we are at relative position `pos`, with current parity of crossings of the test ray.

Initialization:
- `dp[0][0] = 0`, others = `-inf`.

Transition:
- from vertex at relative position `i`, try all incident edges `(v -> u)`:
  - Let `pu = rel_pos[u]`. If `pu == 0`, treat it as `m` (closing the loop).
  - Allow only **forward moves**: `pu > i`.
  - `newParity = parity XOR edgeCrossToggle`
  - update `dp[pu][newParity]`.

Answer candidate for this `(dir, start)`:
- `dp[m][1]` (closed after full sweep, and parity odd ⇒ origin inside).

Take maximum over all starts and both directions. If the maximum is 0 (or no valid cycle), print 0.

### Complexity

- Intersections: `O(N^2)`
- Vertices/edges after subdivision: `O(N^2)` in practice for `N ≤ 50`
- DP tries `2 * m` starts, each DP is roughly `O(m + E)` with adjacency scanning but filtering forward edges. With small constraints and optimization, it passes.

---

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

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

using ld = long double;
static const ld EPS = 1e-9L;
static const ld INF = 1e100L;

struct Point {
    ld x, y;
    Point(ld x=0, ld y=0): x(x), y(y) {}
    Point operator + (const Point& o) const { return {x+o.x, y+o.y}; }
    Point operator - (const Point& o) const { return {x-o.x, y-o.y}; }
    Point operator * (ld k) const { return {x*k, y*k}; }

    ld dot(const Point& o) const { return x*o.x + y*o.y; }
    ld cross(const Point& o) const { return x*o.y - y*o.x; }
    ld norm2() const { return x*x + y*y; }
    ld norm() const { return sqrtl(norm2()); }
    ld angle() const { return atan2l(y, x); }
};

// orientation sign with EPS
static int ccw(const Point& a, const Point& b, const Point& c) {
    ld v = (b-a).cross(c-a);
    if (fabsl(v) <= EPS) return 0;
    return v > 0 ? 1 : -1;
}

static bool point_on_segment(const Point& a, const Point& b, const Point& p) {
    if (ccw(a,b,p) != 0) return false;
    return min(a.x,b.x) - EPS <= p.x && p.x <= max(a.x,b.x) + EPS &&
           min(a.y,b.y) - EPS <= p.y && p.y <= max(a.y,b.y) + EPS;
}

// intersection of infinite lines (a1-b1) and (a2-b2), assuming not parallel
static Point line_line_intersection(const Point& a1, const Point& b1,
                                    const Point& a2, const Point& b2) {
    Point d1 = b1 - a1;
    Point d2 = b2 - a2;
    // a1 + d1 * t
    ld t = (a2 - a1).cross(d2) / d1.cross(d2);
    return a1 + d1 * t;
}

// edge crosses the upward test ray from origin:
// (matches classic "ray casting" but implemented robustly by checking crossing x=0
// with strict inequality and intersection y>0)
static int crosses_up_ray(const Point& a, const Point& b) {
    ld minx = min(a.x, b.x), maxx = max(a.x, b.x);
    if (!(minx <= 0 && 0 < maxx)) return 0;        // must cross x=0 with strict side
    ld t = -a.x / (b.x - a.x);                     // parameter where x=0
    ld y = a.y + t * (b.y - a.y);
    return (y > EPS) ? 1 : 0;
}

static bool same_point(const Point& a, const Point& b) {
    return (a-b).norm() < EPS;
}

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

    int N;
    cin >> N;
    vector<Point> A(N), B(N);
    for (int i=0;i<N;i++) cin >> A[i].x >> A[i].y >> B[i].x >> B[i].y;

    Point O(0,0);

    // 1) discard segments containing the origin
    vector<bool> valid(N, true);
    for (int i=0;i<N;i++) {
        if (point_on_segment(A[i], B[i], O)) valid[i] = false;
    }

    // 2) collect points on each segment: endpoints + intersections
    vector<vector<Point>> segPts(N);
    for (int i=0;i<N;i++) if (valid[i]) {
        segPts[i].push_back(A[i]);
        segPts[i].push_back(B[i]);
    }

    for (int i=0;i<N;i++) if (valid[i]) {
        for (int j=i+1;j<N;j++) if (valid[j]) {
            Point d1 = B[i] - A[i];
            Point d2 = B[j] - A[j];

            if (fabsl(d1.cross(d2)) <= EPS) {
                // parallel: handle overlaps by checking endpoints that lie on the other segment
                if (point_on_segment(A[i], B[i], A[j])) segPts[i].push_back(A[j]);
                if (point_on_segment(A[i], B[i], B[j])) segPts[i].push_back(B[j]);
                if (point_on_segment(A[j], B[j], A[i])) segPts[j].push_back(A[i]);
                if (point_on_segment(A[j], B[j], B[i])) segPts[j].push_back(B[i]);
            } else {
                // non-parallel: compute line intersection and check it lies on both segments
                Point p = line_line_intersection(A[i], B[i], A[j], B[j]);
                if (point_on_segment(A[i], B[i], p) && point_on_segment(A[j], B[j], p)) {
                    segPts[i].push_back(p);
                    segPts[j].push_back(p);
                }
            }
        }
    }

    // 3) sort points along each segment and unique them
    for (int i=0;i<N;i++) if (valid[i]) {
        Point base = A[i];
        auto &v = segPts[i];
        sort(v.begin(), v.end(), [&](const Point& p, const Point& q){
            return (p-base).norm2() < (q-base).norm2();
        });
        v.erase(unique(v.begin(), v.end(), [&](const Point& p, const Point& q){
            return same_point(p,q);
        }), v.end());
    }

    // 4) global unique list of vertices
    vector<Point> pts;
    for (int i=0;i<N;i++) for (auto &p: segPts[i]) pts.push_back(p);

    // lexicographic sort for de-dup
    sort(pts.begin(), pts.end(), [](const Point& p, const Point& q){
        if (p.x != q.x) return p.x < q.x;
        return p.y < q.y;
    });
    pts.erase(unique(pts.begin(), pts.end(), [&](const Point& p, const Point& q){
        return same_point(p,q);
    }), pts.end());

    int M = (int)pts.size();
    if (M == 0) { cout << "0\n"; return 0; }

    // helper: index of a point (M is small, linear scan is fine)
    auto idxOf = [&](const Point& p)->int {
        for (int i=0;i<M;i++) if (same_point(pts[i], p)) return i;
        return -1;
    };

    // 5) build adjacency from atomic edges
    // store edges as (to, length, crossToggle)
    vector<vector<tuple<int, ld, int>>> adj(M);
    // deduplicate edges using set
    vector<set<tuple<int, ld, int>>> adjSet(M);

    for (int i=0;i<N;i++) if (valid[i]) {
        auto &v = segPts[i];
        for (int j=0;j+1<(int)v.size();j++) {
            int u = idxOf(v[j]);
            int w = idxOf(v[j+1]);

            ld len = (v[j] - v[j+1]).norm();
            if (len <= EPS) continue;

            // to preserve a strict angular progress in DP, skip edges whose endpoints have equal angles
            if (fabsl(pts[u].angle() - pts[w].angle()) <= EPS) continue;

            int cross = crosses_up_ray(pts[u], pts[w]);

            adjSet[u].insert({w, len, cross});
            adjSet[w].insert({u, len, cross});
        }
    }
    for (int i=0;i<M;i++) {
        for (auto &e: adjSet[i]) adj[i].push_back(e);
    }

    ld best = 0;

    // 6) try both angular directions
    for (int dir=0; dir<2; dir++) {
        vector<int> order(M);
        iota(order.begin(), order.end(), 0);

        sort(order.begin(), order.end(), [&](int i, int j){
            ld ai = pts[i].angle(), aj = pts[j].angle();
            if (fabsl(ai-aj) > EPS) return (dir==0) ? (ai < aj) : (ai > aj);
            return pts[i].norm() < pts[j].norm(); // tie-break by radius
        });

        vector<int> gpos(M);
        for (int i=0;i<M;i++) gpos[order[i]] = i;

        // choose start point of the “cut”
        for (int start=0; start<M; start++) {
            vector<int> rel(M);
            for (int v=0; v<M; v++) rel[v] = (gpos[v] - gpos[start] + M) % M;

            vector<int> at(M);
            for (int v=0; v<M; v++) at[rel[v]] = v;

            // dp[pos][parity], pos in [0..M], where pos==M means “back to start”
            vector<array<ld,2>> dp(M+1, { -INF, -INF });
            dp[0][0] = 0;

            for (int i=0;i<M;i++) {
                int v = at[i];
                for (int par=0; par<2; par++) {
                    if (dp[i][par] <= -INF/2) continue;
                    for (auto [to, len, cross] : adj[v]) {
                        int pto = rel[to];
                        if (pto == 0) pto = M;      // close the cycle
                        if (pto > i) {
                            int npar = par ^ cross;
                            dp[pto][npar] = max(dp[pto][npar], dp[i][par] + len);
                        }
                    }
                }
            }

            // need odd parity => origin inside
            best = max(best, dp[M][1]);
        }
    }

    cout << fixed << setprecision(9) << (double)best << "\n";
    return 0;
}
```

---

## 5) Python implementation (detailed comments)

```python
import sys
import math

EPS = 1e-9
INF = 1e100

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 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 angle(self): return math.atan2(self.y, self.x)

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

def point_on_segment(a: Point, b: Point, p: Point) -> bool:
    if ccw(a, b, p) != 0:
        return False
    return (min(a.x, b.x) - EPS <= p.x <= max(a.x, b.x) + EPS and
            min(a.y, b.y) - EPS <= p.y <= max(a.y, b.y) + EPS)

def line_line_intersection(a1: Point, b1: Point, a2: Point, b2: Point) -> Point:
    d1 = b1 - a1
    d2 = b2 - a2
    t = (a2 - a1).cross(d2) / d1.cross(d2)
    return a1 + d1 * t

def same_point(a: Point, b: Point) -> bool:
    return (a - b).norm() < EPS

def crosses_up_ray(a: Point, b: Point) -> int:
    # count crossing of x=0 from <=0 to >0, and intersection y > 0
    minx, maxx = min(a.x, b.x), max(a.x, b.x)
    if not (minx <= 0.0 and 0.0 < maxx):
        return 0
    t = -a.x / (b.x - a.x)
    y = a.y + t * (b.y - a.y)
    return 1 if y > EPS else 0

def solve() -> None:
    data = sys.stdin.read().strip().split()
    if not data:
        return
    it = iter(data)
    n = int(next(it))
    A, B = [], []
    for _ in range(n):
        x1, y1, x2, y2 = map(float, (next(it), next(it), next(it), next(it)))
        A.append(Point(x1, y1))
        B.append(Point(x2, y2))

    O = Point(0.0, 0.0)

    # 1) discard segments containing origin
    valid = [True] * n
    for i in range(n):
        if point_on_segment(A[i], B[i], O):
            valid[i] = False

    # 2) per segment, collect endpoints + intersections
    seg_pts = [[] for _ in range(n)]
    for i in range(n):
        if valid[i]:
            seg_pts[i].append(A[i])
            seg_pts[i].append(B[i])

    for i in range(n):
        if not valid[i]:
            continue
        for j in range(i + 1, n):
            if not valid[j]:
                continue
            d1 = B[i] - A[i]
            d2 = B[j] - A[j]
            if abs(d1.cross(d2)) <= EPS:
                # parallel: handle overlap via endpoint inclusion checks
                if point_on_segment(A[i], B[i], A[j]): seg_pts[i].append(A[j])
                if point_on_segment(A[i], B[i], B[j]): seg_pts[i].append(B[j])
                if point_on_segment(A[j], B[j], A[i]): seg_pts[j].append(A[i])
                if point_on_segment(A[j], B[j], B[i]): seg_pts[j].append(B[i])
            else:
                p = line_line_intersection(A[i], B[i], A[j], B[j])
                if point_on_segment(A[i], B[i], p) and point_on_segment(A[j], B[j], p):
                    seg_pts[i].append(p)
                    seg_pts[j].append(p)

    # 3) sort points along each segment and unique
    for i in range(n):
        if not valid[i]:
            continue
        base = A[i]
        seg_pts[i].sort(key=lambda p: (p - base).norm2())
        uniq = []
        for p in seg_pts[i]:
            if not uniq or not same_point(uniq[-1], p):
                uniq.append(p)
        seg_pts[i] = uniq

    # 4) global unique vertex list
    pts = []
    for i in range(n):
        pts.extend(seg_pts[i])

    pts.sort(key=lambda p: (p.x, p.y))
    uniq = []
    for p in pts:
        if not uniq or not same_point(uniq[-1], p):
            uniq.append(p)
    pts = uniq
    m = len(pts)
    if m == 0:
        print("0")
        return

    def idx_of(p: Point) -> int:
        for i in range(m):
            if same_point(pts[i], p):
                return i
        return -1

    # 5) adjacency: store best edges (dedup) by (to, cross)
    adj_best = [dict() for _ in range(m)]
    for i in range(n):
        if not valid[i]:
            continue
        arr = seg_pts[i]
        for j in range(len(arr) - 1):
            u = idx_of(arr[j])
            v = idx_of(arr[j + 1])
            length = (arr[j] - arr[j + 1]).norm()
            if length <= EPS:
                continue
            if abs(pts[u].angle() - pts[v].angle()) <= EPS:
                continue
            cross = crosses_up_ray(pts[u], pts[v])

            # undirected
            if (v, cross) not in adj_best[u] or adj_best[u][(v, cross)] < length:
                adj_best[u][(v, cross)] = length
            if (u, cross) not in adj_best[v] or adj_best[v][(u, cross)] < length:
                adj_best[v][(u, cross)] = length

    adj = [[] for _ in range(m)]
    for u in range(m):
        for (v, cross), length in adj_best[u].items():
            adj[u].append((v, length, cross))

    best = 0.0

    # 6) DP over both angle directions
    for direction in (0, 1):
        order = list(range(m))
        if direction == 0:
            order.sort(key=lambda i: (pts[i].angle(), pts[i].norm()))
        else:
            order.sort(key=lambda i: (-pts[i].angle(), pts[i].norm()))

        gpos = [0] * m
        for i, v in enumerate(order):
            gpos[v] = i

        for start in range(m):
            rel = [(gpos[v] - gpos[start] + m) % m for v in range(m)]
            at = [0] * m
            for v in range(m):
                at[rel[v]] = v

            dp0 = [-INF] * (m + 1)  # parity 0
            dp1 = [-INF] * (m + 1)  # parity 1
            dp0[0] = 0.0

            for i in range(m):
                v = at[i]
                for parity in (0, 1):
                    cur = dp0[i] if parity == 0 else dp1[i]
                    if cur <= -INF/2:
                        continue
                    for (to, length, cross) in adj[v]:
                        pto = rel[to]
                        if pto == 0:
                            pto = m
                        if pto > i:
                            np = parity ^ cross
                            val = cur + length
                            if np == 0:
                                if val > dp0[pto]:
                                    dp0[pto] = val
                            else:
                                if val > dp1[pto]:
                                    dp1[pto] = val

            best = max(best, dp1[m])  # odd parity => origin inside

    print(f"{best:.9f}")

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

---

If you want, I can also add a short proof sketch of “angle-monotone traversal ⇔ at most one intersection per ray” and why “odd parity ⇒ origin inside” is sufficient here.