hdu 4741——Save Labman No.004

题意：求两个空间直线的距离，并求出这两个直线中距离最近两个点的坐标
思路：刚刚复习完高数的立体几何部分。。然后就用了在复习全书上看到的一个方法来求解。
两个向量的叉积表示得到一个垂直于这两个向量的向量，所得向量的膜是这两个向量组成平行四边形的面积。假设题目所给两个直线的方向向量是s1,s2，在两个直线上任选两点，这两点构成向量AB。
则两个直线距离就是(s1 s2 AB)的混合积除以（s1与s2叉积的膜）。混合积的几何意义就是三个向量构成立方体的体积，s1与s2的叉积为这个立方体的底面积。两个相除得到高，就是两个直线的距离。

求两个距离最短点的时候，我的思路是用偏导数来求，不过貌似精度不够。先写在这里。
两个直线用参数方程形式写出来，得到
L1: X1 = at+x1 Y1 = bt+y1 Z1 = ct + z1
L2:X2 = dm +x2 Y2 = em+y2 Z2 = fm+z2
其中，(a,b,c)(d,e,f)是方向向量
距离为D=（X1-X2）^2+(Y1-Y2)^2+(Z1-Z2)^2
当距离最小的时候，D对t的偏导数及D对m的偏导数都是0，最后化成解二元一次方程组。
代码如下：

#include <cstring>
#include <stdio.h>
#include <algorithm>
#include <iostream>
#include <cmath>
#include <map>
#include <string>
#include <queue>
#include <bitset>
#include<assert.h>
using namespace std;

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

struct Point
{
    double x,y,z;
    Point(double _x=0,double _y=0,double _z=0):x(_x),y(_y),z(_z){
        if(dcmp(x)==0)x=0;
        if(dcmp(y)==0)y=0;
        if(dcmp(z)==0)z=0;
    }
    void scan(){
        scanf("%lf%lf%lf",&x,&y,&z);
    }
    void p()
    {
        printf("%.6lf %.6lf %.6lf",x,y,z);
    }
};
typedef Point Vector;

Vector operator + (const Vector &a,const Vector &b){ return Vector(a.x+b.x,a.y+b.y,a.z+b.z);}
Vector operator - (const Vector &a,const Vector &b){ return Vector(a.x - b.x,a.y - b.y, a.z - b.z);}
Vector operator * (const Vector &a,const double &b){ return Vector(a.x*b,a.y*b,a.z*b);}
Vector operator / (const Vector &a,const double &b){ return Vector(a.x/b,a.y/b,a.z/b);}

double dot(const Vector &a,const Vector &b){
    return a.x*b.x+a.y*b.y+a.z*b.z;
}
double Lenth(const Vector &a)
{
    return sqrt((dot(a,a)));
}

Vector det(const Vector &a,const Vector &b)
{
    return Vector(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);
}



int main()
{
//    freopen("data.txt","r",stdin);
    int T;
    scanf("%d",&T);
    while(T--){
        Point a1,a2,b1,b2;
        a1.scan();a2.scan();b1.scan();b2.scan();
        Vector a = a1-a2;
        Vector b = b1-b2;
        Vector AB=b1-a1;
        double ansl=fabs(dot(det(a,b),AB))/fabs(Lenth(det(a,b)));
        printf("%.6lf\n",ansl);
        //======================/这里是我的方法
        double a11=dot(a,a);
        double a12=-(dot(a,b));
        double a21 = -dot(a,b);
        double a22 = dot(b,b);
        if(dcmp(a11)==0)a11=0;
        if(dcmp(a12)==0)a12=0;
        if(dcmp(a21)==0)a21=0;
        if(dcmp(a22)==0)a22=0;
        double c1=dot(a,AB);
        double c2 = -dot(b,AB);
        if(dcmp(c1)==0)c1=0;
        if(dcmp(c2)==0)c2=0;
        double S =a11*a22-a12*a21;

        double t = (-c1*a22+a12*c2)/S;
        double m = (-a11*c2+c1*a21)/S;
        Point ans1=a*(-t)+a1;
        Point ans2 = b*(-m)+b1;
        //=======================/
//++++++++++++++++++++++++++/参考http://blog.sina.com.cn/s/blog_a401a1ea0101ij9z.html
        double t1 = dot(det(AB,b),det(a,b));
        t1/=dot(det(a,b),det(a,b));
        double t2 = dot(det(AB,a),det(a,b));
        t2/=dot(det(a,b),det(a,b));

        Point aa=a1+(a*t1);
        Point bb = b1+(b*t2);
        aa.p();printf(" ");
        bb.p();puts("");
        Vector test1=aa-ans1;
        Vector test2=bb-ans2;
//+++++++++++++++++++     
//下面将两种方法的结果对比了一下，发现当精度在1e-6的时候两个结果是一样的，但是精度到1e-7的时候，就错了   
if(dcmp(test1.x)||dcmp(test1.y)||dcmp(test1.z))cout<<1/0;
        if(dcmp(test2.x)||dcmp(test2.y)||dcmp(test2.z))cout<<1/0;
    }
    return 0;
}