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 *
18 * @run main/timeout=300 PolynomialRoot
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 {
760 if (System.getProperty("os.arch").equals("x86") ||
761 System.getProperty("os.arch").equals("amd64") ||
762 System.getProperty("os.arch").equals("x86_64")){
763 final long t0=System.currentTimeMillis();
764 final double eps=1e-6;
765 //checkRoots();
766 final java.util.Random r=new java.util.Random(-1381923);
767 printSpecialValues();
768
769 final int n_tests=100000;
770 //testRoots(2,n_tests,r,eps);
771 //testRoots(3,n_tests,r,eps);
772 testRoots(4,n_tests,r,eps);
773 final long t1=System.currentTimeMillis();
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 $");
775 System.out.println("PASSED");
776 } else {
777 System.out.println("PASS test for non-x86");
778 }
779 }
780
781
782
783 }