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.complex.Complex;
020import org.apache.commons.math3.complex.ComplexUtils;
021import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
022import org.apache.commons.math3.exception.NoBracketingException;
023import org.apache.commons.math3.exception.NullArgumentException;
024import org.apache.commons.math3.exception.NoDataException;
025import org.apache.commons.math3.exception.TooManyEvaluationsException;
026import org.apache.commons.math3.exception.NumberIsTooLargeException;
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 * <quote>
035 *  <b>A First Course in Numerical Analysis</b>
036 *  ISBN 048641454X, chapter 8.
037 * </quote>
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     * <br/>
180     * Note: This method is not part of the API of {@link BaseUnivariateSolver}.
181     *
182     * @param coefficients Polynomial coefficients.
183     * @param initial Start value.
184     * @return the point at which the function value is zero.
185     * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
186     * if the maximum number of evaluations is exceeded.
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        setup(Integer.MAX_VALUE,
198              new PolynomialFunction(coefficients),
199              Double.NEGATIVE_INFINITY,
200              Double.POSITIVE_INFINITY,
201              initial);
202        return complexSolver.solveAll(ComplexUtils.convertToComplex(coefficients),
203                                      new Complex(initial, 0d));
204    }
205
206    /**
207     * Find a complex root for the polynomial with the given coefficients,
208     * starting from the given initial value.
209     * <br/>
210     * Note: This method is not part of the API of {@link BaseUnivariateSolver}.
211     *
212     * @param coefficients Polynomial coefficients.
213     * @param initial Start value.
214     * @return the point at which the function value is zero.
215     * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
216     * if the maximum number of evaluations is exceeded.
217     * @throws NullArgumentException if the {@code coefficients} is
218     * {@code null}.
219     * @throws NoDataException if the {@code coefficients} array is empty.
220     * @since 3.1
221     */
222    public Complex solveComplex(double[] coefficients,
223                                double initial)
224        throws NullArgumentException,
225               NoDataException,
226               TooManyEvaluationsException {
227        setup(Integer.MAX_VALUE,
228              new PolynomialFunction(coefficients),
229              Double.NEGATIVE_INFINITY,
230              Double.POSITIVE_INFINITY,
231              initial);
232        return complexSolver.solve(ComplexUtils.convertToComplex(coefficients),
233                                   new Complex(initial, 0d));
234    }
235
236    /**
237     * Class for searching all (complex) roots.
238     */
239    private class ComplexSolver {
240        /**
241         * Check whether the given complex root is actually a real zero
242         * in the given interval, within the solver tolerance level.
243         *
244         * @param min Lower bound for the interval.
245         * @param max Upper bound for the interval.
246         * @param z Complex root.
247         * @return {@code true} if z is a real zero.
248         */
249        public boolean isRoot(double min, double max, Complex z) {
250            if (isSequence(min, z.getReal(), max)) {
251                double tolerance = FastMath.max(getRelativeAccuracy() * z.abs(), getAbsoluteAccuracy());
252                return (FastMath.abs(z.getImaginary()) <= tolerance) ||
253                     (z.abs() <= getFunctionValueAccuracy());
254            }
255            return false;
256        }
257
258        /**
259         * Find all complex roots for the polynomial with the given
260         * coefficients, starting from the given initial value.
261         *
262         * @param coefficients Polynomial coefficients.
263         * @param initial Start value.
264         * @return the point at which the function value is zero.
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         */
271        public Complex[] solveAll(Complex coefficients[], Complex initial)
272            throws NullArgumentException,
273                   NoDataException,
274                   TooManyEvaluationsException {
275            if (coefficients == null) {
276                throw new NullArgumentException();
277            }
278            final int n = coefficients.length - 1;
279            if (n == 0) {
280                throw new NoDataException(LocalizedFormats.POLYNOMIAL);
281            }
282            // Coefficients for deflated polynomial.
283            final Complex c[] = new Complex[n + 1];
284            for (int i = 0; i <= n; i++) {
285                c[i] = coefficients[i];
286            }
287
288            // Solve individual roots successively.
289            final Complex root[] = new Complex[n];
290            for (int i = 0; i < n; i++) {
291                final Complex subarray[] = new Complex[n - i + 1];
292                System.arraycopy(c, 0, subarray, 0, subarray.length);
293                root[i] = solve(subarray, initial);
294                // Polynomial deflation using synthetic division.
295                Complex newc = c[n - i];
296                Complex oldc = null;
297                for (int j = n - i - 1; j >= 0; j--) {
298                    oldc = c[j];
299                    c[j] = newc;
300                    newc = oldc.add(newc.multiply(root[i]));
301                }
302            }
303
304            return root;
305        }
306
307        /**
308         * Find a complex root for the polynomial with the given coefficients,
309         * starting from the given initial value.
310         *
311         * @param coefficients Polynomial coefficients.
312         * @param initial Start value.
313         * @return the point at which the function value is zero.
314         * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
315         * if the maximum number of evaluations is exceeded.
316         * @throws NullArgumentException if the {@code coefficients} is
317         * {@code null}.
318         * @throws NoDataException if the {@code coefficients} array is empty.
319         */
320        public Complex solve(Complex coefficients[], Complex initial)
321            throws NullArgumentException,
322                   NoDataException,
323                   TooManyEvaluationsException {
324            if (coefficients == null) {
325                throw new NullArgumentException();
326            }
327
328            final int n = coefficients.length - 1;
329            if (n == 0) {
330                throw new NoDataException(LocalizedFormats.POLYNOMIAL);
331            }
332
333            final double absoluteAccuracy = getAbsoluteAccuracy();
334            final double relativeAccuracy = getRelativeAccuracy();
335            final double functionValueAccuracy = getFunctionValueAccuracy();
336
337            final Complex nC  = new Complex(n, 0);
338            final Complex n1C = new Complex(n - 1, 0);
339
340            Complex z = initial;
341            Complex oldz = new Complex(Double.POSITIVE_INFINITY,
342                                       Double.POSITIVE_INFINITY);
343            while (true) {
344                // Compute pv (polynomial value), dv (derivative value), and
345                // d2v (second derivative value) simultaneously.
346                Complex pv = coefficients[n];
347                Complex dv = Complex.ZERO;
348                Complex d2v = Complex.ZERO;
349                for (int j = n-1; j >= 0; j--) {
350                    d2v = dv.add(z.multiply(d2v));
351                    dv = pv.add(z.multiply(dv));
352                    pv = coefficients[j].add(z.multiply(pv));
353                }
354                d2v = d2v.multiply(new Complex(2.0, 0.0));
355
356                // Check for convergence.
357                final double tolerance = FastMath.max(relativeAccuracy * z.abs(),
358                                                      absoluteAccuracy);
359                if ((z.subtract(oldz)).abs() <= tolerance) {
360                    return z;
361                }
362                if (pv.abs() <= functionValueAccuracy) {
363                    return z;
364                }
365
366                // Now pv != 0, calculate the new approximation.
367                final Complex G = dv.divide(pv);
368                final Complex G2 = G.multiply(G);
369                final Complex H = G2.subtract(d2v.divide(pv));
370                final Complex delta = n1C.multiply((nC.multiply(H)).subtract(G2));
371                // Choose a denominator larger in magnitude.
372                final Complex deltaSqrt = delta.sqrt();
373                final Complex dplus = G.add(deltaSqrt);
374                final Complex dminus = G.subtract(deltaSqrt);
375                final Complex denominator = dplus.abs() > dminus.abs() ? dplus : dminus;
376                // Perturb z if denominator is zero, for instance,
377                // p(x) = x^3 + 1, z = 0.
378                if (denominator.equals(new Complex(0.0, 0.0))) {
379                    z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
380                    oldz = new Complex(Double.POSITIVE_INFINITY,
381                                       Double.POSITIVE_INFINITY);
382                } else {
383                    oldz = z;
384                    z = z.subtract(nC.divide(denominator));
385                }
386                incrementEvaluationCount();
387            }
388        }
389    }
390}