View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.analysis.polynomials;
23  
24  import java.io.Serializable;
25  import java.util.Arrays;
26  
27  import org.hipparchus.CalculusFieldElement;
28  import org.hipparchus.analysis.FieldUnivariateFunction;
29  import org.hipparchus.analysis.ParametricUnivariateFunction;
30  import org.hipparchus.analysis.differentiation.Derivative;
31  import org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction;
32  import org.hipparchus.exception.LocalizedCoreFormats;
33  import org.hipparchus.exception.MathIllegalArgumentException;
34  import org.hipparchus.exception.NullArgumentException;
35  import org.hipparchus.util.FastMath;
36  import org.hipparchus.util.MathUtils;
37  
38  /**
39   * Immutable representation of a real polynomial function with real coefficients.
40   * <p>
41   * <a href="http://mathworld.wolfram.com/HornersMethod.html">Horner's Method</a>
42   * is used to evaluate the function.</p>
43   *
44   */
45  public class PolynomialFunction implements UnivariateDifferentiableFunction, FieldUnivariateFunction, Serializable {
46      /**
47       * Serialization identifier
48       */
49      private static final long serialVersionUID = -7726511984200295583L;
50      /**
51       * The coefficients of the polynomial, ordered by degree -- i.e.,
52       * coefficients[0] is the constant term and coefficients[n] is the
53       * coefficient of x^n where n is the degree of the polynomial.
54       */
55      private final double[] coefficients;
56  
57      /**
58       * Construct a polynomial with the given coefficients.  The first element
59       * of the coefficients array is the constant term.  Higher degree
60       * coefficients follow in sequence.  The degree of the resulting polynomial
61       * is the index of the last non-null element of the array, or 0 if all elements
62       * are null.
63       * <p>
64       * The constructor makes a copy of the input array and assigns the copy to
65       * the coefficients property.</p>
66       *
67       * @param c Polynomial coefficients.
68       * @throws NullArgumentException if {@code c} is {@code null}.
69       * @throws MathIllegalArgumentException if {@code c} is empty.
70       */
71      public PolynomialFunction(double... c)
72          throws MathIllegalArgumentException, NullArgumentException {
73          super();
74          MathUtils.checkNotNull(c);
75          int n = c.length;
76          if (n == 0) {
77              throw new MathIllegalArgumentException(LocalizedCoreFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY);
78          }
79          while ((n > 1) && (c[n - 1] == 0)) {
80              --n;
81          }
82          this.coefficients = new double[n];
83          System.arraycopy(c, 0, this.coefficients, 0, n);
84      }
85  
86      /**
87       * Compute the value of the function for the given argument.
88       * <p>
89       *  The value returned is </p><p>
90       *  {@code coefficients[n] * x^n + ... + coefficients[1] * x  + coefficients[0]}
91       * </p>
92       *
93       * @param x Argument for which the function value should be computed.
94       * @return the value of the polynomial at the given point.
95       *
96       * @see org.hipparchus.analysis.UnivariateFunction#value(double)
97       */
98      @Override
99      public double value(double x) {
100        return evaluate(coefficients, x);
101     }
102 
103     /**
104      * Returns the degree of the polynomial.
105      *
106      * @return the degree of the polynomial.
107      */
108     public int degree() {
109         return coefficients.length - 1;
110     }
111 
112     /**
113      * Returns a copy of the coefficients array.
114      * <p>
115      * Changes made to the returned copy will not affect the coefficients of
116      * the polynomial.</p>
117      *
118      * @return a fresh copy of the coefficients array.
119      */
120     public double[] getCoefficients() {
121         return coefficients.clone();
122     }
123 
124     /**
125      * Uses Horner's Method to evaluate the polynomial with the given coefficients at
126      * the argument.
127      *
128      * @param coefficients Coefficients of the polynomial to evaluate.
129      * @param argument Input value.
130      * @return the value of the polynomial.
131      * @throws MathIllegalArgumentException if {@code coefficients} is empty.
132      * @throws NullArgumentException if {@code coefficients} is {@code null}.
133      */
134     protected static double evaluate(double[] coefficients, double argument)
135         throws MathIllegalArgumentException, NullArgumentException {
136         MathUtils.checkNotNull(coefficients);
137         int n = coefficients.length;
138         if (n == 0) {
139             throw new MathIllegalArgumentException(LocalizedCoreFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY);
140         }
141         double result = coefficients[n - 1];
142         for (int j = n - 2; j >= 0; j--) {
143             result = argument * result + coefficients[j];
144         }
145         return result;
146     }
147 
148 
149     /** {@inheritDoc}
150      * @throws MathIllegalArgumentException if {@code coefficients} is empty.
151      * @throws NullArgumentException if {@code coefficients} is {@code null}.
152      */
153     @Override
154     public <T extends Derivative<T>> T value(final T t)
155         throws MathIllegalArgumentException, NullArgumentException {
156         MathUtils.checkNotNull(coefficients);
157         int n = coefficients.length;
158         if (n == 0) {
159             throw new MathIllegalArgumentException(LocalizedCoreFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY);
160         }
161         T result = t.getField().getZero().add(coefficients[n - 1]);
162         for (int j = n - 2; j >= 0; j--) {
163             result = result.multiply(t).add(coefficients[j]);
164         }
165         return result;
166     }
167 
168     /** {@inheritDoc}
169      * @throws MathIllegalArgumentException if {@code coefficients} is empty.
170      * @throws NullArgumentException if {@code coefficients} is {@code null}.
171      * @since 1.3
172      */
173     @Override
174     public <T extends CalculusFieldElement<T>> T value(final T t)
175         throws MathIllegalArgumentException, NullArgumentException {
176         MathUtils.checkNotNull(coefficients);
177         int n = coefficients.length;
178         if (n == 0) {
179             throw new MathIllegalArgumentException(LocalizedCoreFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY);
180         }
181         T result = t.getField().getZero().add(coefficients[n - 1]);
182         for (int j = n - 2; j >= 0; j--) {
183             result = result.multiply(t).add(coefficients[j]);
184         }
185         return result;
186     }
187 
188     /**
189      * Add a polynomial to the instance.
190      *
191      * @param p Polynomial to add.
192      * @return a new polynomial which is the sum of the instance and {@code p}.
193      */
194     public PolynomialFunction add(final PolynomialFunction p) {
195         // identify the lowest degree polynomial
196         final int lowLength  = FastMath.min(coefficients.length, p.coefficients.length);
197         final 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         System.arraycopy((coefficients.length < p.coefficients.length) ?
205                          p.coefficients : coefficients,
206                          lowLength,
207                          newCoefficients, lowLength,
208                          highLength - lowLength);
209 
210         return new PolynomialFunction(newCoefficients);
211     }
212 
213     /**
214      * Subtract a polynomial from the instance.
215      *
216      * @param p Polynomial to subtract.
217      * @return a new polynomial which is the instance minus {@code p}.
218      */
219     public PolynomialFunction subtract(final PolynomialFunction p) {
220         // identify the lowest degree polynomial
221         int lowLength  = FastMath.min(coefficients.length, p.coefficients.length);
222         int highLength = FastMath.max(coefficients.length, p.coefficients.length);
223 
224         // build the coefficients array
225         double[] newCoefficients = new double[highLength];
226         for (int i = 0; i < lowLength; ++i) {
227             newCoefficients[i] = coefficients[i] - p.coefficients[i];
228         }
229         if (coefficients.length < p.coefficients.length) {
230             for (int i = lowLength; i < highLength; ++i) {
231                 newCoefficients[i] = -p.coefficients[i];
232             }
233         } else {
234             System.arraycopy(coefficients, lowLength, newCoefficients, lowLength,
235                              highLength - lowLength);
236         }
237 
238         return new PolynomialFunction(newCoefficients);
239     }
240 
241     /**
242      * Negate the instance.
243      *
244      * @return a new polynomial with all coefficients negated
245      */
246     public PolynomialFunction negate() {
247         double[] newCoefficients = new double[coefficients.length];
248         for (int i = 0; i < coefficients.length; ++i) {
249             newCoefficients[i] = -coefficients[i];
250         }
251         return new PolynomialFunction(newCoefficients);
252     }
253 
254     /**
255      * Multiply the instance by a polynomial.
256      *
257      * @param p Polynomial to multiply by.
258      * @return a new polynomial equal to this times {@code p}
259      */
260     public PolynomialFunction multiply(final PolynomialFunction p) {
261         double[] newCoefficients = new double[coefficients.length + p.coefficients.length - 1];
262 
263         for (int i = 0; i < newCoefficients.length; ++i) {
264             newCoefficients[i] = 0.0;
265             for (int j = FastMath.max(0, i + 1 - p.coefficients.length);
266                  j < FastMath.min(coefficients.length, i + 1);
267                  ++j) {
268                 newCoefficients[i] += coefficients[j] * p.coefficients[i-j];
269             }
270         }
271 
272         return new PolynomialFunction(newCoefficients);
273     }
274 
275     /**
276      * Returns the coefficients of the derivative of the polynomial with the given coefficients.
277      *
278      * @param coefficients Coefficients of the polynomial to differentiate.
279      * @return the coefficients of the derivative or {@code null} if coefficients has length 1.
280      * @throws MathIllegalArgumentException if {@code coefficients} is empty.
281      * @throws NullArgumentException if {@code coefficients} is {@code null}.
282      */
283     protected static double[] differentiate(double[] coefficients)
284         throws MathIllegalArgumentException, NullArgumentException {
285         MathUtils.checkNotNull(coefficients);
286         int n = coefficients.length;
287         if (n == 0) {
288             throw new MathIllegalArgumentException(LocalizedCoreFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY);
289         }
290         if (n == 1) {
291             return new double[]{0};
292         }
293         double[] result = new double[n - 1];
294         for (int i = n - 1; i > 0; i--) {
295             result[i - 1] = i * coefficients[i];
296         }
297         return result;
298     }
299 
300     /**
301      * Returns an anti-derivative of this polynomial, with 0 constant term.
302      *
303      * @return a polynomial whose derivative has the same coefficients as this polynomial
304      */
305     public PolynomialFunction antiDerivative() {
306         final int d = degree();
307         final double[] anti = new double[d + 2];
308         anti[0] = 0d;
309         for (int i = 1; i <= d + 1; i++) {
310             anti[i] = coefficients[i - 1]  / i;
311         }
312         return new PolynomialFunction(anti);
313     }
314 
315     /**
316      * Returns the definite integral of this polymomial over the given interval.
317      * <p>
318      * [lower, upper] must describe a finite interval (neither can be infinite
319      * and lower must be less than or equal to upper).
320      *
321      * @param lower lower bound for the integration
322      * @param upper upper bound for the integration
323      * @return the integral of this polymomial over the given interval
324      * @throws MathIllegalArgumentException if the bounds do not describe a finite interval
325      */
326     public double integrate(final double lower, final double upper) {
327         if (Double.isInfinite(lower) || Double.isInfinite(upper)) {
328             throw new MathIllegalArgumentException(LocalizedCoreFormats.INFINITE_BOUND);
329         }
330         if (lower > upper) {
331             throw new MathIllegalArgumentException(LocalizedCoreFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND);
332         }
333         final PolynomialFunction anti = antiDerivative();
334         return anti.value(upper) - anti.value(lower);
335     }
336 
337     /**
338      * Returns the derivative as a {@link PolynomialFunction}.
339      *
340      * @return the derivative polynomial.
341      */
342     public PolynomialFunction polynomialDerivative() {
343         return new PolynomialFunction(differentiate(coefficients));
344     }
345 
346     /**
347      * Returns a string representation of the polynomial.
348      *
349      * <p>The representation is user oriented. Terms are displayed lowest
350      * degrees first. The multiplications signs, coefficients equals to
351      * one and null terms are not displayed (except if the polynomial is 0,
352      * in which case the 0 constant term is displayed). Addition of terms
353      * with negative coefficients are replaced by subtraction of terms
354      * with positive coefficients except for the first displayed term
355      * (i.e. we display <code>-3</code> for a constant negative polynomial,
356      * but <code>1 - 3 x + x^2</code> if the negative coefficient is not
357      * the first one displayed).</p>
358      *
359      * @return a string representation of the polynomial.
360      */
361     @Override
362     public String toString() {
363         StringBuilder s = new StringBuilder();
364         if (coefficients[0] == 0.0) {
365             if (coefficients.length == 1) {
366                 return "0";
367             }
368         } else {
369             s.append(toString(coefficients[0]));
370         }
371 
372         for (int i = 1; i < coefficients.length; ++i) {
373             if (coefficients[i] != 0) {
374                 if (s.length() > 0) {
375                     if (coefficients[i] < 0) {
376                         s.append(" - ");
377                     } else {
378                         s.append(" + ");
379                     }
380                 } else {
381                     if (coefficients[i] < 0) {
382                         s.append('-');
383                     }
384                 }
385 
386                 double absAi = FastMath.abs(coefficients[i]);
387                 if ((absAi - 1) != 0) {
388                     s.append(toString(absAi));
389                     s.append(' ');
390                 }
391 
392                 s.append('x');
393                 if (i > 1) {
394                     s.append('^');
395                     s.append(i);
396                 }
397             }
398         }
399 
400         return s.toString();
401     }
402 
403     /**
404      * Creates a string representing a coefficient, removing ".0" endings.
405      *
406      * @param coeff Coefficient.
407      * @return a string representation of {@code coeff}.
408      */
409     private static String toString(double coeff) {
410         final String c = Double.toString(coeff);
411         if (c.endsWith(".0")) {
412             return c.substring(0, c.length() - 2);
413         } else {
414             return c;
415         }
416     }
417 
418     /** {@inheritDoc} */
419     @Override
420     public int hashCode() {
421         final int prime = 31;
422         int result = 1;
423         result = prime * result + Arrays.hashCode(coefficients);
424         return result;
425     }
426 
427     /** {@inheritDoc} */
428     @Override
429     public boolean equals(Object obj) {
430         if (this == obj) {
431             return true;
432         }
433         if (!(obj instanceof PolynomialFunction)) {
434             return false;
435         }
436         PolynomialFunction other = (PolynomialFunction) obj;
437         return Arrays.equals(coefficients, other.coefficients);
438     }
439 
440     /**
441      * Dedicated parametric polynomial class.
442      *
443      */
444     public static class Parametric implements ParametricUnivariateFunction {
445 
446         /** Empty constructor.
447          * <p>
448          * This constructor is not strictly necessary, but it prevents spurious
449          * javadoc warnings with JDK 18 and later.
450          * </p>
451          * @since 3.0
452          */
453         public Parametric() { // NOPMD - unnecessary constructor added intentionally to make javadoc happy
454             // nothing to do
455         }
456 
457         /** {@inheritDoc} */
458         @Override
459         public double[] gradient(double x, double ... parameters) {
460             final double[] gradient = new double[parameters.length];
461             double xn = 1.0;
462             for (int i = 0; i < parameters.length; ++i) {
463                 gradient[i] = xn;
464                 xn *= x;
465             }
466             return gradient;
467         }
468 
469         /** {@inheritDoc} */
470         @Override
471         public double value(final double x, final double ... parameters)
472             throws MathIllegalArgumentException {
473             return PolynomialFunction.evaluate(parameters, x);
474         }
475     }
476 }