sgu 110 Dungeon

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <iostream>
#include <algorithm>
#include <cmath>
#define N 55
#define inf 0x7f7f7f7f
using namespace std;

struct Point3
{
    double x, y, z;
    Point3(double x=0, double y=0, double z=0):x(x), y(y), z(z) { }
};
typedef Point3 Vector3;
const double eps = 1e-9;

int dcmp(double x)
{
    if (fabs(x) < eps) return 0;
    else return x < 0 ? -1 : 1;
}

int n;
double r[N];
Point3 Begin, circle[N];
Vector3 Direct;

Vector3 operator + (Vector3 A, Vector3 B) { return Vector3(A.x+B.x, A.y+B.y, A.z+B.z); }
Vector3 operator - (Point3 A, Point3 B) { return Vector3(A.x-B.x, A.y-B.y, A.z-B.z); }
Vector3 operator * (Vector3 A, double p) { return Vector3(A.x*p, A.y*p, A.z*p); }
Vector3 operator / (Vector3 A, double p) { return Vector3(A.x/p, A.y/p, A.z/p); }

double min(double x, double y) { return x < y ? x : y; }

double Dot (Vector3 A, Vector3 B) { return A.x*B.x + A.y*B.y + A.z*B.z; }
double Length (Vector3 A) { return sqrt(Dot(A, A)); }
double Angle (Vector3 A, Vector3 B) { return acos(Dot(A, B) / Length(A) / Length(B)); }
Vector3 Cross (Vector3 A, Vector3 B) { return Vector3(A.y*B.z - A.z*B.y, A.z*B.x - A.x*B.z, A.y*B.z - A.z*B.y); }
double disPointLine(Point3 P, Point3 A, Point3 B)
{
    Vector3 v1, v2;
    v1 = B - A; v2 = P - A;
    return Length(Cross(v1, v2)) / Length(v1);
}

int main()
{
    scanf("%d", &n);
    double x, y, z;
    for (int i = 1; i <= n; ++i)
    {
        scanf("%lf%lf%lf%lf", &x, &y, &z, &r[i]);
        circle[i] = Point3(x, y, z);
    }
    scanf("%lf%lf%lf", &x, &y, &z); Begin = Point3(x, y ,z);
    scanf("%lf%lf%lf", &x, &y, &z); Direct = Point3(x, y, z) - Begin;
    int nowcircle = 0;
    for (int w = 1; w <= 11; ++w)
    {
        int nextcircle = 0;
        double mindis = inf;
        for (int i = 1; i <= n; ++i)
        {
            if (i == nowcircle) continue;
            double a, b, c;
            Point3 nc = circle[i];
            a = Direct.x*Direct.x + Direct.y*Direct.y + Direct.z*Direct.z;
            b = 2 * ((Begin.x-nc.x) * Direct.x + (Begin.y-nc.y)*Direct.y + (Begin.z-nc.z)*Direct.z);
            c = (Begin.x-nc.x)*(Begin.x-nc.x) + (Begin.y-nc.y)*(Begin.y-nc.y) + (Begin.z-nc.z)*(Begin.z-nc.z) - r[i]*r[i];
            double delta;
            delta = b*b - 4*a*c;
            if (delta < 0) continue;
            double ans1, ans2;
            ans1 = (-b + sqrt(delta)) / (2*a);
            ans2 = (-b - sqrt(delta)) / (2*a);
            if (dcmp(ans1) < 0 && dcmp(ans2) < 0) continue;
            else if (dcmp(ans1) < 0 && ans2 < mindis)
            {
                mindis = ans2;
                nextcircle = i;
            }
            else if (dcmp(ans2) < 0 && ans1 < mindis)
            {
                mindis = ans1;
                nextcircle = i;
            }
            else if (min(ans1, ans2) < mindis)
            {
                mindis = min(ans1, ans2);
                nextcircle = i;
            }
        }
        if (!nextcircle) break;
        else if (w == 11)
        {
            printf(" etc.");
            break;
        }
        else
        {
            if (w == 1) printf("%d", nextcircle);
            else printf(" %d", nextcircle);
            Point3 jiao;
            Vector3 v1, v2, fa;
            jiao = Begin + Direct*mindis; v1 = Direct;
            fa = circle[nextcircle] - jiao;
            v2 = (fa / Length(fa)) * (2 * Dot(v1, fa) / Length(fa));
            nowcircle = nextcircle;
            Begin = jiao; Direct = v1 - v2;
        }
    }
    printf("\n");
}