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

452. Colony Maintenance
Time limit per test: 0.75 second(s)
Memory limit: 262144 kilobytes
input: standard
output: standard

It is the year of 2xxx. Human beings have migrated out of the Earth and resided in space colonies. Because of the severe environment of the space, which comes from cosmic rays and meteorites, those colonies have to be continuously repaired by maintenance robots. In this problem, you are requested to write a program that calculates the shortest distance for a robot to move from the present point to the next repair point on a colony. Each colony can be modeled as a polycube, that is, one or more cubes of the same size all joined by their faces. Note that, as forming a polycube, all cubes of each colony are connected. This implies every cube has at least one face coincident with a face of another cube, except for colonies with a single cube. The figure below illustrates a couple of example colonies.


The maintenance robot can move only on the surface of the polycubes, that is, on faces not in common with other cubes. In addition, due to the structure of the colonies, move of the robot beyond a face is restricted to the following cases: (a) the robot is moving between adjacent faces of the same cube; (b) the robot is moving between adjacent faces of adjacent cubes; and (c) the robot is moving along the inner side of an L-shape, namely, the robot is moving between adjacent faces of two cubes that have the common adjacent cube. Here, adjacent faces denote those with the common edge, and adjacent cubes denote those with the common face.


For the purpose of this problem, we consider an xyz-space (i.e. a three-dimensional space with the Cartesian coordinate system) where all cubes are placed in such a way that each edge is parallel to x-axis, y-axis, or z-axis. Also, the edge length of the cubes is scaled to 100.
Input
The input has the following format:
n
x1 y1 z1
...
xn yn zn
sx sy sz dx dy dz
n is the number of cubes of a colony (1 ≤ n ≤ 16). (xi,yi,zi) represents the coordinates of the center of the i-th cube, where each coordinate value is guaranteed to be a multiple of 100 (i.e. the edge length of the cubes). (sx,sy,sz) and (dx,dy,dz) represent the robot's present point and the next repair point respectively. You may assume that these two points are always different and that neither of them lies on any edge. All coordinate values are integers between -2000 and 2000 inclusive.
Output
For each dataset, print the distance of the shortest route from the present point to the next repair point. The distance may be printed with any number of digits after the decimal point, provided the absolute error does not exceed 10-8.
Example(s)
sample input
sample output
1
0 0 0
0 0 50 30 40 50
50.0

sample input
sample output
1
0 0 0
50 0 0 0 50 0
100.0

sample input
sample output
2
0 0 0
100 0 0
0 0 50 100 0 50
100.0

sample input
sample output
3
0 0 0
100 0 0
0 100 0
100 50 0 50 100 0
100.0

sample input
sample output
7
0 0 0
100 0 0
-100 0 0
0 100 0
0 -100 0
0 0 100
0 0 -100
150 0 0 -150 0 0
416.2277660168

sample input
sample output
6
0 0 0
100 0 0
0 100 0
100 100 0
0 0 100
100 100 100
100 0 50 0 100 50
141.4213562373

sample input
sample output
6
0 0 0
100 0 0
0 100 0
100 100 0
0 0 100
100 100 100
100 100 150 0 0 150
316.2277660168

sample input
sample output
6
0 0 0
100 0 0
0 100 0
100 100 100
100 0 100
0 100 100
0 0 50 100 100 50
341.4213562373

sample input
sample output
8
0 0 0
-100 0 0
-100 100 0
-100 200 0
0 200 0
100 200 0
100 200 100
100 100 100
0 0 50 100 100 50
341.4213562373

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

We are given a connected polycube consisting of at most `16` axis-aligned cubes of side length `100`.

A robot can move only on exposed cube faces. It may cross an exposed face edge according to the physical surface rules: convex folds, flat continuations, and concave L-shaped inner folds.

Given two surface points `S` and `D`, neither lying on an edge, compute the shortest walking distance on the polycube surface.

---

## 2. Key observations

### Observation 1: Surface faces are square patches

Each cube has 6 faces. A face is walkable iff there is no cube directly adjacent through that face.

Represent every exposed face as `(cube_grid_position, outward_direction)` where outward direction is one of `+x, -x, +y, -y, +z, -z`.

### Observation 2: Face adjacency is local

Suppose we are on face `(C, d)` and cross an edge in tangent direction `t`. Let `A = C + t` and `B = C + t + d`. Then:

1. If `A` is missing, fold around the current cube edge: `(C, d) -> (C, t)`.
2. Else if `B` is missing, continue flat onto the adjacent cube: `(C, d) -> (A, d)`.
3. Else both `A` and `B` exist, move along the concave inner L-shape: `(C, d) -> (B, -t)`.

### Observation 3: Shortest surface paths become straight after unfolding

On a polyhedral surface, if we unfold all crossed faces into one plane, a shortest path segment that does not pass through a vertex is a straight line. Therefore, a shortest path can be decomposed into straight unfolded geodesic pieces whose endpoints are among `S`, `D`, and surface vertices.

We build a graph where nodes are `S`, `D`, and surface vertex sheets, and edge weights are shortest straight unfolded geodesic distances. The answer is the shortest path in this graph.

### Observation 4: One geometric vertex may contain multiple sheets

Several face corners can have the same 3D coordinate but may not be connected by walkable surface around that point. The robot cannot teleport through a zero-width touching point. We group face-corners into connected components through surface edges — each component is a **vertex sheet** that becomes one graph node.

### Observation 5: Compute direct geodesic distances by unfolding sweep

For each graph node as source, propagate over faces using unfolding windows. A window stores the current face, entry edge, unfolded source point, visible angular sector, and 2D coordinates of current face corners. When crossing an edge: clip visibility sector, unfold the neighboring face, continue propagation. Whenever a visible graph node is encountered, record its Euclidean distance. Then run Dijkstra on the graph of these direct distances.

---

## 3. Full solution approach

### Step 1: Build exposed faces

For every cube and every direction: if no neighboring cube exists there, store the face as an exposed surface face with its cube grid coordinate, outward direction, 3D face center, and four 3D corners.

### Step 2: Build face adjacency

For each exposed face and each of its four edges, apply the three local rules above to determine which exposed face lies across that edge and which edge index of the neighbor corresponds.

### Step 3: Build vertex sheets using DSU

For every face-corner `(f, k)`, union it with the face-corners at the same 3D position reachable through the two incident surface edges. The DSU components are the vertex sheets.

### Step 4: Locate faces containing S and D

Scan all exposed faces and find which one contains each point.

### Step 5: Compute a coarse upper bound

Run Dijkstra on a graph where each face connects to its four neighbors with cost 100, and S/D connect to their containing faces. This gives a coarse upper bound on the true shortest path.

### Step 6: Run unfolding field sweeps

For each vertex sheet and for S, launch a field sweep:
- From a vertex sheet, seed one window per incident face-corner.
- From S, seed one window in the face containing S.

Each window propagates outward, clipping its angular cone to each new edge and unfolding the neighboring face into the same 2D plane. Record distances whenever visible corners, S, or D fall within the current cone. Prune windows whose nearest edge point exceeds the upper bound.

### Step 7: Dijkstra over the graph nodes

Use the collected direct geodesic distances between all pairs of graph nodes. Run Dijkstra from S to D and output the distance.

---

## 4. C++ implementation

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

---

## 5. Python implementation (detailed comments)

```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()
```
