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


Note
The figure below shows how the robot will move in the last example:

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

```text
(cube_grid_position, outward_direction)
```

where `outward_direction` is one of:

```text
+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:

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

Then:

1. If `A` is missing, we fold around the current cube edge:
   ```text
   (C, d) -> (C, t)
   ```

2. Else if `B` is missing, we continue flat onto the adjacent cube:
   ```text
   (C, d) -> (A, d)
   ```

3. Else both `A` and `B` exist, so we move along the concave inner L-shape:
   ```text
   (C, d) -> (B, -t)
   ```

This gives the neighboring surface face across every surface edge.

---

### 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:

```text
S, D, surface vertices
```

So we build a graph:

```text
nodes = S + D + surface vertex sheets
edge weight = shortest straight unfolded path between two nodes
```

Then the answer is 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.

So instead of using one node per 3D vertex, we group face-corners into connected components through surface edges. Each component is a **vertex sheet**.

---

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

For each graph node as source, we propagate over faces using unfolding windows.

A window stores:

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

When crossing an edge:

1. Clip visibility sector to the crossed edge.
2. Unfold the neighboring face into the same plane.
3. Continue propagation.

Whenever a visible graph node is encountered, record its Euclidean distance in the unfolded plane.

Finally 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 a neighboring cube exists there, the face is internal.
- Otherwise, store the face as an exposed surface face.

For each exposed face, store:

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

---

### Step 2: Build face adjacency

For each exposed face and each of its four edges:

- Determine the tangent direction `t` from the face center to the edge midpoint.
- Use the local rule with cubes `A = C + t`, `B = C + t + d`.
- Find the neighboring exposed face.
- Also record which edge of the neighboring face is the shared edge.

---

### Step 3: Build vertex sheets

Create one DSU element for every face-corner.

For each face-corner, union it with the corresponding corner of the neighboring face across its two incident edges.

After all unions, each DSU component is one valid surface vertex sheet.

---

### Step 4: Locate source and destination faces

For each point `S` and `D`, find the exposed face whose plane contains the point and where the point lies strictly inside the square.

The statement guarantees neither point lies on an edge.

---

### Step 5: Compute a pruning upper bound

Run a rough Dijkstra on the face graph:

- neighboring faces connected with cost `100`,
- source and destination connected to their containing face centers.

This gives an upper bound on the answer.

During unfolding, if a window is already farther than this upper bound, it cannot improve the final answer and is discarded.

---

### Step 6: Run unfolding sweep from every graph node

For every vertex sheet and for `S`:

- Initialize one or more starting windows.
- Propagate across faces.
- Record direct unfolded distances to visible vertex sheets, `S`, and `D`.

This fills a matrix:

```text
field[u][v] = direct geodesic distance from u to v
```

We do not need to run a sweep from `D`, because Dijkstra stops when reaching `D`.

---

### Step 7: Dijkstra on graph nodes

Graph nodes are:

```text
0 ... num_sheets - 1 : vertex sheets
num_sheets           : S
num_sheets + 1       : D
```

Run Dijkstra from `S` using `field[u][v]` as edge weights.

The distance to `D` is the answer.

---

## 4. C++ implementation with detailed comments

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

using coord_t = double;
const coord_t INF = 1e100;

// ------------------------------------------------------------
// 2D point used for unfolded surface geometry.
// ------------------------------------------------------------
struct Point {
    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);
    }

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

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

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

    coord_t norm() const {
        return sqrt(norm2());
    }

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

// ------------------------------------------------------------
// 3D integer vector.
// All coordinates are integral.
// ------------------------------------------------------------
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];
}

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

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

// ------------------------------------------------------------
// Global data.
// ------------------------------------------------------------
int n;
vector<V3> centers;
V3 src_pt, dst_pt;

// Surface face data.
vector<V3> face_center;
vector<int> face_dir;
vector<array<int, 3>> face_grid;
vector<array<V3, 4>> face_corner;

// Maps (cube grid x, y, z, direction) to face id.
map<array<int, 4>, int> face_id;

// Face adjacency.
// nb_face[f][e] is face adjacent to face f through edge e.
// nb_entry[f][e] is the corresponding edge index in that adjacent face.
vector<array<int, 4>> nb_face;
vector<array<int, 4>> nb_entry;

// For each face-corner, stores its vertex sheet id.
vector<int> sheet_of;
int num_sheets;

// Containing faces of source and destination.
int face_of_src, face_of_dst;

// Occupied cube grid cells.
set<array<int, 3>> grid_set;

coord_t ub;

// Canonical 2D square for a newly unfolded face.
const Point canon[4] = {
    {0, 0},
    {100, 0},
    {100, 100},
    {0, 100}
};

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

// ------------------------------------------------------------
// Input.
// ------------------------------------------------------------
void read_input() {
    cin >> n;

    centers.resize(n);

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

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

// ------------------------------------------------------------
// Build all exposed faces.
// ------------------------------------------------------------
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;
        int gy = c[1] / 100;
        int gz = c[2] / 100;

        for (int d = 0; d < 6; d++) {
            // Internal face.
            if (present(gx + dvec[d][0],
                        gy + dvec[d][1],
                        gz + dvec[d][2])) {
                continue;
            }

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

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

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

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

                if (la == -1) 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;
            }

            int id = face_center.size();

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

// ------------------------------------------------------------
// Build adjacency across each surface edge.
// ------------------------------------------------------------
void build_adjacency() {
    int nf = face_center.size();

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

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

        for (int e = 0; e < 4; e++) {
            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
            };

            // Direction from face center to edge midpoint.
            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;
                }
            }

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

            // B = A + outward direction.
            int bx = ax + dvec[d][0];
            int by = ay + dvec[d][1];
            int bz = az + dvec[d][2];

            array<int, 4> key;

            if (!present(ax, ay, az)) {
                // Convex fold around current cube edge.
                key = {g[0], g[1], g[2], t};
            } else if (!present(bx, by, bz)) {
                // Flat continuation onto adjacent cube.
                key = {ax, ay, az, d};
            } else {
                // Concave inner L-shaped fold.
                key = {bx, by, bz, opp(t)};
            }

            int nf2 = face_id[key];

            nb_face[f][e] = nf2;

            // Find which edge in nf2 is the same geometric edge.
            int found = -1;

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

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

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

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

    return -1;
}

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

    vector<int> parent(nf * 4);

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

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

        return x;
    };

    auto unite = [&](int a, int b) {
        parent[find(a)] = find(b);
    };

    for (int f = 0; f < nf; f++) {
        for (int k = 0; k < 4; k++) {
            V3 v = face_corner[f][k];

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

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

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

    map<int, int> remap;
    num_sheets = 0;

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

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

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

// ------------------------------------------------------------
// Find the exposed face containing a given point.
// The point is guaranteed not to lie on an edge.
// ------------------------------------------------------------
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;

        // Must be in the plane of this face.
        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;
}

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

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

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

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

// ------------------------------------------------------------
// Check whether point pt is inside angular sector apex-lo-hi.
// ------------------------------------------------------------
bool in_sector(Point apex, Point lo, Point hi, Point pt) {
    Point dl = lo - apex;
    Point dr = hi - apex;

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

    Point v = pt - apex;

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

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

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

        return;
    }

    coord_t r = -a / b;

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

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

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

    coord_t s0 = 0;
    coord_t s1 = 1;

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

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

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

    return true;
}

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

    coord_t len = ab * ab;

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

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

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

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

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

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

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

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

    Point base = m0 + u * along;

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

    // The new face must lie on the opposite side of the shared edge
    // from the current unfolded source.
    coord_t side = (m1 - m0) ^ (apex - m0);
    coord_t s1 = (m1 - m0) ^ (cand1 - m0);

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

    return cand1;
}

// ------------------------------------------------------------
// Window used by unfolding sweep.
// ------------------------------------------------------------
struct Window {
    int face;
    int entry;
    int depth;

    Point apex;
    Point lo;
    Point hi;

    array<Point, 4> c2;
};

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

    int total = nf + 2;

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

    for (int f = 0; f < nf; f++) {
        for (int e = 0; e < 4; e++) {
            graph[f].push_back({nb_face[f][e], 100.0});
        }
    }

    auto dist3 = [&](V3 a, V3 b) {
        V3 d = sub3(a, b);
        return sqrt((coord_t)dot3(d, d));
    };

    int src_node = nf;
    int dst_node = nf + 1;

    coord_t ws = dist3(src_pt, face_center[face_of_src]);
    coord_t wd = dist3(dst_pt, face_center[face_of_dst]);

    graph[src_node].push_back({face_of_src, ws});
    graph[face_of_src].push_back({src_node, ws});

    graph[dst_node].push_back({face_of_dst, wd});
    graph[face_of_dst].push_back({dst_node, wd});

    vector<coord_t> dist(total, INF);

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

    dist[src_node] = 0;
    pq.push({0, src_node});

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

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

        for (auto [v, w] : graph[u]) {
            coord_t nd = cur + w;

            if (nd + 1e-9 < dist[v]) {
                dist[v] = nd;
                pq.push({nd, v});
            }
        }
    }

    ub = dist[dst_node];
}

// ------------------------------------------------------------
// Run one unfolding sweep from a graph node.
// starts are starting unfolded positions in incident faces.
// ------------------------------------------------------------
void run_field(
    const vector<pair<int, Point>>& starts,
    int source_node,
    vector<coord_t>& dist
) {
    int total = num_sheets + 2;

    dist.assign(total, INF);
    dist[source_node] = 0;

    vector<Window> stack;

    for (auto [face, apex] : starts) {
        Window w;

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

        stack.push_back(w);
    }

    while (!stack.empty()) {
        Window w = stack.back();
        stack.pop_back();

        int f = w.face;

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

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

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

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

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

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

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

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

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

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

            Point q0, q1;

            if (w.entry == -1) {
                // Initial face has full visibility.
                // Avoid degenerate case where source is exactly at a corner.
                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;
                }
            }

            // Prune windows already farther than current upper bound.
            if (point_seg_dist(w.apex, q0, q1) > ub + 1e-6) {
                continue;
            }

            // Safety guard.
            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];

            array<Point, 4> c2g;

            // Unfold neighboring face.
            for (int k = 0; k < 4; k++) {
                V3 cc = face_corner[g][k];

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

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

            stack.push_back(nw);
        }
    }
}

// ------------------------------------------------------------
// Full solver.
// ------------------------------------------------------------
void solve() {
    build_faces();
    build_adjacency();
    build_sheets();

    face_of_src = find_face_of(src_pt);
    face_of_dst = find_face_of(dst_pt);

    compute_ub();

    int nf = face_center.size();

    // Incident face-corners for every vertex sheet.
    vector<vector<pair<int, int>>> incident(num_sheets);

    for (int f = 0; f < nf; f++) {
        for (int k = 0; k < 4; k++) {
            int s = sheet_of[f * 4 + k];

            incident[s].push_back({f, k});
        }
    }

    int total = num_sheets + 2;

    int src_node = num_sheets;
    int dst_node = num_sheets + 1;

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

    // Sweeps from all vertex sheets.
    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]);
    }

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

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

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

    // Dijkstra on graph nodes.
    vector<coord_t> dist(total, INF);

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

    dist[src_node] = 0;
    pq.push({0, src_node});

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

        if (cur > 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 > INF / 2) continue;

            coord_t nd = cur + w;

            if (nd + 1e-12 < dist[v]) {
                dist[v] = nd;
                pq.push({nd, v});
            }
        }
    }

    cout << fixed << setprecision(10) << dist[dst_node] << '\n';
}

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

    read_input();
    solve();

    return 0;
}
```

