## 1. Abridged problem statement

A colony is a connected polycube made of at most `16` axis-aligned cubes of side length `100`; cube centers are given. A robot can move only on the exposed surface faces. Across a surface edge, movement follows the physical surface: convex folds, flat continuation across adjacent cubes, and concave L-shaped inner folds are allowed.

Given two surface points `S` and `D`, neither lying on an edge, compute the shortest distance along the surface from `S` to `D`. Output the answer with absolute error at most `1e-8`.

---

## 2. Detailed editorial

### Geometry of the surface

Each cube has 6 faces. A face is part of the walkable surface iff there is no neighboring cube on the other side of that face.

We represent a surface face by:

```text
(cube grid coordinate, outward direction)
```

There are 6 possible outward directions:

```text
+x, -x, +y, -y, +z, -z
```

For every exposed face, we store:

- its 3D center,
- its four 3D corners,
- its cube grid coordinate,
- its outward direction.

The robot can cross from one surface face to another through an edge. Suppose we are on face `(C, d)` and cross an edge in tangent direction `t`.

Let:

```text
A = C + t
B = C + t + d
```

Then exactly one of the following applies:

1. `A` is absent: we fold around the edge of cube `C`, going to face `(C, t)`.
2. `A` is present but `B` is absent: we continue flat onto face `(A, d)`.
3. Both `A` and `B` are present: we enter the concave corner and go to face `(B, -t)`.

This builds a complete adjacency relation between surface faces.

---

### Why shortest paths can be reduced to vertex-to-vertex geodesics

On a polyhedral surface, a shortest path is locally straight after unfolding the crossed faces into the plane. It only changes direction at surface vertices.

Therefore, any shortest path from `S` to `D` can be decomposed into straight geodesic pieces between:

- `S`,
- surface vertices,
- `D`.

So we build a graph whose nodes are:

```text
surface vertex sheets + S + D
```

The edge weight between two nodes is the length of the shortest straight unfolded geodesic between them that does not bend at another vertex.

Then the final answer is obtained by Dijkstra on this graph.

---

### Why “vertex sheets” are needed

At a geometric cube corner, several surface face-corners may meet. Usually they form one local connected patch, but sometimes multiple patches touch only at the exact point.

The robot cannot teleport through a zero-width pinch between disconnected surface sheets.

So instead of treating a 3D vertex as one graph node, we group incident face-corners by actual surface connectivity through edges. Each connected group is called a **sheet** and becomes one graph node.

This prevents invalid shortcuts through touching-only corners.

---

### Computing edge weights by unfolding windows

For each source graph node, we compute distances to every other graph node by a surface unfolding sweep.

A sweep state, called a window, contains:

```text
current face
entry edge
depth
unfolded source point
visible angular sector
2D coordinates of the current face corners
```

Initially:

- from `S`, start in the face containing `S`;
- from a vertex sheet, start in each incident face-corner of that sheet.

When a face is unfolded into the plane, any shortest path through that face appears as a straight segment from the unfolded source point.

For each current unfolded face:

1. Check all four corners.  
   If a corner lies inside the visible sector, record its Euclidean distance from the source.
2. If this face contains `S` or `D`, check those points similarly.
3. For each outgoing edge except the one we came from:
   - clip the current angular sector to the edge segment;
   - if the clipped interval is non-empty, unfold the neighboring face across that edge;
   - push the new window.

The neighboring face is placed in the same 2D plane by preserving distances to the shared edge.

---

### Bounding the sweep

Without pruning, windows may revisit faces many times. A coarse upper bound is first computed by Dijkstra on face centers, where adjacent face centers are connected by length `100`.

During window propagation, if the closest point of a window edge segment is already farther than this upper bound, the window cannot improve the final answer and is discarded.

This keeps the search finite for `n <= 16`.

---

### Final algorithm

1. Read cubes and points.
2. Build all exposed surface faces.
3. Build face adjacency across edges.
4. Build vertex sheets using DSU.
5. Locate faces containing `S` and `D`.
6. Compute a coarse upper bound.
7. For every sheet and for `S`, run the unfolding field sweep.
8. Build the complete weighted graph implicitly from those distances.
9. Run Dijkstra from `S` to `D`.

---

## 3. Commented C++ solution

Below is the provided solution with explanatory comments.  
Some generic unused geometry helpers from the original template are omitted for readability; the algorithmic structure is unchanged.

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

// All 2D unfolding computations are done in double precision.
using coord_t = double;

// -------------------------
// 2D point for unfolded faces
// -------------------------
struct Point {
    static constexpr coord_t eps = 1e-9;

    coord_t x, y;

    // Construct point.
    Point(coord_t x = 0, coord_t y = 0) : x(x), y(y) {}

    // Vector addition.
    Point operator+(const Point& p) const {
        return Point(x + p.x, y + p.y);
    }

    // Vector subtraction.
    Point operator-(const Point& p) const {
        return Point(x - p.x, y - p.y);
    }

    // Scalar multiplication.
    Point operator*(coord_t c) const {
        return Point(x * c, y * c);
    }

    // Scalar division.
    Point operator/(coord_t c) const {
        return Point(x / c, y / c);
    }

    // Dot product.
    coord_t operator*(const Point& p) const {
        return x * p.x + y * p.y;
    }

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

    // Squared length.
    coord_t norm2() const {
        return x * x + y * y;
    }

    // Euclidean length.
    coord_t norm() const {
        return sqrt(norm2());
    }

