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 *
|
|
18 * @run main 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 final long t0=System.currentTimeMillis();
|
|
761 final double eps=1e-6;
|
|
762 //checkRoots();
|
|
763 final java.util.Random r=new java.util.Random(-1381923);
|
|
764 printSpecialValues();
|
|
765
|
|
766 final int n_tests=10000000;
|
|
767 //testRoots(2,n_tests,r,eps);
|
|
768 //testRoots(3,n_tests,r,eps);
|
|
769 testRoots(4,n_tests,r,eps);
|
|
770 final long t1=System.currentTimeMillis();
|
|
771 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 $");
|
|
772 }
|
|
773
|
|
774
|
|
775
|
|
776 }
|