1 //package com.polytechnik.utils;

     2 /*

     3  * (C) Vladislav Malyshkin 2010

     4  * This file is under GPL version 3.

     5  *

     6  */

     7 

     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  */

    12 

    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 PolynomialRoot

    19 */

    20 

    21 public class PolynomialRoot  {

    22 

    23 

    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 }

    50 

    51 

    52 

    53 static final double SQRT3=Math.sqrt(3.0),SQRT2=Math.sqrt(2.0);

    54 

    55 

    56 private static final boolean PRINT_DEBUG=false;

    57 

    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   }

    69 

    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);

    75 

    76 

    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;

    81 

    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   }

   124 

   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   }

   131 

   132   if(PRINT_DEBUG) System.err.println("G Path c="+c+" d="+d);

   133   final double [] u=new double[3];

   134 

   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;

   147 

   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);

   157 

   158       final double Q = qcub / 9;

   159       final double R = rcub / 54;

   160 

   161       final double Q3 = Q * Q * Q;

   162       final double R2 = R * R;

   163 

   164       final double CR2 = 729 * rcub * rcub;

   165       final double CQ3 = 2916 * qcub * qcub * qcub;

   166 

   167       if(PRINT_DEBUG) System.err.println("CR2="+CR2+" CQ3="+CQ3+" R="+R+" Q="+Q);

   168 

   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);

   199 

   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;

   209 

   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 */

   217 

   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    */

   225 

   226 

   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   }

   331 

   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   }

   352 

   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);

   355 

   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;

   364 

   365   return 4;

   366 }

   367 

   368 

   369 

   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     }

   407 

   408 

   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     }

   444 

   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     }

   455 

   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     }

   465 

   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         }

   492 

   493         findPolynomialRoots(n,p,rex,imx);

   494 

   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     }

   505 

   506 

   507 

   508 

   509     static final double EPS=0;

   510 

   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     }

   522 

   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     }

   560 

   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);

   575 

   576   final double Q = q / 9;

   577   final double R = r / 54;

   578 

   579   final double Q3 = Q * Q * Q;

   580   final double R2 = R * R;

   581 

   582   final double CR2 = 729 * r * r;

   583   final double CQ3 = 2916 * q * q * q;

   584   final double a3=a/3;

   585 

   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 */

   596 

   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. */

   600 

   601       final double sqrtQ = Math.sqrt (Q);

   602 

   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;

   626 

   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;

   638 

   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   }

   647 

   648     }

   649 

   650 

   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     }

   696 

   697 

   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     }

   755 

   756 

   757 

   758     public static void main(final String [] args)

   759     {

   760   final long t0=System.currentTimeMillis();

   761   final double eps=1e-6;

   762   //checkRoots();

   763   final java.util.Random r=new java.util.Random(-1381923);

   764   printSpecialValues();

   765 

   766   final int n_tests=10000000;

   767   //testRoots(2,n_tests,r,eps);

   768   //testRoots(3,n_tests,r,eps);

   769   testRoots(4,n_tests,r,eps);

   770   final long t1=System.currentTimeMillis();

   771   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 $");

   772     }

   773 

   774 

   775 

   776 }