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 才能过，，，，郁闷无极限啊