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

110. Dungeon

time limit per test: 0.25 sec.
memory limit per test: 4096 KB


The mission of space explorers found on planet M the vast dungeon. One of the dungeon halls is fill with the bright spheres. The explorers find out that the light rays reflect from the surface of the spheres according the ordinary law (the incidence angle is equal to the reflectance angle, the incidence ray, the reflected ray and the perpendicular to the sphere surface lay in the one plane). The ancient legend says that if the light ray will reflect from the spheres in the proper order, than the door to the room with very precious ancient knowledge will open. You are not to guess the right sequence; your task is much simpler. You are given the positions and the radii of the spheres, the place where the laser shot was made and the direction of light propagation. And you must find out the sequence in which the light will be reflected from the spheres.


Input

The first line of input contains the single integer n (1≤n≤50) - the amount of the spheres. The next n lines contain the coordinates and the radii of the spheres xi, yi, zi, ri (the integer numbers less or equal to 10000 by absolute value). The last line contains 6 real numbers - the coordinates of two points. The first one gives the coordinates of the place of laser shot, and the second gives the direction in which it was made (the second point is the point on the ray). The starting point of the ray lies strictly outside of any sphere.


Output

Your program must output the sequence of sphere numbers (spheres are numbers from 1 as they was given in input), from which the light ray was reflected. If the ray will reflect more the 10 times, than you must output first 10, then a space and the word 'etc.' (without quotes). Notice: if the light ray goes at a tangent to the sphere you must assume that the ray was reflected by the sphere.


Sample Input 1

1
0 0 2 1
0 0 0 0 0 1
Sample Output 1

1
Sample Input 2

2
0 0 2 1
0 0 -2 1
0 0 0 0 0 100
Sample Output 2

1 2 1 2 1 2 1 2 1 2 etc.

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

You are given up to 50 spheres in 3D space, each with center \((x_i, y_i, z_i)\) and radius \(r_i\). A laser ray starts at point \(P_0\) (outside all spheres) and initially goes in the direction defined by another point \(P_1\) on the ray.

The ray:

- Moves along straight segments.
- Reflects on a sphere using the usual optics law: incidence angle = reflection angle; reflection about the surface normal at the hit point.
- Counts as a hit even if it only touches a sphere tangentially.
- May hit the same sphere multiple times over time.

Simulate the reflections:

- Output the sequence of indices (1-based) of spheres hit, in order.
- If there are more than 10 hits, output only the first 10 indices, followed by `etc.`.


2. Key observations
-------------------

1. **Ray representation**

   A ray from origin \(O\) in direction \(D\) can be written as:
   \[
   P(t) = O + tD,\quad t \ge 0
   \]
   It is convenient and numerically stable to keep \(D\) normalized to unit length.

2. **Sphere equation**

   A sphere with center \(C\) and radius \(r\) is defined as:
   \[
   |P - C|^2 = r^2
   \]

3. **Ray–sphere intersection reduces to a quadratic in \(t\)**

   Substitute the ray into the sphere equation and solve:
   \[
   |O + tD - C|^2 = r^2
   \]
   This yields:
   \[
   a t^2 + b t + c = 0
   \]
   with:
   - \(V = O - C\)
   - \(a = D \cdot D\) (≈ 1 if \(D\) is normalized)
   - \(b = 2 (V \cdot D)\)
   - \(c = V \cdot V - r^2\)

   Discriminant:
   \[
   \Delta = b^2 - 4ac
   \]

   - If \(\Delta < 0\): no intersection.
   - If \(\Delta = 0\): tangent hit (one intersection).
   - If \(\Delta > 0\): two points; take the **smallest positive** \(t\).

   We use a small epsilon (e.g. `1e-6`) and accept hits with \(t > \text{EPS}\).

4. **Choosing the next sphere**

   At any step, given the current ray (origin, direction), the next hit is the sphere (other than the one just hit) with the smallest positive intersection parameter \(t\). If none has such \(t\), the ray escapes, and simulation stops.

   Skipping the last-hit sphere avoids immediately intersecting the same sphere again due to numerical errors (we are starting exactly on the surface).

