hdu 4741 Save Labman No.004 异面直线最短距离及对应的点

#include<stdio.h>
#include<math.h>
const double eps=1e-9;
#define zero(x) (((x)>0?(x):-(x))<eps)
int dbcmp(double a)
{
    if(a>eps) return 1;
    else if(a<(-eps)) return -1;
    return 0;
}
struct point
{
    double x,y,z;
    point(){}
    point (double a,double b,double c)
    {
        x=a;
        y=b;
        z=c;
    }
};
double mo(point a)
{
    double ret=a.x*a.x+a.y*a.y+a.z*a.z;
    ret=sqrt(ret);
    return ret;
}
point jian(point a,point b)
{
    point c;
    c.x=a.x-b.x;
    c.y=a.y-b.y;
    c.z=a.z-b.z;
    return c;
}
point jia(point a,point b)
{
    point c;
    c.x=a.x+b.x;
    c.y=a.y+b.y;
    c.z=a.z+b.z;
    return c;
}
point xmult(point u,point v){
    point 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(point u,point v){
    return u.x*v.x+u.y*v.y+u.z*v.z;
}

//矢量差 U - V
point subt(point u,point v){
    point ret;
    ret.x=u.x-v.x;
    ret.y=u.y-v.y;
    ret.z=u.z-v.z;
    return ret;
}
double vlen(point p){
    return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);
}

double linetoline(point u1,point u2,point v1,point v2){
    point n=xmult(subt(u1,u2),subt(v1,v2));
    return fabs(dmult(subt(u1,v1),n))/vlen(n);
}
point pvec(point s1,point s2,point s3){
    return xmult(subt(s1,s2),subt(s2,s3));
}

int dots_onplane(point a,point b,point c,point d){
    return zero(dmult(pvec(a,b,c),subt(d,a)));
}
point intersection(point u1,point u2,point v1,point v2){
    point ret=u1;
    double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x))
        /((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x));
    ret.x+=(u2.x-u1.x)*t;
    ret.y+=(u2.y-u1.y)*t;
    ret.z+=(u2.z-u1.z)*t;
    return ret;
}
double angle_cos(point u1,point u2,point u3,point v1,point v2,point v3){
    return dmult(pvec(u1,u2,u3),pvec(v1,v2,v3))/vlen(pvec(u1,u2,u3))/vlen(pvec(v1,v2,v3));
}
double ooo(point a)
{
    double ret=a.x+a.y+a.z;
    return ret;
}
point cheng(double a,point b)
{
    b.x*=a;
    b.y*=a;
    b.z*=a;
    return b;
}
double getbi(point a,point b,point c,point d)
{
    point di=jian(d,c);
    point ac=jian(c,a);
    point ab=jian(b,a);
    double sab=mo(xmult(di,ab)),sac=mo(xmult(di,ac));
    
    double coss=angle_cos(a,c,d,a,b,jia(di,a));
    
    double shiji=coss*sac;
    
    
    
    double bi=shiji/sab;

    return bi;
    
}
int main()
{
    point a,b,c,d;
    point c2,d2;
    point fa;
    int t;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%lf%lf%lf%lf%lf%lf",&a.x,&a.y,&a.z,&b.x,&b.y,&b.z);
        scanf("%lf%lf%lf%lf%lf%lf",&c.x,&c.y,&c.z,&d.x,&d.y,&d.z);
        double dist=linetoline(a,b,c,d);
        printf("%.6lf\n",dist);
        
        point j1,j2;

        double bi=getbi(a,b,c,d);
        
        j1=jia(a,cheng(bi,jian(b,a)));
        
        bi=getbi(c,d,a,b);

        j2=jia(c,cheng(bi,jian(d,c)));
        
        printf("%.6lf %.6lf %.6lf %.6lf %.6lf %.6lf\n",j1.x,j1.y,j1.z,j2.x,j2.y,j2.z);
    }
    return 0;
}