Mercurial > hg > truffle
annotate test/compiler/8005956/PolynomialRoot.java @ 18014:a96ae21442f6
Merge
author | asaha |
---|---|
date | Tue, 10 Jun 2014 13:43:58 -0700 |
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 } |