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

469. Ghostbusters
Time limit per test: 0.5 second(s)
Memory limit: 262144 kilobytes
input: standard
output: standard



Now famous Ghostbusters want your help. As you probably know, all ghosts that are caught are put in asylum to stay there forever (since they are already dead and cannot be obliterated). The problem is that though usually ghosts can intersect and even lie inside each other, some new ghosts by yet unknown to Egon reasons are afraid of some old ghosts and cannot intersect with them. That's what your problem is — to find the maximum radius of a new ghost that can be put to the asylum in the worst case (i.e., when it is afraid of all other ghosts in the asylum). You task is simplified by the fact that all ghosts are in fact perfect solid spheres and stay at the same points forever. Also don't forget that new ghost cannot intersect borders of the asylum (but can touch them as well as other ghosts).

Input
The first line of the input file contains three integer numbers W, H and D — sizes of the asylum given in meters (). Asylum is a rectangular parallelepiped. Let us introduce coordinate system in such a way that two opposite corners of the asylum have coordinates (0,0,0) and (W,H,D). Then position of the i-th ghost can be described in this system as sphere centered at point with coordinates (Xi,Yi,Zi) and having radius Ri. Second line of the input file contains integer number N — number of ghosts in the asylum (0 ≤ N ≤ 10). i-th of the next N lines contains four numbers Xi, Yi, Zi and Ri (-1000 < Xi,Yi,Zi <2000, 0 ≤ Ri ≤ 1000). Though new ghost cannot intersect other ghosts and borders of the asylum, old ghosts may violate these rules and cross each other and boundaries of the asylum. However, it is guaranteed that at least some part of a ghost is contained inside the asylum. Ghost of radius 0 is supposed to consist of one point and cannot lie inside your ghost.

Output
To the output file write three numbers — coordinates of a point where ghost of required maximal radius, centered at this point, may stay. If there are several such points, output coordinates of any of them (but only one). Your output will be considered correct if maximal radius of a ghost centered at the point found by your program will be within 1 mm of the optimal one. To avoid any precision problems we recommend to output coordinates with maximal possible accuracy. It is guaranteed that at least one solution with positive radius will exist.

Example(s)
sample input
sample output
9 13 15
2
3 9 6 1
4 8 3 5
5 4 11



In the example, the first ghost is inside the second one. Ghost of maximum radius of 4 can stay at the point (5, 4, 11) — in this case it touches three borders of the asylum and the second ghost.

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

We are given a rectangular box with corners `(0, 0, 0)` and `(W, H, D)`. Inside or around it there are `N ≤ 10` old spherical ghosts, each with center `(Xi, Yi, Zi)` and radius `Ri`.

We need to place a new spherical ghost inside the box so that:

- it does not intersect the box borders,
- it does not intersect any old ghost,
- touching is allowed.

Output any center point of a new ghost whose possible radius is within `1 mm` of the optimal maximum radius.

---

## 2. Key observations needed to solve the problem

### Observation 1: Feasibility for a fixed radius

Suppose we want to check whether radius `r` is possible.

Then the center `p = (x, y, z)` of the new ghost must satisfy:

```text
r ≤ x ≤ W - r
r ≤ y ≤ H - r
r ≤ z ≤ D - r
```

So the center must lie inside the box shrunk by distance `r`.

Also, for every old ghost `i`:

```text
distance(p, center_i) ≥ r + Ri
```

That means the center must lie outside every old ghost sphere inflated by `r`.

So for fixed `r`, we need to know whether this region is non-empty.

---

### Observation 2: The answer is monotonic

If radius `r` is possible, then any smaller radius is also possible.

Therefore we can binary search the maximum radius.

The upper bound is:

```text
min(W, H, D) / 2
```

---

### Observation 3: Only few boundary surfaces matter

For fixed `r`, the boundary surfaces are:

- 6 planes from the shrunken box:

```text
x = r, x = W - r
y = r, y = H - r
z = r, z = D - r
```

- `N` inflated ghost spheres:

```text
|p - c_i| = Ri + r
```

Since `N ≤ 10`, there are at most `16` surfaces.