5. **Reflection direction**

   Let:

   - \(D\) be the **incoming** (normalized) direction.
   - \(P\) the hit point.
   - \(C\) the center of the sphere.
   - Normal at hit point \(N = \frac{P - C}{\|P - C\|}\) (unit vector).

   Reflection formula:
   \[
   R = D - 2(D \cdot N) N
   \]
   Then normalize `R` again.

6. **Simulation limit**

   - We only need to *print* up to 10 hits.
   - To know if there are "more than 10", simulate up to 11 hits:
     - If we get >10 hits, print first 10 and `etc.`.
     - Otherwise print whatever hits we obtained.

7. **Complexity**

   - At each step, check all up to 50 spheres: \(O(n)\).
   - At most 11 steps → \(O(11n) = O(n)\) total.
   - Very small, easily within limits.

8. **Numeric robustness**

   - Use `double` precision.
   - Epsilon (`EPS = 1e-6`) for:
     - ignoring very small or negative `t` near zero,
     - treating slightly negative discriminants (within noise) as 0.
   - Normalize direction vectors consistently.


3. Full solution approach
-------------------------

1. **Parse input**

   - Read `n`.
   - For each sphere: read `xi, yi, zi, ri`.
   - Read 6 real numbers:
     - `x1, y1, z1` — ray starting point \(P_0\).
     - `x2, y2, z2` — another point on the ray.

2. **Initial ray setup**

   - Set ray origin \(O = (x1, y1, z1)\).
   - Set direction:
     \[
     D = (x2 - x1, y2 - y1, z2 - z1)
     \]
     then normalize \(D\) to unit length.

3. **Core functions**

   a) **Intersection function** `ray_sphere_intersect(O, D, sphere)`:

   - Let `V = O - C`.
   - Compute `a`, `b`, `c` as above.
   - Compute `disc = b*b - 4*a*c`.
   - If `disc < -EPS`: return `-1` (no intersection).
   - Else:
     - Set `disc = sqrt(max(0.0, disc))`.
     - Compute `t1 = (-b - disc) / (2*a)`, `t2 = (-b + disc) / (2*a)`.
     - Return the smallest of `t1`, `t2` that exceeds EPS, else `-1`.

   b) **Reflection function** `reflect(D, N)`:

   - Compute: `R = D - 2 * (D·N) * N`.

4. **Simulation loop**

   - Prepare a vector `hits` storing hit sphere indices (1-based) and `last_hit = -1`.
   - For up to 11 steps:
     1. Let `best_t = +∞`, `best_idx = -1`.
     2. For each sphere `i` other than `last_hit`, compute `t`; if `t > 0` and `t < best_t`, update.
     3. If `best_idx == -1`, break.
     4. Append `best_idx + 1` to `hits`, set `last_hit = best_idx`.
     5. Compute hit point `P = O + best_t * D`.
     6. Compute unit normal `N = (P - C_best).normalize()`.
     7. Reflect direction: `D = normalize(D - 2(D·N)N)`.
     8. Set new origin `O = P`.

5. **Output**

   - Print up to the first 10 hits separated by spaces; if there were more than 10, append `etc.`; end with newline.


4. C++ implementation
---------------------

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

const double EPS = 1e-6;

