hdu 4617 计算几何 三维点线面关系

/*
题意：给出n个weapon，每个weapon都是一个镭射光束，相当于一个无限长的实心圆柱，当两个weapon相接触，就会产生爆炸；
问给出n个weapon的某一垂直切面(为圆形)，求这n个weapon之间是否会相接触，会则输出Lucky,否则输出两两之间weapon的最
短距离；

题解：计算几何：两直线的距离+平面法向量+直线与平面交点+平行线段
首先分别求出两个平面的法向量，因为圆柱的中心轴平行于平面法向量，然后求该向量与平面交点to，通过原点o与点to以及
圆柱的中心点center，将线段o-to平移到center-aim(aim为要求的点)，这样center，aim两点就能表示中心轴线，最后求两直
线间的距离，最终减去两个半径即可得出最终答案。

注意：当给出的平面的center为原点时，需要特判，否则会出错
*/
#include <cstdio>
#include <cmath>
#include <algorithm>

#define eps 1e-8
#define zero(x) (((x)>0?(x):-(x))<eps)

struct point3{double x,y,z;};
struct line3{point3 a,b;}linea,lineb;

//计算cross product U x V
point3 xmult(point3 u,point3 v){
    point3 ret;
    ret.x=u.y*v.z-v.y*u.z;
    ret.y=u.z*v.x-u.x*v.z;
    ret.z=u.x*v.y-u.y*v.x;
    return ret;
}

//计算dot product U . V
double dmult(point3 u,point3 v){
    return u.x*v.x+u.y*v.y+u.z*v.z;
}

//矢量差 U - V
point3 subt(point3 u,point3 v){
    point3 ret;
    ret.x=u.x-v.x;
    ret.y=u.y-v.y;
    ret.z=u.z-v.z;
    return ret;
}

// 求平面法向量(原点到返回点即为所求向量)
point3 pvec(point3 s1,point3 s2,point3 s3)
{
    return xmult(subt(s1,s2),subt(s2,s3));
}

//向量大小
double vlen(point3 p){
    return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);
}

// 两直线距离
double linetoline(point3 u1,point3 u2,point3 v1,point3 v2){
    point3 n=xmult(subt(u1,u2),subt(v1,v2));
    return fabs(dmult(subt(u1,v1),n))/vlen(n);
}

// 求直线与平面的交点
point3 intersection(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){
    point3 ret=pvec(s1,s2,s3);
    double t=(ret.x*(s1.x-l1.x)+ret.y*(s1.y-l1.y)+ret.z*(s1.z-l1.z))/
        (ret.x*(l2.x-l1.x)+ret.y*(l2.y-l1.y)+ret.z*(l2.z-l1.z));
    ret.x=l1.x+(l2.x-l1.x)*t;
    ret.y=l1.y+(l2.y-l1.y)*t;
    ret.z=l1.z+(l2.z-l1.z)*t;
    return ret;
}


struct location
{
    point3 a,b,o;
}s[29];

double getdistance(point3 a, point3 b)
{
    return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y) + (a.z-b.z)*(a.z-b.z));
}


int main(void)
{
    int t,n;
    scanf("%d",&t);
    while (t--)
    {
        scanf("%d",&n);
        for(int i=0; i<n; i++)
        {
            scanf("%lf%lf%lf",&s[i].o.x, &s[i].o.y, &s[i].o.z);
            scanf("%lf%lf%lf",&s[i].a.x, &s[i].a.y, &s[i].a.z);
            scanf("%lf%lf%lf",&s[i].b.x, &s[i].b.y, &s[i].b.z);
        }

        double ans = 1e10;
        bool flag = true;
        for(int i=0; i<n-1; i++)
        {
            for(int j=i+1; j<n; j++)
            {
                point3 to1,to2,origin,tmp;
                origin.x = origin.y = origin.z = 0;
                to1 = pvec(s[i].a,s[i].o,s[i].b);
                tmp = intersection(to1,origin,s[i].o,s[i].a,s[i].b);
                // 特判该平面的center是否原点
                if (!(zero(tmp.x-origin.x) && zero(tmp.y-origin.y) && zero(tmp.z-origin.z)))
                {
                    to1.x = s[i].o.x - to1.x;
                    to1.y = s[i].o.y - to1.y;
                    to1.z = s[i].o.z - to1.z;
                }

                to2 = pvec(s[j].a,s[j].o,s[j].b);
                tmp = intersection(to2,origin,s[j].o,s[j].a,s[j].b);
                // 特判该平面的center是否原点
                if (!(zero(tmp.x-origin.x) && zero(tmp.y-origin.y) && zero(tmp.z-origin.z)))
                {
                    to2.x = s[j].o.x - to2.x;
                    to2.y = s[j].o.y - to2.y;
                    to2.z = s[j].o.z - to2.z;
                }
                double radis1 = getdistance(s[i].o,s[i].a);
                double radis2 = getdistance(s[j].o,s[j].a);
                double linedis = linetoline(s[i].o,to1,s[j].o,to2);
                linedis = linedis - radis1 - radis2;
                if (linedis < eps)
                {
                    flag = false;
                    break;
                }
                ans = std::min(ans,linedis);
            }
            if (!flag)
                break;
        }
        if (flag)
            printf("%.2lf\n",ans);
        else
            printf("Lucky\n");

    }
    return 0;
}