If a feasible region is non-empty, one can find a feasible extremal point lying on at least three boundary surfaces. Therefore, we can enumerate all triples of surfaces and test their intersection points.

There are at most:

```text
C(16, 3) = 560
```

triples, which is very small.

---

## 3. Full solution approach based on the observations

We binary search the answer radius `r`.

For each candidate radius `r`, we check feasibility as follows.

There are two types of boundary surfaces:

1. Box planes.
2. Inflated old ghost spheres.

We enumerate every triple of surfaces.

---

### Case 1: Three planes

Three box planes give one point only if they correspond to three different axes. We test whether this corner is valid.

---

### Case 2: At least one sphere

Pick one sphere as the main sphere:

```text
|p - c0| = rho0
```

The other two surfaces are converted to planes.

If the surface is already a box plane, keep it.

If the surface is another sphere, subtract the two sphere equations to get the radical plane:

```text
2(c1 - c0) · p = rho0² - rho1² - |c0|² + |c1|²
```

So the problem becomes: two planes intersect in a line; that line intersects the main sphere in at most two points. We test those candidate points.

---

### Validity check

A candidate point `p` is valid for radius `r` if it lies inside the shrunken box and for every old ghost, `distance(p, center_i) ≥ r + Ri`.

If any candidate point is valid, radius `r` is feasible.

---

### Precision

The required accuracy is `1 mm = 0.001`, so 45 iterations of binary search are more than enough.

Finally, we run the feasibility check once more for the final radius and output the saved witness point.

---

## 4. C++ implementation with comments

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

---

## 5. Python implementation with comments

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


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


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


def cross(a, 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],
    ]


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


witness = [W / 2.0, H / 2.0, D / 2.0]


def valid_center(p, r):
    eps = 1e-7
    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
    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):
    global witness
    u = cross(n1, n2)
    lu2 = dot(u, u)
    if lu2 < 1e-18:
        return False
    t1 = cross(n2, u)
    t2 = cross(u, n1)
    p0 = [(d1 * t1[i] + d2 * t2[i]) / lu2 for i in range(3)]
    w = sub(p0, c)
    b = 2.0 * dot(u, w)
    cc = dot(w, w) - rho * rho
    disc = b * b - 4.0 * lu2 * cc
    if disc < 0.0:
        disc = 0.0
    sq = math.sqrt(disc)
    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):
    global witness
    if 2.0 * r > W or 2.0 * r > H or 2.0 * r > D:
        return False
    m = 6 + n
    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],
    ]
    plane_ds = [r, W - r, r, H - r, r, D - r]

    def sphere_center(sid):
        gx, gy, gz, _ = ghosts[sid - 6]
        return [gx, gy, gz]

    def sphere_radius(sid):
        return ghosts[sid - 6][3] + r

    for a, b, c in combinations(range(m), 3):
        ids = [a, b, c]
        sph = -1
        for sid in ids:
            if sid >= 6:
                sph = sid
        if sph == -1:
            axes = [sid // 2 for sid in ids]
            if len(set(axes)) == 3:
                p = [0.0, 0.0, 0.0]
                for sid in ids:
                    p[sid // 2] = plane_ds[sid]
                if valid_center(p, r):
                    witness = p
                    return True
            continue
        rc = sphere_center(sph)
        rrho = sphere_radius(sph)
        planes_n = []
        planes_d = []
        for sid in ids:
            if sid == sph:
                continue
            if sid < 6:
                planes_n.append(plane_normals[sid])
                planes_d.append(plane_ds[sid])
            else:
                oc = sphere_center(sid)
                orho = sphere_radius(sid)
                normal = [2.0 * (oc[i] - rc[i]) for i in range(3)]
                dval = rrho * rrho - orho * orho - dot(rc, rc) + dot(oc, oc)
                planes_n.append(normal)
                planes_d.append(dval)
        if line_sphere(planes_n[0], planes_d[0], planes_n[1], planes_d[1], rc, rrho, r):
            return True
    return False


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

for _ in range(60):
    mid = (lo + hi) / 2.0
    if feasible(mid):
        lo = mid
    else:
        hi = mid

feasible(lo)

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