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 . Vdouble dmult(point u,point v){    return u.x*v.x+u.y*v.y+u.z*v.z;}//矢量差 U - Vpoint 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;}