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