diff --git a/MathLib/Nonlinear/cubic_roots.hpp b/MathLib/Nonlinear/cubic_roots.hpp
index ebb42e10a2a3f75566e10f9ed114a013966087f9..cae09ec1b498b1ed25400c3b1dba687966771b0d 100644
--- a/MathLib/Nonlinear/cubic_roots.hpp
+++ b/MathLib/Nonlinear/cubic_roots.hpp
@@ -89,7 +89,7 @@ std::array<Real, 3> cubic_roots(Real a, Real b, Real c, Real d)
         // branch, and we have 3 real roots, two are a double root. Take
         // (x+1)^2(x-2) = x^3 - 3x -2 as an example. This clearly has a double
         // root at x = -1, and it gets sent into this branch.
-        Real arg = R * R - Q * Q * Q;
+        Real arg = std::max(0., R * R - Q * Q * Q);
         Real A = (R >= 0 ? -1 : 1) * cbrt(abs(R) + sqrt(arg));
         Real B = 0;
         if (A != 0)