## 1. Abridged problem statement

Given a rectangular box asylum with opposite corners `(0,0,0)` and `(W,H,D)`, and `N ≤ 10` existing spherical ghosts, find a point where the center of a new spherical ghost can be placed so that its radius is as large as possible.

The new ghost:

- must lie completely inside the box,
- must not intersect any existing ghost,
- may touch walls or existing ghosts.

Existing ghosts may intersect each other or the box boundaries.

Output any center point whose achievable radius is within `1 mm` of optimal.

---

## 2. Detailed editorial

### Reformulation

Suppose we want to check whether a new ghost of radius `r` can be placed.

Its center `p = (x,y,z)` must satisfy:

1. It must be at least distance `r` from every wall:

\[
r \le x \le W-r
\]
\[
r \le y \le H-r
\]
\[
r \le z \le D-r
\]

2. It must be outside every old ghost inflated by `r`.

If an old ghost has center `c_i` and radius `R_i`, then:

\[
|p-c_i| \ge R_i+r
\]

So for fixed `r`, we need to know whether the set

\[
F(r)=\{p : p \text{ is inside the shrunken box and outside all inflated ghosts}\}
\]

is non-empty.

The answer radius is monotone:

- if radius `r` is feasible, then every smaller radius is feasible,
- if radius `r` is impossible, every larger radius is impossible.

Therefore we can binary search the maximum radius.

---

### Feasibility checking

For a fixed `r`, the boundary surfaces are:

- six box planes:

\[
x=r,\ x=W-r,\ y=r,\ y=H-r,\ z=r,\ z=D-r
\]

- `N` inflated ghost spheres:

\[
|p-c_i|=R_i+r
\]

There are at most `6 + 10 = 16` surfaces.

If the feasible region is non-empty, there is a feasible point that is locally extremal in some direction. Such an extremal point must lie on at least three boundary surfaces.

So it is enough to enumerate all triples of boundary surfaces and test their intersection points.

There are only:

\[
\binom{16}{3}=560
\]

triples, which is tiny.

---

### Intersections of three surfaces

There are two cases.

---

#### Case 1: Three planes

Three box planes intersect in one point if they correspond to three different axes.

For example:

\[
x=r,\ y=H-r,\ z=D-r
\]

gives one corner of the shrunken box.

We test that point.

---

#### Case 2: At least one sphere

Choose one sphere as the main sphere.

The other two surfaces are converted into planes:

- if the surface is already a box plane, keep it,
- if it is another sphere, subtract its equation from the main sphere equation.

For two spheres:

\[
|p-c_0|^2=\rho_0^2
\]

\[
|p-c_j|^2=\rho_j^2
\]

Subtracting gives a plane:

\[
2(c_j-c_0)\cdot p
=
\rho_0^2-\rho_j^2-|c_0|^2+|c_j|^2
\]

This is the radical plane of the two spheres.

Now the problem becomes:

- intersect two planes, obtaining a line,
- intersect that line with the main sphere.

A line can intersect a sphere in:

- zero points,
- one tangent point,
- two points.

Each candidate point is checked against all constraints.

---

### Binary search

The largest possible radius is at most:

\[
\frac{\min(W,H,D)}{2}
\]

We binary search in `[0, min(W,H,D)/2]`.

The provided solution uses 45 iterations, enough for much better than `1 mm` accuracy.

After the binary search, it calls `feasible(lo)` once more to recover a valid witness point for the final radius.

---

## 3. C++ Solution

