hdu 4741 Save Labman No.004

题目大意：求空间两条直线的最短距离及最近的点坐标

解题思路：http://wenku.baidu.com/view/ea02ff8a6529647d2728520c.html

代码实现1：

#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <iostream>
using namespace std;
struct node
{
    double x;
    double y;
    double z;
}p[10], P, Q;
long double s, t;
long double pf(long double a, long double b)
{
     return (a - b) * (a - b);
}
long double dis(node a, node b)
{
     return sqrt(pf(a.x, b.x) + pf(a.y, b.y) + pf(a.z, b.z));
}
void Cal()
{
     long double t1 = pf(p[2].x, p[1].x) + pf(p[2].y, p[1].y) + pf(p[2].z ,p[1].z);
     long double t2 = (p[2].x - p[1].x) * (p[4].x - p[3].x) + (p[2].y - p[1].y) * (p[4].y - p[3].y)
                 + (p[2].z - p[1].z) * (p[4].z - p[3].z);
     long double t3 = (p[1].x - p[2].x) * (p[1].x - p[3].x) + (p[1].y - p[2].y) * (p[1].y - p[3].y)
                 + (p[1].z - p[2].z) * (p[1].z - p[3].z);

     long double t4 = pf(p[4].x, p[3].x) + pf(p[4].y, p[3].y) + pf(p[4].z ,p[3].z);
     long double t5 = (p[1].x - p[3].x) * (p[4].x - p[3].x) + (p[1].y - p[3].y) * (p[4].y - p[3].y)
                 + (p[1].z - p[3].z) * (p[4].z - p[3].z);
     t = (t1 * t5 + t2 * t3) / (t1 * t4 - t2 * t2);
     s = (t2 * t5 + t3 * t4) / (t1 * t4 - t2 * t2);

}
int main()
{
    int t0;
    scanf("%d", &t0);
    for(int i = 1; i <= t0; i++)
    {
            for(int j = 1; j <= 4; j++)
            {
                    scanf("%lf%lf%lf", &p[j].x, &p[j].y, &p[j].z);
            }
            Cal();
            P.x = p[1].x + s * (p[2].x - p[1].x);
            P.y = p[1].y + s * (p[2].y - p[1].y);
            P.z = p[1].z + s * (p[2].z - p[1].z);

            Q.x = p[3].x + t * (p[4].x - p[3].x);
            Q.y = p[3].y + t * (p[4].y - p[3].y);
            Q.z = p[3].z + t * (p[4].z - p[3].z);

            long double ans = dis(P, Q);
            printf("%.6lf\n", (double)ans);
            printf("%.6lf %.6lf %.6lf %.6lf %.6lf %.6lf\n",(double)P.x, (double)P.y, (double)P.z, (double)Q.x, (double)Q.y, (double)Q.z);

    }
    return 0;
}

http://blog.sina.com.cn/s/blog_a401a1ea0101ij9z.html

容易理解的常规方法：

已知空间中两线段，如果它们无限变粗，判断是否相交。（主要讨论不在同一平面的情况）

线段AB 线段CD

问题的关键是求出这两条任意直线之间的最短距离，以及在这个距离上的两线最接近点坐标，判断该点是否在线段AB和线段CD上。

首先将直线方程化为对称式，得到其方向向量n1=（a1,b1,c1),n2=(a2,b2,c2).

再将两向量叉乘得到其公垂向量N=（x,y,z）,在两直线上分别选取点A,B(任意)，得到向量AB，

求向量AB在向量N方向的投影即为两异面直线间的距离了（就是最短距离啦）。

最短距离的求法：d=|向量N*向量AB|/|向量N|（上面是两向量的数量积，下面是取模）。

设交点为C,D，带入公垂线N的对称式中，又因为C,D两点分别满足一开始的直线方程，所以得到关于C（或D）的两个连等方程，分别解出来就好了！

没有理解的简单方法：

三维向量的叉积公式：

向量a×向量b=| i j k| |a1 b1 c1| |a2 b2 c2|

这是一个三阶行列式

其值为 (b1c2-b2c1, c1a2-a1c2, a1b2-a2b1)

（i、j、k分别为空间中相互垂直的三条坐标轴的单位向量）。

代码实现2：

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;

struct Point
{
     double x, y, z;
     Point(double x = 0, double  y = 0, double z = 0):x(x), y(y), z(z){}
};

typedef Point Vector;

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

Vector Cross(Vector a, 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);
}

double Dot(Vector a, Vector b)
{
     return (a.x * b.x + a.y * b.y + a.z * b.z);
}
double Getlen(Vector a)
{
      return sqrt(a.x*a.x + a.y*a.y + a.z*a.z);
}
int  main()
{
     Point a1, b1, a2, b2;
     int n;
     scanf("%d", &n);
     for(int i = 0; i < n; i ++)
     {
             scanf("%lf %lf %lf", &a1.x, &a1.y, &a1.z);
             scanf("%lf %lf %lf", &b1.x, &b1.y, &b1.z);
             scanf("%lf %lf %lf", &a2.x, &a2.y, &a2.z);
             scanf("%lf %lf %lf", &b2.x, &b2.y, &b2.z);

             Vector n1 = (a1 - b1);
             Vector n2 = (a2 - b2);

             Vector n = Cross(n1, n2);

             Vector ab = (a1 - a2);
             double ans = Dot(n , ab) / Getlen(n);

             Point p1 = a1, p2 = a2;
             Vector d1 = (b1 - a1), d2 = (b2 - a2);

             double t1 = Dot( Cross(p2 - p1, d2), Cross(d1, d2));
             double t2 = Dot( Cross(p2 - p1, d1) , Cross(d1, d2));
             double fm = Getlen(Cross(d1, d2));
             t1 = t1 / (fm * fm);
             t2 = t2 / (fm * fm);

             Point p, q;
             p = (a1 + (b1 - a1) * t1);
             q = (a2 + (b2 - a2) * t2);

             printf("%.6f\n", fabs(ans));
             printf("%.6f %.6f %.6f ",p.x, p.y, p.z);
             printf("%.6f %.6f %.6f\n",q.x, q.y, q.z);


     }
     return 0;
}

有一点我好想吐槽：为什么一定要用long double 才能过，，，，郁闷无极限啊