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