001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math3.analysis.solvers;
018
019import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
020import org.apache.commons.math3.complex.Complex;
021import org.apache.commons.math3.complex.ComplexUtils;
022import org.apache.commons.math3.exception.NoBracketingException;
023import org.apache.commons.math3.exception.NoDataException;
024import org.apache.commons.math3.exception.NullArgumentException;
025import org.apache.commons.math3.exception.NumberIsTooLargeException;
026import org.apache.commons.math3.exception.TooManyEvaluationsException;
027import org.apache.commons.math3.exception.util.LocalizedFormats;
028import org.apache.commons.math3.util.FastMath;
029
030/**
031 * Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html">
032 * Laguerre's Method</a> for root finding of real coefficient polynomials.
033 * For reference, see
034 * <blockquote>
035 *  <b>A First Course in Numerical Analysis</b>,
036 *  ISBN 048641454X, chapter 8.
037 * </blockquote>
038 * Laguerre's method is global in the sense that it can start with any initial
039 * approximation and be able to solve all roots from that point.
040 * The algorithm requires a bracketing condition.
041 *
042 * @since 1.2
043 */
044public class LaguerreSolver extends AbstractPolynomialSolver {
045    /** Default absolute accuracy. */
046    private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
047    /** Complex solver. */
048    private final ComplexSolver complexSolver = new ComplexSolver();
049
050    /**
051     * Construct a solver with default accuracy (1e-6).
052     */
053    public LaguerreSolver() {
054        this(DEFAULT_ABSOLUTE_ACCURACY);
055    }
056    /**
057     * Construct a solver.
058     *
059     * @param absoluteAccuracy Absolute accuracy.
060     */
061    public LaguerreSolver(double absoluteAccuracy) {
062        super(absoluteAccuracy);
063    }
064    /**
065     * Construct a solver.
066     *
067     * @param relativeAccuracy Relative accuracy.
068     * @param absoluteAccuracy Absolute accuracy.
069     */
070    public LaguerreSolver(double relativeAccuracy,
071                          double absoluteAccuracy) {
072        super(relativeAccuracy, absoluteAccuracy);
073    }
074    /**
075     * Construct a solver.
076     *
077     * @param relativeAccuracy Relative accuracy.
078     * @param absoluteAccuracy Absolute accuracy.
079     * @param functionValueAccuracy Function value accuracy.
080     */
081    public LaguerreSolver(double relativeAccuracy,
082                          double absoluteAccuracy,
083                          double functionValueAccuracy) {
084        super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
085    }
086
087    /**
088     * {@inheritDoc}
089     */
090    @Override
091    public double doSolve()
092        throws TooManyEvaluationsException,
093               NumberIsTooLargeException,
094               NoBracketingException {
095        final double min = getMin();
096        final double max = getMax();
097        final double initial = getStartValue();
098        final double functionValueAccuracy = getFunctionValueAccuracy();
099
100        verifySequence(min, initial, max);
101
102        // Return the initial guess if it is good enough.
103        final double yInitial = computeObjectiveValue(initial);
104        if (FastMath.abs(yInitial) <= functionValueAccuracy) {
105            return initial;
106        }
107
108        // Return the first endpoint if it is good enough.
109        final double yMin = computeObjectiveValue(min);
110        if (FastMath.abs(yMin) <= functionValueAccuracy) {
111            return min;
112        }
113
114        // Reduce interval if min and initial bracket the root.
115        if (yInitial * yMin < 0) {
116            return laguerre(min, initial, yMin, yInitial);
117        }
118
119        // Return the second endpoint if it is good enough.
120        final double yMax = computeObjectiveValue(max);
121        if (FastMath.abs(yMax) <= functionValueAccuracy) {
122            return max;
123        }
124
125        // Reduce interval if initial and max bracket the root.
126        if (yInitial * yMax < 0) {
127            return laguerre(initial, max, yInitial, yMax);
128        }
129
130        throw new NoBracketingException(min, max, yMin, yMax);
131    }
132
133    /**
134     * Find a real root in the given interval.
135     *
136     * Despite the bracketing condition, the root returned by
137     * {@link LaguerreSolver.ComplexSolver#solve(Complex[],Complex)} may
138     * not be a real zero inside {@code [min, max]}.
139     * For example, <code> p(x) = x<sup>3</sup> + 1, </code>
140     * with {@code min = -2}, {@code max = 2}, {@code initial = 0}.
141     * When it occurs, this code calls
142     * {@link LaguerreSolver.ComplexSolver#solveAll(Complex[],Complex)}
143     * in order to obtain all roots and picks up one real root.
144     *
145     * @param lo Lower bound of the search interval.
146     * @param hi Higher bound of the search interval.
147     * @param fLo Function value at the lower bound of the search interval.
148     * @param fHi Function value at the higher bound of the search interval.
149     * @return the point at which the function value is zero.
150     * @deprecated This method should not be part of the public API: It will
151     * be made private in version 4.0.
152     */
153    @Deprecated
154    public double laguerre(double lo, double hi,
155                           double fLo, double fHi) {
156        final Complex c[] = ComplexUtils.convertToComplex(getCoefficients());
157
158        final Complex initial = new Complex(0.5 * (lo + hi), 0);
159        final Complex z = complexSolver.solve(c, initial);
160        if (complexSolver.isRoot(lo, hi, z)) {
161            return z.getReal();
162        } else {
163            double r = Double.NaN;
164            // Solve all roots and select the one we are seeking.
165            Complex[] root = complexSolver.solveAll(c, initial);
166            for (int i = 0; i < root.length; i++) {
167                if (complexSolver.isRoot(lo, hi, root[i])) {
168                    r = root[i].getReal();
169                    break;
170                }
171            }
172            return r;
173        }
174    }
175
176    /**
177     * Find all complex roots for the polynomial with the given
178     * coefficients, starting from the given initial value.
179     * <p>
180     * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p>
181     *
182     * @param coefficients Polynomial coefficients.
183     * @param initial Start value.
184     * @return the full set of complex roots of the polynomial
185     * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
186     * if the maximum number of evaluations is exceeded when solving for one of the roots
187     * @throws NullArgumentException if the {@code coefficients} is
188     * {@code null}.
189     * @throws NoDataException if the {@code coefficients} array is empty.
190     * @since 3.1
191     */
192    public Complex[] solveAllComplex(double[] coefficients,
193                                     double initial)
194        throws NullArgumentException,
195               NoDataException,
196               TooManyEvaluationsException {
197       return solveAllComplex(coefficients, initial, Integer.MAX_VALUE);
198    }
199
200    /**
201     * Find all complex roots for the polynomial with the given
202     * coefficients, starting from the given initial value.
203     * <p>
204     * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p>
205     *
206     * @param coefficients polynomial coefficients
207     * @param initial start value
208     * @param maxEval maximum number of evaluations
209     * @return the full set of complex roots of the polynomial
210     * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
211     * if the maximum number of evaluations is exceeded when solving for one of the roots
212     * @throws NullArgumentException if the {@code coefficients} is
213     * {@code null}
214     * @throws NoDataException if the {@code coefficients} array is empty
215     * @since 3.5
216     */
217    public Complex[] solveAllComplex(double[] coefficients,
218                                     double initial, int maxEval)
219        throws NullArgumentException,
220               NoDataException,
221               TooManyEvaluationsException {
222        setup(maxEval,
223              new PolynomialFunction(coefficients),
224              Double.NEGATIVE_INFINITY,
225              Double.POSITIVE_INFINITY,
226              initial);
227        return complexSolver.solveAll(ComplexUtils.convertToComplex(coefficients),
228                                      new Complex(initial, 0d));
229    }
230
231    /**
232     * Find a complex root for the polynomial with the given coefficients,
233     * starting from the given initial value.
234     * <p>
235     * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p>
236     *
237     * @param coefficients Polynomial coefficients.
238     * @param initial Start value.
239     * @return a complex root of the polynomial
240     * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
241     * if the maximum number of evaluations is exceeded.
242     * @throws NullArgumentException if the {@code coefficients} is
243     * {@code null}.
244     * @throws NoDataException if the {@code coefficients} array is empty.
245     * @since 3.1
246     */
247    public Complex solveComplex(double[] coefficients,
248                                double initial)
249        throws NullArgumentException,
250               NoDataException,
251               TooManyEvaluationsException {
252       return solveComplex(coefficients, initial, Integer.MAX_VALUE);
253    }
254
255    /**
256     * Find a complex root for the polynomial with the given coefficients,
257     * starting from the given initial value.
258     * <p>
259     * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p>
260     *
261     * @param coefficients polynomial coefficients
262     * @param initial start value
263     * @param maxEval maximum number of evaluations
264     * @return a complex root of the polynomial
265     * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
266     * if the maximum number of evaluations is exceeded
267     * @throws NullArgumentException if the {@code coefficients} is
268     * {@code null}
269     * @throws NoDataException if the {@code coefficients} array is empty
270     * @since 3.1
271     */
272    public Complex solveComplex(double[] coefficients,
273                                double initial, int maxEval)
274        throws NullArgumentException,
275               NoDataException,
276               TooManyEvaluationsException {
277        setup(maxEval,
278              new PolynomialFunction(coefficients),
279              Double.NEGATIVE_INFINITY,
280              Double.POSITIVE_INFINITY,
281              initial);
282        return complexSolver.solve(ComplexUtils.convertToComplex(coefficients),
283                                   new Complex(initial, 0d));
284    }
285
286    /**
287     * Class for searching all (complex) roots.
288     */
289    private class ComplexSolver {
290        /**
291         * Check whether the given complex root is actually a real zero
292         * in the given interval, within the solver tolerance level.
293         *
294         * @param min Lower bound for the interval.
295         * @param max Upper bound for the interval.
296         * @param z Complex root.
297         * @return {@code true} if z is a real zero.
298         */
299        public boolean isRoot(double min, double max, Complex z) {
300            if (isSequence(min, z.getReal(), max)) {
301                double tolerance = FastMath.max(getRelativeAccuracy() * z.abs(), getAbsoluteAccuracy());
302                return (FastMath.abs(z.getImaginary()) <= tolerance) ||
303                     (z.abs() <= getFunctionValueAccuracy());
304            }
305            return false;
306        }
307
308        /**
309         * Find all complex roots for the polynomial with the given
310         * coefficients, starting from the given initial value.
311         *
312         * @param coefficients Polynomial coefficients.
313         * @param initial Start value.
314         * @return the point at which the function value is zero.
315         * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
316         * if the maximum number of evaluations is exceeded.
317         * @throws NullArgumentException if the {@code coefficients} is
318         * {@code null}.
319         * @throws NoDataException if the {@code coefficients} array is empty.
320         */
321        public Complex[] solveAll(Complex coefficients[], Complex initial)
322            throws NullArgumentException,
323                   NoDataException,
324                   TooManyEvaluationsException {
325            if (coefficients == null) {
326                throw new NullArgumentException();
327            }
328            final int n = coefficients.length - 1;
329            if (n == 0) {
330                throw new NoDataException(LocalizedFormats.POLYNOMIAL);
331            }
332            // Coefficients for deflated polynomial.
333            final Complex c[] = new Complex[n + 1];
334            for (int i = 0; i <= n; i++) {
335                c[i] = coefficients[i];
336            }
337
338            // Solve individual roots successively.
339            final Complex root[] = new Complex[n];
340            for (int i = 0; i < n; i++) {
341                final Complex subarray[] = new Complex[n - i + 1];
342                System.arraycopy(c, 0, subarray, 0, subarray.length);
343                root[i] = solve(subarray, initial);
344                // Polynomial deflation using synthetic division.
345                Complex newc = c[n - i];
346                Complex oldc = null;
347                for (int j = n - i - 1; j >= 0; j--) {
348                    oldc = c[j];
349                    c[j] = newc;
350                    newc = oldc.add(newc.multiply(root[i]));
351                }
352            }
353
354            return root;
355        }
356
357        /**
358         * Find a complex root for the polynomial with the given coefficients,
359         * starting from the given initial value.
360         *
361         * @param coefficients Polynomial coefficients.
362         * @param initial Start value.
363         * @return the point at which the function value is zero.
364         * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
365         * if the maximum number of evaluations is exceeded.
366         * @throws NullArgumentException if the {@code coefficients} is
367         * {@code null}.
368         * @throws NoDataException if the {@code coefficients} array is empty.
369         */
370        public Complex solve(Complex coefficients[], Complex initial)
371            throws NullArgumentException,
372                   NoDataException,
373                   TooManyEvaluationsException {
374            if (coefficients == null) {
375                throw new NullArgumentException();
376            }
377
378            final int n = coefficients.length - 1;
379            if (n == 0) {
380                throw new NoDataException(LocalizedFormats.POLYNOMIAL);
381            }
382
383            final double absoluteAccuracy = getAbsoluteAccuracy();
384            final double relativeAccuracy = getRelativeAccuracy();
385            final double functionValueAccuracy = getFunctionValueAccuracy();
386
387            final Complex nC  = new Complex(n, 0);
388            final Complex n1C = new Complex(n - 1, 0);
389
390            Complex z = initial;
391            Complex oldz = new Complex(Double.POSITIVE_INFINITY,
392                                       Double.POSITIVE_INFINITY);
393            while (true) {
394                // Compute pv (polynomial value), dv (derivative value), and
395                // d2v (second derivative value) simultaneously.
396                Complex pv = coefficients[n];
397                Complex dv = Complex.ZERO;
398                Complex d2v = Complex.ZERO;
399                for (int j = n-1; j >= 0; j--) {
400                    d2v = dv.add(z.multiply(d2v));
401                    dv = pv.add(z.multiply(dv));
402                    pv = coefficients[j].add(z.multiply(pv));
403                }
404                d2v = d2v.multiply(new Complex(2.0, 0.0));
405
406                // Check for convergence.
407                final double tolerance = FastMath.max(relativeAccuracy * z.abs(),
408                                                      absoluteAccuracy);
409                if ((z.subtract(oldz)).abs() <= tolerance) {
410                    return z;
411                }
412                if (pv.abs() <= functionValueAccuracy) {
413                    return z;
414                }
415
416                // Now pv != 0, calculate the new approximation.
417                final Complex G = dv.divide(pv);
418                final Complex G2 = G.multiply(G);
419                final Complex H = G2.subtract(d2v.divide(pv));
420                final Complex delta = n1C.multiply((nC.multiply(H)).subtract(G2));
421                // Choose a denominator larger in magnitude.
422                final Complex deltaSqrt = delta.sqrt();
423                final Complex dplus = G.add(deltaSqrt);
424                final Complex dminus = G.subtract(deltaSqrt);
425                final Complex denominator = dplus.abs() > dminus.abs() ? dplus : dminus;
426                // Perturb z if denominator is zero, for instance,
427                // p(x) = x^3 + 1, z = 0.
428                if (denominator.equals(new Complex(0.0, 0.0))) {
429                    z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
430                    oldz = new Complex(Double.POSITIVE_INFINITY,
431                                       Double.POSITIVE_INFINITY);
432                } else {
433                    oldz = z;
434                    z = z.subtract(nC.divide(denominator));
435                }
436                incrementEvaluationCount();
437            }
438        }
439    }
440}