---

## 5. Python implementation with detailed comments

```python
import sys
import math
import heapq


INF = 10**100


# ------------------------------------------------------------
# 2D point for unfolded geometry.
# ------------------------------------------------------------
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):
        return Point(self.x * value, self.y * value)

    def __truediv__(self, value):
        return Point(self.x / 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 unfolded square.
CANON = [
    Point(0, 0),
    Point(100, 0),
    Point(100, 100),
    Point(0, 100),
]


# ------------------------------------------------------------
# DSU for grouping face-corners into vertex sheets.
# ------------------------------------------------------------
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 exposed surface 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]

                # Not exposed.
                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
                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 adjacency across surface 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]

                bx = ax + DVEC[d][0]
                by = ay + DVEC[d][1]
                bz = az + DVEC[d][2]

                if not self.present(ax, ay, az):
                    # Convex fold.
                    key = (g[0], g[1], g[2], t)
                elif not self.present(bx, by, bz):
                    # Flat continuation.
                    key = (ax, ay, az, d)
                else:
                    # Concave L-shaped fold.
                    key = (bx, by, bz, opp(t))

                gface = self.face_id[key]

                self.nb_face[f][e] = gface

                # Find reverse edge index.
                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 using DSU.
    # ------------------------------------------------------------
    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 of corner k.
                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)

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

        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 face containing a 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

                if abs(p[a] - fc[a]) >= 50:
                    ok = False
                    break

            if ok:
                return f

        return -1

    # ------------------------------------------------------------
    # Convert a 3D point on face f to unfolded 2D coordinates.
    # ------------------------------------------------------------
    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

    # ------------------------------------------------------------
    # Angular sector check.
    # ------------------------------------------------------------
    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

        return (
            cross2(dl, v) >= -tol * dl.norm()
            and cross2(dr, v) <= tol * dr.norm()
        )

    # ------------------------------------------------------------
    # Clip interval [s0, s1] by a + b*s >= 0.
    # ------------------------------------------------------------
    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 c-d 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 non-shared corner after unfolding across edge p0-p1.
    # ------------------------------------------------------------
    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

    # ------------------------------------------------------------
    # Rough 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 unfolding sweep from one 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]

            # State:
            # face, entry edge, depth, source point, sector lo, sector hi,
            # unfolded face corners.
            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 S 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 D 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:
                    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

                # Prune by upper bound.
                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]

                c2g = [None] * 4

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

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

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

        return dist

    # ------------------------------------------------------------
    # Full solver.
    # ------------------------------------------------------------
    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):
                s = self.sheet_of[f * 4 + k]
                incident[s].append((f, k))

        total = self.num_sheets + 2

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

        field = [None] * total

        # Sweeps from all 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)

        # 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 > INF / 2:
                    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 = data[idx]
        y = data[idx + 1]
        z = 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()
```