```cpp
#include <bits/stdc++.h>

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

long double W, H, D;
int n;
vector<array<long double, 4>> ghosts;
array<long double, 3> witness;

void read() {
    cin >> W >> H >> D >> n;
    ghosts.assign(n, {});
    for(auto& g: ghosts) {
        cin >> g[0] >> g[1] >> g[2] >> g[3];
    }
}

using vec = array<long double, 3>;

vec operator-(vec a, vec b) { return {a[0] - b[0], a[1] - b[1], a[2] - b[2]}; }

long double dot(vec a, vec b) {
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}

vec cross(vec a, vec b) {
    return {
        a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2],
        a[0] * b[1] - a[1] * b[0]
    };
}

bool valid_center(vec p, long double r) {
    long double eps = 1e-6L;
    if(p[0] < r - eps || p[0] > W - r + eps || p[1] < r - eps ||
       p[1] > H - r + eps || p[2] < r - eps || p[2] > D - r + eps) {
        return false;
    }
    for(auto& g: ghosts) {
        long double dx = p[0] - g[0], dy = p[1] - g[1], dz = p[2] - g[2];
        long double rho = r + g[3];
        if(sqrtl(dx * dx + dy * dy + dz * dz) < rho - eps) {
            return false;
        }
    }
    return true;
}

bool line_sphere(
    vec n1, long double d1, vec n2, long double d2, vec c, long double rho,
    long double r
) {
    vec u = cross(n1, n2);
    long double lu2 = dot(u, u);
    if(lu2 < 1e-18L) {
        return false;
    }

    vec p0;
    vec t1 = cross(n2, u), t2 = cross(u, n1);
    for(int a = 0; a < 3; a++) {
        p0[a] = (d1 * t1[a] + d2 * t2[a]) / lu2;
    }

    vec w = p0 - c;
    long double b = 2 * dot(u, w);
    long double cc = dot(w, w) - rho * rho;
    long double disc = b * b - 4 * lu2 * cc;
    if(disc < 0) {
        disc = 0;
    }
    long double sq = sqrtl(disc);
    for(int s = -1; s <= 1; s += 2) {
        long double t = (-b + s * sq) / (2 * lu2);
        vec p = {p0[0] + t * u[0], p0[1] + t * u[1], p0[2] + t * u[2]};
        if(valid_center(p, r)) {
            witness = p;
            return true;
        }
    }
    return false;
}

bool feasible(long double r) {
    if(2 * r > W || 2 * r > H || 2 * r > D) {
        return false;
    }

    int m = 6 + n;
    vec pn[6] = {{1, 0, 0}, {1, 0, 0}, {0, 1, 0},
                 {0, 1, 0}, {0, 0, 1}, {0, 0, 1}};
    long double pd[6] = {r, W - r, r, H - r, r, D - r};

    auto sphere_center = [&](int id) -> vec {
        return {ghosts[id - 6][0], ghosts[id - 6][1], ghosts[id - 6][2]};
    };

    auto sphere_rho = [&](int id) -> long double {
        return r + ghosts[id - 6][3];
    };

    for(int a = 0; a < m; a++) {
        for(int b = a + 1; b < m; b++) {
            for(int c = b + 1; c < m; c++) {
                int id[3] = {a, b, c};
                int sph = -1;
                for(int t = 0; t < 3; t++) {
                    if(id[t] >= 6) {
                        sph = id[t];
                    }
                }

                if(sph == -1) {
                    if(id[0] / 2 != id[1] / 2 && id[1] / 2 != id[2] / 2 &&
                       id[0] / 2 != id[2] / 2) {
                        vec p;
                        p[id[0] / 2] = pd[id[0]];
                        p[id[1] / 2] = pd[id[1]];
                        p[id[2] / 2] = pd[id[2]];
                        if(valid_center(p, r)) {
                            witness = p;
                            return true;
                        }
                    }

                    continue;
                }

                vec rc = sphere_center(sph);
                long double rrho = sphere_rho(sph);
                vec pl[2];
                long double pld[2];
                int cnt = 0;
                for(int t = 0; t < 3; t++) {
                    if(id[t] == sph) {
                        continue;
                    }

                    if(id[t] < 6) {
                        pl[cnt] = pn[id[t]];
                        pld[cnt] = pd[id[t]];
                    } else {
                        vec oc = sphere_center(id[t]);
                        long double orho = sphere_rho(id[t]);
                        pl[cnt] = {
                            2 * (oc[0] - rc[0]), 2 * (oc[1] - rc[1]),
                            2 * (oc[2] - rc[2])
                        };

                        pld[cnt] = rrho * rrho - orho * orho - dot(rc, rc) +
                                   dot(oc, oc);
                    }

                    cnt++;
                }

                if(line_sphere(pl[0], pld[0], pl[1], pld[1], rc, rrho, r)) {
                    return true;
                }
            }
        }
    }

    return false;
}

void solve() {
    // The achievable radius is monotone, so we binary search it: a candidate
    // radius r works iff some center can host a sphere of that radius. Fixing r
    // pins the placement region exactly: the center must sit in the shrunken
    // box [r, W - r] x [r, H - r] x [r, D - r] (so the sphere clears the walls)
    // and outside every old ghost inflated by r, i.e. at distance at least
    // R_i + r from ghost i's center. So r is feasible iff this region, an
    // intersection of six half-spaces and n ball-exteriors, is non-empty.
    //
    // With N <= 10 we decide that non-emptiness exactly rather than by
    // sampling. Write the region as F(r) = { p : p in the shrunken box and |p -
    // c_i| >= R_i + r for every i }; as a closed subset of the box it is
    // compact, so a linear functional attains a maximum on it. Maximize one
    // that is generic in the sense of ordering points lexicographically by (x,
    // then y, then z), and let p* be the maximizer. The constraints active at
    // p* are the bounding surfaces passing through it; if at most two were
    // active, their zero sets would meet in a curve (or larger set) through p*
    // along which p* could move a little while staying in F(r), and any such
    // move changes some coordinate, so it strictly raises the lexicographic
    // value, contradicting maximality. Hence p* lies on at least three surfaces
    // at once, and it suffices to test only the points where three surfaces
    // meet.
    //
    // The surfaces are the 6 box planes x = r, x = W - r, ..., z = D - r and
    // the n spheres |p - c_i| = R_i + r, at most 16 in all, so we enumerate all
    // C(16, 3) triples and solve each for its meeting point(s). Three planes
    // meet in a single point only when they use three distinct axes, giving one
    // box corner. A triple with at least one sphere is linearized: keep one
    // sphere S0 as the quadric and turn the other two surfaces into planes. A
    // plane stays itself, and a second sphere Sj is replaced by its radical
    // plane with S0, the plane |p - c0|^2 - (R0 + r)^2 = |p - cj|^2 - (Rj +
    // r)^2 that holds exactly on S0 ∩ Sj. The two planes meet in a line
    // (skipped if parallel), which meets S0 in at most two points through one
    // quadratic. A candidate is accepted when it lies in the box and at
    // distance >= R_i + r from every ghost; if any triple yields such a point
    // then r is feasible, and the center kept at the largest feasible r is the
    // answer.

    long double lo = 0, hi = min({W, H, D}) / 2;
    witness = {W / 2, H / 2, D / 2};
    for(int it = 0; it < 45; it++) {
        long double mid = (lo + hi) / 2;
        if(feasible(mid)) {
            lo = mid;
        } else {
            hi = mid;
        }
    }

    feasible(lo);
    cout << fixed << setprecision(12);
    cout << witness[0] << ' ' << witness[1] << ' ' << witness[2] << '\n';
}

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

    int T = 1;
    // cin >> T;
    for(int test = 1; test <= T; test++) {
        read();
        solve();
    }

    return 0;
}
```
---