    // Rotate vector by 90 degrees counter-clockwise.
    Point perp() const {
        return Point(-y, x);
    }
};

// 3D integer vector.
// All input cube coordinates are integral multiples of 100.
using V3 = array<int64_t, 3>;

// Subtract two 3D vectors.
V3 sub3(V3 a, V3 b) {
    return {a[0] - b[0], a[1] - b[1], a[2] - b[2]};
}

// Dot product of two 3D integer vectors.
int64_t dot3(V3 a, V3 b) {
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}

// Six axis directions.
const int dvec[6][3] = {
    {1, 0, 0},
    {-1, 0, 0},
    {0, 1, 0},
    {0, -1, 0},
    {0, 0, 1},
    {0, 0, -1}
};

// Opposite direction.
// Directions are paired as 0/1, 2/3, 4/5.
int opp(int d) {
    return d ^ 1;
}

// Number of cubes.
int n;

// Cube centers.
vector<V3> centers;

// Source and destination points.
V3 src_pt, dst_pt;

// Per exposed face data.
vector<V3> face_center;              // 3D center of face
vector<int> face_dir;                // outward normal direction
vector<array<int, 3>> face_grid;     // cube grid coordinate
vector<array<V3, 4>> face_corner;    // four 3D corners

// Map from (gx, gy, gz, direction) to face id.
map<array<int, 4>, int> face_id;

// For each face edge, adjacent face and reverse edge index.
vector<array<int, 4>> nb_face;
vector<array<int, 4>> nb_entry;

// For every face-corner, which vertex sheet it belongs to.
vector<int> sheet_of;

// Number of vertex sheets.
int num_sheets;

// Faces containing source and destination.
int face_of_src, face_of_dst;

// Set of occupied cube grid cells.
set<array<int, 3>> grid_set;

// Check whether cube exists at grid coordinate.
bool present(int gx, int gy, int gz) {
    return grid_set.count({gx, gy, gz}) > 0;
}

// Read input.
void read() {
    cin >> n;

    centers.resize(n);

    for (int i = 0; i < n; i++) {
        cin >> centers[i][0] >> centers[i][1] >> centers[i][2];
    }

    cin >> src_pt[0] >> src_pt[1] >> src_pt[2];
    cin >> dst_pt[0] >> dst_pt[1] >> dst_pt[2];
}

// Build all exposed square faces.
void build_faces() {
    // Register occupied grid cells.
    for (auto& c : centers) {
        grid_set.insert({
            int(c[0] / 100),
            int(c[1] / 100),
            int(c[2] / 100)
        });
    }

    // Examine every cube and every direction.
    for (auto& c : centers) {
        int gx = c[0] / 100;
        int gy = c[1] / 100;
        int gz = c[2] / 100;

        for (int d = 0; d < 6; d++) {
            // If another cube touches this face, it is not exposed.
            if (present(gx + dvec[d][0],
                        gy + dvec[d][1],
                        gz + dvec[d][2])) {
                continue;
            }

            // Face center is shifted 50 from cube center.
            V3 fc = {
                c[0] + 50 * dvec[d][0],
                c[1] + 50 * dvec[d][1],
                c[2] + 50 * dvec[d][2]
            };

            // Axis perpendicular to this face.
            int axis = d / 2;

            // The two axes lying in the face.
            int la = -1, lb = -1;

            for (int a = 0; a < 3; a++) {
                if (a == axis) continue;

                if (la < 0) la = a;
                else lb = a;
            }

            // Four corner sign combinations in face-local order.
            const int signs[4][2] = {
                {-1, -1},
                {1, -1},
                {1, 1},
                {-1, 1}
            };

            array<V3, 4> corners;

            // Compute all corners.
            for (int k = 0; k < 4; k++) {
                V3 p = fc;

                p[la] += 50 * signs[k][0];
                p[lb] += 50 * signs[k][1];

                corners[k] = p;
            }

            // Store face.
            face_id[{gx, gy, gz, d}] = face_center.size();
            face_center.push_back(fc);
            face_dir.push_back(d);
            face_grid.push_back({gx, gy, gz});
            face_corner.push_back(corners);
        }
    }
}

// Build adjacency between exposed faces across each edge.
void build_adjacency() {
    int nf = face_center.size();

    nb_face.resize(nf);
    nb_entry.resize(nf);

    for (int f = 0; f < nf; f++) {
        int d = face_dir[f];
        auto g = face_grid[f];

        for (int e = 0; e < 4; e++) {
            // Edge endpoints.
            V3 a = face_corner[f][e];
            V3 b = face_corner[f][(e + 1) % 4];

            // Edge midpoint.
            V3 mid = {
                (a[0] + b[0]) / 2,
                (a[1] + b[1]) / 2,
                (a[2] + b[2]) / 2
            };

            // Direction from face center to edge midpoint.
            V3 delta = sub3(mid, face_center[f]);

            int t = -1;

            // Identify tangent direction t.
            for (int dd = 0; dd < 6; dd++) {
                if (delta[0] == 50 * dvec[dd][0] &&
                    delta[1] == 50 * dvec[dd][1] &&
                    delta[2] == 50 * dvec[dd][2]) {
                    t = dd;
                }
            }

            // A = current cube plus tangent direction.
            int ax = g[0] + dvec[t][0];
            int ay = g[1] + dvec[t][1];
            int az = g[2] + dvec[t][2];

            // B = A plus outward direction.
            int cx = ax + dvec[d][0];
            int cy = ay + dvec[d][1];
            int cz = az + dvec[d][2];

            array<int, 4> key;

            // Convex turn around current cube edge.
            if (!present(ax, ay, az)) {
                key = {g[0], g[1], g[2], t};
            }
            // Flat continuation across adjacent cube.
            else if (!present(cx, cy, cz)) {
                key = {ax, ay, az, d};
            }
            // Concave inner turn.
            else {
                key = {cx, cy, cz, opp(t)};
            }

            // Adjacent face id.
            int gface = face_id[key];

            nb_face[f][e] = gface;

            // Find corresponding edge index in adjacent face.
            int found = -1;

            for (int k = 0; k < 4; k++) {
                V3 p = face_corner[gface][k];
                V3 q = face_corner[gface][(k + 1) % 4];

                if ((p == a && q == b) || (p == b && q == a)) {
                    found = k;
                }
            }

            nb_entry[f][e] = found;
        }
    }
}

