[Toon-members] tag/src quartic.cpp

[Ethan Eade]

CVSROOT:        /cvsroot/toon
Module name:    tag
Changes by:     Ethan Eade <ethaneade>  08/04/01 02:40:53

Modified files:
        src            : quartic.cpp 

Log message:
        Addressed unlikely but possible boundary case in quartic solver.

CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/tag/src/quartic.cpp?cvsroot=toon&r1=1.1&r2=1.2

Patches:
Index: quartic.cpp
===================================================================
RCS file: /cvsroot/toon/tag/src/quartic.cpp,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -b -r1.1 -r1.2
--- quartic.cpp 31 Mar 2008 16:35:41 -0000      1.1
+++ quartic.cpp 1 Apr 2008 02:40:52 -0000       1.2
@@ -1,4 +1,4 @@
-#include <util/quartic.h>
+#include <tag/quartic.h>
 
 #include <cmath>
 #include <cassert>
@@ -46,33 +46,54 @@
     double half_B = 0.5 * B;
     double beta = half_B *(half_B*half_B - C) + D;
     double q_B = 0.25 * B;
-    double gamma = q_B * (q_B *(C - 3 * q_B * q_B) - D) + E;
+    double four_gamma = B * (q_B *(C - 3 * q_B * q_B) - D) + 4*E;
+    
+    if (beta*beta < 1e-18) {
+       double disc = alpha*alpha - four_gamma;
+       if (disc < 0)
+           return 0;
+       double root_disc = sqrt(disc);
+       {
+           double disc2 = -alpha + root_disc;
+           if (disc2 < 0)
+               return 0;
+           double root_disc2 = sqrt(0.5 * disc2);
+           r[0] = -q_B - root_disc2;
+           r[1] = -q_B + root_disc2;
+       }
+       {
+           double disc2 = -alpha - root_disc;
+           if (disc2 < 0)
+               return 2;
+           double root_disc2 = sqrt(0.5 * disc2);
+           r[2] = -q_B - root_disc2;
+           r[3] = -q_B + root_disc2;
+       }
+       return 4;
+    }
+
     double third_alpha = alpha * third;
-    double P = -alpha * third_alpha * 0.25 - gamma;
-    double Q = third_alpha * (-third_alpha * third_alpha * 0.25 + gamma) - 
beta*beta*0.125;
+    double P = -0.25*(alpha * third_alpha + four_gamma);
+    double Q = third_alpha * (0.25*(four_gamma - third_alpha * third_alpha)) - 
beta*beta*0.125;
 
     double dcr[3];
     int ndcr = depressed_cubic_real_roots(P,Q, dcr);
-    double y = dcr[ndcr-1] - third * 2.5 * alpha;
-
-    double disc2 = alpha + 2*y;
-    
+    double disc2 = 2*(dcr[ndcr-1] - third_alpha);
     double W = sqrt(disc2);
 
     double disc_base = 2*alpha + disc2; //3*alpha + 2*y;
     double disc_add = 2*beta / W;
-    double x_base = -0.25 * B;
 
     int count = 0;
     if (disc_base + disc_add <= 0) {
        double sr = sqrt(-disc_base-disc_add);
-       r[count++] = x_base + 0.5 * (W - sr);
-       r[count++] = x_base + 0.5 * (W + sr);
+       r[count++] = -q_B + 0.5 * (W - sr);
+       r[count++] = -q_B + 0.5 * (W + sr);
     }
     if (disc_base - disc_add <= 0) {
        double sr = sqrt(-disc_base + disc_add);
-       r[count++] = x_base - 0.5 * (W + sr);   
-       r[count++] = x_base - 0.5 * (W - sr);
+       r[count++] = -q_B - 0.5 * (W + sr);     
+       r[count++] = -q_B - 0.5 * (W - sr);
     }
     return count;
 }