struct Vec3 {
    double x, y, z;
    Vec3(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {}
    Vec3 operator+(const Vec3& o) const {
        return Vec3(x + o.x, y + o.y, z + o.z);
    }
    Vec3 operator-(const Vec3& o) const {
        return Vec3(x - o.x, y - o.y, z - o.z);
    }
    Vec3 operator*(double k) const { return Vec3(x * k, y * k, z * k); }
    double dot(const Vec3& o) const { return x * o.x + y * o.y + z * o.z; }
    double norm() const { return sqrt(dot(*this)); }
    Vec3 normalize() const {
        double n = norm();
        return Vec3(x / n, y / n, z / n);
    }
};

struct Sphere {
    Vec3 center;
    double r;
};

int n;
vector<Sphere> spheres;
Vec3 origin, dir;

double ray_sphere_intersect(Vec3 O, Vec3 D, const Sphere& s) {
    // Find parameter t where ray P(t) = O + t*D intersects sphere:
    //
    //     |P - C|^2 = r^2
    //     |O + t*D - C|^2 = r^2
    //
    // Let V = O - C, then |V + t*D|^2 = r^2.
    //
    //     (V + t*D) . (V + t*D) = r^2
    //     V . V + 2t(V . D) + t^2(D . D) = r^2
    //
    // This gives quadratic equation:
    //   a*t^2 + b*t + c = 0
    //   where a = D . D, b = 2(V . D), c = V . V - r^2
    //
    // This function returns smallest positive t, or -1 if no valid
    // intersection.

    Vec3 V = O - s.center;
    double a = D.dot(D);
    double b = 2 * V.dot(D);
    double c = V.dot(V) - s.r * s.r;
    double disc = b * b - 4 * a * c;
    if(disc < -EPS) {
        return -1;
    }
    disc = sqrt(max(0.0, disc));
    double t1 = (-b - disc) / (2 * a);
    double t2 = (-b + disc) / (2 * a);
    if(t1 > EPS) {
        return t1;
    }
    if(t2 > EPS) {
        return t2;
    }
    return -1;
}

Vec3 reflect(Vec3 D, Vec3 N) {
    // Reflect incident direction D about normal N.
    // Decompose D into parallel and perpendicular components to N:
    //   D_parallel = (D . N)N, D_perp = D - D_parallel
    // Reflection reverses the parallel component:
    //   R = D_perp - D_parallel = D - 2*D_parallel = D - 2(D . N)N

    return D - N * (2 * D.dot(N));
}

void read() {
    cin >> n;
    spheres.resize(n);
    for(int i = 0; i < n; i++) {
        cin >> spheres[i].center.x >> spheres[i].center.y >>
            spheres[i].center.z >> spheres[i].r;
    }
    double x1, y1, z1, x2, y2, z2;
    cin >> x1 >> y1 >> z1 >> x2 >> y2 >> z2;
    origin = Vec3(x1, y1, z1);
    dir = Vec3(x2 - x1, y2 - y1, z2 - z1).normalize();
}

void solve() {
    // We can solve this problem in a straight forward way - we just need to
    // implement some 3D geometry functions. In particular, we want to be able
    // to get intersections, tangents and reflection vectors, all of which have
    // closed forms. Then we can simulate the first 10 hits.

    vector<int> hits;
    int last_hit = -1;
    for(int step = 0; step < 11; step++) {
        double best_t = 1e18;
        int best_sphere = -1;
        for(int i = 0; i < n; i++) {
            if(i == last_hit) {
                continue;
            }
            double t = ray_sphere_intersect(origin, dir, spheres[i]);
            if(t > 0 && t < best_t) {
                best_t = t;
                best_sphere = i;
            }
        }
        if(best_sphere == -1) {
            break;
        }

        hits.push_back(best_sphere + 1);
        last_hit = best_sphere;

        Vec3 P = origin + dir * best_t;
        Vec3 N = (P - spheres[best_sphere].center).normalize();
        dir = reflect(dir, N).normalize();
        origin = P;
    }

    int output_count = min((int)hits.size(), 10);
    for(int i = 0; i < output_count; i++) {
        if(i > 0) {
            cout << " ";
        }
        cout << hits[i];
    }
    if(hits.size() > 10) {
        cout << " etc.";
    }
    cout << "\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();
        // cout << "Case #" << test << ": ";
        solve();
    }

    return 0;
}
```

5. Python implementation with detailed comments
-----------------------------------------------

```python
import sys
import math

# Small epsilon for floating-point comparisons
EPS = 1e-6


class Vec3:
    """
    Simple 3D vector class with basic operations:
    addition, subtraction, scalar multiplication, dot, norm, normalize.
    """
    __slots__ = ("x", "y", "z")

    def __init__(self, x=0.0, y=0.0, z=0.0):
        self.x = float(x)
        self.y = float(y)
        self.z = float(z)

    def __add__(self, other):
        """Vector addition: self + other"""
        return Vec3(self.x + other.x, self.y + other.y, self.z + other.z)

    def __sub__(self, other):
        """Vector subtraction: self - other"""
        return Vec3(self.x - other.x, self.y - other.y, self.z - other.z)

    def __mul__(self, k):
        """Scalar multiplication: self * k"""
        return Vec3(self.x * k, self.y * k, self.z * k)

    # Allow k * self as well
    __rmul__ = __mul__

    def dot(self, other):
        """Dot product with another Vec3."""
        return self.x * other.x + self.y * other.y + self.z * other.z

    def norm(self):
        """Euclidean norm (length) of the vector."""
        return math.sqrt(self.dot(self))

    def normalize(self):
        """Return a normalized (unit-length) copy of this vector."""
        n = self.norm()
        # The problem guarantees the direction will not be zero-length.
        return Vec3(self.x / n, self.y / n, self.z / n)


