题目链接
本题是2002年ICPC欧洲区域赛东北欧赛区的G题
题意
有 N ( N ≤ 100 ) N(N≤100) N(N≤100)个鬼，每个鬼是中心在 ( X i , Y i , Z i ) ( 1 ≤ X i , Y i , Z i ≤ 10000 ) (X_i,Y_i,Z_i) (1 ≤ X_i,Y_i,Z_i ≤ 10000) (Xi,Yi,Zi)(1≤Xi,Yi,Zi≤10000)半径为 R i ( 1 ≤ R i ≤ m i n ( X i , Y i , Z i ) ) R_i (1 ≤ R_i ≤ min(X_i, Y_i, Z_i)) Ri(1≤Ri≤min(Xi,Yi,Zi))的球。鬼是虚假的，鬼与鬼可以任意相交甚至一个鬼完全另外一个鬼内部。在原点 ( 0 , 0 , 0 ) (0,0,0) (0,0,0)发出一根射线，射线和鬼又接触（相切或者穿过）则能消灭鬼，问最多可以消灭多少个鬼并输出相应的鬼编号(编号从1开始)，有多解时输出任意一组鬼的编号即可。
分析
先分析一下射线何时能消灭一个鬼？答案是射线在原点与鬼外切形成的圆锥内时能消灭鬼。那么要消灭最多的鬼就是要找出这些圆锥相交部分涉及鬼数量最多的情况，但这是一个空间问题，似乎无法直观高效地处理，要想办法转化成平面问题：用一个平面去截这些圆锥就能转化成平面问题了！我们知道圆锥的截面是圆锥曲线，但是圆锥曲线中抛物线和双曲线是不闭合的，只有椭圆（或者退化成圆）才是闭合的。
那么要找一种平面使其和所有圆锥相截都能得到椭圆（或者圆），才能转化成能直观处理的平面问题，幸运的是，这种平面能够找到：法向量为 ( − 1 , − 1 , − 1 ) (-1,-1,-1) (−1,−1,−1)的平面就可行！
根据题目的数值范围，编程实践时推荐选择经过点 ( 10 3 , 10 3 , 10 3 ) (\frac{10}{\sqrt3},\frac{10}{\sqrt3},\frac{10}{\sqrt3}) (310,310,310)的平面（法向量为 ( − 1 , − 1 , − 1 ) (-1,-1,-1) (−1,−1,−1)）此时原点到平面的距离为10。
接下来就是编程实现了，首先要进行坐标转换，将经过点 ( 10 3 , 10 3 , 10 3 ) (\frac{10}{\sqrt3},\frac{10}{\sqrt3},\frac{10}{\sqrt3}) (310,310,310)且法向量为 ( − 1 , − 1 , − 1 ) (-1,-1,-1) (−1,−1,−1)）的面转换成点为 ( 0 , 0 , 0 ) (0,0,0) (0,0,0)，法向量为 ( 0 , 0 , ( 3 ) ) (0,0,\sqrt(3)) (0,0,(3))这个参考UVA-1488的处理方法即可，贴一下英文原文：
Hint: Here are some hints to help you implement rotation in 3-dimensional space.
- Assume that there is a plane a x + b y + c z = 0 ax + by + cz = 0 ax+by+cz=0, and our task is to rotate the whole space so that the plane will be in the position z = 0 z = 0 z=0.
As we all know, the normal vector of the plane is ( a , b , c ) (a, b, c) (a,b,c), so after we rotate the space, the normal vector will be the z-axis. So we can only consider how to rotate the vector to the position
( 0 , 0 , a 2 + b 2 + c 2 ) (0, 0, \sqrt {a^2+b^2+c^2}) (0,0,a2+b2+c2).- We can divide the process into 2 steps. First step, rotate the normal vector to ( 0 , a 2 + b 2 , c ) (0, \sqrt {a^2+b^2}, c) (0,a2+b2,c).
Second step, rotate ( 0 , a 2 + b 2 , c ) (0, \sqrt {a^2+b^2}, c) (0,a2+b2,c) to ( 0 , 0 , a 2 + b 2 + c 2 ) (0, 0, \sqrt {a^2 + b^2 + c^2}) (0,0,a2+b2+c2). Then apply the two rotation steps on
the whole space, we can rotate the plane to a given position.
It is easy to implement the two steps. In first step, the third number in the normal vector does
not change, so this is like rotation in a 2-dimensional space. Similarly, we implement the second
step.
然后要求出每个球被射线外切形成的圆锥被截面截得的椭圆参数：长短半轴、中心点坐标和倾斜角。首先容易得到每个椭圆的底角 α \alpha α， cos α = R i X i 2 + Y i 2 + Z i 2 \cos\alpha=\frac{R_i}{\sqrt {X_i^2 + Y_i^2 + Z_i^2}} cosα=Xi2+Yi2+Zi2Ri，转换后的坐标是 ( x i , y i , 10 ) (x_i,y_i,10) (xi,yi,10),则圆锥的高度 h = x i 2 + y i 2 + 100 h=\sqrt{x_i^2+y_i^2+100} h=xi2+yi2+100,因此截面与圆锥底面的夹角 β \beta β也知道， cos β = 10 h \cos\beta=\frac{10}{h} cosβ=h10,圆锥底面半径 r = h cot α r=h\cot\alpha r=hcotα。首先椭圆倾斜角 t i l t = a t a n 2 ( y i , x i ) tilt=atan2(y_i, x_i) tilt=atan2(yi,xi)，再基于 α \alpha α、 β \beta β和 h h h就可以求出椭圆的长半轴 a a a、短半轴 b b b、中心点 ( c x i , c y i ) (cx_i,cy_i) (cxi,cyi)。
在 △ A O 1 B \triangle AO_1B △AO1B和 △ A O 1 C \triangle AO_1C △AO1C中利用正弦定理得 h sin ∠ A B O 1 = B O 1 sin ( 9 0 ° − α ) , h sin ∠ A C O 1 = C O 1 sin ( 9 0 ° − α ) \frac{h}{\sin \angle ABO_1}=\frac {BO_1}{\sin(90^\degree-\alpha)},\frac{h}{\sin \angle ACO_1}=\frac{CO_1}{\sin(90^\degree-\alpha)} sin∠ABO1h=sin(90°−α)BO1,sin∠ACO1h=sin(90°−α)CO1而 ∠ A C O 1 = α − β , ∠ A B O 1 = 180 ° − ∠ A C O 1 − ∠ B A C = 180 ° − ( α − β ) − ( 180 ° − 2 α ) = α + β \angle ACO_1=\alpha-\beta,\angle ABO_1=180\degree-\angle ACO_1 - \angle BAC = 180\degree - (\alpha-\beta) - (180\degree-2\alpha)=\alpha+\beta ∠ACO1=α−β,∠ABO1=180°−∠ACO1−∠BAC=180°−(α−β)−(180°−2α)=α+β,故 B O 1 = h cos α sin ( α + β ) , C O 1 = h cos α sin ( α − β ) BO_1=\frac{h\cos\alpha}{\sin(\alpha+\beta)},CO_1=\frac{h\cos\alpha}{\sin(\alpha-\beta)} BO1=sin(α+β)hcosα,CO1=sin(α−β)hcosα因此 a = B O 1 + C O 1 2 = h cos α 2 sin ( α + β ) + h cos α 2 sin ( α − β ) a=\frac{BO_1+CO_1}{2}=\frac{h\cos\alpha}{2\sin(\alpha+\beta)}+\frac{h\cos\alpha}{2\sin(\alpha-\beta)} a=2BO1+CO1=2sin(α+β)hcosα+2sin(α−β)hcosα。 R = [ h + ( a − B O 1 ) sin β ] cot α , D = ( a − B O 1 ) cos β R=[h+(a-BO_1)\sin\beta]\cot\alpha,D=(a-BO_1)\cos\beta R=[h+(a−BO1)sinβ]cotα,D=(a−BO1)cosβ,因此由 b = R 2 − D 2 b=\sqrt{R^2-D^2} b=R2−D2可将b求出，最后计算出 k = x i 2 + y i 2 + a − B O 1 x i 2 + y i 2 k=\frac{\sqrt{x_i^2+y_i^2}+a-BO_1}{\sqrt{x_i^2+y_i^2}} k=xi2+yi2xi2+yi2+a−BO1，就能得到 c x i = k x i , c y i = k y i cx_i=kx_i,cy_i=ky_i cxi=kxi,cyi=kyi。
每个椭圆参数都求出来了，接下来找出涉及最多椭圆的相交部分就可以了，这需要对椭圆两两枚举求出所有交点，考虑到某些椭圆可能完全在其他椭圆中，还需要将椭圆的中心点也加入，然后对每个点计算包含此点数量最多的的椭圆就能得出最终答案了。
看来起来确实很直观，然而又一个难题出现了：两椭圆最多可能有4个交点，因此联立两椭圆方程求解的话涉及解四次方程，需要找更好的办法求椭圆的焦点。其实Wikipedia上已经给出了利用线性代数二次型（Quadratic form）的求解理论，详细的编程实践参见我的另外一篇博文：求两个二次曲线交点的理论依据和编程实践。
非常具有挑战性的一道好题，附上一份测试数据。
AC代码
#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;
#define eps 1e-12l
#define N 102
#define DBL long double
DBL h[N][6], px[2*N*N], py[2*N*N], x[4], y[4], q = 0.5773502691896257, s = 0.816496580927726; int g[N], p[N], c, n, t;
void cofactor(const DBL (&a)[3][3], DBL (&c)[3][3]) {
c[0][0] = a[1][1]*a[2][2] - a[1][2]*a[2][1];
c[0][1] = a[1][2]*a[2][0] - a[1][0]*a[2][2];
c[0][2] = a[1][0]*a[2][1] - a[1][1]*a[2][0];
c[1][0] = a[0][2]*a[2][1] - a[0][1]*a[2][2];
c[1][1] = a[0][0]*a[2][2] - a[0][2]*a[2][0];
c[1][2] = a[0][1]*a[2][0] - a[0][0]*a[2][1];
c[2][0] = a[0][1]*a[1][2] - a[0][2]*a[1][1];
c[2][1] = a[0][2]*a[1][0] - a[0][0]*a[1][2];
c[2][2] = a[0][0]*a[1][1] - a[0][1]*a[1][0];
}
void mul(const DBL (&a)[3][3], const DBL (&b)[3][3], DBL (&c)[3][3]) {
c[0][0] = a[0][0]*b[0][0] + a[0][1]*b[1][0] + a[0][2]*b[2][0];
c[0][1] = a[0][0]*b[0][1] + a[0][1]*b[1][1] + a[0][2]*b[2][1];
c[0][2] = a[0][0]*b[0][2] + a[0][1]*b[1][2] + a[0][2]*b[2][2];
c[1][0] = a[1][0]*b[0][0] + a[1][1]*b[1][0] + a[1][2]*b[2][0];
c[1][1] = a[1][0]*b[0][1] + a[1][1]*b[1][1] + a[1][2]*b[2][1];
c[1][2] = a[1][0]*b[0][2] + a[1][1]*b[1][2] + a[1][2]*b[2][2];
c[2][0] = a[2][0]*b[0][0] + a[2][1]*b[1][0] + a[2][2]*b[2][0];
c[2][1] = a[2][0]*b[0][1] + a[2][1]*b[1][1] + a[2][2]*b[2][1];
c[2][2] = a[2][0]*b[0][2] + a[2][1]*b[1][2] + a[2][2]*b[2][2];
}
DBL det(const DBL (&a)[3][3], const DBL (&c)[3][3]) {
return a[0][0]*c[0][0] + a[0][1]*c[0][1] + a[0][2]*c[0][2];
}
DBL trace(const DBL (&a)[3][3], const DBL (&b)[3][3]) {
return a[0][0]*b[0][0] + a[0][1]*b[1][0] + a[0][2]*b[2][0] +
a[1][0]*b[0][1] + a[1][1]*b[1][1] + a[1][2]*b[2][1] + a[2][0]*b[0][2] + a[2][1]*b[1][2] + a[2][2]*b[2][2];
}
DBL cubic_root(DBL x) {
return x<0. ? -pow(-x, 1./3.) : pow(x, 1./3.);
}
DBL find_a_root(DBL a, DBL b, DBL c, DBL d) {
DBL A = b*b - 3*a*c, B = b*c - 9*a*d, C = c*c - 3*b*d, delta = B*B - 4.*A*C;
if (A == 0. && B == 0.) return -c/b;
if (delta > 0.) {
delta = sqrt(delta);
DBL y1 = A*b + 1.5*a*(delta-B), y2 = A*b - 1.5*a*(delta+B);
return -(cubic_root(y1) + cubic_root(y2) + b)/3./a;
} else if (delta > -eps) return -B/2./A;
DBL sA = sqrt(A), t = acos((A*b-1.5*a*B)/A/sA)/3.;
return -(b + 2*sA*cos(t)) / 3. / a;
}
int solve_quadratic_equation(DBL a, DBL b, DBL c, DBL *x) {
if (a == 0.) {
if (b == 0.) return 0;
*x = -c/b;
return 1;
}
DBL delta = b*b - 4.*a*c;
if (delta < 0.) {
if (delta < -eps) return 0;
delta = 0.;
}
a *= 2; delta = sqrt(delta);
if (delta == 0. || delta < a*eps) {
*x = -b/a;
return 1;
}
*x = (-b-delta)/a; *(++x) = (delta-b)/a;
return 2;
}
bool check(const DBL (&q)[3][3], DBL x, DBL y) {
return abs(q[0][0]*x*x + 2.*q[0][1]*x*y + q[1][1]*y*y + 2.*q[0][2]*x + 2.*q[1][2]*y + q[2][2]) < eps;
}
int line_conic_intersection(const DBL (&q)[3][3], DBL a, DBL b, DBL c, DBL *x, DBL *y) {
if (a == 0. && b == 0.) return 0;
if (abs(a) > abs(b)) {
DBL d = q[0][0]*b*b - 2.*q[0][1]*a*b + q[1][1]*a*a, f = q[0][0]*c*c - 2.*q[0][2]*a*c + q[2][2]*a*a;
int m = solve_quadratic_equation(d, 2.*(q[1][2]*a*a + q[0][0]*b*c - q[0][1]*a*c - q[0][2]*a*b), f, y);
for (int i=0; i<m; ++i) *(x+i) = (-c - b * (*(y+i)))/a;
return m;
}
DBL d = q[0][0]*b*b - 2.*q[0][1]*a*b + q[1][1]*a*a, f = q[1][1]*c*c - 2.*q[1][2]*b*c + q[2][2]*b*b;
int m = solve_quadratic_equation(d, 2.*(q[0][2]*b*b + q[1][1]*a*c - q[0][1]*b*c - q[1][2]*a*b), f, x);
for (int i=0; i<m; ++i) *(y+i) = (-c - a * (*(x+i)))/b;
return m;
}
void check(const DBL (&q)[3][3], DBL (&x)[4], DBL (&y)[4], int m, int &n) {
for (int i=n-1; i>=0; --i) {
if (!check(q, x[m+i], y[m+i])) {
if (i < n-1) x[m] = x[m+1], y[m] = y[m+1];
--n; continue;
}
for (int j=0; j<m; ++j) if (abs(x[j]-x[m+i])<eps && abs(y[j]-y[m+i])<eps) {
if (i < n-1) x[m] = x[m+1], y[m] = y[m+1];
--n; break;
}
}
}
int conics_intersection(const DBL (&p)[3][3], const DBL (&q)[3][3], DBL (&x)[4], DBL (&y)[4]) {
DBL c1[3][3], c2[3][3];
cofactor(p, c1); cofactor(q, c2);
DBL r = find_a_root(det(p, c1), trace(c1, q), trace(p, c2), det(q, c2));
DBL a = r*p[0][0] + q[0][0], b = r*p[0][1] + q[0][1], c = r*p[1][1] + q[1][1],
d = r*p[0][2] + q[0][2], e = r*p[1][2] + q[1][2], f = r*p[2][2] + q[2][2], a33 = a*c - b*b;
if (a33 > eps) {
x[0] = (b*e - c*d) / a33; y[0] = (b*d - a*e) / a33;
return check(p, x[0], y[0]) && check(q, x[0], y[0]);
} else if (a33 < -eps) {
a33 = sqrt(-a33); c = (b*d - a*e) / a33;
int m = line_conic_intersection(q, a, b-a33, d-c, x, y);
check(p, x, y, 0, m);
int n = line_conic_intersection(q, a, b+a33, d+c, x+m, y+m);
check(p, x, y, m, n);
return m+n;
} else if (abs(b) < eps) {
if (abs(a) > max(abs(c), eps)) {
int t = solve_quadratic_equation(a, 2.*d, f, x); double s = x[1];
if (!t) return 0;
int m = line_conic_intersection(q, -1., 0., x[0], x, y);
check(p, x, y, 0, m);
int n = t > 1 ? line_conic_intersection(q, -1., 0., s, x+m, y+m) : 0;
check(p, x, y, m, n);
return m+n;
} else if (abs(c) > max(abs(a), eps)) {
int t = solve_quadratic_equation(c, 2.*e, f, y); double s = y[1];
if (!t) return 0;
int m = line_conic_intersection(q, 0., -1., y[0], x, y);
check(p, x, y, 0, m);
int n = t > 1 ? line_conic_intersection(q, 0., -1., s, x+m, y+m) : 0;
check(p, x, y, m, n);
return m+n;
}
int m = line_conic_intersection(q, 2.*d, 2.*e, f, x, y);
check(p, x, y, 0, m);
return m;
}
if (a < 0) a = -a, b = -b, c = -c, d = -d, e = -e, f = -f;
DBL s = d*d - a*f;
if (b*d*e < -eps || abs(a*e*e - c*d*d) > eps || s < -eps) return 0;
if (s < eps) {
int m = line_conic_intersection(q, a, b, d, x, y);
check(p, x, y, 0, m);
return m;
}
s = sqrt(s);
int m = line_conic_intersection(q, a, b, d-s, x, y);
check(p, x, y, 0, m);
int n = line_conic_intersection(q, a, b, d+s, x+m, y+m);
check(p, x, y, m, n);
return m+n;
}
bool check(const DBL (&h)[6], DBL x, DBL y) {
return h[0]*x*x + h[1]*x*y + h[2]*y*y + h[3]*x + h[4]*y + h[5] < eps;
}
void solve() {
cin >> n; t = 0;
for (int i=1; i<=n; ++i) {
DBL x, y, z; int r; cin >> x >> y >> z >> r;
DBL v = x; x = (v-y)*M_SQRT1_2; v = (v+y)*M_SQRT1_2; y = q*v - s*z; z = s*v + q*z;
DBL co = r / sqrt(x*x+y*y+z*z), si = sqrt(1. - co*co);
if (z != 10.) v = 10. / z, x *= v, y *= v, v = sqrt(x*x + y*y);
DBL l = sqrt(v*v + 100.), sinb = v / l, cosb = sqrt(1 - sinb*sinb);
DBL c = 1./(si*cosb+sinb*co), d = 1./(si*cosb-sinb*co);
double a = .5*l*co*(c+d), b = a/si/sqrt(c*d);
if (v > 0.) d = (a-l*co*c+v)/v, x *= d, y *= d;
l = sqrt(x*x + y*y); si = y/l; co = x/l; px[t] = x, py[t++] = y;
h[i][0] = si/b*si/b + co/a*co/a; h[i][1] = 2.*(1./a/a - 1./b/b)*si*co;
h[i][2] = co/b*co/b + si/a*si/a; h[i][3] = -2.*h[i][0]*x - h[i][1]*y;
h[i][4] = -2.*h[i][2]*y - h[i][1]*x; h[i][5] = -(h[i][3]*x + h[i][4]*y)/2. - 1.;
}
for (int i=1; i<n; ++i) {
DBL p[][3] = {{h[i][0], h[i][1]/2., h[i][3]/2.},
{h[i][1]/2., h[i][2], h[i][4]/2.}, {h[i][3]/2., h[i][4]/2., h[i][5]}};
for (int j=i+1; j<=n; ++j) {
DBL q[][3] = {{h[j][0], h[j][1]/2., h[j][3]/2.},
{h[j][1]/2., h[j][2], h[j][4]/2.}, {h[j][3]/2., h[j][4]/2., h[j][5]}};
for (int k = conics_intersection(p, q, x, y)-1; k>=0; --k) px[t] = x[k], py[t++] = y[k];
}
}
c = 0;
for (int i=0; i<t; ++i) {
int cc = 0;
for (int j=1; j<=n; ++j) if (check(h[j], px[i], py[i])) p[cc++] = j;
if (cc > c) {
c = cc;
for (int j=0; j<cc; ++j) g[j] = p[j];
}
}
cout << c;
if (c > 0) cout << endl << g[0];
for (int i=1; i<c; ++i) cout << ' ' << g[i];
cout << endl;
}
int main() {
cout << fixed << setprecision(15);
ios::sync_with_stdio(false); cin.tie(0); cout.tie(0);
int t; cin >> t;
for (int k=0; k<t; ++k) {
if (k) cout << endl;
solve();
}
return 0;
}
