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.polynomials;
018
019import java.io.Serializable;
020import java.util.Arrays;
021
022import org.apache.commons.math3.exception.util.LocalizedFormats;
023import org.apache.commons.math3.exception.NoDataException;
024import org.apache.commons.math3.exception.NullArgumentException;
025import org.apache.commons.math3.analysis.DifferentiableUnivariateFunction;
026import org.apache.commons.math3.analysis.UnivariateFunction;
027import org.apache.commons.math3.analysis.ParametricUnivariateFunction;
028import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
029import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction;
030import org.apache.commons.math3.util.FastMath;
031import 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 1455194 2013-03-11 15:45:54Z luc $
040 */
041public 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            throws NoDataException {
410            return PolynomialFunction.evaluate(parameters, x);
411        }
412    }
413}