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.math4.legacy.analysis.solvers;
018
019import org.apache.commons.numbers.complex.Complex;
020import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialFunction;
021import org.apache.commons.math4.legacy.exception.NoBracketingException;
022import org.apache.commons.math4.legacy.exception.NoDataException;
023import org.apache.commons.math4.legacy.exception.NullArgumentException;
024import org.apache.commons.math4.legacy.exception.NumberIsTooLargeException;
025import org.apache.commons.math4.legacy.exception.TooManyEvaluationsException;
026import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
027import org.apache.commons.math4.core.jdkmath.JdkMath;
028
029/**
030 * Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html">
031 * Laguerre's Method</a> for root finding of real coefficient polynomials.
032 * For reference, see
033 * <blockquote>
034 *  <b>A First Course in Numerical Analysis</b>,
035 *  ISBN 048641454X, chapter 8.
036 * </blockquote>
037 * Laguerre's method is global in the sense that it can start with any initial
038 * approximation and be able to solve all roots from that point.
039 * The algorithm requires a bracketing condition.
040 *
041 * @since 1.2
042 */
043public class LaguerreSolver extends AbstractPolynomialSolver {
044    /** Default absolute accuracy. */
045    private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
046    /** Complex solver. */
047    private final ComplexSolver complexSolver = new ComplexSolver();
048
049    /**
050     * Construct a solver with default accuracy (1e-6).
051     */
052    public LaguerreSolver() {
053        this(DEFAULT_ABSOLUTE_ACCURACY);
054    }
055    /**
056     * Construct a solver.
057     *
058     * @param absoluteAccuracy Absolute accuracy.
059     */
060    public LaguerreSolver(double absoluteAccuracy) {
061        super(absoluteAccuracy);
062    }
063    /**
064     * Construct a solver.
065     *
066     * @param relativeAccuracy Relative accuracy.
067     * @param absoluteAccuracy Absolute accuracy.
068     */
069    public LaguerreSolver(double relativeAccuracy,
070                          double absoluteAccuracy) {
071        super(relativeAccuracy, absoluteAccuracy);
072    }
073    /**
074     * Construct a solver.
075     *
076     * @param relativeAccuracy Relative accuracy.
077     * @param absoluteAccuracy Absolute accuracy.
078     * @param functionValueAccuracy Function value accuracy.
079     */
080    public LaguerreSolver(double relativeAccuracy,
081                          double absoluteAccuracy,
082                          double functionValueAccuracy) {
083        super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
084    }
085
086    /**
087     * {@inheritDoc}
088     */
089    @Override
090    public double doSolve()
091        throws TooManyEvaluationsException,
092               NumberIsTooLargeException,
093               NoBracketingException {
094        final double min = getMin();
095        final double max = getMax();
096        final double initial = getStartValue();
097        final double functionValueAccuracy = getFunctionValueAccuracy();
098
099        verifySequence(min, initial, max);
100
101        // Return the initial guess if it is good enough.
102        final double yInitial = computeObjectiveValue(initial);
103        if (JdkMath.abs(yInitial) <= functionValueAccuracy) {
104            return initial;
105        }
106
107        // Return the first endpoint if it is good enough.
108        final double yMin = computeObjectiveValue(min);
109        if (JdkMath.abs(yMin) <= functionValueAccuracy) {
110            return min;
111        }
112
113        // Reduce interval if min and initial bracket the root.
114        if (yInitial * yMin < 0) {
115            return laguerre(min, initial);
116        }
117
118        // Return the second endpoint if it is good enough.
119        final double yMax = computeObjectiveValue(max);
120        if (JdkMath.abs(yMax) <= functionValueAccuracy) {
121            return max;
122        }
123
124        // Reduce interval if initial and max bracket the root.
125        if (yInitial * yMax < 0) {
126            return laguerre(initial, max);
127        }
128
129        throw new NoBracketingException(min, max, yMin, yMax);
130    }
131
132    /**
133     * Find a real root in the given interval.
134     *
135     * Despite the bracketing condition, the root returned by
136     * {@link LaguerreSolver.ComplexSolver#solve(Complex[],Complex)} may
137     * not be a real zero inside {@code [min, max]}.
138     * For example, <code> p(x) = x<sup>3</sup> + 1, </code>
139     * with {@code min = -2}, {@code max = 2}, {@code initial = 0}.
140     * When it occurs, this code calls
141     * {@link LaguerreSolver.ComplexSolver#solveAll(Complex[],Complex)}
142     * in order to obtain all roots and picks up one real root.
143     *
144     * @param lo Lower bound of the search interval.
145     * @param hi Higher bound of the search interval.
146     * @return the point at which the function value is zero.
147     */
148    private double laguerre(double lo, double hi) {
149        final Complex[] c = real2Complex(getCoefficients());
150
151        final Complex initial = Complex.ofCartesian(0.5 * (lo + hi), 0);
152        final Complex z = complexSolver.solve(c, initial);
153        if (complexSolver.isRoot(lo, hi, z)) {
154            return z.getReal();
155        } else {
156            double r = Double.NaN;
157            // Solve all roots and select the one we are seeking.
158            Complex[] root = complexSolver.solveAll(c, initial);
159            for (int i = 0; i < root.length; i++) {
160                if (complexSolver.isRoot(lo, hi, root[i])) {
161                    r = root[i].getReal();
162                    break;
163                }
164            }
165            return r;
166        }
167    }
168
169    /**
170     * Find all complex roots for the polynomial with the given
171     * coefficients, starting from the given initial value.
172     * <p>
173     * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p>
174     *
175     * @param coefficients Polynomial coefficients.
176     * @param initial Start value.
177     * @return the point at which the function value is zero.
178     * @throws org.apache.commons.math4.legacy.exception.TooManyEvaluationsException
179     * if the maximum number of evaluations is exceeded.
180     * @throws NullArgumentException if the {@code coefficients} is
181     * {@code null}.
182     * @throws NoDataException if the {@code coefficients} array is empty.
183     * @since 3.1
184     */
185    public Complex[] solveAllComplex(double[] coefficients,
186                                     double initial)
187        throws NullArgumentException,
188               NoDataException,
189               TooManyEvaluationsException {
190        setup(Integer.MAX_VALUE,
191              new PolynomialFunction(coefficients),
192              Double.NEGATIVE_INFINITY,
193              Double.POSITIVE_INFINITY,
194              initial);
195        return complexSolver.solveAll(real2Complex(coefficients),
196                                      Complex.ofCartesian(initial, 0d));
197    }
198
199    /**
200     * Find a complex root for the polynomial with the given coefficients,
201     * starting from the given initial value.
202     * <p>
203     * Note: This method is not part of the API of {@link BaseUnivariateSolver}.</p>
204     *
205     * @param coefficients Polynomial coefficients.
206     * @param initial Start value.
207     * @return the point at which the function value is zero.
208     * @throws org.apache.commons.math4.legacy.exception.TooManyEvaluationsException
209     * if the maximum number of evaluations is exceeded.
210     * @throws NullArgumentException if the {@code coefficients} is
211     * {@code null}.
212     * @throws NoDataException if the {@code coefficients} array is empty.
213     * @since 3.1
214     */
215    public Complex solveComplex(double[] coefficients,
216                                double initial)
217        throws NullArgumentException,
218               NoDataException,
219               TooManyEvaluationsException {
220        setup(Integer.MAX_VALUE,
221              new PolynomialFunction(coefficients),
222              Double.NEGATIVE_INFINITY,
223              Double.POSITIVE_INFINITY,
224              initial);
225        return complexSolver.solve(real2Complex(coefficients),
226                                   Complex.ofCartesian(initial, 0d));
227    }
228
229    /**
230     * Class for searching all (complex) roots.
231     */
232    private final class ComplexSolver {
233        /**
234         * Check whether the given complex root is actually a real zero
235         * in the given interval, within the solver tolerance level.
236         *
237         * @param min Lower bound for the interval.
238         * @param max Upper bound for the interval.
239         * @param z Complex root.
240         * @return {@code true} if z is a real zero.
241         */
242        public boolean isRoot(double min, double max, Complex z) {
243            if (isSequence(min, z.getReal(), max)) {
244                double tolerance = JdkMath.max(getRelativeAccuracy() * z.abs(), getAbsoluteAccuracy());
245                return JdkMath.abs(z.getImaginary()) <= tolerance ||
246                     z.abs() <= getFunctionValueAccuracy();
247            }
248            return false;
249        }
250
251        /**
252         * Find all complex roots for the polynomial with the given
253         * coefficients, starting from the given initial value.
254         *
255         * @param coefficients Polynomial coefficients.
256         * @param initial Start value.
257         * @return the point at which the function value is zero.
258         * @throws org.apache.commons.math4.legacy.exception.TooManyEvaluationsException
259         * if the maximum number of evaluations is exceeded.
260         * @throws NullArgumentException if the {@code coefficients} is
261         * {@code null}.
262         * @throws NoDataException if the {@code coefficients} array is empty.
263         */
264        public Complex[] solveAll(Complex[] coefficients, Complex initial)
265            throws NullArgumentException,
266                   NoDataException,
267                   TooManyEvaluationsException {
268            if (coefficients == null) {
269                throw new NullArgumentException();
270            }
271            final int n = coefficients.length - 1;
272            if (n == 0) {
273                throw new NoDataException(LocalizedFormats.POLYNOMIAL);
274            }
275            // Coefficients for deflated polynomial.
276            final Complex[] c = coefficients.clone();
277
278            // Solve individual roots successively.
279            final Complex[] root = new Complex[n];
280            for (int i = 0; i < n; i++) {
281                final Complex[] subarray = new Complex[n - i + 1];
282                System.arraycopy(c, 0, subarray, 0, subarray.length);
283                root[i] = solve(subarray, initial);
284                // Polynomial deflation using synthetic division.
285                Complex newc = c[n - i];
286                Complex oldc = null;
287                for (int j = n - i - 1; j >= 0; j--) {
288                    oldc = c[j];
289                    c[j] = newc;
290                    newc = oldc.add(newc.multiply(root[i]));
291                }
292            }
293
294            return root;
295        }
296
297        /**
298         * Find a complex root for the polynomial with the given coefficients,
299         * starting from the given initial value.
300         *
301         * @param coefficients Polynomial coefficients.
302         * @param initial Start value.
303         * @return the point at which the function value is zero.
304         * @throws org.apache.commons.math4.legacy.exception.TooManyEvaluationsException
305         * if the maximum number of evaluations is exceeded.
306         * @throws NullArgumentException if the {@code coefficients} is
307         * {@code null}.
308         * @throws NoDataException if the {@code coefficients} array is empty.
309         */
310        public Complex solve(Complex[] coefficients, Complex initial)
311            throws NullArgumentException,
312                   NoDataException,
313                   TooManyEvaluationsException {
314            if (coefficients == null) {
315                throw new NullArgumentException();
316            }
317
318            final int n = coefficients.length - 1;
319            if (n == 0) {
320                throw new NoDataException(LocalizedFormats.POLYNOMIAL);
321            }
322
323            final double absoluteAccuracy = getAbsoluteAccuracy();
324            final double relativeAccuracy = getRelativeAccuracy();
325            final double functionValueAccuracy = getFunctionValueAccuracy();
326
327            final Complex nC  = Complex.ofCartesian(n, 0);
328            final Complex n1C = Complex.ofCartesian(n - 1, 0);
329
330            Complex z = initial;
331            Complex oldz = Complex.ofCartesian(Double.POSITIVE_INFINITY,
332                                       Double.POSITIVE_INFINITY);
333            while (true) {
334                // Compute pv (polynomial value), dv (derivative value), and
335                // d2v (second derivative value) simultaneously.
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(2);
345
346                // Check for convergence.
347                final double tolerance = JdkMath.max(relativeAccuracy * z.abs(),
348                                                      absoluteAccuracy);
349                if ((z.subtract(oldz)).abs() <= tolerance) {
350                    return z;
351                }
352                if (pv.abs() <= functionValueAccuracy) {
353                    return z;
354                }
355
356                // Now pv != 0, calculate the new approximation.
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                // Choose a denominator larger in magnitude.
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.abs() > dminus.abs() ? dplus : dminus;
366                // Perturb z if denominator is zero, for instance,
367                // p(x) = x^3 + 1, z = 0.
368                // This uses exact equality to zero. A tolerance may be required here.
369                if (denominator.equals(Complex.ZERO)) {
370                    z = z.add(Complex.ofCartesian(absoluteAccuracy, absoluteAccuracy));
371                    oldz = Complex.ofCartesian(Double.POSITIVE_INFINITY,
372                                       Double.POSITIVE_INFINITY);
373                } else {
374                    oldz = z;
375                    z = z.subtract(nC.divide(denominator));
376                }
377                incrementEvaluationCount();
378            }
379        }
380    }
381
382    /**
383     * Converts a {@code double[]} array to a {@code Complex[]} array.
384     *
385     * @param real array of numbers to be converted to their {@code Complex} equivalent
386     * @return {@code Complex} array
387     */
388    private static Complex[] real2Complex(double[] real) {
389        int index = 0;
390        final Complex[] c = new Complex[real.length];
391        for (final double d : real) {
392            c[index] = Complex.ofCartesian(d, 0);
393            index++;
394        }
395        return c;
396    }
397}