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
[1] Shape Interrogation for Computer Aided Design and Manufacture by Nichola M. Patrikalakis and Takashi Maekawa