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

Example:

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

This gives one corner of the shrunken box.

We test whether this point 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:

```text
|p - c1| = rho1
```

Subtract the two sphere equations:

```text
|p - c0|² - rho0² = |p - c1|² - rho1²
```

This becomes a plane called the radical plane.

So the problem becomes:

```text
plane 1
plane 2
main sphere
```

Two planes intersect in a line. A line intersects a 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,
- for every old ghost:

```text
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 50–60 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.

Complexity:

```text
Binary search iterations: about 60
Triples per check: at most 560
Validation per point: O(N), N ≤ 10
```

This is easily fast enough.

---

## 4. C++ implementation with detailed comments

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

using ld = long double;

struct Vec {
    ld x, y, z;
};

Vec operator+(const Vec& a, const Vec& b) {
    return {a.x + b.x, a.y + b.y, a.z + b.z};
}

Vec operator-(const Vec& a, const Vec& b) {
    return {a.x - b.x, a.y - b.y, a.z - b.z};
}

Vec operator*(const Vec& a, ld k) {
    return {a.x * k, a.y * k, a.z * k};
}

ld dot(const Vec& a, const Vec& b) {
    return a.x * b.x + a.y * b.y + a.z * b.z;
}

Vec cross(const Vec& a, const Vec& b) {
    return {
        a.y * b.z - a.z * b.y,
        a.z * b.x - a.x * b.z,
        a.x * b.y - a.y * b.x
    };
}

struct Ghost {
    Vec c;
    ld r;
};

ld W, H, D;
int N;
vector<Ghost> ghosts;

// A valid center found by the latest feasibility check.
Vec witness;

bool valid_center(const Vec& p, ld r) {
    const ld EPS = 1e-7L;

    // The new sphere must stay inside the box.
    if (p.x < r - EPS || p.x > W - r + EPS) return false;
    if (p.y < r - EPS || p.y > H - r + EPS) return false;
    if (p.z < r - EPS || p.z > D - r + EPS) return false;

    // The new sphere must not intersect any old ghost.
    for (const Ghost& g : ghosts) {
        Vec diff = p - g.c;
        ld dist2 = dot(diff, diff);
        ld need = r + g.r;

        if (dist2 < (need - EPS) * (need - EPS)) {
            return false;
        }
    }

    return true;
}

/*
    Intersect two planes and one sphere.

    Plane 1: n1 · p = d1
    Plane 2: n2 · p = d2
    Sphere : |p - center| = rho

    The intersection of the two planes is a line.
    Then we intersect that line with the sphere.
*/
bool line_sphere(
    Vec n1, ld d1,
    Vec n2, ld d2,
    Vec center, ld rho,
    ld current_radius
) {
    // Direction vector of the line of intersection of two planes.
    Vec u = cross(n1, n2);
    ld lu2 = dot(u, u);

    // Parallel planes do not define a unique line.
    if (lu2 < 1e-18L) {
        return false;
    }

    /*
        Find one point p0 on the intersection line.

        Formula:
        p0 = (d1 * (n2 × u) + d2 * (u × n1)) / |u|²
    */
    Vec t1 = cross(n2, u);
    Vec t2 = cross(u, n1);

    Vec p0 = {
        (d1 * t1.x + d2 * t2.x) / lu2,
        (d1 * t1.y + d2 * t2.y) / lu2,
        (d1 * t1.z + d2 * t2.z) / lu2
    };

    /*
        Any point on the line is:

            p(t) = p0 + t * u

        We solve:

            |p0 + t * u - center|² = rho²
    */
    Vec w = p0 - center;

    ld A = lu2;
    ld B = 2.0L * dot(u, w);
    ld C = dot(w, w) - rho * rho;

    ld disc = B * B - 4.0L * A * C;

    // Clamp small or even negative values; final validity check is decisive.
    if (disc < 0) {
        disc = 0;
    }

    ld sq = sqrtl(disc);

    for (int sign : {-1, 1}) {
        ld t = (-B + sign * sq) / (2.0L * A);

        Vec p = p0 + u * t;

        if (valid_center(p, current_radius)) {
            witness = p;
            return true;
        }
    }

    return false;
}

bool feasible(ld r) {
    // The sphere of radius r cannot fit if diameter exceeds any dimension.
    if (2.0L * r > W || 2.0L * r > H || 2.0L * r > D) {
        return false;
    }

    /*
        Surface IDs:

        0: x = r
        1: x = W - r
        2: y = r
        3: y = H - r
        4: z = r
        5: z = D - r

        6 ... 6 + N - 1: inflated ghost spheres
    */

    int total_surfaces = 6 + N;

    Vec plane_normal[6] = {
        {1, 0, 0}, {1, 0, 0},
        {0, 1, 0}, {0, 1, 0},
        {0, 0, 1}, {0, 0, 1}
    };

    ld plane_d[6] = {
        r, W - r,
        r, H - r,
        r, D - r
    };

    auto sphere_center = [&](int id) -> Vec {
        return ghosts[id - 6].c;
    };

    auto sphere_radius = [&](int id) -> ld {
        return ghosts[id - 6].r + r;
    };

    // Enumerate all triples of boundary surfaces.
    for (int a = 0; a < total_surfaces; a++) {
        for (int b = a + 1; b < total_surfaces; b++) {
            for (int c = b + 1; c < total_surfaces; c++) {
                int ids[3] = {a, b, c};

                // Find one sphere in the triple, if any.
                int main_sphere = -1;
                for (int i = 0; i < 3; i++) {
                    if (ids[i] >= 6) {
                        main_sphere = ids[i];
                    }
                }

                // Case 1: all three surfaces are planes.
                if (main_sphere == -1) {
                    int ax0 = ids[0] / 2;
                    int ax1 = ids[1] / 2;
                    int ax2 = ids[2] / 2;

                    // Need one x-plane, one y-plane, one z-plane.
                    if (ax0 != ax1 && ax0 != ax2 && ax1 != ax2) {
                        Vec p = {0, 0, 0};

                        for (int i = 0; i < 3; i++) {
                            int id = ids[i];
                            int axis = id / 2;

                            if (axis == 0) p.x = plane_d[id];
                            if (axis == 1) p.y = plane_d[id];
                            if (axis == 2) p.z = plane_d[id];
                        }

                        if (valid_center(p, r)) {
                            witness = p;
                            return true;
                        }
                    }

                    continue;
                }

                /*
                    Case 2: at least one sphere.

                    Keep one sphere as the main sphere.
                    Convert the other two surfaces into planes.
                */
                Vec c0 = sphere_center(main_sphere);
                ld rho0 = sphere_radius(main_sphere);

                Vec planes_n[2];
                ld planes_d[2];
                int cnt = 0;

                for (int i = 0; i < 3; i++) {
                    int id = ids[i];

                    if (id == main_sphere) {
                        continue;
                    }

                    if (id < 6) {
                        // Box plane.
                        planes_n[cnt] = plane_normal[id];
                        planes_d[cnt] = plane_d[id];
                    } else {
                        /*
                            Another sphere.

                            Main sphere:
                                |p - c0|² = rho0²

                            Other sphere:
                                |p - c1|² = rho1²

                            Subtracting gives a plane:
                                2(c1 - c0) · p =
                                rho0² - rho1² - |c0|² + |c1|²
                        */
                        Vec c1 = sphere_center(id);
                        ld rho1 = sphere_radius(id);

                        planes_n[cnt] = {
                            2.0L * (c1.x - c0.x),
                            2.0L * (c1.y - c0.y),
                            2.0L * (c1.z - c0.z)
                        };

                        planes_d[cnt] =
                            rho0 * rho0 - rho1 * rho1
                            - dot(c0, c0) + dot(c1, c1);
                    }

                    cnt++;
                }

                if (line_sphere(
                    planes_n[0], planes_d[0],
                    planes_n[1], planes_d[1],
                    c0, rho0,
                    r
                )) {
                    return true;
                }
            }
        }
    }

    return false;
}

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

    cin >> W >> H >> D;
    cin >> N;

    ghosts.resize(N);

    for (int i = 0; i < N; i++) {
        cin >> ghosts[i].c.x >> ghosts[i].c.y >> ghosts[i].c.z >> ghosts[i].r;
    }

    ld lo = 0;
    ld hi = min({W, H, D}) / 2.0L;

    witness = {W / 2.0L, H / 2.0L, D / 2.0L};

    // 60 iterations give far better than millimeter precision.
    for (int it = 0; it < 60; it++) {
        ld mid = (lo + hi) / 2.0L;

        if (feasible(mid)) {
            lo = mid;
        } else {
            hi = mid;
        }
    }

    // Recover a valid center for the final radius.
    feasible(lo);

    cout << fixed << setprecision(12);
    cout << (double)witness.x << ' '
         << (double)witness.y << ' '
         << (double)witness.z << '\n';

    return 0;
}
```

