a) i shall assume that the distance between the two non-parallel lines is not zero (i.e. that the two lines are skew).
the direction vector of the first line is d1 = (a1, b1, c1) and the direction vector of the second line is d1 = (a2, b2, c2). A (x1, y1, z1) is some random point on the first line, and B (x2, y2, z2) is some random point on the second line. now, to find the shortest distance between two skew lines, all we need to do is find the absolute value of the scalar resolute (yay spesh) of vector BA in the direction of d1 x d2. why? picture to yourself two skew lines. picture also the plane on which each skew line lies. now, if you bring the two 'parallel' planes together, then the two direction vectors will lie on the 'same' plane. the cross product of the two direction vectors (d1 x d2) will, therefore, always be in the direction perpendicular to this new plane we created by merging the two planes together (and hence in the direction perpendicular to both original planes). (notice how d1 x d2 = 0 if the two lines are parallel, since d1 would simply be a scalar multiple of the d2.) this is convenient because all we need to do now is take a point on each line, join them up to form a vector, and find the magnitude of the projection of this newly formed vector on the 'normal' vector (d1 x d2) to find the shortest distance between the two lines.
as for the working:
D = (BA.(d1xd2))/abs(d1xd2) = (x1-x2, y1-y2,z1-z2).(d1xd2)/abs(d1xd2), as req. (this formula won't work if the lines are parallel because then you would be dividing by zero.)
b) the area of the parallelogram with edge vectors (a1, b1, c1) and (x2-x1, y2-y1, z2-z1) is given by norm((a1, b1, c1) x (x2-x1, y2-y1, z2-z1)). now we know that norm((a1, b1, c1)) is sqrt(a1^2 + b1^2 + c1^2).
now, realise that the height of the parallelogram corresponds to the shortest distance between the two parallel lines.
h*sqrt(a1^2 + b1^2 + c1^2) = norm((a1, b1, c1) x (x2-x1, y2-y1, z2-z1))
h = norm((a1, b1, c1) x (x2-x1, y2-y1, z2-z1))/sqrt(a1^2 + b1^2 + c1^2)
therefore, D = norm((a1, b1, c1) x (x2-x1, y2-y1, z2-z1))/sqrt(a1^2 + b1^2 + c1^2)
hopefully this is right...haha.