hdu 4741——Save Labman No.004

来源:互联网 发布:js浏览器定位 编辑:程序博客网 时间:2024/06/06 01:54

题意:求两个空间直线的距离,并求出这两个直线中距离最近两个点的坐标
思路:刚刚复习完高数的立体几何部分。。然后就用了在复习全书上看到的一个方法来求解。
两个向量的叉积表示得到一个垂直于这两个向量的向量,所得向量的膜是这两个向量组成平行四边形的面积。假设题目所给两个直线的方向向量是s1,s2,在两个直线上任选两点,这两点构成向量AB。
则两个直线距离就是(s1 s2 AB)的混合积除以(s1s2叉积的膜)。混合积的几何意义就是三个向量构成立方体的体积,s1s2的叉积为这个立方体的底面积。两个相除得到高,就是两个直线的距离。
这里写图片描述
求两个距离最短点的时候,我的思路是用偏导数来求,不过貌似精度不够。先写在这里。
两个直线用参数方程形式写出来,得到
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;}
0 0
原创粉丝点击