题目分析：本题重点就在于对给定入射射线AB计算反射射线BC。
对于空间上的入射射线的起点A，射线与球的交点B，球的圆心C。
因为向量AB总可以表示成|AB|*e，其中e为向量AB的方向向量，|AB|为向量AB的模，所以|AO -AB| = |BO| = R
将A，O，R带进去可以得到一个关于|AB|的一元二次方程，解方程求得x1，x2，取其中正的较小的一个。（可能存在无解，莫忘判断）
得到|AB|后，AB = e*|AB|，同时B = AB + Ａ，BO = O - B，此时BC = AB - AB·BO / |BO| * 2 * (BO / |BO| )，做一个图就可以发现这个公式。
我们就是利用上面的公式，然后将求得的满足|AB|长度最小的交点B，以及相交的球的球心O带入来求得反射射线BC（新的入射射线AB），反射射线的起点B（新的入射射线的起点A）。
需要注意的是，可能起点在一个球上，这时候我们需要判断使得一条射线不能连续在同一个球上反射。
代码如下：
#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std ;
typedef long long LL ;
#define rep( i , a , b ) for ( int i = ( a ) ; i < ( b ) ; ++ i )
#define For( i , a , b ) for ( int i = ( a ) ; i <= ( b ) ; ++ i )
#define rev( i , a , b ) for ( int i = ( a ) ; i >= ( b ) ; -- i )
#define travel( e , H , u ) for ( Edge* e = H[u] ; e ; e = e -> next )
#define clr( a , x ) memset ( a , x , sizeof a )
const int MAXN = 55 ;
const double INF = 1e10 ;
const double eps = 1e-5 ;
struct Point {
double x , y , z ;
Point () {}
Point ( double x , double y , double z ) : x ( x ) , y ( y ) , z ( z ) {}
Point operator + ( const Point& p ) const {
return Point ( x + p.x , y + p.y , z + p.z ) ;
}
Point operator - ( const Point& p ) const {
return Point ( x - p.x , y - p.y , z - p.z ) ;
}
Point operator * ( const double& b ) const {
return Point ( x * b , y * b , z * b ) ;
}
Point operator / ( const double& b ) const {
return Point ( x / b , y / b , z / b ) ;
}
void input () {
scanf ( "%lf%lf%lf" , &x , &y , &z ) ;
}
void show () {
printf ( "( %.5f , %.5f , %.5f )\n" , x , y , z ) ;
}
} ;
typedef Point Vector ;
struct Sphere {
Point p ;
double r ;
void input () {
p.input () ;
scanf ( "%lf" , &r ) ;
}
} ;
Sphere ball[MAXN] ;
int n ;
int dcmp ( double x ) {
return ( x > eps ) - ( x < -eps ) ;
}
double Length ( const Vector& p ) {
return sqrt ( p.x * p.x + p.y * p.y + p.z * p.z ) ;
}
double Dot ( const Vector& a , const Vector& b ) {
return a.x * b.x + a.y * b.y + a.z * b.z ;
}
void solve () {
Point pos ;
Vector dir ;
For ( i , 1 , n ) ball[i].input () ;
pos.input () ;
dir.input () ;
dir = dir - pos ;
int now = 0 ;
For ( i , 1 , 11 ) {
int idx = 0 ;
double len = INF ;
Point A , B , O ;
A = pos ;
For ( j , 1 , n ) {
if ( now == j ) continue ;
Vector AO = ball[j].p - A ;
double a = dir.x * dir.x + dir.y * dir.y + dir.z * dir.z ;
double b = - 2 * ( AO.x * dir.x + AO.y * dir.y + AO.z * dir.z ) ;
double c = AO.x * AO.x + AO.y * AO.y + AO.z * AO.z - ball[j].r * ball[j].r ;
double delta = b * b - 4 * a * c ;
if ( dcmp ( delta ) < 0 ) continue ;
double x1 = ( - b - sqrt ( delta ) ) / ( 2 * a ) ;
double x2 = ( - b + sqrt ( delta ) ) / ( 2 * a ) ;
//printf ( "%.5f %.5f\n" , x1 , x2 ) ;
double x = INF ;
if ( dcmp ( x1 ) > 0 ) x = min ( x , x1 ) ;
if ( dcmp ( x2 ) > 0 ) x = min ( x , x2 ) ;
if ( dcmp ( x - INF ) == 0 ) continue ;
else if ( dcmp ( len - x ) > 0 ) {
O = ball[j].p ;
idx = j ;
len = x ;
}
}
if ( !idx ) break ;
if ( i == 11 ) {
printf ( " etc." ) ;
break ;
}
B = dir * len + pos ;
Vector BO = O - B ;
//dir.show () ;
dir = dir - BO * ( Dot ( dir , BO ) / pow ( Length ( BO ) , 2.0 ) * 2.0 ) ;
//dir.show () ;
pos = B ;
if ( i == 1 ) printf ( "%d" , idx ) ;
else printf ( " %d" , idx ) ;
now = idx ;
}
printf ( "\n" ) ;
}
int main () {
while ( ~scanf ( "%d" , &n ) ) solve () ;
return 0 ;
}