--- zsolve_cubic.c	2008-11-19 10:13:47.000000000 +0100
+++ zsolve_cubic.c.new	2009-04-28 11:55:25.000000000 +0200
@@ -41,8 +41,8 @@
   double Q3 = Q * Q * Q;
   double R2 = R * R;
 
-  double CR2 = 729 * r * r;
-  double CQ3 = 2916 * q * q * q;
+//  double CR2 = 729 * r * r;
+//  double CQ3 = 2916 * q * q * q;
 
   if (R == 0 && Q == 0)
     {
@@ -54,7 +54,7 @@
       GSL_IMAG (*z2) = 0;
       return 3;
     }
-  else if (CR2 == CQ3) 
+  else if (R2 == Q3) 
     {
       /* this test is actually R2 == Q3, written in a form suitable
          for exact computation with integers */
@@ -85,11 +85,17 @@
         }
       return 3;
     }
-  else if (CR2 < CQ3)  /* equivalent to R2 < Q3 */
+  else if (R2 < Q3)  /* equivalent to R2 < Q3 */
     {
       double sqrtQ = sqrt (Q);
       double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ;
-      double theta = acos (R / sqrtQ3);
+      double ctheta = R / sqrtQ3;
+      double theta = 0; 
+      if ( ctheta <= -1.0) 
+         theta = M_PI;
+      else if ( ctheta < 1.0) 
+         theta = acos (R / sqrtQ3);
+      
       double norm = -2 * sqrtQ;
       double r0 = norm * cos (theta / 3) - a / 3;
       double r1 = norm * cos ((theta + 2.0 * M_PI) / 3) - a / 3;