---

## 5. Python implementation with detailed comments

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


# ---------------- Vector operations ----------------

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


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


def mul(a, k):
    return [a[0] * k, a[1] * k, a[2] * k]


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


# ---------------- 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 if p can be the center of a new ghost with radius r.
    """
    eps = 1e-7

    # 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

    # Distance constraints against old ghosts.
    for center, gr in ghosts:
        diff = sub(p, center)
        dist2 = dot(diff, diff)

        need = r + gr

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

    return True


def line_sphere(n1, d1, n2, d2, center, rho, current_radius):
    """
    Intersect:

        n1 · p = d1
        n2 · p = d2
        |p - center| = rho

    The first two equations define a line.
    Then intersect that line with the sphere.
    """
    global witness

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

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

    """
    Find one point p0 on the line.

    Formula:
        p0 = (d1 * (n2 × u) + d2 * (u × n1)) / |u|²
    """
    t1 = cross(n2, u)
    t2 = cross(u, n1)

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

    # Parametrize the line:
    # p(t) = p0 + t * u
    w = sub(p0, center)

    # Solve:
    # |p0 + t*u - center|² = rho²
    A = lu2
    B = 2.0 * dot(u, w)
    C = dot(w, w) - rho * rho

    disc = B * B - 4.0 * A * C

    # Clamp negative values; final validity check ensures correctness.
    if disc < 0.0:
        disc = 0.0

    sq = math.sqrt(disc)

    for sign in (-1.0, 1.0):
        t = (-B + sign * sq) / (2.0 * A)

        p = add(p0, mul(u, t))

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

    return False


def feasible(r):
    """
    Return True if radius r is feasible.
    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

    """
    Surface IDs:

    0: x = r
    1: x = W - r
    2: y = r
    3: y = H - r
    4: z = r
    5: z = D - r

    6 ... 6 + N - 1: inflated ghost spheres
    """

    total_surfaces = 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(surface_id):
        return ghosts[surface_id - 6][0]

    def sphere_radius(surface_id):
        return ghosts[surface_id - 6][1] + r

    # Enumerate every triple of boundary surfaces.
    for a, b, c in combinations(range(total_surfaces), 3):
        ids = [a, b, c]

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

        # Case 1: all three surfaces are box planes.
        if main_sphere == -1:
            axes = [surface_id // 2 for surface_id in ids]

            # Need one x-plane, one y-plane, one z-plane.
            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.

        Keep one sphere as the main sphere.
        Convert the other two surfaces into planes.
        """
        c0 = sphere_center(main_sphere)
        rho0 = sphere_radius(main_sphere)

        planes_n = []
        planes_d = []

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

            if surface_id < 6:
                # Box plane remains unchanged.
                planes_n.append(plane_normals[surface_id])
                planes_d.append(plane_ds[surface_id])
            else:
                """
                Another sphere.

                Main sphere:
                    |p - c0|² = rho0²

                Other sphere:
                    |p - c1|² = rho1²

                Subtracting gives the radical plane.
                """
                c1 = sphere_center(surface_id)
                rho1 = sphere_radius(surface_id)

                normal = [
                    2.0 * (c1[0] - c0[0]),
                    2.0 * (c1[1] - c0[1]),
                    2.0 * (c1[2] - c0[2]),
                ]

                dval = (
                    rho0 * rho0
                    - rho1 * rho1
                    - dot(c0, c0)
                    + dot(c1, c1)
                )

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

        if line_sphere(
            planes_n[0],
            planes_d[0],
            planes_n[1],
            planes_d[1],
            c0,
            rho0,
            r,
        ):
            return True

    return False


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

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

# 60 iterations are much more accurate than 1 millimeter.
for _ in range(60):
    mid = (lo + hi) / 2.0

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

# Recover witness point for the final radius.
feasible(lo)

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