1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 package org.hipparchus.analysis.solvers;
23
24 import org.hipparchus.analysis.polynomials.PolynomialFunction;
25 import org.hipparchus.complex.Complex;
26 import org.hipparchus.complex.ComplexUtils;
27 import org.hipparchus.exception.LocalizedCoreFormats;
28 import org.hipparchus.exception.MathIllegalArgumentException;
29 import org.hipparchus.exception.MathIllegalStateException;
30 import org.hipparchus.exception.NullArgumentException;
31 import org.hipparchus.util.FastMath;
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46 public class LaguerreSolver extends AbstractPolynomialSolver {
47
48 private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
49
50 private final ComplexSolver complexSolver = new ComplexSolver();
51
52
53
54
55 public LaguerreSolver() {
56 this(DEFAULT_ABSOLUTE_ACCURACY);
57 }
58
59
60
61
62
63 public LaguerreSolver(double absoluteAccuracy) {
64 super(absoluteAccuracy);
65 }
66
67
68
69
70
71
72 public LaguerreSolver(double relativeAccuracy,
73 double absoluteAccuracy) {
74 super(relativeAccuracy, absoluteAccuracy);
75 }
76
77
78
79
80
81
82
83 public LaguerreSolver(double relativeAccuracy,
84 double absoluteAccuracy,
85 double functionValueAccuracy) {
86 super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
87 }
88
89
90
91
92 @Override
93 public double doSolve()
94 throws MathIllegalArgumentException, MathIllegalStateException {
95 final double min = getMin();
96 final double max = getMax();
97 final double initial = getStartValue();
98 final double functionValueAccuracy = getFunctionValueAccuracy();
99
100 verifySequence(min, initial, max);
101
102
103 final double yInitial = computeObjectiveValue(initial);
104 if (FastMath.abs(yInitial) <= functionValueAccuracy) {
105 return initial;
106 }
107
108
109 final double yMin = computeObjectiveValue(min);
110 if (FastMath.abs(yMin) <= functionValueAccuracy) {
111 return min;
112 }
113
114
115 if (yInitial * yMin < 0) {
116 return laguerre(min, initial);
117 }
118
119
120 final double yMax = computeObjectiveValue(max);
121 if (FastMath.abs(yMax) <= functionValueAccuracy) {
122 return max;
123 }
124
125
126 if (yInitial * yMax < 0) {
127 return laguerre(initial, max);
128 }
129
130 throw new MathIllegalArgumentException(LocalizedCoreFormats.NOT_BRACKETING_INTERVAL,
131 min, max, yMin, yMax);
132 }
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150 private double laguerre(double lo, double hi) {
151 final Complex c[] = ComplexUtils.convertToComplex(getCoefficients());
152
153 final Complex initial = new Complex(0.5 * (lo + hi), 0);
154 final Complex z = complexSolver.solve(c, initial);
155 if (complexSolver.isRoot(lo, hi, z)) {
156 return z.getReal();
157 } else {
158 double r = Double.NaN;
159
160 Complex[] root = complexSolver.solveAll(c, initial);
161 for (int i = 0; i < root.length; i++) {
162 if (complexSolver.isRoot(lo, hi, root[i])) {
163 r = root[i].getReal();
164 break;
165 }
166 }
167 return r;
168 }
169 }
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186 public Complex[] solveAllComplex(double[] coefficients,
187 double initial)
188 throws MathIllegalArgumentException, NullArgumentException,
189 MathIllegalStateException {
190 setup(Integer.MAX_VALUE,
191 new PolynomialFunction(coefficients),
192 Double.NEGATIVE_INFINITY,
193 Double.POSITIVE_INFINITY,
194 initial);
195 return complexSolver.solveAll(ComplexUtils.convertToComplex(coefficients),
196 new Complex(initial, 0d));
197 }
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214 public Complex solveComplex(double[] coefficients,
215 double initial)
216 throws MathIllegalArgumentException, NullArgumentException,
217 MathIllegalStateException {
218 setup(Integer.MAX_VALUE,
219 new PolynomialFunction(coefficients),
220 Double.NEGATIVE_INFINITY,
221 Double.POSITIVE_INFINITY,
222 initial);
223 return complexSolver.solve(ComplexUtils.convertToComplex(coefficients),
224 new Complex(initial, 0d));
225 }
226
227
228
229
230 private class ComplexSolver {
231
232
233
234
235
236
237
238
239
240 public boolean isRoot(double min, double max, Complex z) {
241 if (isSequence(min, z.getReal(), max)) {
242 final double zAbs = z.norm();
243 double tolerance = FastMath.max(getRelativeAccuracy() * zAbs, getAbsoluteAccuracy());
244 return (FastMath.abs(z.getImaginary()) <= tolerance) ||
245 (zAbs <= getFunctionValueAccuracy());
246 }
247 return false;
248 }
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263 public Complex[] solveAll(Complex coefficients[], Complex initial)
264 throws MathIllegalArgumentException, NullArgumentException,
265 MathIllegalStateException {
266 if (coefficients == null) {
267 throw new NullArgumentException();
268 }
269 final int n = coefficients.length - 1;
270 if (n == 0) {
271 throw new MathIllegalArgumentException(LocalizedCoreFormats.POLYNOMIAL);
272 }
273
274 final Complex c[] = new Complex[n + 1];
275 for (int i = 0; i <= n; i++) {
276 c[i] = coefficients[i];
277 }
278
279
280 final Complex root[] = new Complex[n];
281 for (int i = 0; i < n; i++) {
282 final Complex subarray[] = new Complex[n - i + 1];
283 System.arraycopy(c, 0, subarray, 0, subarray.length);
284 root[i] = solve(subarray, initial);
285
286 Complex newc = c[n - i];
287 Complex oldc;
288 for (int j = n - i - 1; j >= 0; j--) {
289 oldc = c[j];
290 c[j] = newc;
291 newc = oldc.add(newc.multiply(root[i]));
292 }
293 }
294
295 return root;
296 }
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311 public Complex solve(Complex coefficients[], Complex initial)
312 throws MathIllegalArgumentException, NullArgumentException,
313 MathIllegalStateException {
314 if (coefficients == null) {
315 throw new NullArgumentException();
316 }
317
318 final int n = coefficients.length - 1;
319 if (n == 0) {
320 throw new MathIllegalArgumentException(LocalizedCoreFormats.POLYNOMIAL);
321 }
322
323 final double absoluteAccuracy = getAbsoluteAccuracy();
324 final double relativeAccuracy = getRelativeAccuracy();
325 final double functionValueAccuracy = getFunctionValueAccuracy();
326
327 final Complex nC = new Complex(n, 0);
328 final Complex n1C = new Complex(n - 1, 0);
329
330 Complex z = initial;
331 Complex oldz = new Complex(Double.POSITIVE_INFINITY,
332 Double.POSITIVE_INFINITY);
333 while (true) {
334
335
336 Complex pv = coefficients[n];
337 Complex dv = Complex.ZERO;
338 Complex d2v = Complex.ZERO;
339 for (int j = n-1; j >= 0; j--) {
340 d2v = dv.add(z.multiply(d2v));
341 dv = pv.add(z.multiply(dv));
342 pv = coefficients[j].add(z.multiply(pv));
343 }
344 d2v = d2v.multiply(new Complex(2.0, 0.0));
345
346
347 final double tolerance = FastMath.max(relativeAccuracy * z.norm(),
348 absoluteAccuracy);
349 if ((z.subtract(oldz)).norm() <= tolerance) {
350 return z;
351 }
352 if (pv.norm() <= functionValueAccuracy) {
353 return z;
354 }
355
356
357 final Complex G = dv.divide(pv);
358 final Complex G2 = G.multiply(G);
359 final Complex H = G2.subtract(d2v.divide(pv));
360 final Complex delta = n1C.multiply((nC.multiply(H)).subtract(G2));
361
362 final Complex deltaSqrt = delta.sqrt();
363 final Complex dplus = G.add(deltaSqrt);
364 final Complex dminus = G.subtract(deltaSqrt);
365 final Complex denominator = dplus.norm() > dminus.norm() ? dplus : dminus;
366
367
368 if (denominator.isZero()) {
369 z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
370 oldz = new Complex(Double.POSITIVE_INFINITY,
371 Double.POSITIVE_INFINITY);
372 } else {
373 oldz = z;
374 z = z.subtract(nC.divide(denominator));
375 }
376 incrementEvaluationCount();
377 }
378 }
379 }
380 }