这道题是计算几何，这是写的第一道计算几何，主要是难在如何求入射光线的反射光线。
我们可以用入射光线 - 入射光线在法线（交点到圆心的向量）上的投影*2 来计算反射光线，自己画一个图，非常清晰明了。
具体到程序里，我们可以 v2 = v1 - fa / Length(fa) * 2 * ( Dot(v1, fa) / Length(fa)) 来求，简单来说就是用内积（点积）求出入射光线在法线上的长度，然后用法线的单位向量乘这个长度就可以了。
#include#include #include #include #include #include #define N 55#define inf 0x7f7f7f7fusing 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");}