1 /*
   2  * (C) Vladislav Malyshkin 2010
   3  * This file is under GPL version 3.
   4  *
   5  */
   6 
   7 /** Polynomial root.
   8  *  @version $Id: PolynomialRoot.java,v 1.105 2012/08/18 00:00:05 mal Exp $
   9  *  @author Vladislav Malyshkin mal@gromco.com
  10  */
  11 
  12 /**
  13  * @test
  14  * @bug 8005956
  15  * @summary C2: assert(!def_outside->member(r)) failed: Use of external LRG overlaps the same LRG defined in this block
  16  * @library /test/lib
  17  * @modules java.base/jdk.internal.misc
  18  *          java.management
  19  *
  20  * @run main/timeout=300 compiler.c2.PolynomialRoot
  21  */
  22 
  23 package compiler.c2;
  24 
  25 import jdk.test.lib.Utils;
  26 
  27 import java.util.Arrays;
  28 import java.util.Random;
  29 
  30 public class PolynomialRoot  {
  31 
  32 
  33 public static int findPolynomialRoots(final int n,
  34               final double [] p,
  35               final double [] re_root,
  36               final double [] im_root)
  37 {
  38     if(n==4)
  39     {
  40   return root4(p,re_root,im_root);
  41     }
  42     else if(n==3)
  43     {
  44   return root3(p,re_root,im_root);
  45     }
  46     else if(n==2)
  47     {
  48   return root2(p,re_root,im_root);
  49     }
  50     else if(n==1)
  51     {
  52   return root1(p,re_root,im_root);
  53     }
  54     else
  55     {
  56   throw new RuntimeException("n="+n+" is not supported yet");
  57     }
  58 }
  59 
  60 
  61 
  62 static final double SQRT3=Math.sqrt(3.0),SQRT2=Math.sqrt(2.0);
  63 
  64 
  65 private static final boolean PRINT_DEBUG=false;
  66 
  67 public static int root4(final double [] p,final double [] re_root,final double [] im_root)
  68 {
  69   if (PRINT_DEBUG) { System.err.println("=====================root4:p=" + Arrays.toString(p)); }
  70   final double vs=p[4];
  71   if(PRINT_DEBUG) System.err.println("p[4]="+p[4]);
  72   if(!(Math.abs(vs)>EPS))
  73   {
  74       re_root[0]=re_root[1]=re_root[2]=re_root[3]=
  75     im_root[0]=im_root[1]=im_root[2]=im_root[3]=Double.NaN;
  76       return -1;
  77   }
  78 
  79 /* zsolve_quartic.c - finds the complex roots of
  80  *  x^4 + a x^3 + b x^2 + c x + d = 0
  81  */
  82   final double a=p[3]/vs,b=p[2]/vs,c=p[1]/vs,d=p[0]/vs;
  83   if(PRINT_DEBUG) System.err.println("input a="+a+" b="+b+" c="+c+" d="+d);
  84 
  85 
  86   final double r4 = 1.0 / 4.0;
  87   final double q2 = 1.0 / 2.0, q4 = 1.0 / 4.0, q8 = 1.0 / 8.0;
  88   final double q1 = 3.0 / 8.0, q3 = 3.0 / 16.0;
  89   final int mt;
  90 
  91   /* Deal easily with the cases where the quartic is degenerate. The
  92    * ordering of solutions is done explicitly. */
  93   if (0 == b && 0 == c)
  94   {
  95       if (0 == d)
  96       {
  97     re_root[0]=-a;
  98     im_root[0]=im_root[1]=im_root[2]=im_root[3]=0;
  99     re_root[1]=re_root[2]=re_root[3]=0;
 100     return 4;
 101       }
 102       else if (0 == a)
 103       {
 104     if (d > 0)
 105     {
 106         final double sq4 = Math.sqrt(Math.sqrt(d));
 107         re_root[0]=sq4*SQRT2/2;
 108         im_root[0]=re_root[0];
 109         re_root[1]=-re_root[0];
 110         im_root[1]=re_root[0];
 111         re_root[2]=-re_root[0];
 112         im_root[2]=-re_root[0];
 113         re_root[3]=re_root[0];
 114         im_root[3]=-re_root[0];
 115         if(PRINT_DEBUG) System.err.println("Path a=0 d>0");
 116     }
 117     else
 118     {
 119         final double sq4 = Math.sqrt(Math.sqrt(-d));
 120         re_root[0]=sq4;
 121         im_root[0]=0;
 122         re_root[1]=0;
 123         im_root[1]=sq4;
 124         re_root[2]=0;
 125         im_root[2]=-sq4;
 126         re_root[3]=-sq4;
 127         im_root[3]=0;
 128         if(PRINT_DEBUG) System.err.println("Path a=0 d<0");
 129     }
 130     return 4;
 131       }
 132   }
 133 
 134   if (0.0 == c && 0.0 == d)
 135   {
 136       root2(new double []{p[2],p[3],p[4]},re_root,im_root);
 137       re_root[2]=im_root[2]=re_root[3]=im_root[3]=0;
 138       return 4;
 139   }
 140 
 141   if(PRINT_DEBUG) System.err.println("G Path c="+c+" d="+d);
 142   final double [] u=new double[3];
 143 
 144   if(PRINT_DEBUG) System.err.println("Generic Path");
 145   /* For non-degenerate solutions, proceed by constructing and
 146    * solving the resolvent cubic */
 147   final double aa = a * a;
 148   final double pp = b - q1 * aa;
 149   final double qq = c - q2 * a * (b - q4 * aa);
 150   final double rr = d - q4 * a * (c - q4 * a * (b - q3 * aa));
 151   final double rc = q2 * pp , rc3 = rc / 3;
 152   final double sc = q4 * (q4 * pp * pp - rr);
 153   final double tc = -(q8 * qq * q8 * qq);
 154   if(PRINT_DEBUG) System.err.println("aa="+aa+" pp="+pp+" qq="+qq+" rr="+rr+" rc="+rc+" sc="+sc+" tc="+tc);
 155   final boolean flag_realroots;
 156 
 157   /* This code solves the resolvent cubic in a convenient fashion
 158    * for this implementation of the quartic. If there are three real
 159    * roots, then they are placed directly into u[].  If two are
 160    * complex, then the real root is put into u[0] and the real
 161    * and imaginary part of the complex roots are placed into
 162    * u[1] and u[2], respectively. */
 163   {
 164       final double qcub = (rc * rc - 3 * sc);
 165       final double rcub = (rc*(2 * rc * rc - 9 * sc) + 27 * tc);
 166 
 167       final double Q = qcub / 9;
 168       final double R = rcub / 54;
 169 
 170       final double Q3 = Q * Q * Q;
 171       final double R2 = R * R;
 172 
 173       final double CR2 = 729 * rcub * rcub;
 174       final double CQ3 = 2916 * qcub * qcub * qcub;
 175 
 176       if(PRINT_DEBUG) System.err.println("CR2="+CR2+" CQ3="+CQ3+" R="+R+" Q="+Q);
 177 
 178       if (0 == R && 0 == Q)
 179       {
 180     flag_realroots=true;
 181     u[0] = -rc3;
 182     u[1] = -rc3;
 183     u[2] = -rc3;
 184       }
 185       else if (CR2 == CQ3)
 186       {
 187     flag_realroots=true;
 188     final double sqrtQ = Math.sqrt (Q);
 189     if (R > 0)
 190     {
 191         u[0] = -2 * sqrtQ - rc3;
 192         u[1] = sqrtQ - rc3;
 193         u[2] = sqrtQ - rc3;
 194     }
 195     else
 196     {
 197         u[0] = -sqrtQ - rc3;
 198         u[1] = -sqrtQ - rc3;
 199         u[2] = 2 * sqrtQ - rc3;
 200     }
 201       }
 202       else if (R2 < Q3)
 203       {
 204     flag_realroots=true;
 205     final double ratio = (R >= 0?1:-1) * Math.sqrt (R2 / Q3);
 206     final double theta = Math.acos (ratio);
 207     final double norm = -2 * Math.sqrt (Q);
 208 
 209     u[0] = norm * Math.cos (theta / 3) - rc3;
 210     u[1] = norm * Math.cos ((theta + 2.0 * Math.PI) / 3) - rc3;
 211     u[2] = norm * Math.cos ((theta - 2.0 * Math.PI) / 3) - rc3;
 212       }
 213       else
 214       {
 215     flag_realroots=false;
 216     final double A = -(R >= 0?1:-1)*Math.pow(Math.abs(R)+Math.sqrt(R2-Q3),1.0/3.0);
 217     final double B = Q / A;
 218 
 219     u[0] = A + B - rc3;
 220     u[1] = -0.5 * (A + B) - rc3;
 221     u[2] = -(SQRT3*0.5) * Math.abs (A - B);
 222       }
 223       if(PRINT_DEBUG) System.err.println("u[0]="+u[0]+" u[1]="+u[1]+" u[2]="+u[2]+" qq="+qq+" disc="+((CR2 - CQ3) / 2125764.0));
 224   }
 225   /* End of solution to resolvent cubic */
 226 
 227   /* Combine the square roots of the roots of the cubic
 228    * resolvent appropriately. Also, calculate 'mt' which
 229    * designates the nature of the roots:
 230    * mt=1 : 4 real roots
 231    * mt=2 : 0 real roots
 232    * mt=3 : 2 real roots
 233    */
 234 
 235 
 236   final double w1_re,w1_im,w2_re,w2_im,w3_re,w3_im,mod_w1w2,mod_w1w2_squared;
 237   if (flag_realroots)
 238   {
 239       mod_w1w2=-1;
 240       mt = 2;
 241       int jmin=0;
 242       double vmin=Math.abs(u[jmin]);
 243       for(int j=1;j<3;j++)
 244       {
 245     final double vx=Math.abs(u[j]);
 246     if(vx<vmin)
 247     {
 248         vmin=vx;
 249         jmin=j;
 250     }
 251       }
 252       final double u1=u[(jmin+1)%3],u2=u[(jmin+2)%3];
 253       mod_w1w2_squared=Math.abs(u1*u2);
 254       if(u1>=0)
 255       {
 256     w1_re=Math.sqrt(u1);
 257     w1_im=0;
 258       }
 259       else
 260       {
 261     w1_re=0;
 262     w1_im=Math.sqrt(-u1);
 263       }
 264       if(u2>=0)
 265       {
 266     w2_re=Math.sqrt(u2);
 267     w2_im=0;
 268       }
 269       else
 270       {
 271     w2_re=0;
 272     w2_im=Math.sqrt(-u2);
 273       }
 274       if(PRINT_DEBUG) System.err.println("u1="+u1+" u2="+u2+" jmin="+jmin);
 275   }
 276   else
 277   {
 278       mt = 3;
 279       final double w_mod2_sq=u[1]*u[1]+u[2]*u[2],w_mod2=Math.sqrt(w_mod2_sq),w_mod=Math.sqrt(w_mod2);
 280       if(w_mod2_sq<=0)
 281       {
 282     w1_re=w1_im=0;
 283       }
 284       else
 285       {
 286     // calculate square root of a complex number (u[1],u[2])
 287     // the result is in the (w1_re,w1_im)
 288     final double absu1=Math.abs(u[1]),absu2=Math.abs(u[2]),w;
 289     if(absu1>=absu2)
 290     {
 291         final double t=absu2/absu1;
 292         w=Math.sqrt(absu1*0.5 * (1.0 + Math.sqrt(1.0 + t * t)));
 293         if(PRINT_DEBUG) System.err.println(" Path1 ");
 294     }
 295     else
 296     {
 297         final double t=absu1/absu2;
 298         w=Math.sqrt(absu2*0.5 * (t + Math.sqrt(1.0 + t * t)));
 299         if(PRINT_DEBUG) System.err.println(" Path1a ");
 300     }
 301     if(u[1]>=0)
 302     {
 303         w1_re=w;
 304         w1_im=u[2]/(2*w);
 305         if(PRINT_DEBUG) System.err.println(" Path2 ");
 306     }
 307     else
 308     {
 309         final double vi = (u[2] >= 0) ? w : -w;
 310         w1_re=u[2]/(2*vi);
 311         w1_im=vi;
 312         if(PRINT_DEBUG) System.err.println(" Path2a ");
 313     }
 314       }
 315       final double absu0=Math.abs(u[0]);
 316       if(w_mod2>=absu0)
 317       {
 318     mod_w1w2=w_mod2;
 319     mod_w1w2_squared=w_mod2_sq;
 320     w2_re=w1_re;
 321     w2_im=-w1_im;
 322       }
 323       else
 324       {
 325     mod_w1w2=-1;
 326     mod_w1w2_squared=w_mod2*absu0;
 327     if(u[0]>=0)
 328     {
 329         w2_re=Math.sqrt(absu0);
 330         w2_im=0;
 331     }
 332     else
 333     {
 334         w2_re=0;
 335         w2_im=Math.sqrt(absu0);
 336     }
 337       }
 338       if(PRINT_DEBUG) System.err.println("u[0]="+u[0]+"u[1]="+u[1]+" u[2]="+u[2]+" absu0="+absu0+" w_mod="+w_mod+" w_mod2="+w_mod2);
 339   }
 340 
 341   /* Solve the quadratic in order to obtain the roots
 342    * to the quartic */
 343   if(mod_w1w2>0)
 344   {
 345       // a shorcut to reduce rounding error
 346       w3_re=qq/(-8)/mod_w1w2;
 347       w3_im=0;
 348   }
 349   else if(mod_w1w2_squared>0)
 350   {
 351       // regular path
 352       final double mqq8n=qq/(-8)/mod_w1w2_squared;
 353       w3_re=mqq8n*(w1_re*w2_re-w1_im*w2_im);
 354       w3_im=-mqq8n*(w1_re*w2_im+w2_re*w1_im);
 355   }
 356   else
 357   {
 358       // typically occur when qq==0
 359       w3_re=w3_im=0;
 360   }
 361 
 362   final double h = r4 * a;
 363   if(PRINT_DEBUG) System.err.println("w1_re="+w1_re+" w1_im="+w1_im+" w2_re="+w2_re+" w2_im="+w2_im+" w3_re="+w3_re+" w3_im="+w3_im+" h="+h);
 364 
 365   re_root[0]=w1_re+w2_re+w3_re-h;
 366   im_root[0]=w1_im+w2_im+w3_im;
 367   re_root[1]=-(w1_re+w2_re)+w3_re-h;
 368   im_root[1]=-(w1_im+w2_im)+w3_im;
 369   re_root[2]=w2_re-w1_re-w3_re-h;
 370   im_root[2]=w2_im-w1_im-w3_im;
 371   re_root[3]=w1_re-w2_re-w3_re-h;
 372   im_root[3]=w1_im-w2_im-w3_im;
 373 
 374   return 4;
 375 }
 376 
 377 
 378 
 379     static void setRandomP(final double [] p, final int n, Random r)
 380     {
 381   if(r.nextDouble()<0.1)
 382   {
 383       // integer coefficiens
 384       for(int j=0;j<p.length;j++)
 385       {
 386     if(j<=n)
 387     {
 388         p[j]=(r.nextInt(2)<=0?-1:1)*r.nextInt(10);
 389     }
 390     else
 391     {
 392         p[j]=0;
 393     }
 394       }
 395   }
 396   else
 397   {
 398       // real coefficiens
 399       for(int j=0;j<p.length;j++)
 400       {
 401     if(j<=n)
 402     {
 403         p[j]=-1+2*r.nextDouble();
 404     }
 405     else
 406     {
 407         p[j]=0;
 408     }
 409       }
 410   }
 411   if(Math.abs(p[n])<1e-2)
 412   {
 413       p[n]=(r.nextInt(2)<=0?-1:1)*(0.1+r.nextDouble());
 414   }
 415     }
 416 
 417 
 418     static void checkValues(final double [] p,
 419           final int n,
 420           final double rex,
 421           final double imx,
 422           final double eps,
 423           final String txt)
 424     {
 425   double res=0,ims=0,sabs=0;
 426   final double xabs=Math.abs(rex)+Math.abs(imx);
 427   for(int k=n;k>=0;k--)
 428   {
 429       final double res1=(res*rex-ims*imx)+p[k];
 430       final double ims1=(ims*rex+res*imx);
 431       res=res1;
 432       ims=ims1;
 433       sabs+=xabs*sabs+p[k];
 434   }
 435   sabs=Math.abs(sabs);
 436   if(false && sabs>1/eps?
 437      (!(Math.abs(res/sabs)<=eps)||!(Math.abs(ims/sabs)<=eps))
 438      :
 439      (!(Math.abs(res)<=eps)||!(Math.abs(ims)<=eps)))
 440   {
 441       throw new RuntimeException(
 442     getPolinomTXT(p)+"\n"+
 443     "\t x.r="+rex+" x.i="+imx+"\n"+
 444     "res/sabs="+(res/sabs)+" ims/sabs="+(ims/sabs)+
 445     " sabs="+sabs+
 446     "\nres="+res+" ims="+ims+" n="+n+" eps="+eps+" "+
 447     " sabs>1/eps="+(sabs>1/eps)+
 448     " f1="+(!(Math.abs(res/sabs)<=eps)||!(Math.abs(ims/sabs)<=eps))+
 449     " f2="+(!(Math.abs(res)<=eps)||!(Math.abs(ims)<=eps))+
 450     " "+txt);
 451   }
 452     }
 453 
 454     static String getPolinomTXT(final double [] p)
 455     {
 456   final StringBuilder buf=new StringBuilder();
 457   buf.append("order="+(p.length-1)+"\t");
 458   for(int k=0;k<p.length;k++)
 459   {
 460       buf.append("p["+k+"]="+p[k]+";");
 461   }
 462   return buf.toString();
 463     }
 464 
 465     static String getRootsTXT(int nr,final double [] re,final double [] im)
 466     {
 467   final StringBuilder buf=new StringBuilder();
 468   for(int k=0;k<nr;k++)
 469   {
 470       buf.append("x."+k+"("+re[k]+","+im[k]+")\n");
 471   }
 472   return buf.toString();
 473     }
 474 
 475     static void testRoots(final int n,
 476         final int n_tests,
 477         final Random rn,
 478         final double eps)
 479     {
 480   final double [] p=new double [n+1];
 481   final double [] rex=new double [n],imx=new double [n];
 482   for(int i=0;i<n_tests;i++)
 483   {
 484     for(int dg=n;dg-->-1;)
 485     {
 486       for(int dr=3;dr-->0;)
 487       {
 488         setRandomP(p,n,rn);
 489         for(int j=0;j<=dg;j++)
 490         {
 491       p[j]=0;
 492         }
 493         if(dr==0)
 494         {
 495       p[0]=-1+2.0*rn.nextDouble();
 496         }
 497         else if(dr==1)
 498         {
 499       p[0]=p[1]=0;
 500         }
 501 
 502         findPolynomialRoots(n,p,rex,imx);
 503 
 504         for(int j=0;j<n;j++)
 505         {
 506       //System.err.println("j="+j);
 507       checkValues(p,n,rex[j],imx[j],eps," t="+i);
 508         }
 509       }
 510     }
 511   }
 512   System.err.println("testRoots(): n_tests="+n_tests+" OK, dim="+n);
 513     }
 514 
 515 
 516 
 517 
 518     static final double EPS=0;
 519 
 520     public static int root1(final double [] p,final double [] re_root,final double [] im_root)
 521     {
 522   if(!(Math.abs(p[1])>EPS))
 523   {
 524       re_root[0]=im_root[0]=Double.NaN;
 525       return -1;
 526   }
 527   re_root[0]=-p[0]/p[1];
 528   im_root[0]=0;
 529   return 1;
 530     }
 531 
 532     public static int root2(final double [] p,final double [] re_root,final double [] im_root)
 533     {
 534   if(!(Math.abs(p[2])>EPS))
 535   {
 536       re_root[0]=re_root[1]=im_root[0]=im_root[1]=Double.NaN;
 537       return -1;
 538   }
 539   final double b2=0.5*(p[1]/p[2]),c=p[0]/p[2],d=b2*b2-c;
 540   if(d>=0)
 541   {
 542       final double sq=Math.sqrt(d);
 543       if(b2<0)
 544       {
 545     re_root[1]=-b2+sq;
 546     re_root[0]=c/re_root[1];
 547       }
 548       else if(b2>0)
 549       {
 550     re_root[0]=-b2-sq;
 551     re_root[1]=c/re_root[0];
 552       }
 553       else
 554       {
 555     re_root[0]=-b2-sq;
 556     re_root[1]=-b2+sq;
 557       }
 558       im_root[0]=im_root[1]=0;
 559   }
 560   else
 561   {
 562       final double sq=Math.sqrt(-d);
 563       re_root[0]=re_root[1]=-b2;
 564       im_root[0]=sq;
 565       im_root[1]=-sq;
 566   }
 567   return 2;
 568     }
 569 
 570     public static int root3(final double [] p,final double [] re_root,final double [] im_root)
 571     {
 572   final double vs=p[3];
 573   if(!(Math.abs(vs)>EPS))
 574   {
 575       re_root[0]=re_root[1]=re_root[2]=
 576     im_root[0]=im_root[1]=im_root[2]=Double.NaN;
 577       return -1;
 578   }
 579   final double a=p[2]/vs,b=p[1]/vs,c=p[0]/vs;
 580   /* zsolve_cubic.c - finds the complex roots of x^3 + a x^2 + b x + c = 0
 581    */
 582   final double q = (a * a - 3 * b);
 583   final double r = (a*(2 * a * a - 9 * b) + 27 * c);
 584 
 585   final double Q = q / 9;
 586   final double R = r / 54;
 587 
 588   final double Q3 = Q * Q * Q;
 589   final double R2 = R * R;
 590 
 591   final double CR2 = 729 * r * r;
 592   final double CQ3 = 2916 * q * q * q;
 593   final double a3=a/3;
 594 
 595   if (R == 0 && Q == 0)
 596   {
 597       re_root[0]=re_root[1]=re_root[2]=-a3;
 598       im_root[0]=im_root[1]=im_root[2]=0;
 599       return 3;
 600   }
 601   else if (CR2 == CQ3)
 602   {
 603       /* this test is actually R2 == Q3, written in a form suitable
 604          for exact computation with integers */
 605 
 606       /* Due to finite precision some double roots may be missed, and
 607          will be considered to be a pair of complex roots z = x +/-
 608          epsilon i close to the real axis. */
 609 
 610       final double sqrtQ = Math.sqrt (Q);
 611 
 612       if (R > 0)
 613       {
 614     re_root[0] = -2 * sqrtQ - a3;
 615     re_root[1]=re_root[2]=sqrtQ - a3;
 616     im_root[0]=im_root[1]=im_root[2]=0;
 617       }
 618       else
 619       {
 620     re_root[0]=re_root[1] = -sqrtQ - a3;
 621     re_root[2]=2 * sqrtQ - a3;
 622     im_root[0]=im_root[1]=im_root[2]=0;
 623       }
 624       return 3;
 625   }
 626   else if (R2 < Q3)
 627   {
 628       final double sgnR = (R >= 0 ? 1 : -1);
 629       final double ratio = sgnR * Math.sqrt (R2 / Q3);
 630       final double theta = Math.acos (ratio);
 631       final double norm = -2 * Math.sqrt (Q);
 632       final double r0 = norm * Math.cos (theta/3) - a3;
 633       final double r1 = norm * Math.cos ((theta + 2.0 * Math.PI) / 3) - a3;
 634       final double r2 = norm * Math.cos ((theta - 2.0 * Math.PI) / 3) - a3;
 635 
 636       re_root[0]=r0;
 637       re_root[1]=r1;
 638       re_root[2]=r2;
 639       im_root[0]=im_root[1]=im_root[2]=0;
 640       return 3;
 641   }
 642   else
 643   {
 644       final double sgnR = (R >= 0 ? 1 : -1);
 645       final double A = -sgnR * Math.pow (Math.abs (R) + Math.sqrt (R2 - Q3), 1.0 / 3.0);
 646       final double B = Q / A;
 647 
 648       re_root[0]=A + B - a3;
 649       im_root[0]=0;
 650       re_root[1]=-0.5 * (A + B) - a3;
 651       im_root[1]=-(SQRT3*0.5) * Math.abs(A - B);
 652       re_root[2]=re_root[1];
 653       im_root[2]=-im_root[1];
 654       return 3;
 655   }
 656 
 657     }
 658 
 659 
 660     static void root3a(final double [] p,final double [] re_root,final double [] im_root)
 661     {
 662   if(Math.abs(p[3])>EPS)
 663   {
 664       final double v=p[3],
 665     a=p[2]/v,b=p[1]/v,c=p[0]/v,
 666     a3=a/3,a3a=a3*a,
 667     pd3=(b-a3a)/3,
 668     qd2=a3*(a3a/3-0.5*b)+0.5*c,
 669     Q=pd3*pd3*pd3+qd2*qd2;
 670       if(Q<0)
 671       {
 672     // three real roots
 673     final double SQ=Math.sqrt(-Q);
 674     final double th=Math.atan2(SQ,-qd2);
 675     im_root[0]=im_root[1]=im_root[2]=0;
 676     final double f=2*Math.sqrt(-pd3);
 677     re_root[0]=f*Math.cos(th/3)-a3;
 678     re_root[1]=f*Math.cos((th+2*Math.PI)/3)-a3;
 679     re_root[2]=f*Math.cos((th+4*Math.PI)/3)-a3;
 680     //System.err.println("3r");
 681       }
 682       else
 683       {
 684     // one real & two complex roots
 685     final double SQ=Math.sqrt(Q);
 686     final double r1=-qd2+SQ,r2=-qd2-SQ;
 687     final double v1=Math.signum(r1)*Math.pow(Math.abs(r1),1.0/3),
 688         v2=Math.signum(r2)*Math.pow(Math.abs(r2),1.0/3),
 689         sv=v1+v2;
 690     // real root
 691     re_root[0]=sv-a3;
 692     im_root[0]=0;
 693     // complex roots
 694     re_root[1]=re_root[2]=-0.5*sv-a3;
 695     im_root[1]=(v1-v2)*(SQRT3*0.5);
 696     im_root[2]=-im_root[1];
 697     //System.err.println("1r2c");
 698       }
 699   }
 700   else
 701   {
 702       re_root[0]=re_root[1]=re_root[2]=im_root[0]=im_root[1]=im_root[2]=Double.NaN;
 703   }
 704     }
 705 
 706 
 707     static void printSpecialValues()
 708     {
 709   for(int st=0;st<6;st++)
 710   {
 711       //final double [] p=new double []{8,1,3,3.6,1};
 712       final double [] re_root=new double [4],im_root=new double [4];
 713       final double [] p;
 714       final int n;
 715       if(st<=3)
 716       {
 717     if(st<=0)
 718     {
 719         p=new double []{2,-4,6,-4,1};
 720         //p=new double []{-6,6,-6,8,-2};
 721     }
 722     else if(st==1)
 723     {
 724         p=new double []{0,-4,8,3,-9};
 725     }
 726     else if(st==2)
 727     {
 728         p=new double []{-1,0,2,0,-1};
 729     }
 730     else
 731     {
 732         p=new double []{-5,2,8,-2,-3};
 733     }
 734     root4(p,re_root,im_root);
 735     n=4;
 736       }
 737       else
 738       {
 739     p=new double []{0,2,0,1};
 740     if(st==4)
 741     {
 742         p[1]=-p[1];
 743     }
 744     root3(p,re_root,im_root);
 745     n=3;
 746       }
 747       System.err.println("======== n="+n);
 748       for(int i=0;i<=n;i++)
 749       {
 750     if(i<n)
 751     {
 752         System.err.println(String.valueOf(i)+"\t"+
 753                p[i]+"\t"+
 754                re_root[i]+"\t"+
 755                im_root[i]);
 756     }
 757     else
 758     {
 759         System.err.println(String.valueOf(i)+"\t"+p[i]+"\t");
 760     }
 761       }
 762   }
 763     }
 764 
 765 
 766 
 767     public static void main(final String [] args)
 768     {
 769       if (System.getProperty("os.arch").equals("x86") ||
 770          System.getProperty("os.arch").equals("amd64") ||
 771          System.getProperty("os.arch").equals("x86_64")){
 772         final long t0=System.currentTimeMillis();
 773         final double eps=1e-6;
 774         //checkRoots();
 775         final Random r = Utils.getRandomInstance();
 776         printSpecialValues();
 777 
 778         final int n_tests=100000;
 779         //testRoots(2,n_tests,r,eps);
 780         //testRoots(3,n_tests,r,eps);
 781         testRoots(4,n_tests,r,eps);
 782         final long t1=System.currentTimeMillis();
 783         System.err.println("PolynomialRoot.main: "+n_tests+" tests OK done in "+(t1-t0)+" milliseconds. ver=$Id: PolynomialRoot.java,v 1.105 2012/08/18 00:00:05 mal Exp $");
 784         System.out.println("PASSED");
 785      } else {
 786        System.out.println("PASS test for non-x86");
 787      }
 788    }
 789 
 790 
 791 
 792 }