SolveNormalizedCubic fix to return proper real root (#224)

* SolveNormalizedCubic fix to return proper real root Signed-off-by: Piotr Barejko <[email protected]> * Update src/Imath/ImathRoots.h use copy sign Co-authored-by: Nick Porcino <[email protected]> Signed-off-by: Piotr Barejko <[email protected]> Co-authored-by: Nick Porcino <[email protected]>
AcademySoftwareFoundation · Jan 20, 2022 · 7761a3b · 7761a3b
1 parent 4198b33
commit 7761a3b
Show file tree

Hide file tree

Showing 2 changed files with 21 additions and 17 deletions.
diff --git a/src/Imath/ImathRoots.h b/src/Imath/ImathRoots.h
@@ -153,28 +153,31 @@ solveNormalizedCubic (T r, T s, T t, T x[3])
         return 1;
     }
 
-    COMPLEX_NAMESPACE::complex<T> u = COMPLEX_NAMESPACE::pow (
-        -q / 2 + COMPLEX_NAMESPACE::sqrt (COMPLEX_NAMESPACE::complex<T> (D)),
-        T (1) / T (3));
+    if (D > 0)
+    {
+        auto real_root = [] (T a, T x) -> T {
+            T sign = std::copysign(T(1), a);
+            return sign * std::pow (sign * a, T (1) / x);
+        };
 
-    COMPLEX_NAMESPACE::complex<T> v = -p / (T (3) * u);
+        T u = real_root (-q / 2 + std::sqrt (D), 3);
+        T v = -p / (T (3) * u);
 
-    const T sqrt3 = T (1.73205080756887729352744634150587); // enough digits
-                                                            // for long double
-    COMPLEX_NAMESPACE::complex<T> y0 (u + v);
+        x[0] = u + v - r / 3;
+        return 1;
+    }
 
-    COMPLEX_NAMESPACE::complex<T> y1 (-(u + v) / T (2) +
-                                      (u - v) / T (2) * COMPLEX_NAMESPACE::complex<T> (0, sqrt3));
+    namespace CN     = COMPLEX_NAMESPACE;
+    CN::complex<T> u = CN::pow (-q / 2 + CN::sqrt (CN::complex<T> (D)), T (1) / T (3));
+    CN::complex<T> v = -p / (T (3) * u);
 
-    COMPLEX_NAMESPACE::complex<T> y2 (-(u + v) / T (2) -
-                                      (u - v) / T (2) * COMPLEX_NAMESPACE::complex<T> (0, sqrt3));
+    const T sqrt3 = T (1.73205080756887729352744634150587); // enough digits
+                                                            // for long double
+    CN::complex<T> y0 (u + v);
+    CN::complex<T> y1 (-(u + v) / T (2) + (u - v) / T (2) * CN::complex<T> (0, sqrt3));
+    CN::complex<T> y2 (-(u + v) / T (2) - (u - v) / T (2) * CN::complex<T> (0, sqrt3));
 
-    if (D > 0)
-    {
-        x[0] = y0.real() - r / 3;
-        return 1;
-    }
-    else if (D == 0)
+    if (D == 0)
     {
         x[0] = y0.real() - r / 3;
         x[1] = y1.real() - r / 3;

diff --git a/src/ImathTest/testRoots.cpp b/src/ImathTest/testRoots.cpp
@@ -206,6 +206,7 @@ testRoots()
     solve (0, 0, 1, 0, 1, 0, 0, 0);      // real solutions: 0
     solve (0, 0, 0, 1, 0, 0, 0, 0);      // real solutions: none
     solve (0, 0, 0, 0, -1, 0, 0, 0);     // real solutions: [-inf, inf]
+    solve (36, -37, 0, 12, 1, -0.4715155, 0, 0);
 
     cout << endl << "solveQuadratic" << endl;
     // Solve quadratic equations