I don't think there was any trick. One approach is to use the rational quadratic parametrization of the ellipse and substitute that into the squared distance

I don't think there was any trick. One approach is to use the rational quadratic parametrization of the ellipse and substitute that into the squared distance function to get a quartic. >I try to avoid using quadratic solutions in code because of numeric instability issues. Sensible! Back then the quartic solver had been written already and extensively tested by someone else - probably for line X torus intersection. Even so there were numerical problems introduced because the quartic was in power basis and so the coefficients didn't have geometric meaning. I guess phenomena like Wilkinson's polynomial could occur. [0] If I was doing it today, I would probably proceed as follows: 1. special case for the point exactly on an axis. 2. otherwise, choose the quadrant containing the point and represent that as a rational quadratic Bezier curve. This parametrizes the quadrant from 0 to 1 and is guaranteed to contain the global minimum. 3. substitute the quadratic Bezier curve into the square distance function to get a Bernstein basis quartic polynomial. 4. differentiate to get a Bernstein basic cubic. 5. solve the roots of the cubic numerically in the range 0 to 1: either Newton's or clipping or some other hybrid numerical method. A good reference is [1]. [0] https://en.wikipedia.org/wiki/Wilkinson%27s_polynomial https://en.wikipedia.org/wiki/Wilkinson%27s_polynomial [1] Shape Interrogation for Computer Aided Design and Manufacture by Nichola M. Patrikalakis and Takashi Maekawa