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