test/compiler/8005956/PolynomialRoot.java

Mon, 28 Jul 2014 15:06:38 -0700

author
fzhinkin
date
Mon, 28 Jul 2014 15:06:38 -0700
changeset 6997
dbb05f6d93c4
parent 5380
e554162ab094
child 6876
710a3c8b516e
permissions
-rw-r--r--

8051344: JVM crashed in Compile::start() during method parsing w/ UseRTMDeopt turned on
Summary: call rtm_deopt() only if there were no compilation bailouts before.
Reviewed-by: kvn

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

mercurial