Mercurial > hg > truffle
annotate test/compiler/8005956/PolynomialRoot.java @ 19466:cb7e8c0da978
Rename LIRLowTier to PostAllocationOptimizationStage.
| author | Josef Eisl <josef.eisl@jku.at> |
|---|---|
| date | Tue, 17 Feb 2015 21:38:05 +0100 |
| parents | e554162ab094 |
| children |
| rev | line source |
|---|---|
| 11045 | 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 * | |
|
11107
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
18 * @run main/timeout=300 PolynomialRoot |
| 11045 | 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 { | |
|
11107
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
760 if (System.getProperty("os.arch").equals("x86") || |
|
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
761 System.getProperty("os.arch").equals("amd64") || |
|
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
762 System.getProperty("os.arch").equals("x86_64")){ |
|
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
763 final long t0=System.currentTimeMillis(); |
|
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
764 final double eps=1e-6; |
|
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
765 //checkRoots(); |
|
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
766 final java.util.Random r=new java.util.Random(-1381923); |
|
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
767 printSpecialValues(); |
| 11045 | 768 |
|
11107
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
769 final int n_tests=100000; |
|
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
770 //testRoots(2,n_tests,r,eps); |
|
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
771 //testRoots(3,n_tests,r,eps); |
|
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
772 testRoots(4,n_tests,r,eps); |
|
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
773 final long t1=System.currentTimeMillis(); |
|
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
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 $"); |
|
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
775 System.out.println("PASSED"); |
|
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
776 } else { |
|
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
777 System.out.println("PASS test for non-x86"); |
|
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
778 } |
|
e554162ab094
8019625: Test compiler/8005956/PolynomialRoot.java timeouts on Solaris SPARCs
adlertz
parents:
11045
diff
changeset
|
779 } |
| 11045 | 780 |
| 781 | |
| 782 | |
| 783 } |