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

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

mercurial