1 package org.djutils.polynomialroots;
2
3 import org.djutils.complex.Complex;
4 import org.djutils.complex.ComplexMath;
5
6 /**
7 * PolynomialRoots2 implements functions to find all roots of linear, quadratic, cubic and quartic polynomials, as well as
8 * higher order polynomials without restrictions on the order. <br>
9 * <br>
10 * Copyright (c) 2020-2024 Delft University of Technology, Jaffalaan 5, 2628 BX Delft, the Netherlands. All rights reserved. See
11 * for project information <a href="https://djutils.org" target="_blank"> https://djutils.org</a>. The DJUTILS project is
12 * distributed under a three-clause BSD-style license, which can be found at
13 * <a href="https://djutils.org/docs/license.html" target="_blank"> https://djutils.org/docs/license.html</a>. <br>
14 * @author <a href="https://www.tudelft.nl/averbraeck">Alexander Verbraeck</a>
15 * @author <a href="https://www.tudelft.nl/pknoppers">Peter Knoppers</a>
16 */
17 public final class PolynomialRoots2
18 {
19 /**
20 * Do not instantiate.
21 */
22 private PolynomialRoots2()
23 {
24 // Do not instantiate
25 }
26
27 /**
28 * Emulate the F77 sign function.
29 * @param a double; the value to optionally sign invert
30 * @param b double; the sign of which determines what to do
31 * @return double; if b >= 0 then a; else -a
32 */
33 private static double sign(final double a, final double b)
34 {
35 return b >= 0 ? a : -a;
36 }
37
38 /**
39 * LINEAR POLYNOMIAL ROOT SOLVER.
40 * <p>
41 * Calculates the root of the linear polynomial:<br>
42 * q1 * x + q0<br>
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 * @param q0 double; independent coefficient
62 * @return Complex[]; the roots of the equation
63 */
64 public static Complex[] linearRoots(final double q0)
65 {
66 return new Complex[] {new Complex(-q0, 0)};
67 }
68
69 /**
70 * QUADRATIC POLYNOMIAL ROOT SOLVER
71 * <p>
72 * Calculates all real + complex roots of the quadratic polynomial:<br>
73 * q2 * x^2 + q1 * x + q0<br>
74 * The code checks internally if rescaling of the coefficients is needed to avoid overflow.
75 * <p>
76 * The order of the roots is as follows:<br>
77 * 1) For real roots, the order is according to their algebraic value on the number scale (largest positive first, largest
78 * negative last).<br>
79 * 2) Since there can be only one complex conjugate pair root, no order is necessary.<br>
80 * q1 : coefficient of x term q0 : independent coefficient
81 * @param q2 double; coefficient of the quadratic term
82 * @param q1 double; coefficient of the x term
83 * @param q0 double; independent coefficient
84 * @return Complex[]; the roots of the equation
85 */
86 public static Complex[] quadraticRoots(final double q2, final double q1, final double q0)
87 {
88 if (q2 == 0)
89 {
90 return linearRoots(q1, q0);
91 }
92 return quadraticRoots(q1 / q2, q0 / q2);
93 }
94
95 /**
96 * QUADRATIC POLYNOMIAL ROOT SOLVER
97 * <p>
98 * Calculates all real + complex roots of the quadratic polynomial:<br>
99 * x^2 + q1 * x + q0<br>
100 * The code checks internally if rescaling of the coefficients is needed to avoid overflow.
101 * <p>
102 * The order of the roots is as follows:<br>
103 * 1) For real roots, the order is according to their algebraic value on the number scale (largest positive first, largest
104 * negative last).<br>
105 * 2) Since there can be only one complex conjugate pair root, no order is necessary.<br>
106 * q1 : coefficient of x term q0 : independent coefficient
107 * @param q1 double; coefficient of the x term
108 * @param q0 double; independent coefficient
109 * @return Complex[]; the roots of the equation
110 */
111 public static Complex[] quadraticRoots(final double q1, final double q0)
112 {
113 boolean rescale;
114
115 double a0, a1;
116 double k = 0, x, y, z;
117
118 // Handle special cases.
119 if (q0 == 0.0 && q1 == 0.0)
120 {
121 // Two real roots at 0,0
122 return new Complex[] {Complex.ZERO, Complex.ZERO};
123 }
124 else if (q0 == 0.0)
125 {
126 // Two real roots; one of these is 0,0
127 // x^2 + q1 * x == x * (x + q1)
128 Complex nonZeroRoot = new Complex(-q1);
129 return new Complex[] {q1 > 0 ? Complex.ZERO : nonZeroRoot, q1 <= 0 ? nonZeroRoot : Complex.ZERO};
130 }
131 else if (q1 == 0.0)
132 {
133 x = Math.sqrt(Math.abs(q0));
134
135 if (q0 < 0.0)
136 {
137 // Two real roots, symmetrically around 0
138 return new Complex[] {new Complex(x, 0), new Complex(-x, 0)};
139 }
140 else
141 {
142 // Two complex roots, symmetrically around 0
143 return new Complex[] {new Complex(0, x), new Complex(0, -x)};
144 }
145 }
146 else
147 {
148 // The general case. Do rescaling, if either squaring of q1/2 or evaluation of
149 // (q1/2)^2 - q0 will lead to overflow. This is better than to have the solver
150 // crashed. Note, that rescaling might lead to loss of accuracy, so we only
151 // invoke it when absolutely necessary.
152 final double sqrtLPN = Math.sqrt(Double.MAX_VALUE); // Square root of the Largest Positive Number
153 rescale = (q1 > sqrtLPN + sqrtLPN); // this detects overflow of (q1/2)^2
154
155 if (!rescale)
156 {
157 x = q1 * 0.5; // we are sure here that x*x will not overflow
158 rescale = (q0 < x * x - Double.MAX_VALUE); // this detects overflow of (q1/2)^2 - q0
159 }
160
161 if (rescale)
162 {
163 x = Math.abs(q1);
164 y = Math.sqrt(Math.abs(q0));
165
166 if (x > y)
167 {
168 k = x;
169 z = 1.0 / x;
170 a1 = sign(1.0, q1);
171 a0 = (q0 * z) * z;
172 }
173 else
174 {
175 k = y;
176 a1 = q1 / y;
177 a0 = sign(1.0, q0);
178 }
179 }
180 else
181 {
182 a1 = q1;
183 a0 = q0;
184 }
185 // Determine the roots of the quadratic. Note, that either a1 or a0 might
186 // have become equal to zero due to underflow. But both cannot be zero.
187 x = a1 * 0.5;
188 y = x * x - a0;
189
190 if (y >= 0.0)
191 {
192 // Two real roots
193 y = Math.sqrt(y);
194
195 if (x > 0.0)
196 {
197 y = -x - y;
198 }
199 else
200 {
201 y = -x + y;
202 }
203
204 if (rescale)
205 {
206 y = y * k; // very important to convert to original
207 z = q0 / y; // root first, otherwise complete loss of
208 }
209 else // root due to possible a0 = 0 underflow
210 {
211 z = a0 / y;
212 }
213 return new Complex[] {new Complex(Math.max(y, z), 0), new Complex(Math.min(y, z), 0)};
214 }
215 else
216 {
217 // Two complex roots (zero real roots)
218 y = Math.sqrt(-y);
219
220 if (rescale)
221 {
222 x *= k;
223 y *= k;
224 }
225 return new Complex[] {new Complex(-x, y), new Complex(-x, -y)};
226 }
227 }
228 }
229
230 /**
231 * CUBIC POLYNOMIAL ROOT SOLVER.
232 * <p>
233 * Calculates all (real and complex) roots of the cubic polynomial:<br>
234 * a * x^3 + b * x^2 + c * x + d<br>
235 * The roots are found using the Newton-Raphson algorithm for the first (real) root, and then deflate the equation to a
236 * quadratic equation to find the other two roots.
237 * <p>
238 * The order of the roots is as follows: 1) For real roots, the order is according to their algebraic value on the number
239 * scale (largest positive first, largest negative last). 2) Since there can be only one complex conjugate pair root, no
240 * order is necessary. 3) All real roots precede the complex ones.
241 * @param a3 double; coefficient of the cubic term
242 * @param a2 double; coefficient of the quadratic term
243 * @param a1 double; coefficient of the linear term
244 * @param a0 double; coefficient of the independent term
245 * @return Complex[]; array of Complex with all the roots; there can be one, two or three roots.
246 */
247 public static Complex[] cubicRootsNewtonFactor(final double a3, final double a2, final double a1, final double a0)
248 {
249 // corner case: a == 0 --> quadratic equation
250 if (a3 == 0.0)
251 {
252 return quadraticRoots(a2, a1, a0);
253 }
254
255 // corner case: d == 0 --> solve x * (ax^2 + bx + c) = 0
256 if (a0 == 0.0)
257 {
258 Complex[] result = quadraticRoots(a3, a2, a1);
259 if (result.length == 0)
260 {
261 return new Complex[] {Complex.ZERO, Complex.ZERO, Complex.ZERO};
262 }
263 else if (result.length == 1)
264 {
265 return new Complex[] {Complex.ZERO, Complex.ZERO, result[0]};
266 }
267 else
268 {
269 return new Complex[] {Complex.ZERO, result[0], result[1]};
270 }
271 }
272
273 double argmax = maxAbs(a3, a2, a1, a0);
274 double d = a0 / argmax;
275 double c = a1 / argmax;
276 double b = a2 / argmax;
277 double a = a3 / argmax;
278
279 // find the first real root
280 double[] args = new double[] {d, c, b, a};
281 double root1 = rootNewtonRaphson(args, 0.0);
282
283 // check and apply bisection if this did not work
284 if (Double.isNaN(root1) || f(args, root1) > 1E-9)
285 {
286 double min = -64.0;
287 double max = 64.0;
288 int s1, s2;
289 do
290 {
291 min *= 2.0;
292 max *= 2.0;
293 if (max > 1E64)
294 {
295 throw new RuntimeException(
296 String.format("cannot find first root for %fx^3 + %fx^2 + %fx + %f = 0", a, b, c, d));
297 }
298 s1 = (int) Math.signum(f(args, min));
299 s2 = (int) Math.signum(f(args, max));
300 }
301 while (s1 != 0 && s2 != 0 && s1 == s2);
302 root1 = rootBisection(args, min, max, 0.0);
303 // if (Double.isNaN(root1))
304 // {
305 // throw new RuntimeException(
306 // String.format("cannot find first root for %fx^3 + %fx^2 + %fx + %f = 0", a, b, c, d));
307 // }
308 // if (Math.abs(f(args, root1)) > 1E-6)
309 // {
310 // throw new RuntimeException(String.format("f(first root) != 0 for [%fx^3 + %fx^2 + %fx + %f = 0, but %f]", a, b,
311 // c, d, f(args, root1)));
312 // }
313 }
314
315 /*-
316 * Use Factor theory to factor out root 1 (r):
317 * if the equation is ax^3 + bx^2 + cx + d, and there is a root r, the equation equals:
318 * (x-r) (ax^2 + (b+ra)x + (c+r(b+ra)) + (ar^3+br^2+cr+d)/(x-r) where the last term becomes zero
319 * We can then solve find the other two roots by solving ax^2 + (b+ra)x + (c+r(b+ra)) = 0
320 */
321
322 Complex[] rootsQuadaratic = quadraticRoots(a, b + root1 * a, c + root1 * (b + root1 * a));
323 // if (Math.abs(f(args, rootsQuadaratic[0]).norm()) > 1E-6)
324 // {
325 // throw new RuntimeException(String.format("f(second root) != 0 for [%fx^3 + %fx^2 + %fx + %f = 0], but %f", a, b, c,
326 // d, f(args, root1)));
327 // }
328 // if (Math.abs(f(args, rootsQuadaratic[1]).norm()) > 1E-6)
329 // {
330 // throw new RuntimeException(String.format("f(second root) != 0 for [%fx^3 + %fx^2 + %fx + %f = 0], but %f", a, b, c,
331 // d, f(args, root1)));
332 // }
333
334 return new Complex[] {new Complex(root1), rootsQuadaratic[0], rootsQuadaratic[1]};
335 }
336
337 /**
338 * CUBIC POLYNOMIAL ROOT SOLVER.
339 * <p>
340 * Calculates all (real and complex) roots of the cubic polynomial:<br>
341 * a * x^3 + b * x^2 + c * x + d<br>
342 * The roots are found using Cardano's algorithm.
343 * <p>
344 * The order of the roots is as follows: 1) For real roots, the order is according to their algebraic value on the number
345 * scale (largest positive first, largest negative last). 2) Since there can be only one complex conjugate pair root, no
346 * order is necessary. 3) All real roots precede the complex ones.
347 * @param a double; coefficient of the cubic term
348 * @param b double; coefficient of the quadratic term
349 * @param c double; coefficient of the linear term
350 * @param d double; coefficient of the independent term
351 * @return Complex[]; array of Complex with all the roots; there can be one, two or three roots.
352 */
353 @SuppressWarnings("checkstyle:localvariablename")
354 public static Complex[] cubicRootsCardano(final double a, final double b, final double c, final double d)
355 {
356 // corner case: a == 0 --> quadratic equation
357 if (a == 0.0)
358 {
359 return quadraticRoots(b, c, d);
360 }
361
362 // corner case: d == 0 --> solve x * (ax^2 + bx + c) = 0
363 if (d == 0.0)
364 {
365 Complex[] result = quadraticRoots(a, b, c);
366 if (result.length == 0)
367 {
368 return new Complex[] {Complex.ZERO, Complex.ZERO, Complex.ZERO};
369 }
370 else if (result.length == 1)
371 {
372 return new Complex[] {Complex.ZERO, Complex.ZERO, result[0]};
373 }
374 else
375 {
376 return new Complex[] {Complex.ZERO, result[0], result[1]};
377 }
378 }
379
380 // other cases: use Cardano's formula
381 double D0 = b * b - 3.0 * a * c;
382 double D1 = 2.0 * b * b * b - 9.0 * a * b * c + 27.0 * a * a * d;
383 if (D0 == 0 && D1 == 0)
384 {
385 // one real solution
386 Complex root = new Complex(rootNewtonRaphson(new double[] {d, c, b, a}, 0.0));
387 return new Complex[] {root, root, root};
388 }
389 Complex r = ComplexMath.sqrt(new Complex(D1 * D1 - 4.0 * D0 * D0 * D0));
390 Complex s = r.plus(D1).times(0.5);
391 if (s.re == 0.0 && s.im == 0.0)
392 {
393 s = r.times(-1.0).plus(D1).times(0.5);
394 }
395 Complex C = ComplexMath.cbrt(s);
396 Complex x1 = C.plus(b).plus((C.reciprocal().times(D0))).times(-1.0 / (3.0 * a));
397 Complex x2 = C.rotate(Math.toRadians(120.0)).plus(b).plus((C.rotate(Math.toRadians(120.0)).reciprocal().times(D0)))
398 .times(-1.0 / (3.0 * a));
399 Complex x3 = C.rotate(Math.toRadians(-120.0)).plus(b).plus((C.rotate(Math.toRadians(-120.0)).reciprocal().times(D0)))
400 .times(-1.0 / (3.0 * a));
401 return new Complex[] {x1, x2, x3};
402 }
403
404 /**
405 * CUBIC POLYNOMIAL ROOT SOLVER.
406 * <p>
407 * Calculates all (real and complex) roots of the cubic polynomial:<br>
408 * a3 * x^3 + a2 * x^2 + a1 * x + a0<br>
409 * The roots are found using Durand-Kerner's algorithm.
410 * @param a3 double; coefficient of the cubic term
411 * @param a2 double; coefficient of the quadratic term
412 * @param a1 double; coefficient of the linear term
413 * @param a0 double; coefficient of the independent term
414 * @return Complex[]; array of Complex with all the roots; there can be one, two or three roots.
415 */
416 public static Complex[] cubicRootsDurandKerner(final double a3, final double a2, final double a1, final double a0)
417 {
418 // corner case: a3 == 0 --> quadratic equation
419 if (a3 == 0.0)
420 {
421 return quadraticRoots(a2, a1, a0);
422 }
423
424 // other cases: use Durand-Kerner's algorithm
425 return rootsDurandKerner(new Complex[] {new Complex(a0 / a3), new Complex(a1 / a3), new Complex(a2 / a3), Complex.ONE});
426 }
427
428 /**
429 * CUBIC POLYNOMIAL ROOT SOLVER.
430 * <p>
431 * Calculates all (real and complex) roots of the cubic polynomial:<br>
432 * a3 * x^3 + a2 * x^2 + a1 * x + a0<br>
433 * The roots are found using Aberth-Ehrlich's algorithm.
434 * @param a3 double; coefficient of the cubic term
435 * @param a2 double; coefficient of the quadratic term
436 * @param a1 double; coefficient of the linear term
437 * @param a0 double; coefficient of the independent term
438 * @return Complex[]; array of Complex with all the roots; there can be one, two or three roots.
439 */
440 public static Complex[] cubicRootsAberthEhrlich(final double a3, final double a2, final double a1, final double a0)
441 {
442 // corner case: a3 == 0 --> quadratic equation
443 if (a3 == 0.0)
444 {
445 return quadraticRoots(a2, a1, a0);
446 }
447
448 // other cases: use Durand-Kerner's algorithm
449 return rootsAberthEhrlich(
450 new Complex[] {new Complex(a0 / a3), new Complex(a1 / a3), new Complex(a2 / a3), Complex.ONE});
451 }
452
453 /**
454 * QUADRATIC POLYNOMIAL ROOT SOLVER.
455 * <p>
456 * Calculates all (real and complex) roots of the cubic polynomial:<br>
457 * a4 * x^4 + a3 * x^3 + a2 * x^2 + a1 * x + a0<br>
458 * The roots are found using Durand-Kerner's algorithm.
459 * @param a4 double; coefficient of the quartic term
460 * @param a3 double; coefficient of the cubic term
461 * @param a2 double; coefficient of the quadratic term
462 * @param a1 double; coefficient of the linear term
463 * @param a0 double; coefficient of the independent term
464 * @return Complex[]; array of Complex with all the roots; always 4 roots are returned; there can be double roots
465 */
466 public static Complex[] quarticRootsDurandKerner(final double a4, final double a3, final double a2, final double a1,
467 final double a0)
468 {
469 // corner case: a4 == 0 --> cubic equation
470 if (a4 == 0.0)
471 {
472 return cubicRootsDurandKerner(a3, a2, a1, a0);
473 }
474
475 return rootsDurandKerner(new Complex[] {new Complex(a0 / a4), new Complex(a1 / a4), new Complex(a2 / a4),
476 new Complex(a3 / a4), Complex.ONE});
477 }
478
479 /**
480 * QUADRATIC POLYNOMIAL ROOT SOLVER.
481 * <p>
482 * Calculates all (real and complex) roots of the cubic polynomial:<br>
483 * a4 * x^4 + a3 * x^3 + a2 * x^2 + a1 * x + a0<br>
484 * The roots are found using Aberth-Ehrlich's algorithm.
485 * @param a4 double; coefficient of the quartic term
486 * @param a3 double; coefficient of the cubic term
487 * @param a2 double; coefficient of the quadratic term
488 * @param a1 double; coefficient of the linear term
489 * @param a0 double; coefficient of the independent term
490 * @return Complex[]; array of Complex with all the roots; always 4 roots are returned; there can be double roots
491 */
492 public static Complex[] quarticRootsAberthEhrlich(final double a4, final double a3, final double a2, final double a1,
493 final double a0)
494 {
495 // corner case: a4 == 0 --> cubic equation
496 if (a4 == 0.0)
497 {
498 return cubicRootsAberthEhrlich(a3, a2, a1, a0);
499 }
500
501 return rootsAberthEhrlich(new Complex[] {new Complex(a0 / a4), new Complex(a1 / a4), new Complex(a2 / a4),
502 new Complex(a3 / a4), Complex.ONE});
503 }
504
505 /** the number of steps in solving the equation with the Durand-Kerner method. */
506 private static final int MAX_STEPS_DURAND_KERNER = 100;
507
508 /**
509 * Polynomial root finder using the Durand-Kerner method, with complex coefficients for the polynomial equation. Own
510 * implementation. See <a href="https://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method">
511 * https://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method</a> for brief information.
512 * @param a Complex[]; the complex factors of the polynomial, where the index i indicate the factor for x^i, so the
513 * polynomial is<br>
514 * a[n]x^n + a[n-1]x^(n-1) + ... + a[2]a^2 + a[1]x + a[0]
515 * @return Complex[] all roots of the equation, where real roots are coded with Im = 0
516 */
517 public static Complex[] rootsDurandKerner(final Complex[] a)
518 {
519 int n = a.length - 1;
520 double radius = 1 + maxAbs(a);
521
522 // initialize the initial values, not as a real number and not as a root of unity
523 Complex[] p = new Complex[n];
524 p[0] = new Complex(Math.sqrt(radius), Math.cbrt(radius));
525 double rot = 350.123 / n;
526 for (int i = 1; i < n; i++)
527 {
528 p[i] = p[0].rotate(rot * i);
529 }
530
531 double maxError = 1.0;
532 int count = 0;
533 while (maxError > 0 && count < MAX_STEPS_DURAND_KERNER)
534 {
535 maxError = 0.0;
536 count++;
537 for (int i = 0; i < n; i++)
538 {
539 Complex factors = Complex.ONE;
540 for (int j = 0; j < n; j++)
541 {
542 if (i != j)
543 {
544 factors = factors.times(p[i].minus(p[j]));
545 }
546 }
547 Complex newValue = p[i].minus(f(a, p[i]).divideBy(factors));
548 if (!(newValue.equals(p[i])))
549 {
550 double error = newValue.minus(p[i]).norm();
551 if (error > maxError)
552 {
553 maxError = error;
554 }
555 p[i] = newValue;
556 }
557 }
558 }
559
560 // correct for 1 ulp above or below zero; make that value zero instead
561 for (int i = 0; i < n; i++)
562 {
563 if (Math.abs(p[i].im) == Double.MIN_VALUE)
564 {
565 p[i] = new Complex(p[i].re, 0.0);
566 }
567 if (Math.abs(p[i].re) == Double.MIN_VALUE)
568 {
569 p[i] = new Complex(0.0, p[i].im);
570 }
571 }
572
573 return p;
574 }
575
576 /** the number of steps in solving the equation with the Durand-Kerner method. */
577 private static final int MAX_STEPS_ABERTH_EHRLICH = 100;
578
579 /**
580 * Polynomial root finder using the Aberth-Ehrlich method or Aberth method, with complex coefficients for the polynomial
581 * equation. Own implementation. See <a href="https://en.wikipedia.org/wiki/Aberth_method">
582 * https://en.wikipedia.org/wiki/Aberth_method</a> for brief information.
583 * @param a Complex[]; the complex factors of the polynomial, where the index i indicate the factor for x^i, so the
584 * polynomial is<br>
585 * a[n]x^n + a[n-1]x^(n-1) + ... + a[2]a^2 + a[1]x + a[0]
586 * @return Complex[] all roots of the equation, where real roots are coded with Im = 0
587 */
588 public static Complex[] rootsAberthEhrlich(final Complex[] a)
589 {
590 int n = a.length - 1;
591 double radius = 1 + maxAbs(a);
592
593 // initialize the initial values, not as a real number and not as a root of unity
594 Complex[] p = new Complex[n];
595 p[0] = new Complex(Math.sqrt(radius), Math.cbrt(radius));
596 double rot = 350.123 / n;
597 for (int i = 1; i < n; i++)
598 {
599 p[i] = p[0].rotate(rot * i);
600 }
601
602 double maxError = 1.0;
603 int count = 0;
604 while (maxError > 0 && count < MAX_STEPS_ABERTH_EHRLICH)
605 {
606 maxError = 0.0;
607 count++;
608 for (int i = 0; i < n; i++)
609 {
610 Complex sum = Complex.ZERO;
611 for (int j = 0; j < n; j++)
612 {
613 if (i != j)
614 {
615 sum = sum.plus(Complex.ONE.divideBy(p[i].minus(p[j])));
616 }
617 }
618 Complex pDivPDer = f(a, p[i]).divideBy(fDerivative(a, p[i]));
619 Complex w = pDivPDer.divideBy(Complex.ONE.minus(pDivPDer.times(sum)));
620 double error = w.norm();
621 if (error > maxError)
622 {
623 maxError = error;
624 }
625 p[i] = p[i].minus(w);
626 }
627 }
628
629 // correct for 1 ulp above or below zero; make that value zero instead
630 for (int i = 0; i < n; i++)
631 {
632 if (Math.abs(p[i].im) == Double.MIN_VALUE)
633 {
634 p[i] = new Complex(p[i].re, 0.0);
635 }
636 if (Math.abs(p[i].re) == Double.MIN_VALUE)
637 {
638 p[i] = new Complex(0.0, p[i].im);
639 }
640 }
641
642 return p;
643 }
644
645 /** the number of steps in approximating a real root with the Newton-Raphson method. */
646 private static final int MAX_STEPS_NEWTON = 100;
647
648 /**
649 * Polynomial root finder using the Newton-Raphson method. See <a href="https://en.wikipedia.org/wiki/Newton%27s_method">
650 * https://en.wikipedia.org/wiki/Newton%27s_method</a> for brief information about the Newton-Raphson or Newton method for
651 * root finding.
652 * @param a double[]; the factors of the polynomial, where the index i indicate the factor for x^i, so the polynomial is<br>
653 * a[n]x^n + a[n-1]x^(n-1) + ... + a[2]a^2 + a[1]x + a[0]
654 * @param allowedError double; the allowed absolute error in the result
655 * @return double the root of the equation that has been found on the basis of the start value, or NaN if not found within
656 * the allowed error bounds and the allowed number of steps.
657 */
658 public static double rootNewtonRaphson(final double[] a, final double allowedError)
659 {
660 double x = 0.1232232323234; // take a a start point that has an almost zero chance of having f'(x) = 0.
661 double newValue = Double.NaN; // for testing convergence
662 double fx = -1;
663 for (int j = 0; j < MAX_STEPS_NEWTON; j++)
664 {
665 fx = f(a, x);
666 newValue = x - fx / fDerivative(a, x);
667 if (x == newValue || Math.abs(fx) <= allowedError)
668 {
669 return x;
670 }
671 x = newValue;
672 }
673 return Math.abs(fx) <= allowedError ? x : Double.NaN;
674 }
675
676 /** the number of steps in approximating a real root with the bisection method. */
677 private static final int MAX_STEPS_BISECTION = 100;
678
679 /**
680 * Polynomial root finder using the bisection method combined with bisection to avoid cycles and the algorithm going out of
681 * bounds. Implementation based on Numerical Recipes in C, section 9.4, pp 365-366. See
682 * <a href="https://en.wikipedia.org/wiki/Newton%27s_method"> https://en.wikipedia.org/wiki/Newton%27s_method</a> for brief
683 * information about the Newton-Raphson or Newton method for root finding.
684 * @param a double[]; the factors of the polynomial, where the index i indicate the factor for x^i, so the polynomial is<br>
685 * a[n]x^n + a[n-1]x^(n-1) + ... + a[2]a^2 + a[1]x + a[0]
686 * @param startMin double; the lowest initial search value
687 * @param startMax double; the highest initial search value
688 * @param allowedError double; the allowed absolute error in the result
689 * @return double the root of the equation that has been found on the basis of the start values, or NaN if not found within
690 * the allowed error bounds and the allowed number of steps.
691 */
692 public static double rootBisection(final double[] a, final double startMin, final double startMax,
693 final double allowedError)
694 {
695 int n = 0;
696 double xPrev = startMin;
697 double min = startMin;
698 double max = startMax;
699 double fmin = f(a, min);
700 while (n <= MAX_STEPS_BISECTION)
701 {
702 double x = (min + max) / 2.0;
703 double fx = f(a, x);
704 if (x == xPrev || Math.abs(fx) < allowedError)
705 {
706 return x;
707 }
708 if (Math.signum(fx) == Math.signum(fmin))
709 {
710 min = x;
711 fmin = fx;
712 }
713 else
714 {
715 max = x;
716 }
717 xPrev = x;
718 n++;
719 }
720 return Double.NaN;
721 }
722
723 /**
724 * Return the max of the norm of the complex coefficients in an array.
725 * @param array Complex[] the array with complex numbers
726 * @return double; the highest value of the norm of the complex numbers in the array
727 */
728 private static double maxAbs(final Complex[] array)
729 {
730 double m = Double.NEGATIVE_INFINITY;
731 for (Complex c : array)
732 {
733 if (c.norm() > m)
734 {
735 m = c.norm();
736 }
737 }
738 return m;
739 }
740
741 /**
742 * Return the max of the absolute values of the coefficients in an array.
743 * @param values double[] the array with numbers
744 * @return double; the highest absolute value of the norm of the complex numbers in the array
745 */
746 private static double maxAbs(final double... values)
747 {
748 double m = Double.NEGATIVE_INFINITY;
749 for (double d : values)
750 {
751 if (Math.abs(d) > m)
752 {
753 m = Math.abs(d);
754 }
755 }
756 return m;
757 }
758
759 /**
760 * Return the complex value of f(c) where f(x) = a[n]x^n + a[n-1]x^(n-1) + ... + a[2]a^2 + a[1]x + a[0].
761 * @param a Complex[] the complex factors of the equation
762 * @param c Complex; the value for which to calculate f(c)
763 * @return Complex; f(c)
764 */
765 private static Complex f(final Complex[] a, final Complex c)
766 {
767 Complex cPow = Complex.ONE;
768 Complex result = Complex.ZERO;
769 for (int i = 0; i < a.length; i++)
770 {
771 result = result.plus(cPow.times(a[i]));
772 cPow = cPow.times(c);
773 }
774 return result;
775 }
776
777 /**
778 * Return the real value of f(x) where f(x) = a[n]x^n + a[n-1]x^(n-1) + ... + a[2]a^2 + a[1]x + a[0].
779 * @param a double[] the factors of the equation
780 * @param x double; the value for which to calculate f(x)
781 * @return double; f(x)
782 */
783 private static double f(final double[] a, final double x)
784 {
785 double xPow = 1.0;
786 double result = 0.0;
787 for (int i = 0; i < a.length; i++)
788 {
789 result += xPow * a[i];
790 xPow = xPow * x;
791 }
792 return result;
793 }
794
795 /**
796 * Return the complex value of f(c) where f(x) = a[n]x^n + a[n-1]x^(n-1) + ... + a[2]a^2 + a[1]x + a[0].
797 * @param a double[] the complex factors of the equation
798 * @param c Complex; the value for which to calculate f(c)
799 * @return Complex; f(c)
800 */
801 private static Complex f(final double[] a, final Complex c)
802 {
803 Complex cPow = Complex.ONE;
804 Complex result = Complex.ZERO;
805 for (int i = 0; i < a.length; i++)
806 {
807 result = result.plus(cPow.times(a[i]));
808 cPow = cPow.times(c);
809 }
810 return result;
811 }
812
813 /**
814 * Return the real value of f'(x) where f(x) = a[n]x^n + a[n-1]x^(n-1) + ... + a[2]a^2 + a[1]x + a[0].<br>
815 * The derivative function f'(x) = n.a[n]x^(n-1) + (n-1).a[n-1]x^(n-2) + ... + 2.a[2]x^1 + 1.a[1]
816 * @param a double[] the factors of the equation
817 * @param x double; the value for which to calculate f'(x)
818 * @return double; f'(x), the value of the derivative function at point x
819 */
820 private static double fDerivative(final double[] a, final double x)
821 {
822 double xPow = 1.0;
823 double result = 0.0;
824 for (int i = 1; i < a.length; i++)
825 {
826 result += xPow * i * a[i];
827 xPow = xPow * x;
828 }
829 return result;
830 }
831
832 /**
833 * Return the complex value of f'(c) where f(c) = a[n]c^n + a[n-1]c^(n-1) + ... + a[2]a^2 + a[1]c + a[0].<br>
834 * The derivative function f'(c) = n.a[n]c^(n-1) + (n-1).a[n-1]c^(n-2) + ... + 2.a[2]c^1 + 1.a[1]
835 * @param a double[] the factors of the equation
836 * @param c double; the value for which to calculate f'(c)
837 * @return Complex; f'(c), the value of the derivative function at point c
838 */
839 private static Complex fDerivative(final Complex[] a, final Complex c)
840 {
841 Complex cPow = Complex.ONE;
842 Complex result = Complex.ZERO;
843 for (int i = 1; i < a.length; i++)
844 {
845 result = result.plus(cPow.times(a[i]).times(i));
846 cPow = cPow.times(c);
847 }
848 return result;
849 }
850
851 /**
852 * @param args String[] not used
853 */
854 public static void main(final String[] args)
855 {
856 double a = 0.001;
857 double b = 1000;
858 double c = 0;
859 double d = 1;
860 Complex[] roots = cubicRootsNewtonFactor(a, b, c, d);
861 for (Complex root : roots)
862 {
863 System.out.println("NewtonFactor " + root + " f(x) = " + f(new double[] {d, c, b, a}, root));
864 }
865 System.out.println();
866 roots = cubicRootsCardano(a, b, c, d);
867 for (Complex root : roots)
868 {
869 System.out.println("Cardano " + root + " f(x) = " + f(new double[] {d, c, b, a}, root));
870 }
871 System.out.println();
872 roots = cubicRootsAberthEhrlich(a, b, c, d);
873 for (Complex root : roots)
874 {
875 System.out.println("Aberth-Ehrlich " + root + " f(x) = " + f(new double[] {d, c, b, a}, root));
876 }
877 System.out.println();
878 roots = cubicRootsDurandKerner(a, b, c, d);
879 for (Complex root : roots)
880 {
881 System.out.println("Durand-Kerner " + root + " f(x) = " + f(new double[] {d, c, b, a}, root));
882 }
883 System.out.println();
884 roots = PolynomialRoots.cubicRoots(a, b, c, d);
885 for (Complex root : roots)
886 {
887 System.out.println("F77 " + root + " f(x) = " + f(new double[] {d, c, b, a}, root));
888 }
889 }
890 }