Numerical stable way to find the local extremes of a cubic (maximum and/or minimum).

Consider having a cubic of the form

  f(x) = a x^3 + b x^2 + c x + d

and you want to know if this cubic has a local minimum and if so, where that minimum is.

Then do the following in your code (make sure that a isn't equal to zero):

// Calculate the discriminant of the derivative of the cubic.
double one_fourth_D = b * b - 3.0 * a * c;  // Yep, that's a 3 not a 4.

// If the determinant is zero, then the cubic has no local extremes.
if (one_fourth_D <= 0.0)
  return;

// Use a sqrt with the same sign as b.
double const half_Q = std::copysign(std::sqrt(one_fourth_D), b);

// Now the critical point that you are after can be found with
if (looking_for_the_local_maximum == half_Q < 0.0)
  critical_point = -c / (b + half_Q);
else
  critical_point = -(b + half_Q) / (3.0 * a);

Because half_Q has the same sign as b, and we add b to half_Q, this is numerically stable.

For a real-life example, see https://github.com/CarloWood/machine-learning/blob/master/math/AnalyzedCubic.cxx#L7

CarloWood/README_numerical_stable_parabola_roots.md