1 /*
2 * Copyright (C) 2012 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.numerical.roots;
17
18 import com.irurueta.algebra.Complex;
19 import com.irurueta.numerical.LockedException;
20 import com.irurueta.numerical.NotAvailableException;
21 import com.irurueta.numerical.NotReadyException;
22
23 /**
24 * Class to estimate the roots of a third degree polynomial along with other
25 * polynomial properties.
26 * A second degree polynomial is defined by its parameters as p(x) = a * x^3 +
27 * b * x^2 + c * x + d, hence the polynomial can be simply be defined by an
28 * array of length 4 [d, c, b, a]
29 * This class is based on:
30 * <a href="http://en.wikipedia.org/wiki/Cubic_function">http://en.wikipedia.org/wiki/Cubic_function</a>
31 */
32 @SuppressWarnings("DuplicatedCode")
33 public class ThirdDegreePolynomialRootsEstimator extends PolynomialRootsEstimator {
34
35 /**
36 * Constant defining machine precision.
37 */
38 public static final double EPS = 1e-10;
39
40 /**
41 * Constant defining one third.
42 */
43 public static final double THIRD = 1.0 / 3.0;
44
45 /**
46 * Constant defining the squared root of three.
47 */
48 public static final double ROOT_THREE = Math.sqrt(3.0);
49
50 /**
51 * Number of parameters valid for a third degree polynomial.
52 */
53 public static final int VALID_POLY_PARAMS_LENGTH = 4;
54
55 /**
56 * Array containing parameters of a second degree polynomial.
57 */
58 private double[] realPolyParams;
59
60 /**
61 * Empty constructor.
62 */
63 public ThirdDegreePolynomialRootsEstimator() {
64 super();
65 realPolyParams = null;
66 }
67
68 /**
69 * Constructor.
70 *
71 * @param polyParams Array containing polynomial parameters.
72 * @throws IllegalArgumentException Raised if the length of the provided
73 * array is not valid.
74 */
75 public ThirdDegreePolynomialRootsEstimator(final double[] polyParams) {
76 super();
77 internalSetPolynomialParameters(polyParams);
78 }
79
80 /**
81 * Set array of third degree polynomial parameters.
82 * A third degree polynomial is defined by p(x) = a * x^3 + b * x^2 + c * x
83 * + d, and the array must be provided as [d, c, b, a].
84 * Note: This class only supports real polynomial parameters
85 *
86 * @param polyParams Array containing polynomial parameters
87 * @throws LockedException Raised if this instance is locked
88 * @throws IllegalArgumentException Raised if the length of the provided
89 * array is not valid.
90 */
91 public void setPolynomialParameters(final double[] polyParams) throws LockedException {
92 if (isLocked()) {
93 throw new LockedException();
94 }
95 internalSetPolynomialParameters(polyParams);
96 }
97
98 /**
99 * Internal method to set array of third degree polynomial parameters.
100 * A third degree polynomial is defined by p(x) = a * x^3 + b * x^2 + c * d
101 * + d, and the array must be provided as [d, c, b, a].
102 * Note: This class only supports real polynomial parameters
103 * This method does not check if this instance is locked.
104 *
105 * @param polyParams Array containing polynomial parameters
106 * @throws IllegalArgumentException Raised if the length of the provided
107 * array is not valid.
108 */
109 private void internalSetPolynomialParameters(final double[] polyParams) {
110 if (polyParams.length < VALID_POLY_PARAMS_LENGTH) {
111 throw new IllegalArgumentException();
112 }
113 if (!isThirdDegree(polyParams)) {
114 throw new IllegalArgumentException();
115 }
116
117 this.realPolyParams = polyParams;
118 }
119
120 /**
121 * Returns array of third degree polynomial parameters.
122 * A third degree polynomial is defined by p(x) = a * x^3 + b * x^2 + c * x
123 * + d, and the array is returned as [d, c, b, a].
124 * Note: This class only supports real polynomial parameters.
125 *
126 * @return Array of first degree polynomial parameters.
127 * @throws NotAvailableException Raised if polynomial parameter have not yet
128 * been provided.
129 */
130 public double[] getRealPolynomialParameters() throws NotAvailableException {
131 if (!arePolynomialParametersAvailable()) {
132 throw new NotAvailableException();
133 }
134 return realPolyParams;
135 }
136
137 /**
138 * Returns boolean indicating whether REAL polynomial parameters have been
139 * provided and is available for retrieval.
140 * Note: This class only supports real polynomial parameters.
141 *
142 * @return True if available, false otherwise.
143 */
144 @Override
145 public boolean arePolynomialParametersAvailable() {
146 return realPolyParams != null;
147 }
148
149 /**
150 * This method will always raise a NotAvailableException because this class
151 * only supports REAL polynomial parameters.
152 *
153 * @return this method always throws an exception.
154 * @throws com.irurueta.numerical.NotAvailableException always thrown
155 */
156 @Override
157 public Complex[] getPolynomialParameters() throws NotAvailableException {
158 throw new NotAvailableException();
159 }
160
161 /**
162 * Estimates the roots of provided polynomial.
163 *
164 * @throws LockedException Raised if this instance is locked estimating
165 * roots.
166 * @throws NotReadyException Raised if this instance is not ready because
167 * polynomial parameters have not been provided.
168 * @throws RootEstimationException Raised if roots cannot be estimated for
169 * some reason.
170 */
171 @Override
172 public void estimate() throws LockedException, NotReadyException, RootEstimationException {
173 if (isLocked()) {
174 throw new LockedException();
175 }
176 if (!isReady()) {
177 throw new NotReadyException();
178 }
179
180 locked = true;
181
182 roots = new Complex[VALID_POLY_PARAMS_LENGTH - 1];
183
184 final var d = realPolyParams[0];
185 final var c = realPolyParams[1];
186 final var b = realPolyParams[2];
187 final var a = realPolyParams[3];
188
189 final var x1 = new Complex();
190 final var x2 = new Complex();
191 final var x3 = new Complex();
192 solveCubic(a, b, c, d, x1, x2, x3);
193
194 if (Double.isNaN(x1.getReal()) || Double.isNaN(x1.getImaginary()) || Double.isNaN(x2.getReal())
195 || Double.isNaN(x2.getImaginary()) || Double.isNaN(x3.getReal()) || Double.isNaN(x3.getImaginary())) {
196
197 locked = false;
198 throw new RootEstimationException();
199 }
200
201 if (x1.getReal() < x2.getReal() && x1.getReal() < x3.getReal()) {
202 // x1 goes first
203 roots[0] = x1;
204 if (x2.getReal() < x3.getReal()) {
205 // x2 goes second and x3 goes third
206 roots[1] = x2;
207 roots[2] = x3;
208 } else {
209 // x3 goes second and x2 goes third
210 roots[1] = x3;
211 roots[2] = x2;
212 }
213 } else if (x2.getReal() < x1.getReal() && x2.getReal() < x3.getReal()) {
214 // x2 goes first
215 roots[0] = x2;
216 if (x1.getReal() < x3.getReal()) {
217 // x1 goes second and x3 goes third
218 roots[1] = x1;
219 roots[2] = x3;
220 } else {
221 // x3 goes second and x1 goes third
222 roots[1] = x3;
223 roots[2] = x1;
224 }
225 } else {
226 // x3 goes first
227 roots[0] = x3;
228 if (x1.getReal() < x2.getReal()) {
229 // x1 goes second and x2 goes third
230 roots[1] = x1;
231 roots[2] = x2;
232 } else {
233 // x2 goes second and x1 goes second
234 roots[1] = x2;
235 roots[2] = x1;
236 }
237 }
238
239 locked = false;
240 }
241
242 /**
243 * Returns boolean indicating whether provided array of polynomial
244 * parameters correspond to a valid third degree polynomial.
245 * A third degree polynomial is defined by p(x) = a * x^3 + b * x^2 + c *x +
246 * d, and the array is returned as [d, c, b, a].
247 * Note: This class only supports real polynomial parameters
248 *
249 * @param polyParams Array containing polynomial parameters
250 * @return True if is a third degree polynomial, false otherwise
251 */
252 public static boolean isThirdDegree(final double[] polyParams) {
253 final var length = polyParams.length;
254 if (length >= VALID_POLY_PARAMS_LENGTH && Math.abs(polyParams[VALID_POLY_PARAMS_LENGTH - 1]) > EPS) {
255 for (var i = VALID_POLY_PARAMS_LENGTH; i < length; i++) {
256 if (Math.abs(polyParams[i]) > EPS) {
257 return false;
258 }
259 }
260 return true;
261 }
262 return false;
263 }
264
265 /**
266 * Returns boolean indicating whether polynomial parameters provided to this
267 * instance correspond to a valid third degree polynomial.
268 * A third degree polynomial is defined by p(x) = a * x^3 + b * x^2 + c * x
269 * + d, and the array is returned as [d, c, b, a].
270 * Note: This class only supports real polynomial parameters
271 *
272 * @return True if is a second degree polynomial, false otherwise
273 * @throws NotReadyException Raised if this instance is not ready because
274 * an array of polynomial parameters has not yet been provided.
275 */
276 public boolean isThirdDegree() throws NotReadyException {
277 if (!isReady()) {
278 throw new NotReadyException();
279 }
280 return isThirdDegree(realPolyParams);
281 }
282
283 /**
284 * Returns boolean indicating whether the roots of the polynomial are three
285 * distinct and real roots or not.
286 * Because this class only supports polynomials with real parameters, we
287 * know that for third degree polynomials that have three distinct roots,
288 * they must be either real or one real and 2 complex conjugate.
289 *
290 * @param polyParams Array containing polynomial parameters
291 * @return True if roots are distinct and real, false otherwise
292 */
293 public static boolean hasThreeDistinctRealRoots(final double[] polyParams) {
294 if (polyParams.length >= VALID_POLY_PARAMS_LENGTH) {
295 return getDiscriminant(polyParams) > EPS;
296 }
297 return false;
298 }
299
300 /**
301 * Returns boolean indicating whether the roots of the polynomial are three
302 * distinct and real roots or not.
303 * Because this class only supports polynomials with real parameters, we
304 * know that for third degree polynomials that have three distinct roots,
305 * they must be either real or one real and 2 complex conjugate.
306 *
307 * @return True if roots are distinct and real, false otherwise
308 * @throws NotReadyException Raised if polynomial parameters haven't yet
309 * been provided
310 */
311 public boolean hasThreeDistinctRealRoots() throws NotReadyException {
312 if (!isReady()) {
313 throw new NotReadyException();
314 }
315 return hasThreeDistinctRealRoots(realPolyParams);
316 }
317
318 /**
319 * Returns boolean indicating whether the polynomial has two real and equal
320 * roots and a third different one (multiplicity 2), or all three roots are
321 * real and equal (multiplicity 3).
322 *
323 * @param polyParams Array containing polynomial parameters
324 * @return True if there are roots with multiplicity greater than one, false
325 * otherwise
326 */
327 public static boolean hasMultipleRealRoot(final double[] polyParams) {
328 if (polyParams.length >= VALID_POLY_PARAMS_LENGTH) {
329 return Math.abs(getDiscriminant(polyParams)) <= EPS;
330 }
331 return false;
332 }
333
334 /**
335 * Returns boolean indicating whether the polynomial has two real and equal
336 * roots and a third different one (multiplicity 2), or all three roots are
337 * real and equal (multiplicity 3).
338 *
339 * @return True if there are roots with multiplicity greater than one, false
340 * otherwise
341 * @throws NotReadyException Raised if polynomial parameters haven't yet
342 * been provided
343 */
344 public boolean hasMultipleRealRoot() throws NotReadyException {
345 if (!isReady()) {
346 throw new NotReadyException();
347 }
348 return hasMultipleRealRoot(realPolyParams);
349 }
350
351 /**
352 * Returns boolean indicating whether the polynomial has one real root and
353 * two complex conjugate roots.
354 *
355 * @param polyParams Array containing polynomial parameters
356 * @return True if polynomial has 1 real root and 2 complex conjugate roots,
357 * false otherwise
358 */
359 public static boolean hasOneRealRootAndTwoComplexConjugateRoots(final double[] polyParams) {
360 if (polyParams.length >= VALID_POLY_PARAMS_LENGTH) {
361 return getDiscriminant(polyParams) < -EPS;
362 }
363 return false;
364 }
365
366 /**
367 * Returns boolean indicating whether the polynomial has one real root and
368 * two complex conjugate roots.
369 *
370 * @return True if polynomial has 1 real root and 2 complex conjugate roots,
371 * false otherwise
372 * @throws NotReadyException Raised if polynomial parameters haven't yet
373 * been provided
374 */
375 public boolean hasOneRealRootAndTwoComplexConjugateRoots() throws NotReadyException {
376 if (!isReady()) {
377 throw new NotReadyException();
378 }
379 return hasOneRealRootAndTwoComplexConjugateRoots(realPolyParams);
380 }
381
382 /**
383 * This method will always raise an IllegalArgumentException because this
384 * class only supports REAL polynomial parameters.
385 *
386 * @throws IllegalArgumentException always thrown.
387 */
388 @Override
389 protected void internalSetPolynomialParameters(final Complex[] polyParams) {
390 // complex values are not supported
391 throw new IllegalArgumentException();
392 }
393
394 /**
395 * Internal method to compute the discriminant of a 3rd degree polynomial.
396 * Discriminants are helpful to determine properties of a 3rd degree
397 * polynomial
398 *
399 * @param polyParams Array containing polynomial parameters
400 * @return Value of discriminant
401 */
402 private static double getDiscriminant(final double[] polyParams) {
403 final var d = polyParams[0];
404 final var c = polyParams[1];
405 final var b = polyParams[2];
406 final var a = polyParams[3];
407
408 return 18.0 * a * b * c * d
409 - 4.0 * b * b * b * d
410 + b * b * c * c
411 - 4.0 * a * c * c * c
412 - 27.0 * a * a * d * d;
413 }
414
415 /**
416 * Computes the cube root or x^(1/3) of provided value x
417 *
418 * @param x Provided value
419 * @return Cube root
420 */
421 private double cubeRoot(final double x) {
422 if (x < 0.0) {
423 return -Math.pow(-x, THIRD);
424 } else {
425 return Math.pow(x, THIRD);
426 }
427 }
428
429 /**
430 * Finds 3rd degree polynomial roots
431 *
432 * @param a 1st parameter
433 * @param b 2nd parameter
434 * @param c 3rd parameter
435 * @param d 4th parameter
436 * @param x1 1st root (output parameter)
437 * @param x2 2nd root (output parameter)
438 * @param x3 3rd root (output parameter)
439 */
440 private void solveCubic(final double a, final double b, final double c, final double d,
441 final Complex x1, final Complex x2, final Complex x3) {
442
443 // find the discriminant
444 final var f = (3.0 * c / a - Math.pow(b, 2.0) / Math.pow(a, 2.0)) / 3.0;
445 final var g = (2.0 * Math.pow(b, 3.0) / Math.pow(a, 3.0) - 9.0 * b * c
446 / Math.pow(a, 2.0) + 27.0 * d / a) / 27.0;
447 final var h = Math.pow(g, 2.0) / 4.0 + Math.pow(f, 3.0) / 27.0;
448 final var absF = Math.abs(f);
449 final var absG = Math.abs(g);
450 final var absH = Math.abs(h);
451 // evaluate discriminant
452 if (absF <= EPS && absG <= EPS && absH <= EPS) {
453 // 3 equal roots
454
455 // when f, g, and h all equal 0 the roots can be found by the following line
456 final double x = -cubeRoot(d / a);
457 x1.setRealAndImaginary(x, 0.0);
458 x2.setRealAndImaginary(x, 0.0);
459 x3.setRealAndImaginary(x, 0.0);
460 } else if (h <= 0.0) {
461 // 3 real roots
462
463 // complicated maths making use of the method
464 final var i = Math.pow(Math.pow(g, 2.0) / 4 - h, 0.5);
465 final var j = cubeRoot(i);
466 final var k = Math.acos(-(g / (2.0 * i)));
467 final var m = Math.cos(k / 3.0);
468 final var n = ROOT_THREE * Math.sin(k / 3.0);
469 final var p = -(b / (3.0 * a));
470
471 // print solutions
472 x1.setRealAndImaginary(2.0 * j * m + p, 0.0);
473 x2.setRealAndImaginary(-j * (m + n) + p, 0.0);
474 x3.setRealAndImaginary(-j * (m - n) + p, 0.0);
475 } else if (h > 0) {
476 // 1 real root and 2 complex roots
477
478 // complicated maths making use of the method
479 final var r = -(g / 2) + Math.pow(h, 0.5);
480 final var s = cubeRoot(r);
481 final var t = -(g / 2) - Math.pow(h, 0.5);
482 final var u = cubeRoot(t);
483 final var p = -(b / (3 * a));
484
485 // print solutions
486 x1.setRealAndImaginary((s + u) + p, 0.0);
487 final var real = -(s + u) / 2 + p;
488 x2.setReal(real);
489 x3.setReal(real);
490 final var imag = (s - u) * ROOT_THREE / 2;
491 x2.setImaginary(imag);
492 x3.setImaginary(imag);
493 }
494 }
495 }