## 4. Python solution with detailed comments

```python
import sys
import math
from itertools import combinations


# ---------- Basic vector operations ----------

def sub(a, b):
    """Return a - b for 3D vectors."""
    return [a[0] - b[0], a[1] - b[1], a[2] - b[2]]


def dot(a, b):
    """Return dot product of two 3D vectors."""
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]


def cross(a, b):
    """Return cross product of two 3D vectors."""
    return [
        a[1] * b[2] - a[2] * b[1],
        a[2] * b[0] - a[0] * b[2],
        a[0] * b[1] - a[1] * b[0],
    ]


# ---------- Input ----------

data = sys.stdin.read().strip().split()

W = float(data[0])
H = float(data[1])
D = float(data[2])

n = int(data[3])

ghosts = []
ptr = 4

for _ in range(n):
    x = float(data[ptr])
    y = float(data[ptr + 1])
    z = float(data[ptr + 2])
    r = float(data[ptr + 3])
    ptr += 4
    ghosts.append((x, y, z, r))


# Global witness point.
witness = [W / 2.0, H / 2.0, D / 2.0]


def valid_center(p, r):
    """
    Check whether p can be the center of a new ghost with radius r.
    """
    eps = 1e-7

    # Check box constraints.
    if p[0] < r - eps or p[0] > W - r + eps:
        return False
    if p[1] < r - eps or p[1] > H - r + eps:
        return False
    if p[2] < r - eps or p[2] > D - r + eps:
        return False

    # Check distance from every old ghost.
    for gx, gy, gz, gr in ghosts:
        dx = p[0] - gx
        dy = p[1] - gy
        dz = p[2] - gz

        need = r + gr
        dist2 = dx * dx + dy * dy + dz * dz

        if dist2 < (need - eps) * (need - eps):
            return False

    return True


def line_sphere(n1, d1, n2, d2, c, rho, r):
    """
    Intersect:
        n1 · p = d1
        n2 · p = d2
        |p - c| = rho

    Test produced points as possible centers.
    """
    global witness

    # Direction of intersection line of the two planes.
    u = cross(n1, n2)
    lu2 = dot(u, u)

    # Parallel or almost-parallel planes.
    if lu2 < 1e-18:
        return False

    # Find one point p0 on the line of intersection.
    t1 = cross(n2, u)
    t2 = cross(u, n1)

    p0 = [
        (d1 * t1[i] + d2 * t2[i]) / lu2
        for i in range(3)
    ]

    # Parametrize line as p(t) = p0 + t*u.
    w = sub(p0, c)

    # Quadratic equation:
    # |p0 + t*u - c|^2 = rho^2
    b = 2.0 * dot(u, w)
    cc = dot(w, w) - rho * rho

    disc = b * b - 4.0 * lu2 * cc

    # If negative only due to floating error, clamp to zero.
    if disc < 0.0:
        disc = 0.0

    sq = math.sqrt(disc)

    # Try both roots.
    for sign in (-1.0, 1.0):
        t = (-b + sign * sq) / (2.0 * lu2)

        p = [
            p0[0] + t * u[0],
            p0[1] + t * u[1],
            p0[2] + t * u[2],
        ]

        if valid_center(p, r):
            witness = p
            return True

    return False


def feasible(r):
    """
    Return True if a new ghost of radius r can be placed.
    Also stores one valid center in witness.
    """
    global witness

    # Radius too large for the box.
    if 2.0 * r > W or 2.0 * r > H or 2.0 * r > D:
        return False

    # There are six box planes and n ghost spheres.
    m = 6 + n

    # Normals of six box planes.
    plane_normals = [
        [1.0, 0.0, 0.0],
        [1.0, 0.0, 0.0],
        [0.0, 1.0, 0.0],
        [0.0, 1.0, 0.0],
        [0.0, 0.0, 1.0],
        [0.0, 0.0, 1.0],
    ]

    # Constants of six box planes.
    plane_ds = [
        r,
        W - r,
        r,
        H - r,
        r,
        D - r,
    ]

    def sphere_center(surface_id):
        """Return center of ghost sphere surface."""
        gx, gy, gz, _ = ghosts[surface_id - 6]
        return [gx, gy, gz]

    def sphere_radius(surface_id):
        """Return inflated radius of ghost sphere surface."""
        return ghosts[surface_id - 6][3] + r

    # Enumerate all triples of surfaces.
    for a, b, c in combinations(range(m), 3):
        ids = [a, b, c]

        # Pick one sphere from the triple, if any.
        sph = -1
        for surface_id in ids:
            if surface_id >= 6:
                sph = surface_id

        # Case 1: three box planes.
        if sph == -1:
            # They must correspond to three different axes.
            axes = [surface_id // 2 for surface_id in ids]

            if len(set(axes)) == 3:
                p = [0.0, 0.0, 0.0]

                for surface_id in ids:
                    axis = surface_id // 2
                    p[axis] = plane_ds[surface_id]

                if valid_center(p, r):
                    witness = p
                    return True

            continue

        # Case 2: at least one sphere.
        # Use this sphere as the main sphere.
        rc = sphere_center(sph)
        rrho = sphere_radius(sph)

        # Convert the other two surfaces into planes.
        planes_n = []
        planes_d = []

        for surface_id in ids:
            if surface_id == sph:
                continue

            if surface_id < 6:
                # Box plane.
                planes_n.append(plane_normals[surface_id])
                planes_d.append(plane_ds[surface_id])
            else:
                # Sphere converted to radical plane.
                oc = sphere_center(surface_id)
                orho = sphere_radius(surface_id)

                normal = [
                    2.0 * (oc[0] - rc[0]),
                    2.0 * (oc[1] - rc[1]),
                    2.0 * (oc[2] - rc[2]),
                ]

                dval = (
                    rrho * rrho
                    - orho * orho
                    - dot(rc, rc)
                    + dot(oc, oc)
                )

                planes_n.append(normal)
                planes_d.append(dval)

        # Intersect two planes and the main sphere.
        if line_sphere(
            planes_n[0],
            planes_d[0],
            planes_n[1],
            planes_d[1],
            rc,
            rrho,
            r,
        ):
            return True

    return False


# ---------- Binary search ----------

lo = 0.0
hi = min(W, H, D) / 2.0

# 60 iterations are plenty for millimeter accuracy.
for _ in range(60):
    mid = (lo + hi) / 2.0

    if feasible(mid):
        lo = mid
    else:
        hi = mid

# Recover witness for final radius.
feasible(lo)

print("{:.12f} {:.12f} {:.12f}".format(witness[0], witness[1], witness[2]))
```

---

## 5. Compressed editorial

Binary search the answer radius `r`.

For fixed `r`, the center must lie inside the shrunken box:

\[
[r,W-r]\times[r,H-r]\times[r,D-r]
\]

and outside every inflated old ghost sphere:

\[
|p-c_i|\ge R_i+r
\]

The boundary surfaces are six box planes and at most ten spheres. If the feasible region is non-empty, some feasible extremal point lies on at least three boundary surfaces. Therefore enumerate all triples of surfaces.

For three planes, test their unique box-corner intersection if they use three different axes.

For triples containing a sphere, keep one sphere as quadratic. Convert any other sphere into a radical plane by subtracting equations, and keep box planes as planes. Then intersect two planes into a line and intersect that line with the chosen sphere, giving at most two candidate points. Test each candidate against all constraints.

Feasibility is monotone in `r`, so binary search up to `min(W,H,D)/2`. Output the witness point from the final feasible radius.
