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