转载学习：
#include <cstdio> #include <cstdlib> #include <cstring> #include <algorithm> #include <cmath> using namespace std; const double EPS = 1e-9; const int MAXN = 40; struct Point3 //空间点 { double x, y, z; Point3( double x=0, double y=0, double z=0 ): x(x), y(y), z(z) { } Point3( const Point3& a ) { x = a.x; y = a.y; z = a.z; return; } void showP() { printf("%f %f %f ", x, y, z); } Point3 operator+( Point3& rhs ) { return Point3( x+rhs.x, y+rhs.y, z+rhs.z ); } }; struct Line3 //空间直线 { Point3 a, b; }; struct plane3 //空间平面 { Point3 a, b, c; plane3() {} plane3( Point3 a, Point3 b, Point3 c ): a(a), b(b), c(c) { } void showPlane() { a.showP(); b.showP(); c.showP(); return; } }; double dcmp( double a ) { if ( fabs( a ) < EPS ) return 0; return a < 0 ? -1 : 1; } //三维叉积 Point3 Cross3( 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; } //三维点积 double Dot3( Point3 u, Point3 v ) { return u.x * v.x + u.y * v.y + u.z * v.z; } //矢量差 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; } //两点距离 double TwoPointDistance( Point3 p1, Point3 p2 ) { return sqrt( (p1.x - p2.x)*(p1.x - p2.x) + (p1.y - p2.y)*(p1.y - p2.y) + (p1.z - p2.z)*(p1.z - p2.z) ); } //向量的模 double VectorLenth( Point3 p ) { return sqrt( p.x*p.x + p.y*p.y + p.z*p.z ); } //空间直线距离 double LineToLine( Line3 u, Line3 v, Point3& tmp ) { tmp = Cross3( Subt( u.a, u.b ), Subt( v.a, v.b ) ); return fabs( Dot3( Subt(u.a, v.a), tmp ) ) / VectorLenth(tmp); } //取平面法向量 Point3 pvec( plane3 s ) { return Cross3( Subt( s.a, s.b ), Subt( s.b, s.c ) ); } //空间平面与直线的交点 Point3 Intersection( Line3 l, plane3 s ) { Point3 ret = pvec(s); double t = ( ret.x*(s.a.x-l.a.x)+ret.y*(s.a.y-l.a.y)+ret.z*(s.a.z-l.a.z) )/( ret.x*(l.b.x-l.a.x)+ret.y*(l.b.y-l.a.y)+ret.z*(l.b.z-l.a.z) ); ret.x = l.a.x + ( l.b.x - l.a.x ) * t; ret.y = l.a.y + ( l.b.y - l.a.y ) * t; ret.z = l.a.z + ( l.b.z - l.a.z ) * t; return ret; } /************以上模板*************/ void solved( Line3 A, Line3 B ) { Point3 normal; double dis = LineToLine( A, B, normal ); printf( "%.6f ", dis ); plane3 pla; pla = plane3( A.a, A.b, A.a + normal ); Point3 u = Intersection( B, pla ); pla = plane3( B.a, B.b, B.a + normal ); Point3 v = Intersection( A, pla ); printf("%.6f %.6f %.6f %.6f %.6f %.6f ", v.x, v.y, v.z, u.x, u.y, u.z ); return; } int main() { int T; scanf( "%d", &T ); while ( T-- ) { Line3 A, B; scanf("%lf%lf%lf", &A.a.x, &A.a.y, &A.a.z ); scanf("%lf%lf%lf", &A.b.x, &A.b.y, &A.b.z ); scanf("%lf%lf%lf", &B.a.x, &B.a.y, &B.a.z ); scanf("%lf%lf%lf", &B.b.x, &B.b.y, &B.b.z ); solved( A, B ); } return 0; }
不知精度误差的WA
#include<stdio.h> #include<math.h> #define eps 1e-12 double myfabs(double x) { if(x<0)x=-x; return x; } int main() { int _case; /*double xa,xb,xc,xd; double ya,yb,yc,yd; double za,zb,zc,zd; */ double Xa,Xb,Xc,Xd,Ya,Yb,Yc,Yd,Za,Zb,Zc,Zd; scanf("%d",&_case); while(_case--) { /*scanf("%lf%lf%lf%lf%lf%lf",&xa,&ya,&za,&xb,&yb,&zb); scanf("%lf%lf%lf%lf%lf%lf",&xc,&yc,&zc,&xd,&yd,&zd);*/ scanf("%lf%lf%lf%lf%lf%lf",&Xa,&Ya,&Za,&Xb,&Yb,&Zb); scanf("%lf%lf%lf%lf%lf%lf",&Xc,&Yc,&Zc,&Xd,&Yd,&Zd); /*long double Xa=(long double)xa; long double Xb=(long double)xb; long double Xc=(long double)xc; long double Xd=(long double)xd; long double Ya=(long double)ya; long double Yb=(long double)yb; long double Yc=(long double)yc; long double Yd=(long double)yd; long double Za=(long double)za; long double Zb=(long double)zb; long double Zc=(long double)zc; long double Zd=(long double)zd;*/ double F11=(Xb-Xa)*(Xb-Xa)+(Yb-Ya)*(Yb-Ya)+(Zb-Za)*(Zb-Za); double F12= (Xd-Xc)*(Xd-Xc)+(Yd-Yc)*(Yd-Yc)+(Zd-Zc)*(Zd-Zc); double F2=(Xb-Xa)*(Xd-Xc)+(Yb-Ya)*(Yd-Yc)+(Zb-Za)*(Zd-Zc); double F31=(Xb-Xa)*(Xc-Xa)+(Yb-Ya)*(Yc-Ya)+(Zb-Za)*(Zc-Za); double F32=(Xd-Xc)*(Xc-Xa)+(Yd-Yc)*(Yc-Ya)+(Zd-Zc)*(Zc-Za); double y=F11*F12-F2*F2; //if(myfabs(y)<eps)y=eps; double t1=(F31*F12-F32*F2)/y; double t2=(F32*F11-F2*F31)/(-y); double Xm=t1*(Xb-Xa)+Xa;//=(Xb-Xa)*[F31*F12-F32*F2]/[F11*F12-F2*F2]+Xa; double Ym=t1*(Yb-Ya)+Ya;//=(Yb-Ya)*[F31*F12-F32*F2]/[F11*F12-F2*F2]+Ya; double Zm=t1*(Zb-Za)+Za;//=(Zb-Za)*[F31*F12-F32*F2]/[F11*F12-F2*F2]+Za; double Xn=t2*(Xd-Xc)+Xc;//=(Xd-Xc)*[F3(c,d)*F1(a,b)-F3(a,b)*F2()]/[F2()*F2()-F1(a,b)*F1(c,d)]+Xc; double Yn=t2*(Yd-Yc)+Yc;//=(Yd-Yc)*[F3(c,d)*F1(a,b)-F3(a,b)*F2()]/[F2()*F2()-F1(a,b)*F1(c,d)]+Yc; double Zn=t2*(Zd-Zc)+Zc; double s=sqrt((Xn-Xm)*(Xn-Xm)+(Yn-Ym)*(Yn-Ym)+(Zn-Zm)*(Zn-Zm)); /*double xm=(double)Xm;//=t1*(Xb-Xa)+Xa;//=(Xb-Xa)*[F31*F12-F32*F2]/[F11*F12-F2*F2]+Xa; double ym=(double)Ym;//=t1*(Yb-Ya)+Ya;//=(Yb-Ya)*[F31*F12-F32*F2]/[F11*F12-F2*F2]+Ya; double zm=(double)Zm;//=t1*(Zb-Za)+Za;//=(Zb-Za)*[F31*F12-F32*F2]/[F11*F12-F2*F2]+Za; double xn=(double)Xn;//=t2*(Xd-Xc)+Xc;//=(Xd-Xc)*[F3(c,d)*F1(a,b)-F3(a,b)*F2()]/[F2()*F2()-F1(a,b)*F1(c,d)]+Xc; double yn=(double)Yn;//=t2*(Yd-Yc)+Yc;//=(Yd-Yc)*[F3(c,d)*F1(a,b)-F3(a,b)*F2()]/[F2()*F2()-F1(a,b)*F1(c,d)]+Yc; double zn=(double)Zn; //printf("%lf ",eps);*/ printf("%.6lf %.6lf %.6lf %.6lf %.6lf %.6lf %.6lf ",s,Xm,Ym,Zm,Xn,Yn,Zn); //printf("%.6lf %.6lf %.6lf %.6lf %.6lf %.6lf %.6lf ",s,xm,ym,zm,xn,yn,zn); } return 0; }