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 * @version $Id: LaguerreSolver.java 1422195 2012-12-15 06:45:18Z psteitz $
043 * @since 1.2
044 */
045public class LaguerreSolver extends AbstractPolynomialSolver {
046    /** Default absolute accuracy. */
047    private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
048    /** Complex solver. */
049    private final ComplexSolver complexSolver = new ComplexSolver();
050
051    /**
052     * Construct a solver with default accuracy (1e-6).
053     */
054    public LaguerreSolver() {
055        this(DEFAULT_ABSOLUTE_ACCURACY);
056    }
057    /**
058     * Construct a solver.
059     *
060     * @param absoluteAccuracy Absolute accuracy.
061     */
062    public LaguerreSolver(double absoluteAccuracy) {
063        super(absoluteAccuracy);
064    }
065    /**
066     * Construct a solver.
067     *
068     * @param relativeAccuracy Relative accuracy.
069     * @param absoluteAccuracy Absolute accuracy.
070     */
071    public LaguerreSolver(double relativeAccuracy,
072                          double absoluteAccuracy) {
073        super(relativeAccuracy, absoluteAccuracy);
074    }
075    /**
076     * Construct a solver.
077     *
078     * @param relativeAccuracy Relative accuracy.
079     * @param absoluteAccuracy Absolute accuracy.
080     * @param functionValueAccuracy Function value accuracy.
081     */
082    public LaguerreSolver(double relativeAccuracy,
083                          double absoluteAccuracy,
084                          double functionValueAccuracy) {
085        super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
086    }
087
088    /**
089     * {@inheritDoc}
090     */
091    @Override
092    public double doSolve()
093        throws TooManyEvaluationsException,
094               NumberIsTooLargeException,
095               NoBracketingException {
096        final double min = getMin();
097        final double max = getMax();
098        final double initial = getStartValue();
099        final double functionValueAccuracy = getFunctionValueAccuracy();
100
101        verifySequence(min, initial, max);
102
103        // Return the initial guess if it is good enough.
104        final double yInitial = computeObjectiveValue(initial);
105        if (FastMath.abs(yInitial) <= functionValueAccuracy) {
106            return initial;
107        }
108
109        // Return the first endpoint if it is good enough.
110        final double yMin = computeObjectiveValue(min);
111        if (FastMath.abs(yMin) <= functionValueAccuracy) {
112            return min;
113        }
114
115        // Reduce interval if min and initial bracket the root.
116        if (yInitial * yMin < 0) {
117            return laguerre(min, initial, yMin, yInitial);
118        }
119
120        // Return the second endpoint if it is good enough.
121        final double yMax = computeObjectiveValue(max);
122        if (FastMath.abs(yMax) <= functionValueAccuracy) {
123            return max;
124        }
125
126        // Reduce interval if initial and max bracket the root.
127        if (yInitial * yMax < 0) {
128            return laguerre(initial, max, yInitial, yMax);
129        }
130
131        throw new NoBracketingException(min, max, yMin, yMax);
132    }
133
134    /**
135     * Find a real root in the given interval.
136     *
137     * Despite the bracketing condition, the root returned by
138     * {@link LaguerreSolver.ComplexSolver#solve(Complex[],Complex)} may
139     * not be a real zero inside {@code [min, max]}.
140     * For example, <code>p(x) = x<sup>3</sup> + 1,</code>
141     * with {@code min = -2}, {@code max = 2}, {@code initial = 0}.
142     * When it occurs, this code calls
143     * {@link LaguerreSolver.ComplexSolver#solveAll(Complex[],Complex)}
144     * in order to obtain all roots and picks up one real root.
145     *
146     * @param lo Lower bound of the search interval.
147     * @param hi Higher bound of the search interval.
148     * @param fLo Function value at the lower bound of the search interval.
149     * @param fHi Function value at the higher bound of the search interval.
150     * @return the point at which the function value is zero.
151     * @deprecated This method should not be part of the public API: It will
152     * be made private in version 4.0.
153     */
154    @Deprecated
155    public double laguerre(double lo, double hi,
156                           double fLo, double fHi) {
157        final Complex c[] = ComplexUtils.convertToComplex(getCoefficients());
158
159        final Complex initial = new Complex(0.5 * (lo + hi), 0);
160        final Complex z = complexSolver.solve(c, initial);
161        if (complexSolver.isRoot(lo, hi, z)) {
162            return z.getReal();
163        } else {
164            double r = Double.NaN;
165            // Solve all roots and select the one we are seeking.
166            Complex[] root = complexSolver.solveAll(c, initial);
167            for (int i = 0; i < root.length; i++) {
168                if (complexSolver.isRoot(lo, hi, root[i])) {
169                    r = root[i].getReal();
170                    break;
171                }
172            }
173            return r;
174        }
175    }
176
177    /**
178     * Find all complex roots for the polynomial with the given
179     * coefficients, starting from the given initial value.
180     * <br/>
181     * Note: This method is not part of the API of {@link BaseUnivariateSolver}.
182     *
183     * @param coefficients Polynomial coefficients.
184     * @param initial Start value.
185     * @return the point at which the function value is zero.
186     * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
187     * if the maximum number of evaluations is exceeded.
188     * @throws NullArgumentException if the {@code coefficients} is
189     * {@code null}.
190     * @throws NoDataException if the {@code coefficients} array is empty.
191     * @since 3.1
192     */
193    public Complex[] solveAllComplex(double[] coefficients,
194                                     double initial)
195        throws NullArgumentException,
196               NoDataException,
197               TooManyEvaluationsException {
198        setup(Integer.MAX_VALUE,
199              new PolynomialFunction(coefficients),
200              Double.NEGATIVE_INFINITY,
201              Double.POSITIVE_INFINITY,
202              initial);
203        return complexSolver.solveAll(ComplexUtils.convertToComplex(coefficients),
204                                      new Complex(initial, 0d));
205    }
206
207    /**
208     * Find a complex root for the polynomial with the given coefficients,
209     * starting from the given initial value.
210     * <br/>
211     * Note: This method is not part of the API of {@link BaseUnivariateSolver}.
212     *
213     * @param coefficients Polynomial coefficients.
214     * @param initial Start value.
215     * @return the point at which the function value is zero.
216     * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
217     * if the maximum number of evaluations is exceeded.
218     * @throws NullArgumentException if the {@code coefficients} is
219     * {@code null}.
220     * @throws NoDataException if the {@code coefficients} array is empty.
221     * @since 3.1
222     */
223    public Complex solveComplex(double[] coefficients,
224                                double initial)
225        throws NullArgumentException,
226               NoDataException,
227               TooManyEvaluationsException {
228        setup(Integer.MAX_VALUE,
229              new PolynomialFunction(coefficients),
230              Double.NEGATIVE_INFINITY,
231              Double.POSITIVE_INFINITY,
232              initial);
233        return complexSolver.solve(ComplexUtils.convertToComplex(coefficients),
234                                   new Complex(initial, 0d));
235    }
236
237    /**
238     * Class for searching all (complex) roots.
239     */
240    private class ComplexSolver {
241        /**
242         * Check whether the given complex root is actually a real zero
243         * in the given interval, within the solver tolerance level.
244         *
245         * @param min Lower bound for the interval.
246         * @param max Upper bound for the interval.
247         * @param z Complex root.
248         * @return {@code true} if z is a real zero.
249         */
250        public boolean isRoot(double min, double max, Complex z) {
251            if (isSequence(min, z.getReal(), max)) {
252                double tolerance = FastMath.max(getRelativeAccuracy() * z.abs(), getAbsoluteAccuracy());
253                return (FastMath.abs(z.getImaginary()) <= tolerance) ||
254                     (z.abs() <= getFunctionValueAccuracy());
255            }
256            return false;
257        }
258
259        /**
260         * Find all complex roots for the polynomial with the given
261         * coefficients, starting from the given initial value.
262         *
263         * @param coefficients Polynomial coefficients.
264         * @param initial Start value.
265         * @return the point at which the function value is zero.
266         * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
267         * if the maximum number of evaluations is exceeded.
268         * @throws NullArgumentException if the {@code coefficients} is
269         * {@code null}.
270         * @throws NoDataException if the {@code coefficients} array is empty.
271         */
272        public Complex[] solveAll(Complex coefficients[], Complex initial)
273            throws NullArgumentException,
274                   NoDataException,
275                   TooManyEvaluationsException {
276            if (coefficients == null) {
277                throw new NullArgumentException();
278            }
279            final int n = coefficients.length - 1;
280            if (n == 0) {
281                throw new NoDataException(LocalizedFormats.POLYNOMIAL);
282            }
283            // Coefficients for deflated polynomial.
284            final Complex c[] = new Complex[n + 1];
285            for (int i = 0; i <= n; i++) {
286                c[i] = coefficients[i];
287            }
288
289            // Solve individual roots successively.
290            final Complex root[] = new Complex[n];
291            for (int i = 0; i < n; i++) {
292                final Complex subarray[] = new Complex[n - i + 1];
293                System.arraycopy(c, 0, subarray, 0, subarray.length);
294                root[i] = solve(subarray, initial);
295                // Polynomial deflation using synthetic division.
296                Complex newc = c[n - i];
297                Complex oldc = null;
298                for (int j = n - i - 1; j >= 0; j--) {
299                    oldc = c[j];
300                    c[j] = newc;
301                    newc = oldc.add(newc.multiply(root[i]));
302                }
303            }
304
305            return root;
306        }
307
308        /**
309         * Find a complex root for the polynomial with the given coefficients,
310         * 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 solve(Complex coefficients[], Complex initial)
322            throws NullArgumentException,
323                   NoDataException,
324                   TooManyEvaluationsException {
325            if (coefficients == null) {
326                throw new NullArgumentException();
327            }
328
329            final int n = coefficients.length - 1;
330            if (n == 0) {
331                throw new NoDataException(LocalizedFormats.POLYNOMIAL);
332            }
333
334            final double absoluteAccuracy = getAbsoluteAccuracy();
335            final double relativeAccuracy = getRelativeAccuracy();
336            final double functionValueAccuracy = getFunctionValueAccuracy();
337
338            final Complex nC  = new Complex(n, 0);
339            final Complex n1C = new Complex(n - 1, 0);
340
341            Complex z = initial;
342            Complex oldz = new Complex(Double.POSITIVE_INFINITY,
343                                       Double.POSITIVE_INFINITY);
344            while (true) {
345                // Compute pv (polynomial value), dv (derivative value), and
346                // d2v (second derivative value) simultaneously.
347                Complex pv = coefficients[n];
348                Complex dv = Complex.ZERO;
349                Complex d2v = Complex.ZERO;
350                for (int j = n-1; j >= 0; j--) {
351                    d2v = dv.add(z.multiply(d2v));
352                    dv = pv.add(z.multiply(dv));
353                    pv = coefficients[j].add(z.multiply(pv));
354                }
355                d2v = d2v.multiply(new Complex(2.0, 0.0));
356
357                // Check for convergence.
358                final double tolerance = FastMath.max(relativeAccuracy * z.abs(),
359                                                      absoluteAccuracy);
360                if ((z.subtract(oldz)).abs() <= tolerance) {
361                    return z;
362                }
363                if (pv.abs() <= functionValueAccuracy) {
364                    return z;
365                }
366
367                // Now pv != 0, calculate the new approximation.
368                final Complex G = dv.divide(pv);
369                final Complex G2 = G.multiply(G);
370                final Complex H = G2.subtract(d2v.divide(pv));
371                final Complex delta = n1C.multiply((nC.multiply(H)).subtract(G2));
372                // Choose a denominator larger in magnitude.
373                final Complex deltaSqrt = delta.sqrt();
374                final Complex dplus = G.add(deltaSqrt);
375                final Complex dminus = G.subtract(deltaSqrt);
376                final Complex denominator = dplus.abs() > dminus.abs() ? dplus : dminus;
377                // Perturb z if denominator is zero, for instance,
378                // p(x) = x^3 + 1, z = 0.
379                if (denominator.equals(new Complex(0.0, 0.0))) {
380                    z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
381                    oldz = new Complex(Double.POSITIVE_INFINITY,
382                                       Double.POSITIVE_INFINITY);
383                } else {
384                    oldz = z;
385                    z = z.subtract(nC.divide(denominator));
386                }
387                incrementEvaluationCount();
388            }
389        }
390    }
391}