## 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. C++ Solution

```cpp
#include <bits/stdc++.h>
// #include <coding_library/geometry/point.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;

struct Point {
    static constexpr coord_t eps = 1e-9;
    static inline const coord_t PI = acos((coord_t)-1.0);

    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
        );
    }

    friend coord_t arc_area(
        const Point& center, coord_t r, const Point& p1, const Point& p2
    ) {
        coord_t theta1 = (p1 - center).angle();
        coord_t theta2 = (p2 - center).angle();
        if(theta2 < theta1 - eps) {
            theta2 += 2 * PI;
        }

        coord_t d_theta = theta2 - theta1;
        coord_t cx = center.x, cy = center.y;
        coord_t area = r * cx * (sin(theta2) - sin(theta1)) -
                       r * cy * (cos(theta2) - cos(theta1)) + r * r * d_theta;
        return area / 2.0;
    }

    friend vector<Point> intersect_circles(
        const Point& c1, coord_t r1, const Point& c2, coord_t r2
    ) {
        Point d = c2 - c1;
        coord_t dist = d.norm();

        if(dist > r1 + r2 + eps || dist < abs(r1 - r2) - eps || dist < eps) {
            return {};
        }

        coord_t a = (r1 * r1 - r2 * r2 + dist * dist) / (2 * dist);
        coord_t h_sq = r1 * r1 - a * a;
        if(h_sq < -eps) {
            return {};
        }
        if(h_sq < 0) {
            h_sq = 0;
        }
        coord_t h = sqrt(h_sq);

        Point mid = c1 + d.unit() * a;
        Point perp_dir = d.perp().unit();

        if(h < eps) {
            return {mid};
        }
        return {mid + perp_dir * h, mid - perp_dir * h};
    }

    friend optional<Point> intersect_ray_segment(
        const Point& ray_start, const Point& ray_through, const Point& seg_a,
        const Point& seg_b
    ) {
        Point ray_dir = ray_through - ray_start;
        if(ray_dir.norm2() < Point::eps) {
            return {};
        }
        Point seg_dir = seg_b - seg_a;
        coord_t denom = ray_dir ^ seg_dir;
        if(fabs(denom) < eps) {
            return {};
        }
        coord_t t = ((seg_a - ray_start) ^ seg_dir) / denom;
        if(t < eps) {
            return {};
        }
        coord_t s = ((seg_a - ray_start) ^ ray_dir) / denom;
        if(s < eps || s > 1 - eps) {
            return {};
        }
        return ray_start + ray_dir * t;
    }
};

using V3 = array<int64_t, 3>;

V3 sub3(V3 a, V3 b) { return {a[0] - b[0], a[1] - b[1], a[2] - b[2]}; }
int64_t dot3(V3 a, V3 b) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; }

const int dvec[6][3] = {{1, 0, 0},  {-1, 0, 0}, {0, 1, 0},
                        {0, -1, 0}, {0, 0, 1},  {0, 0, -1}};

int opp(int d) { return d ^ 1; }

int n;
vector<V3> centers;
V3 src_pt, dst_pt;

vector<V3> face_center;
vector<int> face_dir;
vector<array<int, 3>> face_grid;
vector<array<V3, 4>> face_corner;
map<array<int, 4>, int> face_id;
vector<array<int, 4>> nb_face;
vector<array<int, 4>> nb_entry;
vector<int> sheet_of;
int num_sheets;
int face_of_src, face_of_dst;

set<array<int, 3>> grid_set;

bool present(int gx, int gy, int gz) {
    return grid_set.count({gx, gy, gz}) > 0;
}

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];
}

void build_faces() {
    for(auto& c: centers) {
        grid_set.insert(
            {(int)(c[0] / 100), (int)(c[1] / 100), (int)(c[2] / 100)}
        );
    }

    for(auto& c: centers) {
        int gx = c[0] / 100, gy = c[1] / 100, gz = c[2] / 100;
        for(int d = 0; d < 6; d++) {
            if(present(gx + dvec[d][0], gy + dvec[d][1], gz + dvec[d][2])) {
                continue;
            }

            V3 fc = {
                c[0] + 50 * dvec[d][0], c[1] + 50 * dvec[d][1],
                c[2] + 50 * dvec[d][2]
            };
            int axis = d / 2;
            int la = -1, lb = -1;
            for(int a = 0; a < 3; a++) {
                if(a == axis) {
                    continue;
                }

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

            const int signs[4][2] = {{-1, -1}, {1, -1}, {1, 1}, {-1, 1}};
            array<V3, 4> 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;
            }

            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);
        }
    }
}

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++) {
            V3 a = face_corner[f][e];
            V3 b = face_corner[f][(e + 1) % 4];
            V3 mid = {(a[0] + b[0]) / 2, (a[1] + b[1]) / 2, (a[2] + b[2]) / 2};
            V3 delta = sub3(mid, face_center[f]);
            int t = -1;
            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;
                }
            }

            int ax = g[0] + dvec[t][0], ay = g[1] + dvec[t][1],
                az = g[2] + dvec[t][2];
            int cx = ax + dvec[d][0], cy = ay + dvec[d][1],
                cz = az + dvec[d][2];

            array<int, 4> key;
            if(!present(ax, ay, az)) {
                key = {g[0], g[1], g[2], t};
            } else if(!present(cx, cy, cz)) {
                key = {ax, ay, az, d};
            } else {
                key = {cx, cy, cz, opp(t)};
            }

            int gface = face_id[key];
            nb_face[f][e] = gface;

            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;
        }
    }
}

int corner_in(int g, V3 v) {
    for(int k = 0; k < 4; k++) {
        if(face_corner[g][k] == v) {
            return k;
        }
    }

    return -1;
}

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); };

    for(int f = 0; f < nf; f++) {
        for(int k = 0; k < 4; k++) {
            V3 v = face_corner[f][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];
    }
}

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;
        if(p[axis] != fc[axis]) {
            continue;
        }

        bool ok = true;
        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;
}

const Point canon[4] = {{0, 0}, {100, 0}, {100, 100}, {0, 100}};

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;
}

bool in_sector(Point apex, Point lo, Point hi, Point pt) {
    Point dl = lo - apex, 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;
}

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);
    }
}

bool clip_sector(
    Point apex, Point lo, Point hi, Point c, Point d, Point& q0, Point& q1
) {
    Point dl = lo - apex, dr = hi - apex;
    if((dl ^ dr) < 0) {
        swap(dl, dr);
    }

    coord_t s0 = 0, 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;
}

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();
}

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;
}

struct Win {
    int face, entry, depth;
    Point apex, lo, hi;
    array<Point, 4> c2;
};

coord_t ub;

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

        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);
            }
        }

        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());
            }
        }

        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());
            }
        }

        for(int e = 0; e < 4; e++) {
            if(e == w.entry) {
                continue;
            }

            Point c = w.c2[e], d = w.c2[(e + 1) % 4];
            Point q0, q1;
            if(w.entry == -1) {
                if((c - w.apex).norm() < 1e-9 || (d - w.apex).norm() < 1e-9) {
                    continue;
                }
                q0 = c;
                q1 = d;
            } else if(!clip_sector(w.apex, w.lo, w.hi, c, d, q0, q1)) {
                continue;
            }

            if(point_seg_dist(w.apex, q0, q1) > ub + 1e-6) {
                continue;
            }

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

            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, m1 = d;

            array<Point, 4> c2g;
            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);
        }
    }
}

void compute_ub() {
    int nf = face_center.size();
    int total = nf + 2;
    vector<vector<pair<int, coord_t>>> g(total);
    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));
    };

    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];
}

void solve() {
    // The robot lives on the surface of the polycube, which is a closed
    // manifold tiled by square faces of side 100 - exactly the cube faces that
    // are not glued to another cube. Two faces that share an edge are walkable
    // across that edge, and which face you arrive on is decided locally by the
    // three cubes around the edge. For a face on cube C with outward normal d,
    // step across the edge lying towards a tangent direction t and look at the
    // side cube A = C + t and the diagonal cube B = C + t + d:
    //
    //     - if A is missing you fold over C's own edge onto face (C, t), a
    //       convex 90 degree turn (rule a);
    //
    //     - else if B is missing you cross flat onto the coplanar face (A, d),
    //       a straight 180 degree step (rule b);
    //
    //     - else you turn into the concave corner onto face (B, -t), a 270
    //       degree step (rule c).
    //
    // Every surface edge falls into exactly one of these cases, so the gluing
    // is unambiguous; in particular two cubes that touch only along an edge are
    // never joined, because each of their faces folds back onto its own cube.
    //
    // A shortest route on such a surface is a geodesic: unfold the strip of
    // faces it crosses into the plane and it becomes a straight segment, and it
    // can change direction only at a vertex of the polycube (a cube corner). We
    // therefore turn the problem into a graph whose nodes are the surface
    // vertices together with the start S and the goal D, and whose edges carry
    // the length of the straight geodesic that joins two nodes while passing
    // through face interiors only. Dijkstra over this graph reconstructs any
    // geodesic by chaining these straight pieces between consecutive bends.
    //
    // A vertex is not always a single node. Around a corner the incident faces
    // may split into several sheets that meet only at the point itself - the
    // touching tips of two cubes glued along an edge are the simplest example -
    // and the robot may walk within one sheet but cannot cross between sheets
    // through the bare point. We group the incident face-corners into connected
    // sheets using the edge gluing above and make each sheet its own node,
    // which stops such a pinch from leaking a free shortcut across the gap.
    //
    // The edge lengths come from a single-source sweep that unfolds the surface
    // on the fly. Starting from a point (or from a vertex sheet, which radiates
    // into each of its incident faces) we carry the unfolded image of the
    // source in one common plane and walk outward face by face, keeping for the
    // current face the visible angular cone seen from that image. Crossing into
    // a neighbour we place its far corners from their known distances to the
    // shared edge, then clip the cone against each onward edge; whenever the
    // cone covers a vertex or the target we record the Euclidean distance to
    // its unfolded image. The cone only narrows as it advances, so a window
    // whose nearest point already lies past an upper bound on the answer -
    // taken from a coarse path through the face centres - can never improve any
    // node and is dropped, which is what keeps the otherwise unbounded sweep
    // finite.

    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();
    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, dst_node = num_sheets + 1;

    vector<vector<coord_t>> field(total);
    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]);
    }

    {
        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]);
    }

    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);

    int T = 1;
    for(int test = 1; test <= T; test++) {
        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.
