    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    }