This would, of course, work. And if you were trying to be fancy, you could even do a binary search instead of a linear one.

```
# I'm assuming some functions are defined:
# project_distance() projects a vector onto a unit vector, and gets the distance
# between the vector and it's projection.
# Could be optimized to return the squared distance
# project_point() projects a vector onto a unit vector, and returns the projected point
# distance() finds the euclidean distance between two points.
# Could be optimized to return the squared distance
# addition and subtraction are defined between points/vectors
# multiplication and division are defined between points/vectors and scalars
def closest_on_skew(p1, p2, p3, p4, tol=0.00001):
p1_p2_normal = (p2-p1).normal()
p3_p4_normal = (p4-p3).normal()
# Subtract p3 out from the check points. Because if we do this relative to p3
# then it simplifies a lot of the math, and we can just add it back in later
# This is an easy optimization, so long as it's commented
c1 = p1-p3
c2 = p2-p3
c1_dist = project_distance(c1, p3_p4_normal)
c2_dist = project_distance(c2, p3_p4_normal)
while distance(c1, c2) > tol:
midpoint = (c1 + c2) / 2
mid_dist = project_distance(midpoint, p3_p4_normal)
if c1_dist > c2_dist:
c1 = midpoint
c1_dist = mid_dist
elif c2_dist > c1_dist:
c2 = midpoint
c2_dist = mid_dist
else:
# They're equal. So the only *possible* cases are that
# the lines are parallel, or the midpoint is the exact
# minimum.
break
# We've either hit tolerance, found parallel lines, or an exact min
# so take one last midpoint and return the points as requested
midpoint = (c1 + c2) / 2
# Don't forget to add p3 back into the mix
return midpoint + p3, project_point(midpoint, p3_p4_normal) + p3
```

However, since thereâ€™s a direct algebraic solution that you can find by solving for this property

`dot(n1-n0, p1_p2_normal) == 0 and dot(n1-n0, p3_p4_normal) == 0`

(Meaning the shortest segment that connects two skew lines is perpendicular to both of the skew lines)

I would just do the algebra once, and write code for the exact solution