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.NotReadyException;
21
22 import java.util.Arrays;
23
24 /**
25 * This class estimates the roots of a polynomial of degree n.
26 * p(x) = a0 * x^n + a1 * x^(n - 1) + ... a(n-1) * x + an
27 * then the array of parameters is [an, a(n-1), ... a1, a0]
28 * This class supports polynomials having either real or complex parameters.
29 */
30 public class LaguerrePolynomialRootsEstimator extends PolynomialRootsEstimator {
31
32 // In this implementation we have increased MR and MT to increase accuracy
33 // by iterating a larger but finite number of times
34
35 /**
36 * Constant that affects the number of iterations.
37 */
38 public static final int MR = 80;
39
40 /**
41 * Constant that affects the number of iterations.
42 */
43 public static final int MT = 100;
44
45 /**
46 * Maximum number of iterations.
47 */
48 public static final int MAXIT = MT * MR;
49
50 /**
51 * Constant considered as machine precision for Laguerre method.
52 */
53 public static final double LAGUER_EPS = 1e-10;
54
55 /**
56 * Constant considered as machine precision.
57 */
58 public static final double EPS = 1e-14;
59
60 /**
61 * Constant indicating whether roots will be refined.
62 */
63 public static final boolean DEFAULT_POLISH_ROOTS = true;
64
65 /**
66 * Minimum allowed length in polynomial parameters.
67 */
68 public static final int MIN_VALID_POLY_PARAMS_LENGTH = 2;
69
70 /**
71 * Array containing values for Laguerre method.
72 */
73 private static final double[] frac =
74 {0.0, 0.5, 0.25, 0.75, 0.13, 0.38, 0.62, 0.88, 1.0};
75
76
77 /**
78 * Indicates if roots should be refined.
79 */
80 private boolean polishRoots;
81
82 /**
83 * Constructor.
84 *
85 * @param polishRoots Boolean to determine whether roots should be refined.
86 */
87 public LaguerrePolynomialRootsEstimator(final boolean polishRoots) {
88 super();
89 this.polishRoots = polishRoots;
90 }
91
92 /**
93 * Empty constructor.
94 */
95 public LaguerrePolynomialRootsEstimator() {
96 super();
97 this.polishRoots = DEFAULT_POLISH_ROOTS;
98 }
99
100 /**
101 * Constructor.
102 *
103 * @param polyParams Array containing polynomial parameters.
104 * @param polishRoots Boolean indicating whether roots will be refined.
105 * @throws IllegalArgumentException Raised if length of provided parameters
106 * is not valid. It has to be greater or equal than 2.
107 */
108 public LaguerrePolynomialRootsEstimator(final Complex[] polyParams, final boolean polishRoots) {
109 super();
110 this.polishRoots = polishRoots;
111 internalSetPolynomialParameters(polyParams);
112 }
113
114 /**
115 * Constructor.
116 *
117 * @param polyParams Array containing polynomial parameters.
118 * @throws IllegalArgumentException Raised if length of provided parameters
119 * is not valid. It has to be greater or equal than 2.
120 */
121 public LaguerrePolynomialRootsEstimator(final Complex[] polyParams) {
122 super();
123 this.polishRoots = DEFAULT_POLISH_ROOTS;
124 internalSetPolynomialParameters(polyParams);
125 }
126
127 /**
128 * Estimates the roots of provided polynomial.
129 *
130 * @throws LockedException Raised if this instance is locked estimating a
131 * root.
132 * @throws NotReadyException Raised if this instance is not ready because
133 * polynomial parameters have not been provided.
134 * @throws RootEstimationException Raised if roots cannot be estimated for
135 * some reason (lack of convergence, etc.).
136 */
137 @Override
138 public void estimate() throws LockedException, NotReadyException, RootEstimationException {
139
140 if (isLocked()) {
141 throw new LockedException();
142 }
143 if (!isReady()) {
144 throw new NotReadyException();
145 }
146
147 // polynomial must be at least degree 1
148 if (polyParams.length < MIN_VALID_POLY_PARAMS_LENGTH) {
149 throw new RootEstimationException();
150 }
151
152 locked = true;
153
154 final var a = polyParams;
155 roots = new Complex[a.length - 1];
156
157 int i;
158 final var its = new int[1];
159 var x = new Complex();
160 Complex b;
161 Complex c;
162 final var m = a.length - 1;
163
164 final var ad = Arrays.copyOf(a, a.length);
165 for (var j = m - 1; j >= 0; j--) {
166 x.setRealAndImaginary(0.0, 0.0);
167 final var adV = Arrays.copyOf(ad, j + 2);
168 internalLaguer(adV, x, its);
169 if (Math.abs(x.getImaginary()) <= 2.0 * EPS * Math.abs(x.getReal())) {
170 x = new Complex(x.getReal(), 0.0);
171 }
172 roots[j] = new Complex(x);
173 b = new Complex(ad[j + 1]);
174 for (var jj = j; jj >= 0; jj--) {
175 c = new Complex(ad[jj]);
176 ad[jj] = new Complex(b);
177 b.multiply(x);
178 b.add(c);
179 }
180 }
181 if (polishRoots) {
182 for (var j = 0; j < m; j++) {
183 internalLaguer(a, roots[j], its);
184 }
185 }
186 for (var j = 1; j < m; j++) {
187 x = new Complex(roots[j]);
188 for (i = j - 1; i >= 0; i--) {
189 if (roots[i].getReal() <= x.getReal()) {
190 break;
191 }
192 roots[i + 1] = new Complex(roots[i]);
193 }
194 roots[i + 1] = new Complex(x);
195 }
196
197 locked = false;
198 }
199
200 /**
201 * Returns boolean indicating whether roots are refined after an initial
202 * estimation.
203 *
204 * @return True if roots are refined, false otherwise.
205 */
206 public boolean areRootsPolished() {
207 return polishRoots;
208 }
209
210 /**
211 * Sets boolean indicating whether roots will be refined after an initial
212 * estimation.
213 *
214 * @param enable True if roots will be refined, false otherwise.
215 * @throws LockedException Raised if this instance is locked.
216 */
217 public void setPolishRootsEnabled(final boolean enable) throws LockedException {
218 if (isLocked()) {
219 throw new LockedException();
220 }
221 polishRoots = enable;
222 }
223
224 /**
225 * Internal method to set parameters of a polynomial, taking into account
226 * that a polynomial of degree n is defined as:
227 * p(x) = a0 * x^n + a1 * x^(n - 1) + ... a(n-1) * x + an
228 * then the array of parameters is [an, a(n - 1), ... a1, a0]
229 * Polynomial parameters can be either real or complex values
230 * This method does not check if this class is locked.
231 *
232 * @param polyParams Polynomial parameters.
233 * @throws IllegalArgumentException Raised if the length of the array is not
234 * valid.
235 */
236 @Override
237 protected final void internalSetPolynomialParameters(final Complex[] polyParams) {
238 if (polyParams.length < MIN_VALID_POLY_PARAMS_LENGTH) {
239 throw new IllegalArgumentException();
240 }
241 this.polyParams = polyParams;
242 }
243
244 /**
245 * Internal method to compute a root after decomposing and decreasing the
246 * degree of the polynomial.
247 *
248 * @param a Remaining polynomial parameters (on 1st iteration, the whole
249 * polynomial is provided, on subsequent iterations, the polynomial is
250 * deflated and the degree is reduced).
251 * @param x Estimated root.
252 * @param its number of iterations needed to achieve the estimation.
253 * @throws RootEstimationException Raised if root couldn't be estimated
254 * because of lack of convergence.
255 */
256 private void internalLaguer(final Complex[] a, final Complex x, final int[] its) throws RootEstimationException {
257
258 Complex x1;
259 Complex b;
260 Complex g;
261 Complex g2;
262 final var dx = new Complex();
263 final var d = new Complex();
264 final var f = new Complex();
265 final var h = new Complex();
266 final var sq = new Complex();
267 var gp = new Complex();
268 final var gm = new Complex();
269 final var m = a.length - 1;
270 for (var iter = 1; iter <= MAXIT; iter++) {
271 its[0] = iter;
272 b = new Complex(a[m]);
273 var err = b.getModulus();
274 d.setRealAndImaginary(0.0, 0.0);
275 f.setRealAndImaginary(0.0, 0.0);
276 final var abx = x.getModulus();
277 for (var j = m - 1; j >= 0; j--) {
278 f.multiply(x);
279 f.add(d);
280
281 d.multiply(x);
282 d.add(b);
283
284 b.multiply(x);
285 b.add(a[j]);
286
287 err = b.getModulus() + abx * err;
288 }
289 err *= LAGUER_EPS;
290 if (b.getModulus() <= err) {
291 return;
292 }
293 g = d.divideAndReturnNew(b);
294 g2 = g.powAndReturnNew(2.0);
295 f.divide(b, h);
296 h.multiplyByScalar(-2.0);
297 h.add(g2);
298
299 h.multiplyByScalar(m, sq);
300 sq.subtract(g2);
301 sq.multiplyByScalar(m - 1.0);
302 sq.sqrt();
303
304 g.add(sq, gp);
305 g.subtract(sq, gm);
306
307 final var abp = gp.getModulus();
308 final var abm = gm.getModulus();
309 if (abp < abm) {
310 gp = gm;
311 }
312 if (Math.max(abp, abm) > 0.0) {
313 dx.setRealAndImaginary(m, 0.0);
314 dx.divide(gp);
315 } else {
316 dx.setModulusAndPhase(1.0 + abx, iter);
317 }
318 x1 = x.subtractAndReturnNew(dx);
319 if (x.equals(x1)) {
320 return;
321 }
322 if (iter % MT != 0) {
323 x.copyFrom(x1);
324 } else {
325 int pos = Math.min(iter / MT, frac.length - 1);
326 x.subtract(dx.multiplyByScalarAndReturnNew(frac[pos]));
327 }
328 }
329 // too many iterations in Laguerre
330 locked = false;
331 throw new RootEstimationException();
332 }
333 }