/* * (C) Vladislav Malyshkin 2010 * This file is under GPL version 3. * */ /** Polynomial root. * @version $Id: PolynomialRoot.java,v 1.105 2012/08/18 00:00:05 mal Exp $ * @author Vladislav Malyshkin mal@gromco.com */ /** * @test * @bug 8005956 * @summary C2: assert(!def_outside->member(r)) failed: Use of external LRG overlaps the same LRG defined in this block * @library /testlibrary * @build com.oracle.java.testlibrary.* * @run main/timeout=300 PolynomialRoot */ import com.oracle.java.testlibrary.Utils; import java.util.Arrays; import java.util.Random; public class PolynomialRoot { public static int findPolynomialRoots(final int n, final double [] p, final double [] re_root, final double [] im_root) { if(n==4) { return root4(p,re_root,im_root); } else if(n==3) { return root3(p,re_root,im_root); } else if(n==2) { return root2(p,re_root,im_root); } else if(n==1) { return root1(p,re_root,im_root); } else { throw new RuntimeException("n="+n+" is not supported yet"); } } static final double SQRT3=Math.sqrt(3.0),SQRT2=Math.sqrt(2.0); private static final boolean PRINT_DEBUG=false; public static int root4(final double [] p,final double [] re_root,final double [] im_root) { if (PRINT_DEBUG) { System.err.println("=====================root4:p=" + Arrays.toString(p)); } final double vs=p[4]; if(PRINT_DEBUG) System.err.println("p[4]="+p[4]); if(!(Math.abs(vs)>EPS)) { re_root[0]=re_root[1]=re_root[2]=re_root[3]= im_root[0]=im_root[1]=im_root[2]=im_root[3]=Double.NaN; return -1; } /* zsolve_quartic.c - finds the complex roots of * x^4 + a x^3 + b x^2 + c x + d = 0 */ final double a=p[3]/vs,b=p[2]/vs,c=p[1]/vs,d=p[0]/vs; if(PRINT_DEBUG) System.err.println("input a="+a+" b="+b+" c="+c+" d="+d); final double r4 = 1.0 / 4.0; final double q2 = 1.0 / 2.0, q4 = 1.0 / 4.0, q8 = 1.0 / 8.0; final double q1 = 3.0 / 8.0, q3 = 3.0 / 16.0; final int mt; /* Deal easily with the cases where the quartic is degenerate. The * ordering of solutions is done explicitly. */ if (0 == b && 0 == c) { if (0 == d) { re_root[0]=-a; im_root[0]=im_root[1]=im_root[2]=im_root[3]=0; re_root[1]=re_root[2]=re_root[3]=0; return 4; } else if (0 == a) { if (d > 0) { final double sq4 = Math.sqrt(Math.sqrt(d)); re_root[0]=sq4*SQRT2/2; im_root[0]=re_root[0]; re_root[1]=-re_root[0]; im_root[1]=re_root[0]; re_root[2]=-re_root[0]; im_root[2]=-re_root[0]; re_root[3]=re_root[0]; im_root[3]=-re_root[0]; if(PRINT_DEBUG) System.err.println("Path a=0 d>0"); } else { final double sq4 = Math.sqrt(Math.sqrt(-d)); re_root[0]=sq4; im_root[0]=0; re_root[1]=0; im_root[1]=sq4; re_root[2]=0; im_root[2]=-sq4; re_root[3]=-sq4; im_root[3]=0; if(PRINT_DEBUG) System.err.println("Path a=0 d<0"); } return 4; } } if (0.0 == c && 0.0 == d) { root2(new double []{p[2],p[3],p[4]},re_root,im_root); re_root[2]=im_root[2]=re_root[3]=im_root[3]=0; return 4; } if(PRINT_DEBUG) System.err.println("G Path c="+c+" d="+d); final double [] u=new double[3]; if(PRINT_DEBUG) System.err.println("Generic Path"); /* For non-degenerate solutions, proceed by constructing and * solving the resolvent cubic */ final double aa = a * a; final double pp = b - q1 * aa; final double qq = c - q2 * a * (b - q4 * aa); final double rr = d - q4 * a * (c - q4 * a * (b - q3 * aa)); final double rc = q2 * pp , rc3 = rc / 3; final double sc = q4 * (q4 * pp * pp - rr); final double tc = -(q8 * qq * q8 * qq); if(PRINT_DEBUG) System.err.println("aa="+aa+" pp="+pp+" qq="+qq+" rr="+rr+" rc="+rc+" sc="+sc+" tc="+tc); final boolean flag_realroots; /* This code solves the resolvent cubic in a convenient fashion * for this implementation of the quartic. If there are three real * roots, then they are placed directly into u[]. If two are * complex, then the real root is put into u[0] and the real * and imaginary part of the complex roots are placed into * u[1] and u[2], respectively. */ { final double qcub = (rc * rc - 3 * sc); final double rcub = (rc*(2 * rc * rc - 9 * sc) + 27 * tc); final double Q = qcub / 9; final double R = rcub / 54; final double Q3 = Q * Q * Q; final double R2 = R * R; final double CR2 = 729 * rcub * rcub; final double CQ3 = 2916 * qcub * qcub * qcub; if(PRINT_DEBUG) System.err.println("CR2="+CR2+" CQ3="+CQ3+" R="+R+" Q="+Q); if (0 == R && 0 == Q) { flag_realroots=true; u[0] = -rc3; u[1] = -rc3; u[2] = -rc3; } else if (CR2 == CQ3) { flag_realroots=true; final double sqrtQ = Math.sqrt (Q); if (R > 0) { u[0] = -2 * sqrtQ - rc3; u[1] = sqrtQ - rc3; u[2] = sqrtQ - rc3; } else { u[0] = -sqrtQ - rc3; u[1] = -sqrtQ - rc3; u[2] = 2 * sqrtQ - rc3; } } else if (R2 < Q3) { flag_realroots=true; final double ratio = (R >= 0?1:-1) * Math.sqrt (R2 / Q3); final double theta = Math.acos (ratio); final double norm = -2 * Math.sqrt (Q); u[0] = norm * Math.cos (theta / 3) - rc3; u[1] = norm * Math.cos ((theta + 2.0 * Math.PI) / 3) - rc3; u[2] = norm * Math.cos ((theta - 2.0 * Math.PI) / 3) - rc3; } else { flag_realroots=false; final double A = -(R >= 0?1:-1)*Math.pow(Math.abs(R)+Math.sqrt(R2-Q3),1.0/3.0); final double B = Q / A; u[0] = A + B - rc3; u[1] = -0.5 * (A + B) - rc3; u[2] = -(SQRT3*0.5) * Math.abs (A - B); } 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)); } /* End of solution to resolvent cubic */ /* Combine the square roots of the roots of the cubic * resolvent appropriately. Also, calculate 'mt' which * designates the nature of the roots: * mt=1 : 4 real roots * mt=2 : 0 real roots * mt=3 : 2 real roots */ final double w1_re,w1_im,w2_re,w2_im,w3_re,w3_im,mod_w1w2,mod_w1w2_squared; if (flag_realroots) { mod_w1w2=-1; mt = 2; int jmin=0; double vmin=Math.abs(u[jmin]); for(int j=1;j<3;j++) { final double vx=Math.abs(u[j]); if(vx=0) { w1_re=Math.sqrt(u1); w1_im=0; } else { w1_re=0; w1_im=Math.sqrt(-u1); } if(u2>=0) { w2_re=Math.sqrt(u2); w2_im=0; } else { w2_re=0; w2_im=Math.sqrt(-u2); } if(PRINT_DEBUG) System.err.println("u1="+u1+" u2="+u2+" jmin="+jmin); } else { mt = 3; 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); if(w_mod2_sq<=0) { w1_re=w1_im=0; } else { // calculate square root of a complex number (u[1],u[2]) // the result is in the (w1_re,w1_im) final double absu1=Math.abs(u[1]),absu2=Math.abs(u[2]),w; if(absu1>=absu2) { final double t=absu2/absu1; w=Math.sqrt(absu1*0.5 * (1.0 + Math.sqrt(1.0 + t * t))); if(PRINT_DEBUG) System.err.println(" Path1 "); } else { final double t=absu1/absu2; w=Math.sqrt(absu2*0.5 * (t + Math.sqrt(1.0 + t * t))); if(PRINT_DEBUG) System.err.println(" Path1a "); } if(u[1]>=0) { w1_re=w; w1_im=u[2]/(2*w); if(PRINT_DEBUG) System.err.println(" Path2 "); } else { final double vi = (u[2] >= 0) ? w : -w; w1_re=u[2]/(2*vi); w1_im=vi; if(PRINT_DEBUG) System.err.println(" Path2a "); } } final double absu0=Math.abs(u[0]); if(w_mod2>=absu0) { mod_w1w2=w_mod2; mod_w1w2_squared=w_mod2_sq; w2_re=w1_re; w2_im=-w1_im; } else { mod_w1w2=-1; mod_w1w2_squared=w_mod2*absu0; if(u[0]>=0) { w2_re=Math.sqrt(absu0); w2_im=0; } else { w2_re=0; w2_im=Math.sqrt(absu0); } } 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); } /* Solve the quadratic in order to obtain the roots * to the quartic */ if(mod_w1w2>0) { // a shorcut to reduce rounding error w3_re=qq/(-8)/mod_w1w2; w3_im=0; } else if(mod_w1w2_squared>0) { // regular path final double mqq8n=qq/(-8)/mod_w1w2_squared; w3_re=mqq8n*(w1_re*w2_re-w1_im*w2_im); w3_im=-mqq8n*(w1_re*w2_im+w2_re*w1_im); } else { // typically occur when qq==0 w3_re=w3_im=0; } final double h = r4 * a; 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); re_root[0]=w1_re+w2_re+w3_re-h; im_root[0]=w1_im+w2_im+w3_im; re_root[1]=-(w1_re+w2_re)+w3_re-h; im_root[1]=-(w1_im+w2_im)+w3_im; re_root[2]=w2_re-w1_re-w3_re-h; im_root[2]=w2_im-w1_im-w3_im; re_root[3]=w1_re-w2_re-w3_re-h; im_root[3]=w1_im-w2_im-w3_im; return 4; } static void setRandomP(final double [] p, final int n, Random r) { if(r.nextDouble()<0.1) { // integer coefficiens for(int j=0;j=0;k--) { final double res1=(res*rex-ims*imx)+p[k]; final double ims1=(ims*rex+res*imx); res=res1; ims=ims1; sabs+=xabs*sabs+p[k]; } sabs=Math.abs(sabs); if(false && sabs>1/eps? (!(Math.abs(res/sabs)<=eps)||!(Math.abs(ims/sabs)<=eps)) : (!(Math.abs(res)<=eps)||!(Math.abs(ims)<=eps))) { throw new RuntimeException( getPolinomTXT(p)+"\n"+ "\t x.r="+rex+" x.i="+imx+"\n"+ "res/sabs="+(res/sabs)+" ims/sabs="+(ims/sabs)+ " sabs="+sabs+ "\nres="+res+" ims="+ims+" n="+n+" eps="+eps+" "+ " sabs>1/eps="+(sabs>1/eps)+ " f1="+(!(Math.abs(res/sabs)<=eps)||!(Math.abs(ims/sabs)<=eps))+ " f2="+(!(Math.abs(res)<=eps)||!(Math.abs(ims)<=eps))+ " "+txt); } } static String getPolinomTXT(final double [] p) { final StringBuilder buf=new StringBuilder(); buf.append("order="+(p.length-1)+"\t"); for(int k=0;k-1;) { for(int dr=3;dr-->0;) { setRandomP(p,n,rn); for(int j=0;j<=dg;j++) { p[j]=0; } if(dr==0) { p[0]=-1+2.0*rn.nextDouble(); } else if(dr==1) { p[0]=p[1]=0; } findPolynomialRoots(n,p,rex,imx); for(int j=0;j