class Sphere:
    """Sphere with center (Vec3) and radius r."""
    __slots__ = ("center", "r")

    def __init__(self, center, r):
        self.center = center
        self.r = float(r)


def ray_sphere_intersect(O, D, sphere):
    """
    Compute intersection of ray P(t) = O + t*D with given sphere.
    Returns the smallest t > EPS if intersection exists; otherwise -1.0.

    O: Vec3, ray origin
    D: Vec3, ray direction (unit or not, doesn't matter)
    sphere: Sphere
    """
    V = O - sphere.center

    a = D.dot(D)                # Usually 1.0 if D is normalized
    b = 2.0 * V.dot(D)
    c = V.dot(V) - sphere.r * sphere.r

    # Discriminant of quadratic a t^2 + b t + c = 0
    disc = b * b - 4.0 * a * c

    # If disc < 0 by more than EPS, no intersection
    if disc < -EPS:
        return -1.0

    # Clamp negative disc to 0 for numerical safety, then sqrt
    disc = math.sqrt(max(0.0, disc))

    t1 = (-b - disc) / (2.0 * a)
    t2 = (-b + disc) / (2.0 * a)

    # We want the smallest t that is certainly in front (t > EPS)
    if t1 > EPS:
        return t1
    if t2 > EPS:
        return t2

    return -1.0


def reflect(D, N):
    """
    Reflect direction D about normal N.
    Assumes N is a unit vector.
    Reflection formula: R = D - 2(D·N)N
    """
    k = 2.0 * D.dot(N)
    return D - k * N


def main():
    data = sys.stdin.read().strip().split()
    if not data:
        return

    it = iter(data)

    try:
        n = int(next(it))
    except StopIteration:
        return

    spheres = []
    for _ in range(n):
        x = float(next(it))
        y = float(next(it))
        z = float(next(it))
        r = float(next(it))
        spheres.append(Sphere(Vec3(x, y, z), r))

    # Two points defining the initial ray
    x1 = float(next(it))
    y1 = float(next(it))
    z1 = float(next(it))
    x2 = float(next(it))
    y2 = float(next(it))
    z2 = float(next(it))

    # Initial ray origin
    origin = Vec3(x1, y1, z1)

    # Initial direction: from (x1,y1,z1) to (x2,y2,z2), then normalized
    dir_vec = Vec3(x2 - x1, y2 - y1, z2 - z1).normalize()

    hits = []       # 1-based sphere indices
    last_hit = -1   # 0-based index of last-hit sphere, or -1 if none

    # Simulate up to 11 reflections (one extra used to detect >10 hits)
    for step in range(11):
        best_t = float('inf')
        best_idx = -1

        # Check intersection with all spheres
        for i, s in enumerate(spheres):
            # Avoid immediately re-hitting the same sphere due to precision
            if i == last_hit:
                continue

            t = ray_sphere_intersect(origin, dir_vec, s)

            if t > 0.0 and t < best_t:
                best_t = t
                best_idx = i

        # If no sphere is hit, break
        if best_idx == -1:
            break

        hits.append(best_idx + 1)    # store 1-based index
        last_hit = best_idx          # update last-hit sphere

        # Compute hit point
        P = origin + dir_vec * best_t

        # Normal at hit point (outward from center)
        N = (P - spheres[best_idx].center).normalize()

        # Reflect direction and normalize
        dir_vec = reflect(dir_vec, N).normalize()

        # Move origin to the hit point
        origin = P

    # Output: up to first 10 hits
    out_hits = hits[:10]
    if out_hits:
        sys.stdout.write(" ".join(str(x) for x in out_hits))

    # If more than 10 hits, append " etc."
    if len(hits) > 10:
        if out_hits:
            sys.stdout.write(" ")
        sys.stdout.write("etc.")

    sys.stdout.write("\n")


if __name__ == "__main__":
    main()
```
