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     */
017    package org.apache.commons.math3.analysis.polynomials;
018    
019    import java.io.Serializable;
020    import java.util.Arrays;
021    
022    import org.apache.commons.math3.exception.util.LocalizedFormats;
023    import org.apache.commons.math3.exception.NoDataException;
024    import org.apache.commons.math3.exception.NullArgumentException;
025    import org.apache.commons.math3.analysis.DifferentiableUnivariateFunction;
026    import org.apache.commons.math3.analysis.UnivariateFunction;
027    import org.apache.commons.math3.analysis.ParametricUnivariateFunction;
028    import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
029    import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction;
030    import org.apache.commons.math3.util.FastMath;
031    import org.apache.commons.math3.util.MathUtils;
032    
033    /**
034     * Immutable representation of a real polynomial function with real coefficients.
035     * <p>
036     * <a href="http://mathworld.wolfram.com/HornersMethod.html">Horner's Method</a>
037     * is used to evaluate the function.</p>
038     *
039     * @version $Id: PolynomialFunction.java 1383441 2012-09-11 14:56:39Z luc $
040     */
041    public class PolynomialFunction implements UnivariateDifferentiableFunction, DifferentiableUnivariateFunction, Serializable {
042        /**
043         * Serialization identifier
044         */
045        private static final long serialVersionUID = -7726511984200295583L;
046        /**
047         * The coefficients of the polynomial, ordered by degree -- i.e.,
048         * coefficients[0] is the constant term and coefficients[n] is the
049         * coefficient of x^n where n is the degree of the polynomial.
050         */
051        private final double coefficients[];
052    
053        /**
054         * Construct a polynomial with the given coefficients.  The first element
055         * of the coefficients array is the constant term.  Higher degree
056         * coefficients follow in sequence.  The degree of the resulting polynomial
057         * is the index of the last non-null element of the array, or 0 if all elements
058         * are null.
059         * <p>
060         * The constructor makes a copy of the input array and assigns the copy to
061         * the coefficients property.</p>
062         *
063         * @param c Polynomial coefficients.
064         * @throws NullArgumentException if {@code c} is {@code null}.
065         * @throws NoDataException if {@code c} is empty.
066         */
067        public PolynomialFunction(double c[])
068            throws NullArgumentException, NoDataException {
069            super();
070            MathUtils.checkNotNull(c);
071            int n = c.length;
072            if (n == 0) {
073                throw new NoDataException(LocalizedFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY);
074            }
075            while ((n > 1) && (c[n - 1] == 0)) {
076                --n;
077            }
078            this.coefficients = new double[n];
079            System.arraycopy(c, 0, this.coefficients, 0, n);
080        }
081    
082        /**
083         * Compute the value of the function for the given argument.
084         * <p>
085         *  The value returned is <br/>
086         *  <code>coefficients[n] * x^n + ... + coefficients[1] * x  + coefficients[0]</code>
087         * </p>
088         *
089         * @param x Argument for which the function value should be computed.
090         * @return the value of the polynomial at the given point.
091         * @see UnivariateFunction#value(double)
092         */
093        public double value(double x) {
094           return evaluate(coefficients, x);
095        }
096    
097        /**
098         * Returns the degree of the polynomial.
099         *
100         * @return the degree of the polynomial.
101         */
102        public int degree() {
103            return coefficients.length - 1;
104        }
105    
106        /**
107         * Returns a copy of the coefficients array.
108         * <p>
109         * Changes made to the returned copy will not affect the coefficients of
110         * the polynomial.</p>
111         *
112         * @return a fresh copy of the coefficients array.
113         */
114        public double[] getCoefficients() {
115            return coefficients.clone();
116        }
117    
118        /**
119         * Uses Horner's Method to evaluate the polynomial with the given coefficients at
120         * the argument.
121         *
122         * @param coefficients Coefficients of the polynomial to evaluate.
123         * @param argument Input value.
124         * @return the value of the polynomial.
125         * @throws NoDataException if {@code coefficients} is empty.
126         * @throws NullArgumentException if {@code coefficients} is {@code null}.
127         */
128        protected static double evaluate(double[] coefficients, double argument)
129            throws NullArgumentException, NoDataException {
130            MathUtils.checkNotNull(coefficients);
131            int n = coefficients.length;
132            if (n == 0) {
133                throw new NoDataException(LocalizedFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY);
134            }
135            double result = coefficients[n - 1];
136            for (int j = n - 2; j >= 0; j--) {
137                result = argument * result + coefficients[j];
138            }
139            return result;
140        }
141    
142    
143        /** {@inheritDoc}
144         * @since 3.1
145         * @throws NoDataException if {@code coefficients} is empty.
146         * @throws NullArgumentException if {@code coefficients} is {@code null}.
147         */
148        public DerivativeStructure value(final DerivativeStructure t)
149            throws NullArgumentException, NoDataException {
150            MathUtils.checkNotNull(coefficients);
151            int n = coefficients.length;
152            if (n == 0) {
153                throw new NoDataException(LocalizedFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY);
154            }
155            DerivativeStructure result =
156                    new DerivativeStructure(t.getFreeParameters(), t.getOrder(), coefficients[n - 1]);
157            for (int j = n - 2; j >= 0; j--) {
158                result = result.multiply(t).add(coefficients[j]);
159            }
160            return result;
161        }
162    
163        /**
164         * Add a polynomial to the instance.
165         *
166         * @param p Polynomial to add.
167         * @return a new polynomial which is the sum of the instance and {@code p}.
168         */
169        public PolynomialFunction add(final PolynomialFunction p) {
170            // identify the lowest degree polynomial
171            final int lowLength  = FastMath.min(coefficients.length, p.coefficients.length);
172            final int highLength = FastMath.max(coefficients.length, p.coefficients.length);
173    
174            // build the coefficients array
175            double[] newCoefficients = new double[highLength];
176            for (int i = 0; i < lowLength; ++i) {
177                newCoefficients[i] = coefficients[i] + p.coefficients[i];
178            }
179            System.arraycopy((coefficients.length < p.coefficients.length) ?
180                             p.coefficients : coefficients,
181                             lowLength,
182                             newCoefficients, lowLength,
183                             highLength - lowLength);
184    
185            return new PolynomialFunction(newCoefficients);
186        }
187    
188        /**
189         * Subtract a polynomial from the instance.
190         *
191         * @param p Polynomial to subtract.
192         * @return a new polynomial which is the difference the instance minus {@code p}.
193         */
194        public PolynomialFunction subtract(final PolynomialFunction p) {
195            // identify the lowest degree polynomial
196            int lowLength  = FastMath.min(coefficients.length, p.coefficients.length);
197            int highLength = FastMath.max(coefficients.length, p.coefficients.length);
198    
199            // build the coefficients array
200            double[] newCoefficients = new double[highLength];
201            for (int i = 0; i < lowLength; ++i) {
202                newCoefficients[i] = coefficients[i] - p.coefficients[i];
203            }
204            if (coefficients.length < p.coefficients.length) {
205                for (int i = lowLength; i < highLength; ++i) {
206                    newCoefficients[i] = -p.coefficients[i];
207                }
208            } else {
209                System.arraycopy(coefficients, lowLength, newCoefficients, lowLength,
210                                 highLength - lowLength);
211            }
212    
213            return new PolynomialFunction(newCoefficients);
214        }
215    
216        /**
217         * Negate the instance.
218         *
219         * @return a new polynomial.
220         */
221        public PolynomialFunction negate() {
222            double[] newCoefficients = new double[coefficients.length];
223            for (int i = 0; i < coefficients.length; ++i) {
224                newCoefficients[i] = -coefficients[i];
225            }
226            return new PolynomialFunction(newCoefficients);
227        }
228    
229        /**
230         * Multiply the instance by a polynomial.
231         *
232         * @param p Polynomial to multiply by.
233         * @return a new polynomial.
234         */
235        public PolynomialFunction multiply(final PolynomialFunction p) {
236            double[] newCoefficients = new double[coefficients.length + p.coefficients.length - 1];
237    
238            for (int i = 0; i < newCoefficients.length; ++i) {
239                newCoefficients[i] = 0.0;
240                for (int j = FastMath.max(0, i + 1 - p.coefficients.length);
241                     j < FastMath.min(coefficients.length, i + 1);
242                     ++j) {
243                    newCoefficients[i] += coefficients[j] * p.coefficients[i-j];
244                }
245            }
246    
247            return new PolynomialFunction(newCoefficients);
248        }
249    
250        /**
251         * Returns the coefficients of the derivative of the polynomial with the given coefficients.
252         *
253         * @param coefficients Coefficients of the polynomial to differentiate.
254         * @return the coefficients of the derivative or {@code null} if coefficients has length 1.
255         * @throws NoDataException if {@code coefficients} is empty.
256         * @throws NullArgumentException if {@code coefficients} is {@code null}.
257         */
258        protected static double[] differentiate(double[] coefficients)
259            throws NullArgumentException, NoDataException {
260            MathUtils.checkNotNull(coefficients);
261            int n = coefficients.length;
262            if (n == 0) {
263                throw new NoDataException(LocalizedFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY);
264            }
265            if (n == 1) {
266                return new double[]{0};
267            }
268            double[] result = new double[n - 1];
269            for (int i = n - 1; i > 0; i--) {
270                result[i - 1] = i * coefficients[i];
271            }
272            return result;
273        }
274    
275        /**
276         * Returns the derivative as a {@link PolynomialFunction}.
277         *
278         * @return the derivative polynomial.
279         */
280        public PolynomialFunction polynomialDerivative() {
281            return new PolynomialFunction(differentiate(coefficients));
282        }
283    
284        /**
285         * Returns the derivative as a {@link UnivariateFunction}.
286         *
287         * @return the derivative function.
288         */
289        public UnivariateFunction derivative() {
290            return polynomialDerivative();
291        }
292    
293        /**
294         * Returns a string representation of the polynomial.
295         *
296         * <p>The representation is user oriented. Terms are displayed lowest
297         * degrees first. The multiplications signs, coefficients equals to
298         * one and null terms are not displayed (except if the polynomial is 0,
299         * in which case the 0 constant term is displayed). Addition of terms
300         * with negative coefficients are replaced by subtraction of terms
301         * with positive coefficients except for the first displayed term
302         * (i.e. we display <code>-3</code> for a constant negative polynomial,
303         * but <code>1 - 3 x + x^2</code> if the negative coefficient is not
304         * the first one displayed).</p>
305         *
306         * @return a string representation of the polynomial.
307         */
308        @Override
309        public String toString() {
310            StringBuilder s = new StringBuilder();
311            if (coefficients[0] == 0.0) {
312                if (coefficients.length == 1) {
313                    return "0";
314                }
315            } else {
316                s.append(toString(coefficients[0]));
317            }
318    
319            for (int i = 1; i < coefficients.length; ++i) {
320                if (coefficients[i] != 0) {
321                    if (s.length() > 0) {
322                        if (coefficients[i] < 0) {
323                            s.append(" - ");
324                        } else {
325                            s.append(" + ");
326                        }
327                    } else {
328                        if (coefficients[i] < 0) {
329                            s.append("-");
330                        }
331                    }
332    
333                    double absAi = FastMath.abs(coefficients[i]);
334                    if ((absAi - 1) != 0) {
335                        s.append(toString(absAi));
336                        s.append(' ');
337                    }
338    
339                    s.append("x");
340                    if (i > 1) {
341                        s.append('^');
342                        s.append(Integer.toString(i));
343                    }
344                }
345            }
346    
347            return s.toString();
348        }
349    
350        /**
351         * Creates a string representing a coefficient, removing ".0" endings.
352         *
353         * @param coeff Coefficient.
354         * @return a string representation of {@code coeff}.
355         */
356        private static String toString(double coeff) {
357            final String c = Double.toString(coeff);
358            if (c.endsWith(".0")) {
359                return c.substring(0, c.length() - 2);
360            } else {
361                return c;
362            }
363        }
364    
365        /** {@inheritDoc} */
366        @Override
367        public int hashCode() {
368            final int prime = 31;
369            int result = 1;
370            result = prime * result + Arrays.hashCode(coefficients);
371            return result;
372        }
373    
374        /** {@inheritDoc} */
375        @Override
376        public boolean equals(Object obj) {
377            if (this == obj) {
378                return true;
379            }
380            if (!(obj instanceof PolynomialFunction)) {
381                return false;
382            }
383            PolynomialFunction other = (PolynomialFunction) obj;
384            if (!Arrays.equals(coefficients, other.coefficients)) {
385                return false;
386            }
387            return true;
388        }
389    
390        /**
391         * Dedicated parametric polynomial class.
392         *
393         * @since 3.0
394         */
395        public static class Parametric implements ParametricUnivariateFunction {
396            /** {@inheritDoc} */
397            public double[] gradient(double x, double ... parameters) {
398                final double[] gradient = new double[parameters.length];
399                double xn = 1.0;
400                for (int i = 0; i < parameters.length; ++i) {
401                    gradient[i] = xn;
402                    xn *= x;
403                }
404                return gradient;
405            }
406    
407            /** {@inheritDoc} */
408            public double value(final double x, final double ... parameters) {
409                return PolynomialFunction.evaluate(parameters, x);
410            }
411        }
412    }