// Find index of corner v in face g.
int corner_in(int g, V3 v) {
    for (int k = 0; k < 4; k++) {
        if (face_corner[g][k] == v) {
            return k;
        }
    }

    return -1;
}

// Build vertex sheets using DSU over all face-corners.
void build_sheets() {
    int nf = face_center.size();

    vector<int> par(nf * 4);

    iota(par.begin(), par.end(), 0);

    function<int(int)> find = [&](int x) {
        while (par[x] != x) {
            x = par[x] = par[par[x]];
        }

        return x;
    };

    auto uni = [&](int a, int b) {
        par[find(a)] = find(b);
    };

    // Connect two face-corners if they are connected through a surface edge.
    for (int f = 0; f < nf; f++) {
        for (int k = 0; k < 4; k++) {
            V3 v = face_corner[f][k];

            // The two edges incident to corner k.
            for (int e : {k, (k + 3) % 4}) {
                int g = nb_face[f][e];

                int kk = corner_in(g, v);

                uni(f * 4 + k, g * 4 + kk);
            }
        }
    }

    sheet_of.assign(nf * 4, -1);

    num_sheets = 0;

    map<int, int> remap;

    for (int i = 0; i < nf * 4; i++) {
        int r = find(i);

        if (!remap.count(r)) {
            remap[r] = num_sheets++;
        }

        sheet_of[i] = remap[r];
    }
}

// Find exposed face containing point p.
int find_face_of(V3 p) {
    int nf = face_center.size();

    for (int f = 0; f < nf; f++) {
        V3 fc = face_center[f];

        int d = face_dir[f];

        int axis = d / 2;

        // Point must lie in the face plane.
        if (p[axis] != fc[axis]) {
            continue;
        }

        bool ok = true;

        // Other two coordinates must be inside the square.
        for (int a = 0; a < 3; a++) {
            if (a == axis) continue;

            if (llabs(p[a] - fc[a]) >= 50) {
                ok = false;
            }
        }

        if (ok) return f;
    }

    return -1;
}

// Canonical square coordinates used when a face is first unfolded.
const Point canon[4] = {
    {0, 0},
    {100, 0},
    {100, 100},
    {0, 100}
};

// Convert a 3D point lying on face f into current unfolded 2D frame.
Point to_frame(V3 p, int f, const array<Point, 4>& c2) {
    V3 o = face_corner[f][0];

    coord_t alpha =
        (coord_t)dot3(sub3(p, o), sub3(face_corner[f][1], o)) / 10000.0;

    coord_t beta =
        (coord_t)dot3(sub3(p, o), sub3(face_corner[f][3], o)) / 10000.0;

    return c2[0] + (c2[1] - c2[0]) * alpha + (c2[3] - c2[0]) * beta;
}

// Test if point pt lies inside angular sector apex-lo-hi.
bool in_sector(Point apex, Point lo, Point hi, Point pt) {
    Point dl = lo - apex;
    Point dr = hi - apex;

    if ((dl ^ dr) < 0) {
        swap(dl, dr);
    }

    Point v = pt - apex;

    coord_t tol = 1e-9 * v.norm() + 1e-12;

    coord_t tl = tol * dl.norm();
    coord_t tr = tol * dr.norm();

    return (dl ^ v) >= -tl && (dr ^ v) <= tr;
}

// Clip parameter interval [s0,s1] by half-plane a + b*s >= 0.
void clip_half(coord_t a, coord_t b, coord_t& s0, coord_t& s1) {
    if (fabs(b) < 1e-15) {
        if (a < -1e-12) {
            s0 = 1;
            s1 = 0;
        }

        return;
    }

    coord_t r = -a / b;

    if (b > 0) {
        s0 = max(s0, r);
    } else {
        s1 = min(s1, r);
    }
}

// Clip an edge segment cd to the current angular sector.
bool clip_sector(
    Point apex,
    Point lo,
    Point hi,
    Point c,
    Point d,
    Point& q0,
    Point& q1
) {
    Point dl = lo - apex;
    Point dr = hi - apex;

    if ((dl ^ dr) < 0) {
        swap(dl, dr);
    }

    coord_t s0 = 0;
    coord_t s1 = 1;

    clip_half(dl ^ (c - apex), dl ^ (d - c), s0, s1);
    clip_half(-(dr ^ (c - apex)), -(dr ^ (d - c)), s0, s1);

    if (s0 > s1 - 1e-12) {
        return false;
    }

    q0 = c + (d - c) * s0;
    q1 = c + (d - c) * s1;

    return true;
}

