Rounding issues in the quadratic formula have bitten me several times: in particular, the inverse of a 2D bilinear mapping requires solving a quadratic equation which (if things are skewed) can be very badly conditioned. On the other hand, if one uses a bilinear map on a parallelogram then the quadratic term is gone and one root is zero: also annoying.

In practice I just check the coefficients and then use appropriate versions of the quadratic formula. I guess that I am luckier than Mr. Cook here since I only need one root.

code (some of it was written by me, some improvements are by other people): https://github.com/dealii/dealii/blob/master/source/fe/mappi...

This set of cases (which can probably be further improved) is equivalent to his first answer. I don't think its possible to do better in this test case due to the catastrophic cancelation in subtracting the discriminant from b.