OT: I had the problem as part of my 1st day job in 1990 on Parasolid.

IIRC, there was more than one way of converting the problem into solving a quartic polynomial, which is solvable by radicals and so doesn't need Newton's method.

If the point is inside the evolute of the ellipse then the quartic has four real solutions, which gives two local minima and two local maxima of distance.

If the point is outside or on the evolute of the ellipse then the quartic has two real solutions, which gives one local minima and one local maxima of distance.

Having said that, all methods using quartic polynomials could be unstable if the ellipse was pathological, for example very long and thin. Maybe that would be a time to revert to Newton's method.

I try to avoid using quadratic solutions in code because of numeric instability issues. Granted I never implemented the direct solution, but the number of opportunities for catastrophic cancellation gave me shivers.

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

How the heck did "" get in there?

> 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.

Is there a reason not to use the cubic formula in this situation?

Something like

   (E(t)-P).(E(t)-P)

>Is there a reason not to use the cubic formula in this situation?

Numerical stability in edge cases, which will occur for code exercised as extensively as Parasolid's.

I guess you could a Bernstein version of the cubic formula, though I've not seen one derived in the wild. It would be nice to have as well as Bernstein quartic formula. The Bernstein quadratic is easy.