// Distance from point p to segment ab.
coord_t point_seg_dist(Point p, Point a, Point b) {
    Point ab = b - a;

    coord_t len = ab * ab;

    if (len < 1e-18) {
        return (p - a).norm();
    }

    coord_t t = ((p - a) * ab) / len;

    t = max((coord_t)0.0, min((coord_t)1.0, t));

    return (p - (a + ab * t)).norm();
}

// Place a non-shared corner of neighboring face after unfolding across edge p0-p1.
Point place_corner(V3 p0, V3 p1, V3 q, Point m0, Point m1, Point apex) {
    coord_t along =
        (coord_t)dot3(sub3(q, p0), sub3(p1, p0)) / 100.0;

    coord_t dsq =
        (coord_t)dot3(sub3(q, p0), sub3(q, p0));

    coord_t h = sqrt(max((coord_t)0.0, dsq - along * along));

    Point u = (m1 - m0) * (1.0 / 100.0);

    Point base = m0 + u * along;

    Point cand1 = base + u.perp() * h;
    Point cand2 = base - u.perp() * h;

    coord_t side = (m1 - m0) ^ (apex - m0);
    coord_t s1 = (m1 - m0) ^ (cand1 - m0);

    if ((side >= 0) == (s1 >= 0)) {
        return cand2;
    }

    return cand1;
}

// Window state for unfolding sweep.
struct Win {
    int face;              // current face
    int entry;             // edge through which this face was entered
    int depth;             // unfolding depth
    Point apex;            // unfolded source point
    Point lo, hi;          // visible angular sector limits
    array<Point, 4> c2;    // unfolded positions of current face corners
};

// Global upper bound used for pruning.
coord_t ub;

// Sweep from one source node and compute direct geodesic distances.
void run_field(
    const vector<pair<int, Point>>& starts,
    int source_node,
    vector<coord_t>& dist
) {
    int total = num_sheets + 2;

    dist.assign(total, 1e18);

    dist[source_node] = 0;

    vector<Win> stk;

    // Initial windows.
    for (auto& st : starts) {
        Win w;

        w.face = st.first;
        w.entry = -1;
        w.depth = 0;
        w.apex = st.second;
        w.lo = w.hi = Point(0, 0);
        w.c2 = {canon[0], canon[1], canon[2], canon[3]};

        stk.push_back(w);
    }

    while (!stk.empty()) {
        Win w = stk.back();

        stk.pop_back();

        int f = w.face;

        // Record visible face corners.
        for (int k = 0; k < 4; k++) {
            Point pt = w.c2[k];

            if (w.entry == -1 || in_sector(w.apex, w.lo, w.hi, pt)) {
                coord_t dd = (pt - w.apex).norm();

                int node = sheet_of[f * 4 + k];

                dist[node] = min(dist[node], dd);
            }
        }

        // Record source point if visible.
        if (face_of_src == f) {
            Point sp = to_frame(src_pt, f, w.c2);

            if (w.entry == -1 || in_sector(w.apex, w.lo, w.hi, sp)) {
                dist[num_sheets] = min(dist[num_sheets], (sp - w.apex).norm());
            }
        }

        // Record destination point if visible.
        if (face_of_dst == f) {
            Point dp = to_frame(dst_pt, f, w.c2);

            if (w.entry == -1 || in_sector(w.apex, w.lo, w.hi, dp)) {
                dist[num_sheets + 1] =
                    min(dist[num_sheets + 1], (dp - w.apex).norm());
            }
        }

        // Try crossing every edge except the entry edge.
        for (int e = 0; e < 4; e++) {
            if (e == w.entry) continue;

            Point c = w.c2[e];
            Point d = w.c2[(e + 1) % 4];

            Point q0, q1;

            // First face has full visibility.
            if (w.entry == -1) {
                if ((c - w.apex).norm() < 1e-9 ||
                    (d - w.apex).norm() < 1e-9) {
                    continue;
                }

                q0 = c;
                q1 = d;
            }
            // Otherwise clip edge by current sector.
            else if (!clip_sector(w.apex, w.lo, w.hi, c, d, q0, q1)) {
                continue;
            }

            // Prune if even the closest point is already beyond known answer.
            if (point_seg_dist(w.apex, q0, q1) > ub + 1e-6) {
                continue;
            }

            // Safety limit.
            if (w.depth + 1 > 256) {
                continue;
            }

            // Neighboring face across this edge.
            int g = nb_face[f][e];

            int eg = nb_entry[f][e];

            V3 p0 = face_corner[f][e];
            V3 p1 = face_corner[f][(e + 1) % 4];

            Point m0 = c;
            Point m1 = d;

            array<Point, 4> c2g;

            // Place all four corners of the neighboring face in the unfolded plane.
            for (int k = 0; k < 4; k++) {
                V3 cc = face_corner[g][k];

                if (cc == p0) {
                    c2g[k] = m0;
                } else if (cc == p1) {
                    c2g[k] = m1;
                } else {
                    c2g[k] = place_corner(p0, p1, cc, m0, m1, w.apex);
                }
            }

            Win nw;

            nw.face = g;
            nw.entry = eg;
            nw.depth = w.depth + 1;
            nw.apex = w.apex;
            nw.lo = q0;
            nw.hi = q1;
            nw.c2 = c2g;

            stk.push_back(nw);
        }
    }
}

