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