1 package org.djutils.polynomialroots;
2
3 import static org.junit.Assert.assertArrayEquals;
4 import static org.junit.Assert.assertEquals;
5
6 import org.djutils.complex.Complex;
7 import org.junit.Test;
8
9 /**
10 * TestPolynomialRoots.java.
11 * <p>
12 * Copyright (c) 2020-2020 Delft University of Technology, PO Box 5, 2600 AA, Delft, the Netherlands. All rights reserved. <br>
13 * BSD-style license. See <a href="https://djutils.org/docs/current/djutils/licenses.html">DJUTILS License</a>.
14 * </p>
15 * @author <a href="https://www.tudelft.nl/averbraeck">Alexander Verbraeck</a>
16 * @author <a href="https://www.tudelft.nl/pknoppers">Peter Knoppers</a>
17 */
18 public class TestPolynomialRoots
19 {
20
21 /**
22 * Test the solver for quadratic (and simpler) cases and ensure that the cubic solver falls back to the quadratic solver
23 * when appropriate.
24 */
25 @Test
26 public void quadraticTest()
27 {
28 double[] paramValues = new double[] { 0, 1, 2, 3, 4, Math.PI, -Math.E, 0.0001, 10000 };
29 for (double a : paramValues)
30 {
31 for (double b : paramValues)
32 {
33 for (double c : paramValues)
34 {
35 // System.out.print(a + " x^2 + " + b + " x " + c + " = 0");
36 Complex[] roots = PolynomialRoots.quadraticRoots(a, b, c);
37 int expectedNumber = a != 0 ? 2 : b != 0 ? 1 : 0;
38 assertEquals("number of roots", expectedNumber, roots.length);
39 // System.out.print(": " + expectedNumber + " solutions");
40 // Check that each root is indeed a root
41 for (Complex root : roots)
42 {
43 // System.out.print(" " + root);
44 Complex result = root.times(root).times(a).plus(root.times(b)).plus(c);
45 double margin = root.norm() / 1e6;
46 assertEquals("root is a root of the equation; re", 0, result.re, margin);
47 assertEquals("root is a root of the equation; im", 0, result.im, margin);
48 }
49 // System.out.println();
50 Complex[] alternate = PolynomialRoots.cubicRoots(0, a, b, c);
51 assertArrayEquals("cubic solver returns same results as quadratic solver if a is 0", roots, alternate);
52 }
53 }
54 }
55 double q1 = 3 * Math.sqrt(Double.MAX_VALUE);
56 Complex[] roots = PolynomialRoots.quadraticRoots(q1, 1);
57 assertEquals("two roots", 2, roots.length);
58 assertEquals("first root is 0: re", 0, roots[0].re, 0.0001);
59 assertEquals("first root is 0: im", 0, roots[0].im, 0.0001);
60 // System.out.println(roots[1]);
61 assertEquals("second root is 0: re", -q1, roots[1].re, q1 / 1e7);
62 assertEquals("second root is 0: im", 0, roots[1].im, 0.0001);
63 }
64
65 /**
66 * Test the cubic solver.
67 */
68 @Test
69 public void cubicTest()
70 {
71 double[] paramValues = new double[] { 0, 1, 2, 3, 4, Math.PI, -Math.E, 0.001, 1000 };
72 for (double a : paramValues)
73 {
74 for (double b : paramValues)
75 {
76 for (double c : paramValues)
77 {
78 for (double d : paramValues)
79 {
80 // System.out.print(a + " x^3 " + b + " x^2 + " + c + " x + " + d + " = 0");
81 Complex[] roots = PolynomialRoots.cubicRoots(a, b, c, d);
82 int expectedNumber = a != 0 ? 3 : b != 0 ? 2 : c != 0 ? 1 : 0;
83 assertEquals("number of roots", expectedNumber, roots.length);
84 // System.out.print(": " + expectedNumber + " solutions");
85 // Check that each root is indeed a root
86 for (Complex root : roots)
87 {
88 // System.out.print(" " + root);
89 Complex result = root.times(root).times(root).times(a).plus(root.times(root).times(b))
90 .plus(root.times(c)).plus(d);
91 double margin = root.norm() / 1e6;
92 assertEquals("root is a root of the equation; re", 0, result.re, margin);
93 assertEquals("root is a root of the equation; im", 0, result.im, margin);
94 }
95 // System.out.println();
96 }
97 }
98 }
99 }
100 }
101
102 /**
103 * Test the quartic solver.
104 */
105 @Test
106 public void quarticTest()
107 {
108 double[] paramValues = new double[] { 0, 1, 2, 3, 4, Math.PI, -Math.E, 0.01, 100 };
109 for (double a : paramValues)
110 {
111 for (double b : paramValues)
112 {
113 for (double c : paramValues)
114 {
115 for (double d : paramValues)
116 {
117 for (double e : paramValues)
118 {
119 // System.out.print(a + " x^4 " + b + " x^3 + " + c + " x^2 + " + d + " x + " + e + " = 0");
120 Complex[] roots = PolynomialRoots.quarticRoots(a, b, c, d, e);
121 int expectedNumber = a != 0 ? 4 : b != 0 ? 3 : c != 0 ? 2 : d != 0 ? 1 : 0;
122 assertEquals("number of roots", expectedNumber, roots.length);
123 // System.out.print(": " + expectedNumber + " solutions");
124 // Check that each root is indeed a root
125 for (Complex root : roots)
126 {
127 // System.out.print(" " + root);
128 // TODO compute result with less chance of loss of precision
129 Complex result = root.times(root).times(root).times(root).times(a)
130 .plus(root.times(root).times(root).times(b)).plus(root.times(root).times(c))
131 .plus(root.times(d)).plus(e);
132 double margin = root.norm() / 2e5;
133 assertEquals("root is a root of the equation; re", 0, result.re, margin);
134 assertEquals("root is a root of the equation; im", 0, result.im, margin);
135 }
136 // System.out.println();
137 }
138 }
139 }
140 }
141 }
142 }
143 }