空间直线L0, L1的参数方程为
L0:  P = a+u*s;
L1: Q = b+v*t;
其中a为L0的原点，u为单位方向向量, 参数s从-inf到+inf

则 L0上的点P和L1的点Q的距离的平方为
J = (P-Q)*(P-Q) = (P-Q)^2
展开步骤为

=>  (a+u * s - b - v * t)^2
=>  u^2 * s^2 + v^2 * t^2  - 2*u*v *s*t + 2*(a-b)*u*s - 2*(a-b)*v*t + (a-b)^2

J其实是个z = f(x,y)空间曲面, s, t则是x,y轴, z是距离的平方值

根据微分几何，这个曲面z值的最低点就是L0和L1的最短距离，即极值

这个最低点的对s,t各自的偏导数为零，即过该最低点的x,y两个方向的切线水平

即: 

ro_J / ro_s = 0

ro_J / ro_t = 0

写出这两个偏导的表达式为

2*u^2*s - 2*u*v*t + 2*(a-b)*u = 0

2*v^2*t - 2*u*v*s + 2*(a-b)*v = 0

继续整理得

u^2*s - u*v*t = (b-a)*u

u*v*s - v^2*t = (b-a)*v

因为u,v是单位向量，即u^2 = 1,则

s - u*v*t = (b-a)*u

u*v*s - t = (b-a)*v

这里令  C = u*v, D0 = (b-a)*u, D1 = (b-a)*v

s - C*t = D0

C*s - t = D1  

则

s = (D0-C*D1) / (1-C^2)

t = (C*D0-D1) / (1-C^2)

当C为接近1的数, 即u和v几乎平行，则s和t无穷大，编写代码时应该避免.

计算出了s,t则带入J,  则最短距离即得到.