// Compute coarse upper bound using face-center graph.
void compute_ub() {
    int nf = face_center.size();

    int total = nf + 2;

    vector<vector<pair<int, coord_t>>> g(total);

    // Adjacent faces have approximate cost 100.
    for (int f = 0; f < nf; f++) {
        for (int e = 0; e < 4; e++) {
            g[f].push_back({nb_face[f][e], 100.0});
        }
    }

    auto dist3 = [&](V3 a, V3 b) {
        V3 dd = sub3(a, b);

        return (coord_t)sqrt((coord_t)dot3(dd, dd));
    };

    // Connect source and destination to their containing face centers.
    g[nf].push_back({face_of_src, dist3(src_pt, face_center[face_of_src])});
    g[face_of_src].push_back({nf, dist3(src_pt, face_center[face_of_src])});

    g[nf + 1].push_back({face_of_dst, dist3(dst_pt, face_center[face_of_dst])});
    g[face_of_dst].push_back({nf + 1, dist3(dst_pt, face_center[face_of_dst])});

    vector<coord_t> dist(total, 1e18);

    priority_queue<
        pair<coord_t, int>,
        vector<pair<coord_t, int>>,
        greater<>
    > pq;

    dist[nf] = 0;

    pq.push({0, nf});

    while (!pq.empty()) {
        auto [c, u] = pq.top();

        pq.pop();

        if (c > dist[u] + 1e-9) continue;

        for (auto& [v, w] : g[u]) {
            if (dist[v] > dist[u] + w + 1e-9) {
                dist[v] = dist[u] + w;

                pq.push({dist[v], v});
            }
        }
    }

    ub = dist[nf + 1];
}

// Main solving routine.
void solve() {
    build_faces();

    build_adjacency();

    build_sheets();

    face_of_src = find_face_of(src_pt);

    face_of_dst = find_face_of(dst_pt);

    compute_ub();

    int nf = face_center.size();

    // For each vertex sheet, list all incident face-corners.
    vector<vector<pair<int, int>>> incident(num_sheets);

    for (int f = 0; f < nf; f++) {
        for (int k = 0; k < 4; k++) {
            incident[sheet_of[f * 4 + k]].push_back({f, k});
        }
    }

    int total = num_sheets + 2;

    int src_node = num_sheets;

    int dst_node = num_sheets + 1;

    // field[u][v] stores direct unfolded geodesic distance from u to v.
    vector<vector<coord_t>> field(total);

    // Run field sweep from every vertex sheet.
    for (int v = 0; v < num_sheets; v++) {
        vector<pair<int, Point>> starts;

        for (auto& [f, k] : incident[v]) {
            starts.push_back({f, canon[k]});
        }

        run_field(starts, v, field[v]);
    }

    // Run field sweep from source point.
    {
        array<Point, 4> c2 = {
            canon[0],
            canon[1],
            canon[2],
            canon[3]
        };

        Point sp = to_frame(src_pt, face_of_src, c2);

        run_field({{face_of_src, sp}}, src_node, field[src_node]);
    }

    // Dijkstra over graph nodes.
    vector<coord_t> dist(total, 1e18);

    priority_queue<
        pair<coord_t, int>,
        vector<pair<coord_t, int>>,
        greater<>
    > pq;

    dist[src_node] = 0;

    pq.push({0, src_node});

    while (!pq.empty()) {
        auto [c, u] = pq.top();

        pq.pop();

        if (c > dist[u] + 1e-12) continue;

        if (u == dst_node) continue;

        for (int v = 0; v < total; v++) {
            coord_t w = field[u][v];

            if (w > 1e17) continue;

            if (dist[v] > dist[u] + w + 1e-12) {
                dist[v] = dist[u] + w;

                pq.push({dist[v], v});
            }
        }
    }

    printf("%.10f\n", dist[dst_node]);
}

int main() {
    ios_base::sync_with_stdio(false);

    cin.tie(nullptr);

    read();

    solve();

    return 0;
}
```

---

## 4. Python solution

```python
import sys
import math
import heapq


EPS = 1e-9
INF = 10**100


# ------------------------------------------------------------
# 2D point class used for unfolded faces.
# ------------------------------------------------------------
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, value):
        # Multiplication by scalar.
        return Point(self.x * value, self.y * value)

    def __truediv__(self, value):
        return Point(self.x / value, self.y / value)

    def norm2(self):
        return self.x * self.x + self.y * self.y

    def norm(self):
        return math.sqrt(self.norm2())

    def perp(self):
        return Point(-self.y, self.x)


def dot2(a, b):
    return a.x * b.x + a.y * b.y


def cross2(a, b):
    return a.x * b.y - a.y * b.x


def sub3(a, b):
    return (a[0] - b[0], a[1] - b[1], a[2] - b[2])


def dot3(a, b):
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]


# Six axis directions.
DVEC = [
    (1, 0, 0),
    (-1, 0, 0),
    (0, 1, 0),
    (0, -1, 0),
    (0, 0, 1),
    (0, 0, -1),
]


def opp(d):
    return d ^ 1


# Canonical 2D square.
CANON = [
    Point(0, 0),
    Point(100, 0),
    Point(100, 100),
    Point(0, 100),
]


