Added workaround for bug signaled by Emanuele Siego

ebertolazzi · Jan 15, 2021 · 8e8fbe3 · 8e8fbe3
1 parent 3d0259d
commit 8e8fbe3
Show file tree

Hide file tree

Showing 3 changed files with 24 additions and 4 deletions.
diff --git a/...codeproj/project.xcworkspace/xcuserdata/enrico.xcuserdatad/UserInterfaceState.xcuserstate b/...codeproj/project.xcworkspace/xcuserdata/enrico.xcuserdatad/UserInterfaceState.xcuserstate
diff --git a/src/PolynomialRoots-3-Quartic.cc b/src/PolynomialRoots-3-Quartic.cc
@@ -495,6 +495,7 @@ namespace PolynomialRoots {
     valueType t = qsolve.real_root1();
     valueType s = qsolve.real_root2();
 
+    bool must_refine_r3 = true;
     if ( !qsolve.complexRoots() ) {
       valueType Qs = evalMonicQuartic( s, q3, q2, q1, q0 );
       valueType Qu = evalMonicQuartic( u, q3, q2, q1, q0 );
@@ -521,9 +522,13 @@ namespace PolynomialRoots {
         else if ( Qu <= epsi ) r3 = u;
         else                   nreal = 0;
         */
-        if      ( isZero(Qs) ) r3    = s;
-        else if ( isZero(Qu) ) r3    = u;
-        else                   nreal = 0;
+        if ( isZero(Qs) ) {
+          must_refine_r3 = false; r3 = s;
+        } else if ( isZero(Qu) ) {
+          must_refine_r3 = false; r3 = u;
+        } else {
+          nreal = 0;
+        }
       }
     } else {
       // one single real root (only 1 minimum)
@@ -548,7 +553,8 @@ namespace PolynomialRoots {
     ..  oscillation brackets.
     */
     if ( nreal > 0 ) {
-      iter += zeroQuarticByNewtonBisection( q3, q2, q1, q0, r3 );
+      if ( must_refine_r3 )
+        iter += zeroQuarticByNewtonBisection( q3, q2, q1, q0, r3 );
       r3 *= scale;
 
       /*

diff --git a/src_tests/check_3_quartic.cc b/src_tests/check_3_quartic.cc
@@ -115,6 +115,17 @@ static double quartic19[] = {
   1.0 * ( -23.818259047289210 )
 };
 
+// Emanuele Siego example #2
+static double rootQuarticReal20[] = { 0.14271266144792,   0.14271266144792,  0.36842201250793, 0.36842202205209 };
+static double rootQuarticImag20[] = { 0.060284865351686, -0.060284865351686, 0, 0 };
+static double quartic20[] = {
+  5.0 * ( -4.1254841444253390e+02 ),
+  4.0 * ( 5.9401042551598744e+02 ),
+  3.0 * ( -3.1743361677460911e+02 ),
+  2.0 * ( 7.7477065082592929e+01 ),
+  1.0 * ( -8.9469482641160880e+00 )
+};
+
 #define DO_TEST( N ) \
   cout << "\n\nText N." << N << '\n'; \
   do_test( quartic##N, rootQuarticReal##N, rootQuarticImag##N )
@@ -143,6 +154,7 @@ do_test(
 int
 main() {
   cout.precision(14);
+  #if 1
   DO_TEST(0);
   DO_TEST(1);
   DO_TEST(2);
@@ -163,6 +175,8 @@ main() {
   DO_TEST(17);
   DO_TEST(18);
   DO_TEST(19);
+  #endif
+  DO_TEST(20);
   cout << "\n\nALL DONE!\n";
   return 0;
 }