/*
    * Copyright (C) 2016 Alberto Irurueta Carro (alberto@irurueta.com)
    *
    * Licensed under the Apache License, Version 2.0 (the "License");
    * you may not use this file except in compliance with the License.
    * You may obtain a copy of the License at
    *
    *         http://www.apache.org/licenses/LICENSE-2.0
    *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  package com.irurueta.numerical.polynomials;
  
  import com.irurueta.algebra.ArrayUtils;
  import com.irurueta.algebra.Complex;
  import com.irurueta.numerical.NumericalException;
  import com.irurueta.numerical.roots.FirstDegreePolynomialRootsEstimator;
  import com.irurueta.numerical.roots.LaguerrePolynomialRootsEstimator;
  import com.irurueta.numerical.roots.PolynomialRootsEstimator;
  import com.irurueta.numerical.roots.SecondDegreePolynomialRootsEstimator;
  import com.irurueta.numerical.roots.ThirdDegreePolynomialRootsEstimator;
  import com.irurueta.numerical.signal.processing.Convolver1D;
  
  import java.io.Serializable;
  import java.util.ArrayList;
  
  /**
   * Contains a polynomial and common operations done with polynomials.
   * This implementation only supports polynomials with real parameters.
   */
  public class Polynomial implements Serializable {
  
      /**
       * Minimum derivative / integration order.
       */
      private static final int MIN_ORDER = 1;
  
      /**
       * Minimum allowed length in polynomial parameters.
       */
      public static final int MIN_VALID_POLY_PARAMS_LENGTH = 1;
  
      /**
       * Constant defining machine precision
       */
      public static final double EPS = 1e-10;
  
      /**
       * Array containing parameters defining a polynomial.
       * For a polynomial having the expression p(x) = a + b*x + c*x^2 + ...
       * provided array must be [a, b, c, ...]
       */
      private double[] polyParams;
  
      /**
       * Constructor.
       * Creates a polynomial initialized to zero.
       */
      public Polynomial() {
          polyParams = new double[MIN_VALID_POLY_PARAMS_LENGTH];
      }
  
      /**
       * Constructor.
       *
       * @param numberOfParameters number of parameters of polynomial to create.
       * @throws IllegalArgumentException if number of parameters is less than 1.
       */
      public Polynomial(final int numberOfParameters) {
          if (numberOfParameters < MIN_VALID_POLY_PARAMS_LENGTH) {
              throw new IllegalArgumentException("at least 1 parameter is required");
          }
          polyParams = new double[numberOfParameters];
      }
  
      /**
       * Constructor.
       * For a polynomial having the expression p(x) = a + b*x + c*x^2 + ...
       * provided array must be [a, b, c, ...]
       *
       * @param polyParams parameters defining a polynomial.
       * @throws IllegalArgumentException if provided array does not have at least
       *                                  length 2.
       */
      public Polynomial(final double... polyParams) {
          setPolyParams(polyParams);
      }
  
      /**
       * Gets array defining parameters of polynomial.
       * A polynomial having the expression p(x) = a + b*x + c*x^2 + ...
       * has an array of the form [a, b, c, ...].
       *
       * @return parameters defining a polynomial.
       */
     public double[] getPolyParams() {
         return polyParams;
     }
 
     /**
      * Sets array defining parameters of polynomial.
      * A polynomial having the expression p(x) = a + b*x + c*x^2 + ...
      * has an array of the form [a, b, c, ...].
      *
      * @param polyParams array defining parameters of polynomial. Must have at
      *                   least length 2.
      * @throws IllegalArgumentException if provided array does not have at least
      *                                  length 2.
      */
     public final void setPolyParams(final double... polyParams) {
         if (polyParams.length < MIN_VALID_POLY_PARAMS_LENGTH) {
             throw new IllegalArgumentException("must have at least length 1");
         }
 
         this.polyParams = polyParams;
     }
 
     /**
      * Gets degree of polynomial.
      *
      * @return degree of polynomial.
      */
     public int getDegree() {
         for (var i = polyParams.length - 1; i >= 1; i--) {
             if (Math.abs(polyParams[i]) > EPS) {
                 return i;
             }
         }
 
         return 0;
     }
 
     /**
      * Adds this polynomial to another one and stores the result into provided
      * instance.
      *
      * @param other  other polynomial to be added.
      * @param result instance where result will be stored.
      */
     @SuppressWarnings("Duplicates")
     public void add(final Polynomial other, final Polynomial result) {
         final var maxLength = Math.max(polyParams.length, other.polyParams.length);
         final var minLength = Math.min(polyParams.length, other.polyParams.length);
 
         var resultPolyParams = result.polyParams;
         if (resultPolyParams.length != maxLength) {
             resultPolyParams = new double[maxLength];
         }
 
         for (var i = 0; i < minLength; i++) {
             resultPolyParams[i] = polyParams[i] + other.polyParams[i];
         }
 
         if (polyParams.length > other.polyParams.length) {
             // this is longer than other
             System.arraycopy(polyParams, minLength, resultPolyParams, minLength, maxLength - minLength);
         } else {
             // other is longer than this
             System.arraycopy(other.polyParams, minLength, resultPolyParams, minLength, maxLength - minLength);
         }
 
         result.polyParams = resultPolyParams;
     }
 
     /**
      * Adds another polynomial to this polynomial.
      *
      * @param other other polynomial to be added.
      */
     public void add(final Polynomial other) {
         add(other, this);
     }
 
     /**
      * Adds this polynomial to another one and returns a new polynomial as a
      * result.
      *
      * @param other other polynomial to be added.
      * @return a new polynomial containing the sum.
      */
     public Polynomial addAndReturnNew(final Polynomial other) {
         final var length = Math.max(polyParams.length, other.polyParams.length);
         final var result = new Polynomial(length);
         add(other, result);
 
         return result;
     }
 
     /**
      * Subtract other polynomial from this one and stores the result into
      * provided instance.
      *
      * @param other  other polynomial to be subtracted from this one.
      * @param result instance where result will be stored.
      */
     @SuppressWarnings("Duplicates")
     public void subtract(final Polynomial other, final Polynomial result) {
         final var maxLength = Math.max(polyParams.length, other.polyParams.length);
         final var minLength = Math.min(polyParams.length, other.polyParams.length);
 
         var resultPolyParams = result.polyParams;
         if (resultPolyParams.length != maxLength) {
             resultPolyParams = new double[maxLength];
         }
 
         for (var i = 0; i < minLength; i++) {
             resultPolyParams[i] = polyParams[i] - other.polyParams[i];
         }
 
         if (polyParams.length > other.polyParams.length) {
             // this is longer than other
             System.arraycopy(polyParams, minLength, resultPolyParams, minLength, maxLength - minLength);
         } else {
             // other is longer than this
             for (var i = minLength; i < maxLength; i++) {
                 resultPolyParams[i] = -other.polyParams[i];
             }
         }
 
         result.polyParams = resultPolyParams;
     }
 
     /**
      * Subtracts another polynomial form this one.
      *
      * @param other other polynomial to be subtracted from this one.
      */
     public void subtract(final Polynomial other) {
         subtract(other, this);
     }
 
     /**
      * Subtract other polynomial from this one and returns a new polynomial as a
      * result.
      *
      * @param other other polynomial to be subtracted from this one.
      * @return a new polynomial containing result of subtraction.
      */
     public Polynomial subtractAndReturnNew(final Polynomial other) {
         final var length = Math.max(polyParams.length, other.polyParams.length);
         final var result = new Polynomial(length);
         subtract(other, result);
 
         return result;
     }
 
     /**
      * Multiplies two polynomials.
      *
      * @param other  other polynomial to multiply with.
      * @param result instance where resulting polynomial will be stored.
      */
     public void multiply(final Polynomial other, final Polynomial result) {
         final var thisLength = polyParams.length;
         final var otherLength = other.polyParams.length;
         final var resultLength = thisLength + otherLength - 1;
         if (result.polyParams.length != resultLength || result == this) {
             // if length does not match or result is stored in this polynomial,
             // create new polynomial array of parameters
             result.polyParams = Convolver1D.convolve(polyParams, other.polyParams);
         } else {
             // if length is the same, overwrite values
             Convolver1D.convolve(polyParams, other.polyParams, result.polyParams);
         }
     }
 
     /**
      * Multiplies this polynomial with another one.
      *
      * @param other other polynomial to multiply with.
      */
     public void multiply(final Polynomial other) {
         multiply(other, this);
     }
 
     /**
      * Multiplies two polynomials and returns a new instance containing result.
      *
      * @param other other polynomial to multiply with.
      * @return a new polynomial containing result of multiplication.
      */
     public Polynomial multiplyAndReturnNew(final Polynomial other) {
         final var thisLength = polyParams.length;
         final var otherLength = other.polyParams.length;
         final var resultLength = thisLength + otherLength - 1;
         final var result = new Polynomial(resultLength);
         Convolver1D.convolve(polyParams, other.polyParams, result.polyParams);
 
         return result;
     }
 
     /**
      * Multiplies all parameters of this polynomial by a scalar and stores the
      * result into provided polynomial instance.
      *
      * @param scalar scalar to multiply parameters with.
      * @param result instance where result will be stored.
      */
     public void multiplyByScalar(final double scalar, final Polynomial result) {
         var resultPolyParams = result.polyParams;
         if (resultPolyParams.length != polyParams.length || result == this) {
             resultPolyParams = new double[polyParams.length];
         }
         ArrayUtils.multiplyByScalar(polyParams, scalar, resultPolyParams);
         result.polyParams = resultPolyParams;
     }
 
     /**
      * Multiplies all parameters of this polynomial by provided scalar.
      *
      * @param scalar scalar to multiply parameters with.
      */
     public void multiplyByScalar(final double scalar) {
         multiplyByScalar(scalar, this);
     }
 
     /**
      * Multiplies all parameters of this polynomial by a scalar and returns a
      * new polynomial containing the result.
      *
      * @param scalar scalar to multiply parameters with.
      * @return a new polynomial containing the result of the operation.
      */
     public Polynomial multiplyByScalarAndReturnNew(final double scalar) {
         final var result = new Polynomial(polyParams.length);
         multiplyByScalar(scalar, result);
         return result;
     }
 
     /**
      * Gets roots of polynomial.
      *
      * @return estimated roots of this polynomial
      * @throws NumericalException if roots estimation fails.
      */
     public Complex[] getRoots() throws NumericalException {
         final var degree = getDegree();
 
         final PolynomialRootsEstimator estimator;
         switch (degree) {
             case 0:
                 // no roots
                 return null;
             case 1:
                 // first degree
                 estimator = new FirstDegreePolynomialRootsEstimator(polyParams);
                 break;
             case 2:
                 // second degree
                 estimator = new SecondDegreePolynomialRootsEstimator(polyParams);
                 break;
             case 3:
                 // third degree
                 estimator = new ThirdDegreePolynomialRootsEstimator(polyParams);
                 break;
             default:
                 // greater degree
 
                 // copy real parameters into complex values
                 final var params = new Complex[this.polyParams.length];
                 for (int i = 0; i < this.polyParams.length; i++) {
                     params[i] = new Complex(this.polyParams[i]);
                 }
                 estimator = new LaguerrePolynomialRootsEstimator(params);
                 break;
         }
 
         estimator.estimate();
         return estimator.getRoots();
     }
 
     /**
      * Evaluates polynomial at provided value.
      *
      * @param x value to evaluate polynomial at.
      * @return result of polynomial evaluation.
      */
     public double evaluate(final double x) {
         var result = 0.0;
         var powX = 1.0;
         for (var polyParam : polyParams) {
             result += polyParam * powX;
             powX *= x;
         }
 
         return result;
     }
 
     /**
      * Computes derivative of polynomial.
      *
      * @param result instance where derivative will be stored.
      */
     @SuppressWarnings("Duplicates")
     public void derivative(final Polynomial result) {
         final var resultLength = polyParams.length - 1;
         final var resultLength2 = Math.max(resultLength, 1);
 
         var resultPolyParams = result.polyParams;
         if (resultPolyParams.length != resultLength2 || result == this) {
             resultPolyParams = new double[resultLength2];
         }
         if (resultLength == 0) {
             resultPolyParams[0] = 0.0;
         }
 
         for (int i = 0, j = 1; i < resultLength; i++, j++) {
             resultPolyParams[i] = j * polyParams[j];
         }
 
         result.polyParams = resultPolyParams;
     }
 
     /**
      * Replaces this instance by its derivative.
      */
     public void derivative() {
         derivative(this);
     }
 
     /**
      * Computes derivative of polynomial.
      *
      * @return a new instance containing derivative.
      */
     public Polynomial derivativeAndReturnNew() {
         final var resultLength = Math.max(polyParams.length - 1, 1);
         final var result = new Polynomial(resultLength);
         derivative(result);
         return result;
     }
 
     /**
      * Evaluates derivative of polynomial at provided value.
      *
      * @param x value to evaluate derivative of polynomial at.
      * @return result of evaluation of derivative.
      */
     public double evaluateDerivative(final double x) {
         var result = 0.0;
         var powX = 1.0;
         for (var j = 1; j < polyParams.length; j++) {
             result += j * polyParams[j] * powX;
             powX *= x;
         }
 
         return result;
     }
 
     /**
      * Computes second derivative of polynomial.
      *
      * @param result instance where second derivative will be stored.
      */
     @SuppressWarnings("Duplicates")
     public void secondDerivative(final Polynomial result) {
         final var resultLength = polyParams.length - 2;
         final var resultLength2 = Math.max(resultLength, 1);
 
         var resultPolyParams = result.polyParams;
         if (resultPolyParams.length != resultLength2 || result == this) {
             resultPolyParams = new double[resultLength2];
         }
         if (resultLength == 0) {
             resultPolyParams[0] = 0.0;
         }
 
         for (int i = 0, j = 2, k = 1; i < resultLength; i++, j++, k++) {
             resultPolyParams[i] = j * k * polyParams[j];
         }
 
         result.polyParams = resultPolyParams;
     }
 
     /**
      * Replaces this instance by its second derivative.
      */
     public void secondDerivative() {
         secondDerivative(this);
     }
 
     /**
      * Computes second derivative of polynomial.
      *
      * @return a new instance containing second derivative.
      */
     public Polynomial secondDerivativeAndReturnNew() {
         final var resultLength = Math.max(polyParams.length - 2, 1);
         final var result = new Polynomial(resultLength);
         secondDerivative(result);
         return result;
     }
 
     /**
      * Evaluates second derivative of polynomial at provided value.
      *
      * @param x value to evaluate second derivative of polynomial at.
      * @return result of evaluation of second derivative.
      */
     public double evaluateSecondDerivative(final double x) {
         var result = 0.0;
         var powX = 1.0;
         for (int j = 2, k = 1; j < polyParams.length; j++, k++) {
             result += j * k * polyParams[j] * powX;
             powX *= x;
         }
 
         return result;
     }
 
     /**
      * Computes nth-order derivative of polynomial.
      *
      * @param order  order of derivative to compute. Must be at least 1.
      * @param result instance where nth-order derivative will be stored.
      * @throws IllegalArgumentException if provided order is less than 1.
      */
     @SuppressWarnings("Duplicates")
     public void nthDerivative(final int order, final Polynomial result) {
         if (order < MIN_ORDER) {
             throw new IllegalArgumentException();
         }
 
         final var resultLength = polyParams.length - order;
         final var resultLength2 = Math.max(resultLength, 1);
 
         var resultPolyParams = result.polyParams;
         if (resultPolyParams.length != resultLength2 || result == this) {
             resultPolyParams = new double[resultLength2];
         }
         if (resultLength == 0) {
             resultPolyParams[0] = 0.0;
         }
 
         for (int i = 0, j = order; i < resultLength; i++, j++) {
             var param = j;
             for (var k = 1; k < order; k++) {
                 param *= j - k;
             }
             resultPolyParams[i] = param * polyParams[j];
         }
 
         result.polyParams = resultPolyParams;
     }
 
     /**
      * Replaces this instance by its nth-order derivative.
      *
      * @param order order of derivative to compute. Must be at least 1.
      * @throws IllegalArgumentException if provided order is less than 1.
      */
     public void nthDerivative(final int order) {
         nthDerivative(order, this);
     }
 
     /**
      * Computes nth-order derivative of polynomial.
      *
      * @param order order of derivative to compute. Must be at least 1.
      * @return a new instance containing nth-order derivative.
      * @throws IllegalArgumentException if provided order is less than 1.
      */
     public Polynomial nthDerivativeAndReturnNew(final int order) {
         final var resultLength = Math.max(polyParams.length - order, 1);
         final var result = new Polynomial(resultLength);
         nthDerivative(order, result);
         return result;
     }
 
     /**
      * Evaluates nth-derivative of polynomial at provided value.
      *
      * @param x     value to evaluate nth-derivative of polynomial at.
      * @param order order of derivative to evaluate. Must be at least 1.
      * @return result of evaluation of nth-derivative.
      * @throws IllegalArgumentException if provided order is less than 1.
      */
     public double evaluateNthDerivative(final double x, final int order) {
         if (order < MIN_ORDER) {
             throw new IllegalArgumentException("order must be at least 1");
         }
 
         var result = 0.0;
         var powX = 1.0;
         for (var i = order; i < polyParams.length; i++) {
             var param = i;
             for (var j = 1; j < order; j++) {
                 param *= i - j;
             }
             result += param * polyParams[i] * powX;
             powX *= x;
         }
 
         return result;
     }
 
     /**
      * Computes polynomial containing the integration of current one.
      * Because infinite polynomials exist with different constant values,
      * constant term can be provided as well.
      *
      * @param result   instance where resulting polynomial will be stored.
      * @param constant constant term.
      */
     public void integration(final Polynomial result, final double constant) {
         final var resultLength = polyParams.length + 1;
         var resultPolyParams = result.polyParams;
         if (resultPolyParams.length != resultLength || result == this) {
             resultPolyParams = new double[resultLength];
         }
 
         resultPolyParams[0] = constant;
         for (int i = 0, j = 1; i < polyParams.length; i++, j++) {
             resultPolyParams[j] = polyParams[i] / j;
         }
 
         result.polyParams = resultPolyParams;
     }
 
     /**
      * Computes polynomial containing the integration of current one and
      * assuming a zero constant term.
      *
      * @param result instance where resulting polynomial will be stored.
      */
     public void integration(final Polynomial result) {
         integration(result, 0.0);
     }
 
     /**
      * Updates this instance to contain its integration.
      *
      * @param constant constant term.
      */
     public void integration(final double constant) {
         integration(this, constant);
     }
 
     /**
      * Updates this instance to contain its integration using a zero constant
      * term.
      */
     public void integration() {
         integration(this);
     }
 
     /**
      * Computes polynomial containing the integration of current one.
      * Because infinite polynomials exist with different constant values,
      * constant term can be provided as well.
      *
      * @param constant constant term.
      * @return a new instance containing integration polynomial.
      */
     public Polynomial integrationAndReturnNew(final double constant) {
         final var result = new Polynomial(polyParams.length + 1);
         integration(result, constant);
         return result;
     }
 
     /**
      * Computes polynomial containing the integration of current one and
      * assuming a zero constant term.
      *
      * @return a new instance containing integration polynomial.
      */
     public Polynomial integrationAndReturnNew() {
         return integrationAndReturnNew(0.0);
     }
 
     /**
      * Integrate polynomial within provided interval.
      *
      * @param startX start of integration interval.
      * @param endX   end of integration interval.
      * @return result of integration.
      */
     public double integrateInterval(final double startX, final double endX) {
 
         var resultStart = 0.0;
         var resultEnd = 0.0;
         var powStartX = startX;
         var powEndX = endX;
         double polyParam;
         for (int i = 0, j = 1; i < polyParams.length; i++, j++) {
             polyParam = polyParams[i] / j;
             resultStart += polyParam * powStartX;
             powStartX *= startX;
 
             resultEnd += polyParam * powEndX;
             powEndX *= endX;
         }
 
         return resultEnd - resultStart;
     }
 
     /**
      * Computes polynomial containing the nth-order integration of current one.
      * Because infinite polynomials exist with different constant values,
      * constant terms for each integration order can be provided as well.
      *
      * @param order     order of integration to compute. Must be at least 1.
      * @param result    instance where resulting polynomial will be stored.
      * @param constants constant terms for each integration order. Must have a
      *                  length equal to order if provided.
      * @throws IllegalArgumentException if provided order is less than 1 or if
      *                                  constants does not have length equal to order.
      */
     public void nthIntegration(final int order, final Polynomial result, final double[] constants) {
         if (order < MIN_ORDER) {
             throw new IllegalArgumentException("order must be at least 1");
         }
         if (constants != null && constants.length != order) {
             throw new IllegalArgumentException("length of constants must be order");
         }
         final var resultLength = polyParams.length + order;
         var resultPolyParams = result.polyParams;
         if (resultPolyParams.length != resultLength || result == this) {
             resultPolyParams = new double[resultLength];
         }
 
         for (var i = 0; i < order; i++) {
             if (constants != null) {
                 var param = 1;
                 for (var k = 1; k <= i; k++) {
                     param *= k;
                 }
                 resultPolyParams[i] = constants[i] / param;
             } else {
                 resultPolyParams[i] = 0.0;
             }
         }
         for (int i = 0, j = order; i < polyParams.length; i++, j++) {
             var param = j;
             for (var k = 1; k < order; k++) {
                 param *= j - k;
             }
             resultPolyParams[j] = polyParams[i] / param;
         }
 
         result.polyParams = resultPolyParams;
     }
 
     /**
      * Computes polynomial containing the nth-order integration of current one.
      *
      * @param order  order of integration to compute. Must be at least 1.
      * @param result instance where resulting polynomial will be stored.
      * @throws IllegalArgumentException if provided order is less than 1.
      */
     public void nthIntegration(final int order, final Polynomial result) {
         nthIntegration(order, result, null);
     }
 
     /**
      * Computes polynomial containing the nth-order integration of current one.
      * Because infinite polynomials exist with different constant values,
      * constant terms for each integration order can be provided as well.
      *
      * @param order     order of integration to compute. Must be at least 1.
      * @param constants constant terms for each integration order. Must have a
      *                  length equal to order if provided.
      * @throws IllegalArgumentException if provided order is less than 1 or if
      *                                  constants does not have length equal to order.
      */
     public void nthIntegration(final int order, final double[] constants) {
         nthIntegration(order, this, constants);
     }
 
     /**
      * Computes polynomial containing the nth-order integration of current one.
      *
      * @param order order of integration to compute. Must be at least 1.
      */
     public void nthIntegration(final int order) {
         nthIntegration(order, (double[]) null);
     }
 
     /**
      * Computes polynomial containing the nth-order integration of current one.
      * Because infinite polynomials exist with different constant values,
      * constant terms for each integration order can be provided as well.
      *
      * @param order     order of integration to compute. Must be at least 1.
      * @param constants constant terms for each integration order. Must have a
      *                  length equal to order if provided.
      * @return a new polynomial containing the nth-order integration.
      * @throws IllegalArgumentException if provided order is less than 1 or if
      *                                  constants does not have length equal to order.
      */
     public Polynomial nthIntegrationAndReturnNew(final int order, final double[] constants) {
         final var result = new Polynomial();
         nthIntegration(order, result, constants);
         return result;
     }
 
     /**
      * Computes polynomial containing the nth-order integration of current one.
      *
      * @param order order of integration to compute. Must be at least 1.
      * @return a new polynomial containing the nth-order integration.
      * @throws IllegalArgumentException if provided order is less than 1 or if
      *                                  constants does not have length equal to order.
      */
     public Polynomial nthIntegrationAndReturnNew(final int order) {
         return nthIntegrationAndReturnNew(order, null);
     }
 
     /**
      * Computes nth-integration over provided interval.
      *
      * @param startX    start of integration interval.
      * @param endX      end of integration interval.
      * @param order     order of integration. Must be at least 1.
      * @param constants constant terms for each integration order. Must have a
      *                  length equal to order if provided.
      * @return result of integration.
      * @throws IllegalArgumentException if provided order is less than 1 or if
      *                                  constants does not have length equal to order.
      */
     public double nthOrderIntegrateInterval(
             final double startX, final double endX, final int order, final double[] constants) {
         if (order < MIN_ORDER) {
             throw new IllegalArgumentException();
         }
         if (constants != null && constants.length != order) {
             throw new IllegalArgumentException();
         }
 
         var resultStart = 0.0;
         var resultEnd = 0.0;
         var powStartX = 1.0;
         var powEndX = 1.0;
         double polyParam;
         for (var i = 0; i < order; i++) {
             if (constants != null) {
                 var param = 1;
                 for (var k = 1; k <= i; k++) {
                     param *= k;
                 }
                 polyParam = constants[i] / param;
                 resultStart += polyParam * powStartX;
                 resultEnd += polyParam * powEndX;
             }
             powStartX *= startX;
             powEndX *= endX;
         }
 
         for (int i = 0, j = order; i < polyParams.length; i++, j++) {
             var param = j;
             for (var k = 1; k < order; k++) {
                 param *= j - k;
             }
             polyParam = polyParams[i] / param;
             resultStart += polyParam * powStartX;
             powStartX *= startX;
 
             resultEnd += polyParam * powEndX;
             powEndX *= endX;
         }
 
         return resultEnd - resultStart;
     }
 
     /**
      * Computes nth-integration over provided interval.
      *
      * @param startX start of integration interval.
      * @param endX   end of integration interval.
      * @param order  order of integration. Must be at least 1.
      * @return result of integration.
      * @throws IllegalArgumentException if provided order is less than 1.
      */
     public double nthOrderIntegrateInterval(final double startX, final double endX, final int order) {
         return nthOrderIntegrateInterval(startX, endX, order, null);
     }
 
     /**
      * Trims polynomial to remove all terms above degree that can be neglected.
      *
      * @param result instance where result will be stored.
      */
     public void trim(final Polynomial result) {
         final var degree = getDegree();
         final var resultLength = degree + 1;
 
         final double[] resultPolyParams;
         if (result.polyParams.length != resultLength) {
             resultPolyParams = new double[resultLength];
         } else {
             resultPolyParams = result.polyParams;
         }
         System.arraycopy(polyParams, 0, resultPolyParams, 0, resultLength);
 
         result.polyParams = resultPolyParams;
     }
 
     /**
      * Trims this polynomial to remove all terms above degree that can be
      * neglected.
      */
     public void trim() {
         trim(this);
     }
 
     /**
      * Trims this polynomial to remove all terms above degree that can be
      * neglected and returns the result as a new polynomial.
      *
      * @return a new trimmed polynomial.
      */
     public Polynomial trimAndReturnNew() {
         final var result = new Polynomial();
         trim(result);
         return result;
     }
 
     /**
      * Normalizes parameters of this polynomial so that the array of parameters
      * has unitary norm and stores result into provided instance.
      * Normalization keeps location of real roots, but other roots or
      * properties of polynomials might change.
      *
      * @param result instance where normalized polynomial will be stored.
      */
     public void normalize(final Polynomial result) {
         var resultPolyParams = result.polyParams;
         if (resultPolyParams.length != polyParams.length) {
             resultPolyParams = new double[polyParams.length];
         }
         ArrayUtils.normalize(polyParams, resultPolyParams);
         result.polyParams = resultPolyParams;
     }
 
     /**
      * Normalizes this polynomial so that the array of parameters has unitary
      * norm.
      * Normalization keeps location of real roots, but other roots or
      * properties of polynomials might change.
      */
     public void normalize() {
         normalize(this);
     }
 
     /**
      * Normalizes parameters of this polynomial so that the array of parameters
      * has unitary norm and returns result as a new polynomial instance.
      * Normalization keeps location of real roots, but other roots or
      * properties of polynomials might change.
      *
      * @return a new normalized polynomial instance.
      */
     public Polynomial normalizeAndReturnNew() {
         final var result = new Polynomial(polyParams.length);
         normalize(result);
         return result;
     }
 
     /**
      * Normalizes parameters of this polynomial so that the highest degree term
      * becomes 1.0 and stores result into provided instance.
      *
      * @param result instance where result of normalization will be stored.
      */
     public void normalizeHighestDegreeTerm(final Polynomial result) {
         final var degree = getDegree();
         final var term = polyParams[degree];
         var resultPolyParams = result.polyParams;
         if (resultPolyParams.length != polyParams.length) {
             resultPolyParams = new double[polyParams.length];
         }
         ArrayUtils.multiplyByScalar(polyParams, 1.0 / term, resultPolyParams);
         result.polyParams = resultPolyParams;
     }
 
     /**
      * Normalizes parameters of this polynomial so that the highest degree term
      * becomes 1.0.
      */
     public void normalizeHighestDegreeTerm() {
         normalizeHighestDegreeTerm(this);
     }
 
     /**
      * Normalizes parameters of this polynomial so that the highest degree term
      * becomes 1.0 and returns the result as a new instance.
      *
      * @return a new normalized polynomial.
      */
     public Polynomial normalizeHighestDegreeTermAndReturnNew() {
         final var result = new Polynomial(polyParams.length);
         normalizeHighestDegreeTerm(result);
         return result;
     }
 
     /**
      * Gets location of maxima in this polynomial.
      *
      * @return location of maxima or null if polynomial has no maxima.
      * @throws NumericalException if maxima cannot be determined due to
      *                            numerical instabilities.
      */
     public double[] getMaxima() throws NumericalException {
         return getMaxima(EPS);
     }
 
     /**
      * Gets location of maxima in this polynomial.
      *
      * @param threshold threshold to allow possible small deviations in first
      *                  derivative respect to pure real roots. This should be a very small
      *                  positive value.
      * @return location of maxima or null if polynomial has no maxima.
      * @throws NumericalException       if maxima cannot be determined due to
      *                                  numerical instabilities.
      * @throws IllegalArgumentException if provided threshold is negative.
      */
     @SuppressWarnings("Duplicates")
     public double[] getMaxima(final double threshold) throws NumericalException {
         if (threshold < 0.0) {
             throw new IllegalArgumentException();
         }
 
         final var derivative = derivativeAndReturnNew();
 
         // roots of derivative contains either minima or maxima.
         final var derivativeRoots = derivative.getRoots();
         final var maxima = new ArrayList<Complex>();
         if (derivativeRoots != null) {
             for (var derivativeRoot : derivativeRoots) {
                 if (Math.abs(derivativeRoot.getImaginary()) > threshold) {
                     // root is imaginary (not allowed)
                     continue;
                 }
 
                 final var x = derivativeRoot.getReal();
                 final var secondDerivativeEval = evaluateSecondDerivative(x);
                 if (secondDerivativeEval < 0.0) {
                     // is maxima
                     maxima.add(derivativeRoot);
                 }
             }
         }
 
         // return real parts of maxima, since we only allow real roots of first
         // derivative
         if (maxima.isEmpty()) {
             return null;
         }
 
         final var result = new double[maxima.size()];
         int i = 0;
         for (final var m : maxima) {
             result[i] = m.getReal();
             i++;
         }
 
         return result;
     }
 
     /**
      * Gets location of minima in this polynomial.
      *
      * @return location of minima or null if polynomial has no minima.
      * @throws NumericalException if minima cannot be determined due to
      *                            numerical instabilities.
      */
     public double[] getMinima() throws NumericalException {
         return getMinima(EPS);
     }
 
     /**
      * Gets location of minima in this polynomial.
      *
      * @param threshold threshold to allow possible small deviations in first
      *                  derivative respect to pure real roots. This should be a very small
      *                  positive value.
      * @return location of minima or null if polynomial has no minima.
      * @throws NumericalException       if minima cannot be determined due to
      *                                  numerical instabilities.
      * @throws IllegalArgumentException if provided threshold is negative.
      */
     @SuppressWarnings("Duplicates")
     public double[] getMinima(double threshold) throws NumericalException {
         if (threshold < 0.0) {
             throw new IllegalArgumentException();
         }
 
         final var derivative = derivativeAndReturnNew();
 
         // roots of derivative contains either minima or maxima.
         final var derivativeRoots = derivative.getRoots();
         final var minima = new ArrayList<Complex>();
         if (derivativeRoots != null) {
             for (final var derivativeRoot : derivativeRoots) {
                 if (Math.abs(derivativeRoot.getImaginary()) > threshold) {
                     //root is imaginary (not allowed)
                     continue;
                 }
 
                 final var x = derivativeRoot.getReal();
                 final var secondDerivativeEval = evaluateSecondDerivative(x);
                 if (secondDerivativeEval >= 0.0) {
                     // is minima
                     minima.add(derivativeRoot);
                 }
             }
         }
 
         // return real parts of minima, since we only allow real roots of first
         // derivative
         if (minima.isEmpty()) {
             return null;
         }
 
         final var result = new double[minima.size()];
         var i = 0;
         for (final var m : minima) {
             result[i] = m.getReal();
             i++;
         }
 
         return result;
     }
 
     /**
      * Gets location of minima or maxima (i.e. extrema) in this polynomial.
      *
      * @return location of minima or maxima, or null if polynomial has no
      * minima or maxima.
      * @throws NumericalException if minima or maxima cannot be determined due
      *                            to numerical instabilities.
      */
     public double[] getExtrema() throws NumericalException {
         return getExtrema(EPS);
     }
 
     /**
      * Gets location of minima or maxima (i.e. extrema) in this polynomial.
      *
      * @param threshold threshold to allow possible small deviations in first
      *                  derivative respect to pure real roots. This should be a very small
      *                  positive value.
      * @return location of minima or maxima, or null if polynomial has no minima
      * or maxima.
      * @throws NumericalException       if minima or maxima cannot be determined due
      *                                  to numerical instabilities.
      * @throws IllegalArgumentException if provided threshold is negative.
      */
     @SuppressWarnings("DuplicatedCode")
     public double[] getExtrema(final double threshold) throws NumericalException {
         if (threshold < 0.0) {
             throw new IllegalArgumentException("threshold must be positive");
         }
 
         final var derivative = derivativeAndReturnNew();
 
         // roots of derivative contains either minima or maxima.
         final var derivativeRoots = derivative.getRoots();
         final var minimaOrMaxima = new ArrayList<Complex>();
         if (derivativeRoots != null) {
             for (final var derivativeRoot : derivativeRoots) {
                 if (Math.abs(derivativeRoot.getImaginary()) > threshold) {
                     // root is imaginary (not allowed)
                     continue;
                 }
 
                 minimaOrMaxima.add(derivativeRoot);
             }
         }
 
         // return real parts of roots, since we only allow real roots of first
         // derivative
         if (minimaOrMaxima.isEmpty()) {
             return null;
         }
 
         final var result = new double[minimaOrMaxima.size()];
         int i = 0;
         for (final var m : minimaOrMaxima) {
             result[i] = m.getReal();
             i++;
         }
 
         return result;
     }
 }