class DSU:
    def __init__(self, n):
        self.parent = list(range(n))

    def find(self, x):
        while self.parent[x] != x:
            self.parent[x] = self.parent[self.parent[x]]
            x = self.parent[x]
        return x

    def union(self, a, b):
        ra = self.find(a)
        rb = self.find(b)
        if ra != rb:
            self.parent[ra] = rb


class Solver:
    def __init__(self, centers, src_pt, dst_pt):
        self.centers = centers
        self.src_pt = src_pt
        self.dst_pt = dst_pt

        self.grid_set = set()

        self.face_center = []
        self.face_dir = []
        self.face_grid = []
        self.face_corner = []
        self.face_id = {}

        self.nb_face = []
        self.nb_entry = []

        self.sheet_of = []
        self.num_sheets = 0

        self.face_of_src = -1
        self.face_of_dst = -1

        self.ub = INF

    def present(self, gx, gy, gz):
        return (gx, gy, gz) in self.grid_set

    # ------------------------------------------------------------
    # Build all exposed faces.
    # ------------------------------------------------------------
    def build_faces(self):
        for c in self.centers:
            self.grid_set.add((c[0] // 100, c[1] // 100, c[2] // 100))

        for c in self.centers:
            gx = c[0] // 100
            gy = c[1] // 100
            gz = c[2] // 100

            for d in range(6):
                nx = gx + DVEC[d][0]
                ny = gy + DVEC[d][1]
                nz = gz + DVEC[d][2]

                # Hidden face.
                if self.present(nx, ny, nz):
                    continue

                # Face center.
                fc = (
                    c[0] + 50 * DVEC[d][0],
                    c[1] + 50 * DVEC[d][1],
                    c[2] + 50 * DVEC[d][2],
                )

                axis = d // 2

                # The two axes lying in this face.
                axes = [a for a in range(3) if a != axis]
                la, lb = axes[0], axes[1]

                signs = [(-1, -1), (1, -1), (1, 1), (-1, 1)]

                corners = []

                for sa, sb in signs:
                    p = [fc[0], fc[1], fc[2]]
                    p[la] += 50 * sa
                    p[lb] += 50 * sb
                    corners.append(tuple(p))

                fid = len(self.face_center)

                self.face_id[(gx, gy, gz, d)] = fid
                self.face_center.append(fc)
                self.face_dir.append(d)
                self.face_grid.append((gx, gy, gz))
                self.face_corner.append(corners)

    # ------------------------------------------------------------
    # Build surface adjacency across edges.
    # ------------------------------------------------------------
    def build_adjacency(self):
        nf = len(self.face_center)

        self.nb_face = [[-1] * 4 for _ in range(nf)]
        self.nb_entry = [[-1] * 4 for _ in range(nf)]

        for f in range(nf):
            d = self.face_dir[f]
            g = self.face_grid[f]

            for e in range(4):
                a = self.face_corner[f][e]
                b = self.face_corner[f][(e + 1) % 4]

                mid = (
                    (a[0] + b[0]) // 2,
                    (a[1] + b[1]) // 2,
                    (a[2] + b[2]) // 2,
                )

                delta = sub3(mid, self.face_center[f])

                t = -1

                for dd in range(6):
                    if delta == (
                        50 * DVEC[dd][0],
                        50 * DVEC[dd][1],
                        50 * DVEC[dd][2],
                    ):
                        t = dd
                        break

                ax = g[0] + DVEC[t][0]
                ay = g[1] + DVEC[t][1]
                az = g[2] + DVEC[t][2]

                cx = ax + DVEC[d][0]
                cy = ay + DVEC[d][1]
                cz = az + DVEC[d][2]

                if not self.present(ax, ay, az):
                    key = (g[0], g[1], g[2], t)
                elif not self.present(cx, cy, cz):
                    key = (ax, ay, az, d)
                else:
                    key = (cx, cy, cz, opp(t))

                gface = self.face_id[key]

                self.nb_face[f][e] = gface

                found = -1

                for k in range(4):
                    p = self.face_corner[gface][k]
                    q = self.face_corner[gface][(k + 1) % 4]

                    if (p == a and q == b) or (p == b and q == a):
                        found = k
                        break

                self.nb_entry[f][e] = found

    def corner_in(self, face, vertex):
        for k in range(4):
            if self.face_corner[face][k] == vertex:
                return k
        return -1

    # ------------------------------------------------------------
    # Build vertex sheets.
    # ------------------------------------------------------------
    def build_sheets(self):
        nf = len(self.face_center)

        dsu = DSU(nf * 4)

        for f in range(nf):
            for k in range(4):
                vertex = self.face_corner[f][k]

                # Two incident edges at this face corner.
                for e in (k, (k + 3) % 4):
                    g = self.nb_face[f][e]
                    kk = self.corner_in(g, vertex)

                    dsu.union(f * 4 + k, g * 4 + kk)

        remap = {}
        self.sheet_of = [-1] * (nf * 4)
        self.num_sheets = 0

        for i in range(nf * 4):
            r = dsu.find(i)

            if r not in remap:
                remap[r] = self.num_sheets
                self.num_sheets += 1

            self.sheet_of[i] = remap[r]

    # ------------------------------------------------------------
    # Find exposed face containing a given point.
    # ------------------------------------------------------------
    def find_face_of(self, p):
        for f, fc in enumerate(self.face_center):
            d = self.face_dir[f]
            axis = d // 2

            if p[axis] != fc[axis]:
                continue

            ok = True

            for a in range(3):
                if a == axis:
                    continue

                # Input guarantees the point is not on an edge.
                if abs(p[a] - fc[a]) >= 50:
                    ok = False
                    break

            if ok:
                return f

        return -1

    # ------------------------------------------------------------
    # Convert a 3D point on a face into the current unfolded 2D frame.
    # ------------------------------------------------------------
    def to_frame(self, p, f, c2):
        o = self.face_corner[f][0]

        edge01 = sub3(self.face_corner[f][1], o)
        edge03 = sub3(self.face_corner[f][3], o)

        po = sub3(p, o)

        alpha = dot3(po, edge01) / 10000.0
        beta = dot3(po, edge03) / 10000.0

        return c2[0] + (c2[1] - c2[0]) * alpha + (c2[3] - c2[0]) * beta

    def in_sector(self, apex, lo, hi, pt):
        dl = lo - apex
        dr = hi - apex

        if cross2(dl, dr) < 0:
            dl, dr = dr, dl

        v = pt - apex

        tol = 1e-9 * v.norm() + 1e-12

        tl = tol * dl.norm()
        tr = tol * dr.norm()

        return cross2(dl, v) >= -tl and cross2(dr, v) <= tr

    def clip_half(self, a, b, s0, s1):
        if abs(b) < 1e-15:
            if a < -1e-12:
                return 1.0, 0.0
            return s0, s1

        r = -a / b

        if b > 0:
            s0 = max(s0, r)
        else:
            s1 = min(s1, r)

        return s0, s1

    # ------------------------------------------------------------
    # Clip segment cd to angular sector.
    # ------------------------------------------------------------
    def clip_sector(self, apex, lo, hi, c, d):
        dl = lo - apex
        dr = hi - apex

        if cross2(dl, dr) < 0:
            dl, dr = dr, dl

        s0 = 0.0
        s1 = 1.0

        dc = d - c

        s0, s1 = self.clip_half(cross2(dl, c - apex), cross2(dl, dc), s0, s1)
        s0, s1 = self.clip_half(-cross2(dr, c - apex), -cross2(dr, dc), s0, s1)

        if s0 > s1 - 1e-12:
            return None

        q0 = c + dc * s0
        q1 = c + dc * s1

        return q0, q1

    def point_seg_dist(self, p, a, b):
        ab = b - a

        length2 = dot2(ab, ab)

        if length2 < 1e-18:
            return (p - a).norm()

        t = dot2(p - a, ab) / length2
        t = max(0.0, min(1.0, t))

        proj = a + ab * t

        return (p - proj).norm()

    # ------------------------------------------------------------
    # Place a far corner of the neighboring face in the unfolded plane.
    # ------------------------------------------------------------
    def place_corner(self, p0, p1, q, m0, m1, apex):
        along = dot3(sub3(q, p0), sub3(p1, p0)) / 100.0

        dsq = dot3(sub3(q, p0), sub3(q, p0))

        h = math.sqrt(max(0.0, dsq - along * along))

        u = (m1 - m0) * (1.0 / 100.0)

        base = m0 + u * along

        cand1 = base + u.perp() * h
        cand2 = base - u.perp() * h

        side = cross2(m1 - m0, apex - m0)
        s1 = cross2(m1 - m0, cand1 - m0)

        if (side >= 0) == (s1 >= 0):
            return cand2

        return cand1

    # ------------------------------------------------------------
    # Coarse upper bound using face-center graph.
    # ------------------------------------------------------------
    def compute_ub(self):
        nf = len(self.face_center)

        total = nf + 2

        graph = [[] for _ in range(total)]

        for f in range(nf):
            for e in range(4):
                graph[f].append((self.nb_face[f][e], 100.0))

        def dist3(a, b):
            d = sub3(a, b)
            return math.sqrt(dot3(d, d))

        src_node = nf
        dst_node = nf + 1

        ws = dist3(self.src_pt, self.face_center[self.face_of_src])
        wd = dist3(self.dst_pt, self.face_center[self.face_of_dst])

        graph[src_node].append((self.face_of_src, ws))
        graph[self.face_of_src].append((src_node, ws))

        graph[dst_node].append((self.face_of_dst, wd))
        graph[self.face_of_dst].append((dst_node, wd))

        dist = [INF] * total
        dist[src_node] = 0.0

        pq = [(0.0, src_node)]

        while pq:
            cur, u = heapq.heappop(pq)

            if cur > dist[u] + 1e-9:
                continue

            for v, w in graph[u]:
                nd = cur + w

                if nd + 1e-9 < dist[v]:
                    dist[v] = nd
                    heapq.heappush(pq, (nd, v))

        self.ub = dist[dst_node]

    # ------------------------------------------------------------
    # Run one unfolding sweep from a source graph node.
    # ------------------------------------------------------------
    def run_field(self, starts, source_node):
        total = self.num_sheets + 2

        dist = [INF] * total
        dist[source_node] = 0.0

        stack = []

        for face, apex in starts:
            c2 = [Point(p.x, p.y) for p in CANON]

            stack.append((face, -1, 0, apex, Point(0, 0), Point(0, 0), c2))

        while stack:
            face, entry, depth, apex, lo, hi, c2 = stack.pop()

            # Record visible corners.
            for k in range(4):
                pt = c2[k]

                if entry == -1 or self.in_sector(apex, lo, hi, pt):
                    node = self.sheet_of[face * 4 + k]
                    dd = (pt - apex).norm()

                    if dd < dist[node]:
                        dist[node] = dd

            # Record source point if visible.
            if self.face_of_src == face:
                sp = self.to_frame(self.src_pt, face, c2)

                if entry == -1 or self.in_sector(apex, lo, hi, sp):
                    dd = (sp - apex).norm()

                    if dd < dist[self.num_sheets]:
                        dist[self.num_sheets] = dd

            # Record destination point if visible.
            if self.face_of_dst == face:
                dp = self.to_frame(self.dst_pt, face, c2)

                if entry == -1 or self.in_sector(apex, lo, hi, dp):
                    dd = (dp - apex).norm()

                    if dd < dist[self.num_sheets + 1]:
                        dist[self.num_sheets + 1] = dd

            # Propagate through all edges except the entry edge.
            for e in range(4):
                if e == entry:
                    continue

                c = c2[e]
                d = c2[(e + 1) % 4]

                if entry == -1:
                    # Avoid degenerate windows starting exactly along an incident edge.
                    if (c - apex).norm() < 1e-9 or (d - apex).norm() < 1e-9:
                        continue

                    q0, q1 = c, d
                else:
                    clipped = self.clip_sector(apex, lo, hi, c, d)

                    if clipped is None:
                        continue

                    q0, q1 = clipped

                # If this window is already farther than known answer, discard.
                if self.point_seg_dist(apex, q0, q1) > self.ub + 1e-6:
                    continue

                if depth + 1 > 256:
                    continue

                g = self.nb_face[face][e]
                eg = self.nb_entry[face][e]

                p0 = self.face_corner[face][e]
                p1 = self.face_corner[face][(e + 1) % 4]

                m0 = c
                m1 = d

                c2g = [None] * 4

                # Place neighboring face corners.
                for k in range(4):
                    cc = self.face_corner[g][k]

                    if cc == p0:
                        c2g[k] = m0
                    elif cc == p1:
                        c2g[k] = m1
                    else:
                        c2g[k] = self.place_corner(p0, p1, cc, m0, m1, apex)

                stack.append((g, eg, depth + 1, apex, q0, q1, c2g))

        return dist

    # ------------------------------------------------------------
    # Full solve routine.
    # ------------------------------------------------------------
    def solve(self):
        self.build_faces()
        self.build_adjacency()
        self.build_sheets()

        self.face_of_src = self.find_face_of(self.src_pt)
        self.face_of_dst = self.find_face_of(self.dst_pt)

        self.compute_ub()

        nf = len(self.face_center)

        incident = [[] for _ in range(self.num_sheets)]

        for f in range(nf):
            for k in range(4):
                sheet = self.sheet_of[f * 4 + k]
                incident[sheet].append((f, k))

        total = self.num_sheets + 2

        src_node = self.num_sheets
        dst_node = self.num_sheets + 1

        field = [None] * total

        # Field sweeps from vertex sheets.
        for v in range(self.num_sheets):
            starts = []

            for f, k in incident[v]:
                starts.append((f, Point(CANON[k].x, CANON[k].y)))

            field[v] = self.run_field(starts, v)

        # Field sweep from source point.
        c2 = [Point(p.x, p.y) for p in CANON]
        sp = self.to_frame(self.src_pt, self.face_of_src, c2)

        field[src_node] = self.run_field([(self.face_of_src, sp)], src_node)

        # Dijkstra over graph nodes.
        dist = [INF] * total
        dist[src_node] = 0.0

        pq = [(0.0, src_node)]

        while pq:
            cur, u = heapq.heappop(pq)

            if cur > dist[u] + 1e-12:
                continue

            if u == dst_node:
                continue

            row = field[u]

            for v in range(total):
                w = row[v]

                if w > 1e90:
                    continue

                nd = cur + w

                if nd + 1e-12 < dist[v]:
                    dist[v] = nd
                    heapq.heappush(pq, (nd, v))

        return dist[dst_node]


def main():
    data = list(map(int, sys.stdin.read().split()))

    if not data:
        return

    idx = 0

    n = data[idx]
    idx += 1

    centers = []

    for _ in range(n):
        x, y, z = data[idx], data[idx + 1], data[idx + 2]
        idx += 3
        centers.append((x, y, z))

    sx, sy, sz, dx, dy, dz = data[idx:idx + 6]

    src_pt = (sx, sy, sz)
    dst_pt = (dx, dy, dz)

    solver = Solver(centers, src_pt, dst_pt)

    ans = solver.solve()

    print(f"{ans:.10f}")


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

---

## 5. Compressed editorial

Build the exposed square faces of the polycube. For each face edge, determine the adjacent surface face using the local cube configuration: convex fold, flat continuation, or concave fold.

A shortest path on a polyhedral surface is straight when the crossed faces are unfolded and may bend only at vertices. However, one geometric vertex can contain several disconnected surface sheets, so group face-corners into connected sheets using DSU. The graph nodes are all vertex sheets plus `S` and `D`.

For every node, run an unfolding sweep. A sweep window stores the current face, its unfolded 2D corners, the unfolded source point, and the angular interval of visibility. When crossing an edge, clip the visibility cone to that edge and unfold the neighboring face into the same plane. Whenever a visible corner, `S`, or `D` is encountered, record its Euclidean distance from the source. A coarse face-center Dijkstra gives an upper bound used to prune windows.

These sweeps produce direct geodesic distances between graph nodes. Finally, run Dijkstra from `S` to `D` on this graph and output the distance.