+ Real delta0 = (b*b) - Real(3) * a * c;
+ Real delta1 = Real(2) * (b * b * b) - Real(9) * a * b * c + Real(27) * (a * a) * d;
+
+ std::vector<Real> roots;
+
+ Real C = pow((delta1 + Sqrt((delta1 * delta1) - 4 * (delta0 * delta0 * delta0)) ) / Real(2), 1/3);
+
+ if (false && discriminant < 0)
+ {
+ Real real_root = (Real(-1) / (Real(3) * a)) * (b + C + delta0 / C);
+
+ roots.push_back(real_root);
+
+ return roots;
+
+ }
+