1 package org.djutils.polynomialroots;
2
3 import org.djutils.complex.Complex;
4
5 /**
6 * PolynomialRoots.java. Polynomial234RootSolvers - final all roots of linear, quadratic, cubic and quartic polynomials. Derived
7 * from <a href="https://netlib.org/toms/954.zip">https://dl.acm.org/doi/10.1145/2699468</a>. Manual translation from Fortran90
8 * to java by Peter Knoppers.
9 * <p>
10 * Copyright (c) 2020-2024 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved. <br>
11 * BSD-style license. See <a href="https://djutils.org/docs/current/djutils/licenses.html">DJUTILS License</a>.
12 * </p>
13 * @author <a href="https://www.tudelft.nl/averbraeck">Alexander Verbraeck</a>
14 * @author <a href="https://www.tudelft.nl/pknoppers">Peter Knoppers</a>
15 */
16 public final class PolynomialRoots
17 {
18 /**
19 * Do not instantiate.
20 */
21 private PolynomialRoots()
22 {
23 // Do not instantiate
24 }
25
26 /**
27 * Emulate the F77 sign function.
28 * @param a double; the value to optionally sign invert
29 * @param b double; the sign of which determines what to do
30 * @return double; if b >= 0 then a; else -a
31 */
32 private static double sign(final double a, final double b)
33 {
34 return b >= 0 ? a : -a;
35 }
36
37 /**
38 * LINEAR POLYNOMIAL ROOT SOLVER.
39 * <p>
40 * Calculates the root of the linear polynomial:<br>
41 * q1 * x + q0<br>
42 * Unlike the quadratic, cubic and quartic code, this is NOT derived from that Fortran90 code; it was added for completenes.
43 * @param q1 double; coefficient of the x term
44 * @param q0 double; independent coefficient
45 * @return Complex[]; the roots of the equation
46 */
47 public static Complex[] linearRoots(final double q1, final double q0)
48 {
49 if (q1 == 0)
50 {
51 return new Complex[] {}; // No roots; return empty array
52 }
53 return linearRoots(q0 / q1);
54 }
55
56 /**
57 * LINEAR POLYNOMIAL ROOT SOLVER.
58 * <p>
59 * Calculates the root of the linear polynomial:<br>
60 * x + q0<br>
61 * Unlike the quadratic, cubic and quartic code, this is NOT derived from that Fortran90 code; it was added for completenes.
62 * @param q0 double; independent coefficient
63 * @return Complex[]; the roots of the equation
64 */
65 public static Complex[] linearRoots(final double q0)
66 {
67 return new Complex[] { new Complex(-q0, 0) };
68 }
69
70 /**
71 * QUADRATIC POLYNOMIAL ROOT SOLVER
72 * <p>
73 * Calculates all real + complex roots of the quadratic polynomial:<br>
74 * q2 * x^2 + q1 * x + q0<br>
75 * The code checks internally if rescaling of the coefficients is needed to avoid overflow.
76 * <p>
77 * The order of the roots is as follows:<br>
78 * 1) For real roots, the order is according to their algebraic value on the number scale (largest positive first, largest
79 * negative last).<br>
80 * 2) Since there can be only one complex conjugate pair root, no order is necessary.<br>
81 * q1 : coefficient of x term q0 : independent coefficient
82 * @param q2 double; coefficient of the quadratic term
83 * @param q1 double; coefficient of the x term
84 * @param q0 double; independent coefficient
85 * @return Complex[]; the roots of the equation
86 */
87 public static Complex[] quadraticRoots(final double q2, final double q1, final double q0)
88 {
89 if (q2 == 0)
90 {
91 return linearRoots(q1, q0);
92 }
93 return quadraticRoots(q1 / q2, q0 / q2);
94 }
95
96 /**
97 * QUADRATIC POLYNOMIAL ROOT SOLVER
98 * <p>
99 * Calculates all real + complex roots of the quadratic polynomial:<br>
100 * x^2 + q1 * x + q0<br>
101 * The code checks internally if rescaling of the coefficients is needed to avoid overflow.
102 * <p>
103 * The order of the roots is as follows:<br>
104 * 1) For real roots, the order is according to their algebraic value on the number scale (largest positive first, largest
105 * negative last).<br>
106 * 2) Since there can be only one complex conjugate pair root, no order is necessary.<br>
107 * q1 : coefficient of x term q0 : independent coefficient
108 * @param q1 double; coefficient of the x term
109 * @param q0 double; independent coefficient
110 * @return Complex[]; the roots of the equation
111 */
112 public static Complex[] quadraticRoots(final double q1, final double q0)
113 {
114 boolean rescale;
115
116 double a0, a1;
117 double k = 0, x, y, z;
118
119 // Handle special cases.
120 if (q0 == 0.0 && q1 == 0.0)
121 {
122 // Two real roots at 0,0
123 return new Complex[] { Complex.ZERO, Complex.ZERO };
124 }
125 else if (q0 == 0.0)
126 {
127 // Two real roots; one of these is 0,0
128 // x^2 + q1 * x == x * (x + q1)
129 Complex nonZeroRoot = new Complex(-q1);
130 return new Complex[] { q1 > 0 ? Complex.ZERO : nonZeroRoot, q1 <= 0 ? nonZeroRoot : Complex.ZERO };
131 }
132 else if (q1 == 0.0)
133 {
134 x = Math.sqrt(Math.abs(q0));
135
136 if (q0 < 0.0)
137 {
138 // Two real roots, symmetrically around 0
139 return new Complex[] { new Complex(x, 0), new Complex(-x, 0) };
140 }
141 else
142 {
143 // Two complex roots, symmetrically around 0
144 return new Complex[] { new Complex(0, x), new Complex(0, -x) };
145 }
146 }
147 else
148 {
149 // The general case. Do rescaling, if either squaring of q1/2 or evaluation of
150 // (q1/2)^2 - q0 will lead to overflow. This is better than to have the solver
151 // crashed. Note, that rescaling might lead to loss of accuracy, so we only
152 // invoke it when absolutely necessary.
153 final double sqrtLPN = Math.sqrt(Double.MAX_VALUE); // Square root of the Largest Positive Number
154 rescale = (q1 > sqrtLPN + sqrtLPN); // this detects overflow of (q1/2)^2
155
156 if (!rescale)
157 {
158 x = q1 * 0.5; // we are sure here that x*x will not overflow
159 rescale = (q0 < x * x - Double.MAX_VALUE); // this detects overflow of (q1/2)^2 - q0
160 }
161
162 if (rescale)
163 {
164 x = Math.abs(q1);
165 y = Math.sqrt(Math.abs(q0));
166
167 if (x > y)
168 {
169 k = x;
170 z = 1.0 / x;
171 a1 = sign(1.0, q1);
172 a0 = (q0 * z) * z;
173 }
174 else
175 {
176 k = y;
177 a1 = q1 / y;
178 a0 = sign(1.0, q0);
179 }
180 }
181 else
182 {
183 a1 = q1;
184 a0 = q0;
185 }
186 // Determine the roots of the quadratic. Note, that either a1 or a0 might
187 // have become equal to zero due to underflow. But both cannot be zero.
188 x = a1 * 0.5;
189 y = x * x - a0;
190
191 if (y >= 0.0)
192 {
193 // Two real roots
194 y = Math.sqrt(y);
195
196 if (x > 0.0)
197 {
198 y = -x - y;
199 }
200 else
201 {
202 y = -x + y;
203 }
204
205 if (rescale)
206 {
207 y = y * k; // very important to convert to original
208 z = q0 / y; // root first, otherwise complete loss of
209 }
210 else // root due to possible a0 = 0 underflow
211 {
212 z = a0 / y;
213 }
214 return new Complex[] { new Complex(Math.max(y, z), 0), new Complex(Math.min(y, z), 0) };
215 }
216 else
217 {
218 // Two complex roots (zero real roots)
219 y = Math.sqrt(-y);
220
221 if (rescale)
222 {
223 x *= k;
224 y *= k;
225 }
226 return new Complex[] { new Complex(-x, y), new Complex(-x, -y) };
227 }
228 }
229 }
230
231 /**
232 * CUBIC POLYNOMIAL ROOT SOLVER.
233 * <p>
234 * Calculates all (real and complex) roots of the cubic polynomial:<br>
235 * c3 * x^3 + c2 * x^2 + c1 * x + c0<br>
236 * The first real root (which always exists) is obtained using an optimized Newton-Raphson scheme. The other remaining roots
237 * are obtained through composite deflation into a quadratic.
238 * <P>
239 * The cubic root solver can handle any size of cubic coefficients and there is no danger of overflow due to proper
240 * rescaling of the cubic polynomial. The order of the roots is as follows: 1) For real roots, the order is according to
241 * their algebraic value on the number scale (largest positive first, largest negative last). 2) Since there can be only one
242 * complex conjugate pair root, no order is necessary. 3) All real roots precede the complex ones.
243 * @param c3 double; coefficient of the cubic term
244 * @param c2 double; coefficient of the quadratic term
245 * @param c1 double; coefficient of the linear term
246 * @param c0 double; coefficient of the independent term
247 * @return Complex[]; array of Complex with all the roots
248 */
249 public static Complex[] cubicRoots(final double c3, final double c2, final double c1, final double c0)
250 {
251 if (c3 == 0)
252 {
253 return quadraticRoots(c2, c1, c0);
254 }
255 return cubicRoots(c2 / c3, c1 / c3, c0 / c3, false);
256 }
257
258 /**
259 * CUBIC POLYNOMIAL ROOT SOLVER.
260 * <p>
261 * Calculates all (real and complex) roots of the cubic polynomial:<br>
262 * x^3 + c2 * x^2 + c1 * x + c0<br>
263 * The first real root (which always exists) is obtained using an optimized Newton-Raphson scheme. The other remaining roots
264 * are obtained through composite deflation into a quadratic.
265 * <P>
266 * The cubic root solver can handle any size of cubic coefficients and there is no danger of overflow due to proper
267 * rescaling of the cubic polynomial. The order of the roots is as follows: 1) For real roots, the order is according to
268 * their algebraic value on the number scale (largest positive first, largest negative last). 2) Since there can be only one
269 * complex conjugate pair root, no order is necessary. 3) All real roots precede the complex ones.
270 * @param c2 double; coefficient of the quadratic term
271 * @param c1 double; coefficient of the linear term
272 * @param c0 double; coefficient of the independent term
273 * @return Complex[]; array of Complex with all the roots
274 */
275 public static Complex[] cubicRoots(final double c2, final double c1, final double c0)
276 {
277 return cubicRoots(c2, c1, c0, false);
278 }
279
280 /**
281 * CUBIC POLYNOMIAL ROOT SOLVER.
282 * <p>
283 * Calculates all (real and complex) roots of the cubic polynomial:<br>
284 * x^3 + c2 * x^2 + c1 * x + c0<br>
285 * The first real root (which always exists) is obtained using an optimized Newton-Raphson scheme. The other remaining roots
286 * are obtained through composite deflation into a quadratic. An option for printing detailed info about the intermediate
287 * stages in solving the cubic is available.
288 * <P>
289 * The cubic root solver can handle any size of cubic coefficients and there is no danger of overflow due to proper
290 * rescaling of the cubic polynomial. The order of the roots is as follows: 1) For real roots, the order is according to
291 * their algebraic value on the number scale (largest positive first, largest negative last). 2) Since there can be only one
292 * complex conjugate pair root, no order is necessary. 3) All real roots precede the complex ones.
293 * @param c2 double; coefficient of the quadratic term
294 * @param c1 double; coefficient of the linear term
295 * @param c0 double; coefficient of the independent term
296 * @param verbose boolean; if true; produce debugging output; if false; do not produce debugging output
297 * @return Complex[]; array of Complex with all the roots
298 */
299 public static Complex[] cubicRoots(final double c2, final double c1, final double c0, final boolean verbose)
300 {
301 final int cubicType;
302 int deflateCase;
303
304 final double macheps = Math.ulp(1.0);
305 final double one27th = 1.0 / 27.0;
306 final double two27th = 2.0 / 27.0;
307 final double third = 1.0 / 3.0;
308
309 // Newton-Raphson coeffs for class 5 and 6
310 final double p51 = 8.78558e-1;
311 final double p52 = 1.92823e-1;
312 final double p53 = 1.19748;
313 final double p54 = 3.45219e-1;
314 final double q51 = 5.71888e-1;
315 final double q52 = 5.66324e-1;
316 final double q53 = 2.83772e-1;
317 final double q54 = 4.01231e-1;
318 final double r51 = 7.11154e-1;
319 final double r52 = 5.05734e-1;
320 final double r53 = 8.37476e-1;
321 final double r54 = 2.07216e-1;
322 final double s51 = 3.22313e-1;
323 final double s52 = 2.64881e-1;
324 final double s53 = 3.56228e-1;
325 final double s54 = 4.45532e-3;
326
327 final int allzero = 0;
328 final int linear = 1;
329 final int quadratic = 2;
330 final int general = 3;
331
332 double a0 = 0, a1 = 0, a2 = 0;
333 double a = 0, b = 0, c = 0, k = 0, s = 0, t, u = 0, x = 0, y, z;
334 double xShift = 0;
335
336 if (verbose)
337 {
338 System.out.println("initial cubic c2 = " + c2);
339 System.out.println("initial cubic c1 = " + c1);
340 System.out.println("initial cubic c0 = " + c0);
341 System.out.println("------------------------------------------------");
342 }
343 // Handle special cases.
344 //
345 // 1) all terms zero
346 // 2) only quadratic term is nonzero -> linear equation.
347 // 3) only independent term is zero -> quadratic equation.
348 if (c0 == 0.0 && c1 == 0.0 && c2 == 0.0)
349 {
350 cubicType = allzero;
351 }
352 else if (c0 == 0.0 && c1 == 0.0)
353 {
354 k = 1.0;
355 a2 = c2;
356 cubicType = linear;
357 }
358 else if (c0 == 0.0)
359 {
360 k = 1.0;
361 a2 = c2;
362 a1 = c1;
363 cubicType = quadratic;
364 }
365 else
366 {
367 // The general case. Rescale cubic polynomial, such that largest absolute coefficient
368 // is (exactly!) equal to 1. Honor the presence of a special cubic case that might have
369 // been obtained during the rescaling process (due to underflow in the coefficients).
370 x = Math.abs(c2);
371 y = Math.sqrt(Math.abs(c1));
372 z = Math.cbrt(Math.abs(c0));
373 u = Math.max(Math.max(x, y), z);
374
375 if (u == x)
376 {
377 k = 1.0 / x;
378 a2 = sign(1.0, c2);
379 a1 = (c1 * k) * k;
380 a0 = ((c0 * k) * k) * k;
381 }
382 else if (u == y)
383 {
384 k = 1.0 / y;
385 a2 = c2 * k;
386 a1 = sign(1.0, c1);
387 a0 = ((c0 * k) * k) * k;
388 }
389 else
390 {
391 k = 1.0 / z;
392 a2 = c2 * k;
393 a1 = (c1 * k) * k;
394 a0 = sign(1.0, c0);
395 }
396
397 if (verbose)
398 {
399 System.out.println("rescaling factor = " + k);
400 System.out.println("------------------------------------------------");
401 System.out.println("rescaled cubic c2 = " + a2);
402 System.out.println("rescaled cubic c1 = " + a1);
403 System.out.println("rescaled cubic c0 = " + a0);
404 System.out.println("------------------------------------------------");
405 }
406
407 k = 1.0 / k;
408
409 if (a0 == 0.0 && a1 == 0.0 && a2 == 0.0)
410 {
411 cubicType = allzero;
412 }
413 else if (a0 == 0.0 && a1 == 0.0)
414 {
415 cubicType = linear;
416 }
417 else if (a0 == 0.0)
418 {
419 cubicType = quadratic;
420 }
421 else
422 {
423 cubicType = general;
424 }
425 }
426
427 switch (cubicType)
428 {
429 case allzero: // 1) Only zero roots
430 return new Complex[] { Complex.ZERO, Complex.ZERO, Complex.ZERO };
431
432 case linear: // 2) The linear equation case -> additional 2 zeros.
433 {
434 double root = -a2 * k;
435 return new Complex[] { new Complex(Math.max(0.0, root), 0), Complex.ZERO, new Complex(Math.min(0.0, root), 0) };
436 }
437
438 case quadratic: // 3) The quadratic equation case -> additional 1 zero.
439 {
440 Complex[] otherRoots = quadraticRoots(a2, a1);
441
442 if (otherRoots[0].isReal()) // There are guaranteed to be two roots of the quadratic sub case
443 {
444 // Three real roots
445 double xx = otherRoots[0].re * k;
446 double yy = otherRoots[1].re * k;
447 return new Complex[] { new Complex(Math.max(xx, 0.0), 0), new Complex(Math.max(yy, Math.min(xx, 0.0)), 0),
448 new Complex(Math.min(yy, 0.0), 0) };
449 }
450 else
451 {
452 // One real root and two complex roots
453 return new Complex[] { Complex.ZERO, otherRoots[0].times(k), otherRoots[1].times(k) };
454 }
455 }
456 case general:
457 {
458 // 3) The general cubic case. Set the best Newton-Raphson root estimates for the cubic.
459 // The easiest and most robust conditions are checked first. The most complicated
460 // ones are last and only done when absolutely necessary.
461 // Newton-Raphson coefficients for class 1 and 2
462 final double p1 = 1.09574;
463 final double q1 = 3.23900e-1;
464 final double r1 = 3.23900e-1;
465 final double s1 = 9.57439e-2;
466
467 if (a0 == 1.0)
468 {
469 x = -p1 + q1 * a1 - a2 * (r1 - s1 * a1);
470
471 a = a2;
472 b = a1;
473 c = a0;
474 xShift = 0.0;
475 }
476 else if (a0 == -1.0)
477 {
478 x = p1 - q1 * a1 - a2 * (r1 - s1 * a1);
479
480 a = a2;
481 b = a1;
482 c = a0;
483 xShift = 0.0;
484 }
485 else if (a1 == 1.0)
486 {
487 // Newton-Raphson coeffs for class 4
488 final double q4 = 7.71845e-1;
489 final double s4 = 2.28155e-1;
490
491 if (a0 > 0.0)
492 {
493 x = a0 * (-q4 - s4 * a2);
494 }
495 else
496 {
497 x = a0 * (-q4 + s4 * a2);
498 }
499
500 a = a2;
501 b = a1;
502 c = a0;
503 xShift = 0.0;
504 }
505 else if (a1 == -1.0)
506 {
507 // Newton-Raphson coeffs for class 3
508 final double p3 = 1.14413;
509 final double q3 = 2.75509e-1;
510 final double r3 = 4.45578e-1;
511 final double s3 = 2.59342e-2;
512
513 y = -two27th;
514 y = y * a2;
515 y = y * a2 - third;
516 y = y * a2;
517
518 if (a0 < y)
519 {
520 x = +p3 - q3 * a0 - a2 * (r3 + s3 * a0); // + guess
521 }
522 else
523 {
524 x = -p3 - q3 * a0 - a2 * (r3 - s3 * a0); // - guess
525 }
526
527 a = a2;
528 b = a1;
529 c = a0;
530 xShift = 0.0;
531 }
532 else if (a2 == 1.0)
533 {
534 b = a1 - third;
535 c = a0 - one27th;
536
537 if (Math.abs(b) < macheps && Math.abs(c) < macheps) // triple -1/3 root
538 {
539 x = -third * k;
540 Complex root = new Complex(x, 0);
541 return new Complex[] { root, root, root };
542 }
543 else
544 {
545 y = third * a1 - two27th;
546
547 if (a1 <= third)
548 {
549 if (a0 > y)
550 {
551 x = -p51 - q51 * a0 + a1 * (r51 - s51 * a0); // - guess
552 }
553 else
554 {
555 x = +p52 - q52 * a0 - a1 * (r52 + s52 * a0); // + guess
556 }
557 }
558 else if (a0 > y)
559 {
560 x = -p53 - q53 * a0 + a1 * (r53 - s53 * a0); // <-1/3 guess
561 }
562 else
563 {
564 x = +p54 - q54 * a0 - a1 * (r54 + s54 * a0); // >-1/3 guess
565 }
566 }
567 if (Math.abs(b) < 1.e-2 && Math.abs(c) < 1.e-2) // use shifted root
568 {
569 c = -third * b + c;
570 if (Math.abs(c) < macheps)
571 {
572 c = 0.0; // prevent random noise
573 }
574 a = 0.0;
575 xShift = third;
576 x = x + xShift;
577 }
578 else
579 {
580 a = a2;
581 b = a1;
582 c = a0;
583 xShift = 0.0;
584 }
585 }
586 else if (a2 == -1.0)
587 {
588 b = a1 - third;
589 c = a0 + one27th;
590
591 if (Math.abs(b) < macheps && Math.abs(c) < macheps) // triple 1/3 root
592 {
593 x = third * k;
594 Complex root = new Complex(x, 0);
595 return new Complex[] { root, root, root };
596 }
597 else
598 {
599 y = two27th - third * a1;
600
601 if (a1 <= third)
602 {
603 if (a0 < y)
604 {
605 x = +p51 - q51 * a0 - a1 * (r51 + s51 * a0); // +1 guess
606 }
607 else
608 {
609 x = -p52 - q52 * a0 + a1 * (r52 - s52 * a0); // -1 guess
610 }
611 }
612 else if (a0 < y)
613 {
614 x = +p53 - q53 * a0 - a1 * (r53 + s53 * a0); // >1/3 guess
615 }
616 else
617 {
618 x = -p54 - q54 * a0 + a1 * (r54 - s54 * a0); // <1/3 guess
619 }
620 }
621
622 if (Math.abs(b) < 1.e-2 && Math.abs(c) < 1.e-2) // use shifted root
623 {
624 c = third * b + c;
625 if (Math.abs(c) < macheps)
626 {
627 c = 0.0; // prevent random noise
628 }
629 a = 0.0;
630 xShift = -third;
631 x = x + xShift;
632 }
633 else
634 {
635 a = a2;
636 b = a1;
637 c = a0;
638 xShift = 0.0;
639 }
640 }
641
642 // Perform Newton/Bisection iterations on x^3 + ax^2 + bx + c.
643 z = x + a;
644 y = x + z;
645 z = z * x + b;
646 y = y * x + z; // C'(x)
647 z = z * x + c; // C(x)
648 t = z; // save C(x) for sign comparison
649 x = x - z / y; // 1st improved root
650
651 int oscillate = 0;
652 boolean bisection = false;
653 boolean converged = false;
654
655 while ((!converged) && (!bisection)) // Newton-Raphson iterates
656 {
657 z = x + a;
658 y = x + z;
659 z = z * x + b;
660 y = y * x + z;
661 z = z * x + c;
662
663 if (z * t < 0.0) // does Newton start oscillating ?
664 {
665 if (z < 0.0)
666 {
667 oscillate = oscillate + 1; // increment oscillation counter
668 s = x; // save lower bisection bound
669 }
670 else
671 {
672 u = x; // save upper bisection bound
673 }
674 t = z; // save current C(x)
675 }
676
677 y = z / y; // Newton correction
678 x = x - y; // new Newton root
679
680 bisection = oscillate > 2; // activate bisection
681 converged = Math.abs(y) <= Math.abs(x) * macheps; // Newton convergence indicator
682
683 if (verbose)
684 {
685 System.out.println("Newton root = " + x);
686 }
687 }
688
689 if (bisection)
690 {
691 t = u - s; // initial bisection interval
692 while (Math.abs(t) > Math.abs(x) * macheps) // bisection iterates
693 {
694 z = x + a;
695 z = z * x + b; // C (x)
696 z = z * x + c;
697
698 if (z < 0.0)
699 {
700 s = x;
701 }
702 else
703 {
704 u = x; // keep bracket on root
705 }
706
707 t = 0.5 * (u - s); // new bisection interval
708 x = s + t; // new bisection root
709
710 if (verbose)
711 {
712 System.out.println("Bisection root = " + x);
713 }
714 }
715 }
716
717 if (verbose)
718 {
719 System.out.println("------------------------------------------------");
720 }
721
722 x = x - xShift; // unshift root
723 // Forward / backward deflate rescaled cubic (if needed) to check for other real roots.
724 // The deflation analysis is performed on the rescaled cubic. The actual deflation must
725 // be performed on the original cubic, not the rescaled one. Otherwise deflation errors
726 // will be enhanced when undoing the rescaling on the extra roots.
727 z = Math.abs(x);
728 s = Math.abs(a2);
729 t = Math.abs(a1);
730 u = Math.abs(a0);
731
732 y = z * Math.max(s, z); // take maximum between |x^2|,|a2 * x|
733
734 deflateCase = 1; // up to now, the maximum is |x^3| or |a2 * x^2|
735
736 if (y < t) // check maximum between |x^2|,|a2 * x|,|a1|
737 {
738 y = t * z; // the maximum is |a1 * x|
739 deflateCase = 2; // up to now, the maximum is |a1 * x|
740 }
741 else
742 {
743 y = y * z; // the maximum is |x^3| or |a2 * x^2|
744 }
745
746 if (y < u) // check maximum between |x^3|,|a2 * x^2|,|a1 * x|,|a0|
747 {
748 deflateCase = 3; // the maximum is |a0|
749 }
750
751 y = x * k; // real root of original cubic
752
753 switch (deflateCase)
754 {
755 case 1:
756 x = 1.0 / y;
757 t = -c0 * x; // t -> backward deflation on unscaled cubic
758 s = (t - c1) * x; // s -> backward deflation on unscaled cubic
759 break;
760 case 2:
761 s = c2 + y; // s -> forward deflation on unscaled cubic
762 t = -c0 / y; // t -> backward deflation on unscaled cubic
763 break;
764 case 3:
765 s = c2 + y; // s -> forward deflation on unscaled cubic
766 t = c1 + s * y; // t -> forward deflation on unscaled cubic
767 break;
768
769 default:
770 throw new RuntimeException("Bad switch; cannot happen");
771 }
772
773 if (verbose)
774 {
775 System.out.println("Residual quadratic q1 = " + s);
776 System.out.println("Residual quadratic q0 = " + t);
777 System.out.println("------------------------------------------------");
778 }
779
780 Complex[] quadraticRoots = quadraticRoots(s, t);
781 // call quadraticRoots (s, t, nReal, root (1:2,1:2))
782
783 if (quadraticRoots[0].isReal())
784 {
785 // Three real roots
786 return new Complex[] { new Complex(Math.max(quadraticRoots[0].re, y), 0),
787 new Complex(Math.max(quadraticRoots[1].re, Math.min(quadraticRoots[0].re, y)), 0),
788 new Complex(Math.min(quadraticRoots[1].re, y), 0) };
789 }
790 else
791 {
792 // One real root and two complex roots
793 return new Complex[] { new Complex(y, 0), quadraticRoots[0], quadraticRoots[1] };
794 }
795 }
796
797 default:
798 throw new RuntimeException("Bad switch; cannot happen");
799 }
800 }
801
802 /**
803 * QUARTIC POLYNOMIAL ROOT SOLVER
804 * <p>
805 * Calculates all real + complex roots of the quartic polynomial:<br>
806 * q4 * x^4 + q3 * x^3 + q2 * x^2 + q1 * x + q0
807 * <p>
808 * The quartic root solver can handle any size of quartic coefficients and there is no danger of overflow, due to proper
809 * rescaling of the quartic polynomial.
810 * <p>
811 * The order of the roots is as follows:<br>
812 * 1) For real roots, the order is according to their algebraic value on the number scale (largest positive first, largest
813 * negative last).<br>
814 * 2) For complex conjugate pair roots, the order is according to the algebraic value of their real parts (largest positive
815 * first). If the real parts are equal, the order is according to the algebraic value of their imaginary parts (largest
816 * first).<br>
817 * 3) All real roots precede the complex ones.
818 * @param q4 double; coefficient of the quartic term
819 * @param q3 double; coefficient of the cubic term
820 * @param q2 double; coefficient of the quadratic term
821 * @param q1 double; coefficient of the linear term
822 * @param q0 double; independent coefficient
823 * @return Complex[]; array of Complex with all the roots
824 */
825 public static Complex[] quarticRoots(final double q4, final double q3, final double q2, final double q1, final double q0)
826 {
827 if (q4 == 0)
828 {
829 return cubicRoots(q3, q2, q1, q0);
830 }
831 return quarticRoots(q3 / q4, q2 / q4, q1 / q4, q0 / q4);
832 }
833
834 /**
835 * QUARTIC POLYNOMIAL ROOT SOLVER
836 * <p>
837 * Calculates all real + complex roots of the quartic polynomial:<br>
838 * x^4 + q3 * x^3 + q2 * x^2 + q1 * x + q0
839 * <p>
840 * The quartic root solver can handle any size of quartic coefficients and there is no danger of overflow, due to proper
841 * rescaling of the quartic polynomial.
842 * <p>
843 * The order of the roots is as follows:<br>
844 * 1) For real roots, the order is according to their algebraic value on the number scale (largest positive first, largest
845 * negative last).<br>
846 * 2) For complex conjugate pair roots, the order is according to the algebraic value of their real parts (largest positive
847 * first). If the real parts are equal, the order is according to the algebraic value of their imaginary parts (largest
848 * first).<br>
849 * 3) All real roots precede the complex ones.
850 * @param q3 double; coefficient of the cubic term
851 * @param q2 double; coefficient of the quadratic term
852 * @param q1 double; coefficient of the linear term
853 * @param q0 double; independent coefficient
854 * @return Complex[]; array of Complex with all the roots
855 */
856 public static Complex[] quarticRoots(final double q3, final double q2, final double q1, final double q0)
857 {
858 return quarticRoots(q3, q2, q1, q0, false);
859 }
860
861 /**
862 * QUARTIC POLYNOMIAL ROOT SOLVER
863 * <p>
864 * Calculates all real + complex roots of the quartic polynomial:<br>
865 * x^4 + q3 * x^3 + q2 * x^2 + q1 * x + q0
866 * <p>
867 * An option for printing detailed info about the intermediate stages in solving the quartic is available. This enables a
868 * detailed check in case something went wrong and the roots obtained are not proper.<br>
869 * The quartic root solver can handle any size of quartic coefficients and there is no danger of overflow, due to proper
870 * rescaling of the quartic polynomial.
871 * <p>
872 * The order of the roots is as follows:<br>
873 * 1) For real roots, the order is according to their algebraic value on the number scale (largest positive first, largest
874 * negative last).<br>
875 * 2) For complex conjugate pair roots, the order is according to the algebraic value of their real parts (largest positive
876 * first). If the real parts are equal, the order is according to the algebraic value of their imaginary parts (largest
877 * first).<br>
878 * 3) All real roots precede the complex ones.
879 * @param q3 double; coefficient of the cubic term
880 * @param q2 double; coefficient of the quadratic term
881 * @param q1 double; coefficient of the linear term
882 * @param q0 double; independent coefficient
883 * @param verbose boolean; if true; produce debugging output; if false; do not produce debugging output
884 * @return Complex[]; array of Complex with all the roots
885 */
886 public static Complex[] quarticRoots(final double q3, final double q2, final double q1, final double q0,
887 final boolean verbose)
888 {
889 int quarticType;
890
891 final int biquadratic = 2;
892 final int cubic = 3;
893 final int general = 4;
894
895 double a0 = 0, a1 = 0, a2, a3 = 0;
896 double a, b, c, d, k, s, t, u = 0, x, y, z;
897
898 final double macheps = Math.ulp(1.0);
899
900 if (verbose)
901 {
902 System.out.println("initial quartic q3 = " + q3);
903 System.out.println("initial quartic q2 = " + q2);
904 System.out.println("initial quartic q1 = " + q1);
905 System.out.println("initial quartic q0 = " + q0);
906 System.out.println("------------------------------------------------");
907 }
908 /*
909 * Handle special cases. Since the cubic solver handles all its special cases by itself, we need to check only for two
910 * cases:<br> 1) independent term is zero -> solve cubic and include the zero root <br> 2) the biquadratic case.
911 */
912 if (q0 == 0.0)
913 {
914 k = 1.0;
915 a3 = q3;
916 a2 = q2;
917 a1 = q1;
918 quarticType = cubic;
919 }
920 else if (q3 == 0.0 && q1 == 0.0)
921 {
922 k = 1.0;
923 a2 = q2;
924 a0 = q0;
925 quarticType = biquadratic;
926 }
927 else
928 {
929 /*
930 * The general case. Rescale quartic polynomial, such that largest absolute coefficient is (exactly!) equal to 1.
931 * Honor the presence of a special quartic case that might have been obtained during the rescaling process (due to
932 * underflow in the coefficients).
933 */
934 s = Math.abs(q3);
935 t = Math.sqrt(Math.abs(q2));
936 u = Math.cbrt(Math.abs(q1));
937 x = Math.sqrt(Math.sqrt(Math.abs(q0)));
938 y = Math.max(Math.max(Math.max(s, t), u), x);
939
940 if (y == s)
941 {
942 k = 1.0 / s;
943 a3 = sign(1.0, q3);
944 a2 = (q2 * k) * k;
945 a1 = ((q1 * k) * k) * k;
946 a0 = (((q0 * k) * k) * k) * k;
947 }
948 else if (y == t)
949 {
950 k = 1.0 / t;
951 a3 = q3 * k;
952 a2 = sign(1.0, q2);
953 a1 = ((q1 * k) * k) * k;
954 a0 = (((q0 * k) * k) * k) * k;
955 }
956 else if (y == u)
957 {
958 k = 1.0 / u;
959 a3 = q3 * k;
960 a2 = (q2 * k) * k;
961 a1 = sign(1.0, q1);
962 a0 = (((q0 * k) * k) * k) * k;
963 }
964 else
965 {
966 k = 1.0 / x;
967 a3 = q3 * k;
968 a2 = (q2 * k) * k;
969 a1 = ((q1 * k) * k) * k;
970 a0 = sign(1.0, q0);
971 }
972
973 k = 1.0 / k;
974
975 if (verbose)
976 {
977 System.out.println("rescaling factor = " + k);
978 System.out.println("------------------------------------------------");
979 System.out.println("rescaled quartic q3 = " + a3);
980 System.out.println("rescaled quartic q2 = " + a2);
981 System.out.println("rescaled quartic q1 = " + a1);
982 System.out.println("rescaled quartic q0 = " + a0);
983 System.out.println("------------------------------------------------");
984 }
985
986 if (a0 == 0.0)
987 {
988 quarticType = cubic;
989 }
990 else if (a3 == 0.0 && a1 == 0.0)
991 {
992 quarticType = biquadratic;
993 }
994 else
995 {
996 quarticType = general;
997 }
998 }
999
1000 switch (quarticType)
1001 {
1002 case cubic: // 1) The quartic with independent term = 0 -> solve cubic and add a zero root.
1003 {
1004 // x^4 + q3 * x^3 + q2 * x^2 + q1 * x + 0 == x * (x^3 + q3 * x^2 + q2 * x + q1)
1005 Complex[] cubicRoots = cubicRoots(a3, a2, a1, verbose);
1006
1007 if (cubicRoots.length == 3 && cubicRoots[0].isReal() && cubicRoots[1].isReal())
1008 {
1009 // Three real roots of the cubic; ordered x >= y >= z
1010 x = cubicRoots[0].re * k;
1011 y = cubicRoots[1].re * k;
1012 z = cubicRoots[2].re * k;
1013
1014 return new Complex[] { new Complex(Math.max(x, 0), 0), new Complex(Math.max(y, Math.min(x, 0)), 0),
1015 new Complex(Math.max(z, Math.min(y, 0)), 0), new Complex(Math.min(z, 0), 0) };
1016 }
1017 else
1018 {
1019 // One real cubic root; should be in the first entry of the array
1020 if (!cubicRoots[0].isReal())
1021 {
1022 throw new RuntimeException("Cannot happen");
1023 }
1024 x = cubicRoots[0].re * k;
1025
1026 return new Complex[] { new Complex(Math.max(0, x), 0), new Complex(Math.min(0, x), 0), cubicRoots[1],
1027 cubicRoots[2] };
1028 }
1029 }
1030
1031 case biquadratic: // The quartic with x^3 and x terms = 0 -> solve biquadratic.
1032 {
1033 Complex[] quadraticRoots = quadraticRoots(q2, q0);
1034 if (quadraticRoots.length == 2 && quadraticRoots[0].isReal() && quadraticRoots[1].isReal())
1035 {
1036 x = quadraticRoots[0].re; // real roots of quadratic are ordered x >= y
1037 y = quadraticRoots[1].re;
1038
1039 if (y >= 0.0)
1040 {
1041 x = Math.sqrt(x) * k;
1042 y = Math.sqrt(y) * k;
1043 return new Complex[] { new Complex(x, 0), new Complex(y, 0), new Complex(-y, 0), new Complex(-x, 0) };
1044 }
1045 else if (x >= 0.0 && y < 0.0)
1046 {
1047 x = Math.sqrt(x) * k;
1048 y = Math.sqrt(Math.abs(y)) * k;
1049 return new Complex[] { new Complex(x, 0), new Complex(-x, 0), new Complex(0, y), new Complex(0, -y) };
1050 }
1051 else if (x < 0.0)
1052 {
1053 x = Math.sqrt(Math.abs(x)) * k;
1054 y = Math.sqrt(Math.abs(y)) * k;
1055 return new Complex[] { new Complex(0, y), new Complex(0, x), new Complex(0, -x), new Complex(0, -y) };
1056 }
1057 }
1058 else
1059 {
1060 // complex conjugate pair biquadratic roots x +/- iy.
1061
1062 x = quadraticRoots[0].re * 0.5;
1063 y = quadraticRoots[0].im * 0.5;
1064 z = Math.sqrt(x * x + y * y);
1065 y = Math.sqrt(z - x) * k;
1066 x = Math.sqrt(z + x) * k;
1067 return new Complex[] { new Complex(x, y), new Complex(x, -y), new Complex(-x, y), new Complex(-x, -y) };
1068
1069 }
1070 }
1071
1072 case general:
1073 {
1074 int nReal;
1075 /*
1076 * 3) The general quartic case. Search for stationary points. Set the first derivative polynomial (cubic) equal
1077 * to zero and find its roots. Check, if any minimum point of Q(x) is below zero, in which case we must have
1078 * real roots for Q(x). Hunt down only the real root, which will potentially converge fastest during Newton
1079 * iterates. The remaining roots will be determined by deflation Q(x) -> cubic.<p> The best roots for the Newton
1080 * iterations are the two on the opposite ends, i.e. those closest to the +2 and -2. Which of these two roots to
1081 * take, depends on the location of the Q(x) minima x = s and x = u, with s > u. There are three cases:<br> 1)
1082 * both Q(s) and Q(u) < 0<br> The best root is the one that corresponds to the lowest of these minima. If Q(s)
1083 * is lowest -> start Newton from +2 downwards (or zero, if s < 0 and a0 > 0). If Q(u) is lowest -> start Newton
1084 * from -2 upwards (or zero, if u > 0 and a0 > 0).<br> 2) only Q(s) < 0>br? With both sides +2 and -2 possible
1085 * as a Newton starting point, we have to avoid the area in the Q(x) graph, where inflection points are present.
1086 * Solving Q''(x) = 0, leads to solutions x = -a3/4 +/- discriminant, i.e. they are centered around -a3/4. Since
1087 * both inflection points must be either on the r.h.s or l.h.s. from x = s, a simple test where s is in relation
1088 * to -a3/4 allows us to avoid the inflection point area.<br> 3) only Q(u) < 0<br> Same of what has been said
1089 * under 2) but with x = u.
1090 */
1091
1092 x = 0.75 * a3;
1093 y = 0.50 * a2;
1094 z = 0.25 * a1;
1095
1096 if (verbose)
1097 {
1098 System.out.println("dQ(x)/dx cubic c2 = " + x);
1099 System.out.println("dQ(x)/dx cubic c1 = " + y);
1100 System.out.println("dQ(x)/dx cubic c0 = " + z);
1101 System.out.println("------------------------------------------------");
1102 }
1103
1104 Complex[] cubicRoots = cubicRoots(x, y, z);
1105
1106 s = cubicRoots[0].re; // Q'(x) root s (real for sure)
1107 x = s + a3;
1108 x = x * s + a2;
1109 x = x * s + a1;
1110 x = x * s + a0; // Q(s)
1111
1112 y = 1.0; // dual info: Q'(x) has more real roots, and if so, is Q(u) < 0 ?
1113
1114 if (cubicRoots[1].isReal()) // then they're all real
1115 {
1116 u = cubicRoots[2].re; // Q'(x) root u
1117 y = u + a3;
1118 y = y * u + a2;
1119 y = y * u + a1;
1120 y = y * u + a0; // Q(u)
1121 }
1122
1123 if (verbose)
1124 {
1125 System.out.println("dQ(x)/dx root s = " + s);
1126 System.out.println("Q(s) = " + x);
1127 System.out.println("dQ(x)/dx root u = " + u);
1128 System.out.println("Q(u) = " + y);
1129 System.out.println("------------------------------------------------");
1130 }
1131
1132 if (x < 0.0 && y < 0.0)
1133 {
1134 if (x < y)
1135 {
1136 if (s < 0.0)
1137 {
1138 x = 1.0 - sign(1.0, a0);
1139 }
1140 else
1141 {
1142 x = 2.0;
1143 }
1144 }
1145 else if (u > 0.0)
1146 {
1147 x = -1.0 + sign(1.0, a0);
1148 }
1149 else
1150 {
1151 x = -2.0;
1152 }
1153
1154 nReal = 1;
1155 }
1156 else if (x < 0.0)
1157 {
1158 if (s < -a3 * 0.25)
1159 {
1160 if (s > 0.0)
1161 {
1162 x = -1.0 + sign(1.0, a0);
1163 }
1164 else
1165 {
1166 x = -2.0;
1167 }
1168 }
1169 else if (s < 0.0)
1170 {
1171 x = 1.0 - sign(1.0, a0);
1172 }
1173 else
1174 {
1175 x = 2.0;
1176 }
1177 nReal = 1;
1178 }
1179 else if (y < 0.0)
1180 {
1181 if (u < -a3 * 0.25)
1182 {
1183 if (u > 0.0)
1184 {
1185 x = -1.0 + sign(1.0, a0);
1186 }
1187 else
1188 {
1189 x = -2.0;
1190 }
1191 }
1192 else if (u < 0.0)
1193 {
1194 x = 1.0 - sign(1.0, a0);
1195 }
1196 else
1197 {
1198 x = 2.0;
1199 }
1200 nReal = 1;
1201 }
1202 else
1203 {
1204 nReal = 0;
1205 }
1206 /*
1207 * Do all necessary Newton iterations. In case we have more than 2 oscillations, exit the Newton iterations and
1208 * switch to bisection. Note, that from the definition of the Newton starting point, we always have Q(x) > 0 and
1209 * Q'(x) starts (-ve/+ve) for the (-2/+2) starting points and (increase/decrease) smoothly and staying (< 0 / >
1210 * 0). In practice, for extremely shallow Q(x) curves near the root, the Newton procedure can overshoot slightly
1211 * due to rounding errors when approaching the root. The result are tiny oscillations around the root. If such a
1212 * situation happens, the Newton iterations are abandoned after 3 oscillations and further location of the root
1213 * is done using bisection starting with the oscillation brackets.
1214 */
1215 if (nReal > 0)
1216 {
1217 int oscillate = 0;
1218 boolean bisection = false;
1219 boolean converged = false;
1220 int deflateCase;
1221
1222 while ((!converged) && (!bisection)) // Newton-Raphson iterates
1223 {
1224 y = x + a3;
1225 z = x + y;
1226 y = y * x + a2; // y = Q(x)
1227 z = z * x + y;
1228 y = y * x + a1; // z = Q'(x)
1229 z = z * x + y;
1230 y = y * x + a0;
1231
1232 if (y < 0.0) // does Newton start oscillating ?
1233 {
1234 oscillate = oscillate + 1; // increment oscillation counter
1235 s = x; // save lower bisection bound
1236 }
1237 else
1238 {
1239 u = x; // save upper bisection bound
1240 }
1241
1242 y = y / z; // Newton correction
1243 x = x - y; // new Newton root
1244
1245 bisection = oscillate > 2; // activate bisection
1246 converged = Math.abs(y) <= Math.abs(x) * macheps; // Newton convergence indicator
1247
1248 if (verbose)
1249 {
1250 System.out.println("Newton root = " + x);
1251 }
1252 }
1253
1254 if (bisection)
1255 {
1256 t = u - s; // initial bisection interval
1257 while (Math.abs(t) > Math.abs(x) * macheps) // bisection iterates
1258 {
1259 y = x + a3;
1260 y = y * x + a2; // y = Q(x)
1261 y = y * x + a1;
1262 y = y * x + a0;
1263
1264 if (y < 0.0)
1265 {
1266 s = x;
1267 }
1268 else // keep bracket on root
1269 {
1270 u = x;
1271 }
1272
1273 t = 0.5 * (u - s); // new bisection interval
1274 x = s + t; // new bisection root
1275
1276 if (verbose)
1277 {
1278 System.out.println("Bisection root = " + x);
1279 }
1280 }
1281 }
1282
1283 if (verbose)
1284 {
1285 System.out.println("------------------------------------------------");
1286 }
1287 /*
1288 * Find remaining roots -> reduce to cubic. The reduction to a cubic polynomial is done using composite
1289 * deflation to minimize rounding errors. Also, while the composite deflation analysis is done on the
1290 * reduced quartic, the actual deflation is being performed on the original quartic again to avoid enhanced
1291 * propagation of root errors.
1292 */
1293 z = Math.abs(x);
1294 a = Math.abs(a3);
1295 b = Math.abs(a2); // prepare for composite deflation
1296 c = Math.abs(a1);
1297 d = Math.abs(a0);
1298
1299 y = z * Math.max(a, z); // take maximum between |x^2|,|a3 * x|
1300
1301 deflateCase = 1; // up to now, the maximum is |x^4| or |a3 * x^3|
1302
1303 if (y < b) // check maximum between |x^2|,|a3 * x|,|a2|
1304 {
1305 y = b * z; // the maximum is |a2| -> form |a2 * x|
1306 deflateCase = 2; // up to now, the maximum is |a2 * x^2|
1307 }
1308 else
1309 {
1310 y = y * z; // the maximum is |x^3| or |a3 * x^2|
1311 }
1312
1313 if (y < c) // check maximum between |x^3|,|a3 * x^2|,|a2 * x|,|a1|
1314 {
1315 y = c * z; // the maximum is |a1| -> form |a1 * x|
1316 deflateCase = 3; // up to now, the maximum is |a1 * x|
1317 }
1318 else
1319 {
1320 y = y * z; // the maximum is |x^4|,|a3 * x^3| or |a2 * x^2|
1321 }
1322
1323 if (y < d) // check maximum between |x^4|,|a3 * x^3|,|a2 * x^2|,|a1 * x|,|a0|
1324 {
1325 deflateCase = 4; // the maximum is |a0|
1326 }
1327
1328 x = x * k; // 1st real root of original Q(x)
1329
1330 switch (deflateCase)
1331 {
1332 case 1:
1333 z = 1.0 / x;
1334 u = -q0 * z; // u -> backward deflation on original Q(x)
1335 t = (u - q1) * z; // t -> backward deflation on original Q(x)
1336 s = (t - q2) * z; // s -> backward deflation on original Q(x)
1337 break;
1338 case 2:
1339 z = 1.0 / x;
1340 u = -q0 * z; // u -> backward deflation on original Q(x)
1341 t = (u - q1) * z; // t -> backward deflation on original Q(x)
1342 s = q3 + x; // s -> forward deflation on original Q(x)
1343 break;
1344 case 3:
1345 s = q3 + x; // s -> forward deflation on original Q(x)
1346 t = q2 + s * x; // t -> forward deflation on original Q(x)
1347 u = -q0 / x; // u -> backward deflation on original Q(x)
1348 break;
1349 case 4:
1350 s = q3 + x; // s -> forward deflation on original Q(x)
1351 t = q2 + s * x; // t -> forward deflation on original Q(x)
1352 u = q1 + t * x; // u -> forward deflation on original Q(x)
1353 break;
1354 default:
1355 throw new RuntimeException("Bad switch; cannot happen");
1356 }
1357
1358 if (verbose)
1359 {
1360 System.out.println("Residual cubic c2 = " + s);
1361 System.out.println("Residual cubic c1 = " + t);
1362 System.out.println("Residual cubic c0 = " + u);
1363 System.out.println("------------------------------------------------");
1364 }
1365
1366 cubicRoots = cubicRoots(s, t, u, verbose);
1367 nReal = 0;
1368 for (Complex complex : cubicRoots)
1369 {
1370 if (complex.isReal())
1371 {
1372 nReal++;
1373 }
1374 }
1375 if (nReal == 3)
1376 {
1377 s = cubicRoots[0].re;
1378 t = cubicRoots[1].re; // real roots of cubic are ordered s >= t >= u
1379 u = cubicRoots[2].re;
1380 // Construct a new array and insert x at the appropriate place
1381 return new Complex[] { new Complex(Math.max(s, x), 0), new Complex(Math.max(t, Math.min(s, x)), 0),
1382 new Complex(Math.max(u, Math.min(t, x)), 0), new Complex(Math.min(u, x), 0) };
1383 }
1384 else // there is only one real cubic root here
1385 {
1386 s = cubicRoots[0].re;
1387 return new Complex[] { new Complex(Math.max(s, x), 0), new Complex(Math.min(s, x), 0), cubicRoots[1],
1388 cubicRoots[2] };
1389 }
1390 }
1391 else
1392 {
1393 /*
1394 * As no real roots have been found by now, only complex roots are possible. Find real parts of roots first,
1395 * followed by imaginary components.
1396 */
1397 s = a3 * 0.5;
1398 t = s * s - a2;
1399 u = s * t + a1; // value of Q'(-a3/4) at stationary point -a3/4
1400
1401 boolean notZero = (Math.abs(u) >= macheps); // H(-a3/4) is considered > 0 at stationary point
1402
1403 if (verbose)
1404 {
1405 System.out.println("dQ/dx (-a3/4) value = " + u);
1406 System.out.println("------------------------------------------------");
1407 }
1408
1409 boolean minimum;
1410 if (a3 != 0.0)
1411 {
1412 s = a1 / a3;
1413 minimum = (a0 > s * s); // H''(-a3/4) > 0 -> minimum
1414 }
1415 else
1416 {
1417 minimum = (4 * a0 > a2 * a2); // H''(-a3/4) > 0 -> minimum
1418 }
1419
1420 boolean iterate = notZero || (!notZero && minimum);
1421
1422 if (iterate)
1423 {
1424 x = sign(2.0, a3); // initial root -> target = smaller mag root
1425
1426 int oscillate = 0;
1427 boolean bisection = false;
1428 boolean converged = false;
1429
1430 while (!converged && !bisection) // ! Newton-Raphson iterates
1431 {
1432 a = x + a3;
1433 b = x + a; // a = Q(x)
1434 c = x + b;
1435 d = x + c; // b = Q'(x)
1436 a = a * x + a2;
1437 b = b * x + a; // c = Q''(x) / 2
1438 c = c * x + b;
1439 a = a * x + a1; // d = Q'''(x) / 6
1440 b = b * x + a;
1441 a = a * x + a0;
1442 y = a * d * d - b * c * d + b * b; // y = H(x), usually < 0
1443 z = 2 * d * (4 * a - b * d - c * c); // z = H'(x)
1444
1445 if (y > 0.0) // does Newton start oscillating ?
1446 {
1447 oscillate = oscillate + 1; // increment oscillation counter
1448 s = x; // save upper bisection bound
1449 }
1450 else
1451 {
1452 u = x; // save lower bisection bound
1453 }
1454
1455 y = y / z; // Newton correction
1456 x = x - y; // new Newton root
1457
1458 bisection = oscillate > 2; // activate bisection
1459 converged = Math.abs(y) <= Math.abs(x) * macheps; // Newton convergence criterion
1460
1461 if (verbose)
1462 {
1463 System.out.println("Newton H(x) root = " + x);
1464 }
1465
1466 }
1467
1468 if (bisection)
1469 {
1470 t = u - s; // initial bisection interval
1471 while (Math.abs(t) > Math.abs(x * macheps)) // bisection iterates
1472 {
1473 a = x + a3;
1474 b = x + a; // a = Q(x)
1475 c = x + b;
1476 d = x + c; // b = Q'(x)
1477 a = a * x + a2;
1478 b = b * x + a; // c = Q''(x) / 2
1479 c = c * x + b;
1480 a = a * x + a1; // d = Q'''(x) / 6
1481 b = b * x + a;
1482 a = a * x + a0;
1483 y = a * d * d - b * c * d + b * b; // y = H(x)
1484
1485 if (y > 0.0)
1486 {
1487 s = x;
1488 }
1489 else // keep bracket on root
1490 {
1491 u = x;
1492 }
1493
1494 t = 0.5 * (u - s); // new bisection interval
1495 x = s + t; // new bisection root
1496
1497 if (verbose)
1498 {
1499 System.out.println("Bisection H(x) root = " + x);
1500 }
1501
1502 }
1503 }
1504
1505 if (verbose)
1506 {
1507 System.out.println("------------------------------------------------");
1508 }
1509
1510 a = x * k; // 1st real component -> a
1511 b = -0.5 * q3 - a; // 2nd real component -> b
1512
1513 x = 4 * a + q3; // Q'''(a)
1514 y = x + q3 + q3;
1515 y = y * a + q2 + q2; // Q'(a)
1516 y = y * a + q1;
1517 y = Math.max(y / x, 0.0); // ensure Q'(a) / Q'''(a) >= 0
1518 x = 4 * b + q3; // Q'''(b)
1519 z = x + q3 + q3;
1520 z = z * b + q2 + q2; // Q'(b)
1521 z = z * b + q1;
1522 z = Math.max(z / x, 0.0); // ensure Q'(b) / Q'''(b) >= 0
1523 c = a * a; // store a^2 for later
1524 d = b * b; // store b^2 for later
1525 s = c + y; // magnitude^2 of (a + iy) root
1526 t = d + z; // magnitude^2 of (b + iz) root
1527
1528 if (s > t) // minimize imaginary error
1529 {
1530 c = Math.sqrt(y); // 1st imaginary component -> c
1531 d = Math.sqrt(q0 / s - d); // 2nd imaginary component -> d
1532 }
1533 else
1534 {
1535 c = Math.sqrt(q0 / t - c); // 1st imaginary component -> c
1536 d = Math.sqrt(z); // 2nd imaginary component -> d
1537 }
1538 }
1539 else // no bisection -> real components equal
1540 {
1541
1542 a = -0.25 * q3; // 1st real component -> a
1543 b = a; // 2nd real component -> b = a
1544
1545 x = a + q3;
1546 x = x * a + q2; // Q(a)
1547 x = x * a + q1;
1548 x = x * a + q0;
1549 y = -0.1875 * q3 * q3 + 0.5 * q2; // Q''(a) / 2
1550 z = Math.max(y * y - x, 0.0); // force discriminant to be >= 0
1551 z = Math.sqrt(z); // square root of discriminant
1552 y = y + sign(z, y); // larger magnitude root
1553 x = x / y; // smaller magnitude root
1554 c = Math.max(y, 0.0); // ensure root of biquadratic > 0
1555 d = Math.max(x, 0.0); // ensure root of biquadratic > 0
1556 c = Math.sqrt(c); // large magnitude imaginary component
1557 d = Math.sqrt(d); // small magnitude imaginary component
1558 }
1559
1560 if (a > b)
1561 {
1562 return new Complex[] { new Complex(a, c), new Complex(a, -c), new Complex(b, d), new Complex(b, -d) };
1563 }
1564 else if (a < b)
1565 {
1566 return new Complex[] { new Complex(b, d), new Complex(b, -d), new Complex(a, c), new Complex(a, -c) };
1567 }
1568 else
1569 {
1570 return new Complex[] { new Complex(a, c), new Complex(a, -c), new Complex(a, d), new Complex(a, -d) };
1571 }
1572
1573 } // # of real roots 'if'
1574 }
1575
1576 default:
1577 throw new RuntimeException("Bad switch; cannot happen");
1578 } // end select ! quartic type select
1579 }
1580 }