/*
 * Open Source Physics software is free software as described near the bottom of this code file.
 *
 * For additional information and documentation on Open Source Physics please see:
 * <http://www.opensourcephysics.org/>
 */

package org.opensourcephysics.numerics;

/**
 * Polynomial implements a mathematical polynomial:
 * c[0] + c[1] * x + c[2] * x^2 + ....
 *
 * This class is based on code published in Object-Oriented Implementation of
 * Numerical Methods by Didier H. Besset.
 * The code has been adapted to the OSP framework by Wolfgang Christian.
 */
public class Polynomial implements Function {
  /**
   * Polynomial coefficients.
   */
  protected double[] coefficients;

  /**
   * Constructs a polynomial with the given coefficients.
   * @param coef polynomial coefficients.
   */
  public Polynomial(double[] coef) {
    coefficients = coef;
  }

  /**
   * Gets a clone of the polynomial coefficients c:
   * c[0] + c[1] * x + c[2] * x^2 + ....
   *
   * @return double[]
   */
  public double[] getCoefficients() {
    return coefficients.clone();
  }

  /**
   * Constructs a polynomial with the given coefficients.
   * @param coef polynomial coefficients.
   */
  public Polynomial(String[] coef) {
    coefficients = new double[coef.length];
    for(int i = 0, n = coef.length; i<n; i++) {
      try {
        coefficients[i] = Double.parseDouble(coef[i]);
      } catch(NumberFormatException ex) {
        coefficients[i] = 0;
      }
    }
  }

  /**
   * Evaluates a polynomial using the given coefficients.
   *
   * @param x
   * @param coeff the polynomial coefficients.
   */
  public static double evalPolynomial(final double x, final double[] coeff) {
    int n = coeff.length-1;
    double y = coeff[n];
    for(int i = n-1; i>=0; i--) {
      y = coeff[i]+(y*x);
    }
    return y;
  }

  /**
   *
   * @param r double    number added to the polynomial.
   * @return Polynomial
   */
  public Polynomial add(double r) {
    int n = coefficients.length;
    double coef[] = new double[n];
    coef[0] = coefficients[0]+r;
    for(int i = 1; i<n; i++) {
      coef[i] = coefficients[i];
    }
    return new Polynomial(coef);
  }

  /**
   * Adds the given polynomial to this polynomial.
   * @param p Polynomial
   * @return Polynomial
   */
  public Polynomial add(Polynomial p) {
    int n = Math.max(p.degree(), degree())+1;
    double[] coef = new double[n];
    for(int i = 0; i<n; i++) {
      coef[i] = coefficient(i)+p.coefficient(i);
    }
    return new Polynomial(coef);
  }

  /**
   * Gets the coefficient value at the desired position
   * @param n int    the position of the coefficient to be returned
   * @return double the coefficient value
   * @version 1.2
   */
  public double coefficient(int n) {
    return(n<coefficients.length) ? coefficients[n] : 0;
  }

  /**
   * Deflates the polynomial by removing the root.
   *
   * @param r double     a root of the polynomial (no check made).
   * @return Polynomial the receiver divided by polynomial (x - r).
   */
  public Polynomial deflate(double r) {
    int n = degree();
    double remainder = coefficients[n];
    double[] coef = new double[n];
    for(int k = n-1; k>=0; k--) {
      coef[k] = remainder;
      remainder = remainder*r+coefficients[k];
    }
    return new Polynomial(coef);
  }

  /**
   * Gets the degree of this polynomial function.
   * @return int degree of this polynomial function
   */
  public int degree() {
    return coefficients.length-1;
  }

  /**
   * Gets the derivative of this polynomial.
   * @return Polynomial the derivative.
   */
  public Polynomial derivative() {
    int n = degree();
    if(n==0) {
      double coef[] = {0};
      return new Polynomial(coef);
    }
    double coef[] = new double[n];
    for(int i = 1; i<=n; i++) {
      coef[i-1] = coefficients[i]*i;
    }
    return new Polynomial(coef);
  }

  /**
   * Divides this polynomial by a constant.
   * @param r double
   * @return Polynomial
   */
  public Polynomial divide(double r) {
    return multiply(1/r);
  }

  /**
   * Divides this polynomial by another polynomial.
   *
   * The remainder is dropped.
   *
   * @param p Polynomial
   * @return Polynomial
   */
  public Polynomial divide(Polynomial p) {
    return divideWithRemainder(p)[0];
  }

  /**
   * Divides this polynomial by another polynomial.
   *
   * @param p polynomial
   * @return polynomial array containing the answer and remainder
   */
  public Polynomial[] divideWithRemainder(Polynomial p) {
    Polynomial[] answer = new Polynomial[2];
    int m = degree();
    int n = p.degree();
    if(m<n) {
      double[] q = {0};
      answer[0] = new Polynomial(q);
      answer[1] = p;
      return answer;
    }
    double[] quotient = new double[m-n+1];
    double[] coef = new double[m+1];
    for(int k = 0; k<=m; k++) {
      coef[k] = coefficients[k];
    }
    double norm = 1/p.coefficient(n);
    for(int k = m-n; k>=0; k--) {
      quotient[k] = coef[n+k]*norm;
      for(int j = n+k-1; j>=k; j--) {
        coef[j] -= quotient[k]*p.coefficient(j-k);
      }
    }
    double[] remainder = new double[n];
    for(int k = 0; k<n; k++) {
      remainder[k] = coef[k];
    }
    answer[0] = new Polynomial(quotient);
    answer[1] = new Polynomial(remainder);
    return answer;
  }

  /**
   * Integrates this polynomial.
   * The integral has the value 0 at x = 0.
   *
   * @return Polynomial the integral
   */
  public Polynomial integral() {
    return integral(0);
  }

  /**
   * Integrates this polynomial having the specified value for x = 0.
   * @param value double    value of the integral at x=0
   * @return Polynomial the integral.
   */
  public Polynomial integral(double value) {
    int n = coefficients.length+1;
    double coef[] = new double[n];
    coef[0] = value;
    for(int i = 1; i<n; i++) {
      coef[i] = coefficients[i-1]/i;
    }
    return new Polynomial(coef);
  }

  /**
   * Multiplies this polynomial by a constant.
   * @param r double
   * @return Polynomial
   */
  public Polynomial multiply(double r) {
    int n = coefficients.length;
    double coef[] = new double[n];
    for(int i = 0; i<n; i++) {
      coef[i] = coefficients[i]*r;
    }
    return new Polynomial(coef);
  }

  /**
   * Multiplies this polynomial by another polynomial.
   * @param p Polynomial
   * @return Polynomial
   */
  public Polynomial multiply(Polynomial p) {
    int n = p.degree()+degree();
    double[] coef = new double[n+1];
    for(int i = 0; i<=n; i++) {
      coef[i] = 0;
      for(int k = 0; k<=i; k++) {
        coef[i] += p.coefficient(k)*coefficient(i-k);
      }
    }
    return new Polynomial(coef);
  }

  /**
   * Gets the complex roots of this polynomial.
   * @return double[]
   */
  public double[][] roots() {
    int highest = coefficients.length-1; // start with the degree
    for(int i = highest; i>0; i--) {
      if(coefficients[highest]!=0) {
        break;
      }
      highest = i;
    }
    if(highest==0)return new double[2][0];
    double ch = coefficients[highest]; //highest nonzero coefficient
    double[][] companion = new double[highest][highest];
    companion[0][highest-1] = -coefficients[0]/ch;
    for(int i = 0; i<highest-1; i++) {
      companion[0][i] = -coefficients[highest-i-1]/ch;
      companion[i+1][i] = 1;
    }
    EigenvalueDecomposition eigen = new EigenvalueDecomposition(companion);
    double[][] roots = new double[2][];
    roots[0] = eigen.getRealEigenvalues();
    roots[1] = eigen.getImagEigenvalues();
    return roots;
  }

  /**
   * Gets the real roots of this polynomial.
   * @return double[]
   */
  public double[] rootsReal() {
    double[][] roots = roots();
    int n = roots[0].length;
    double[] temp = new double[n];
    int count = 0;
    for(int i = 0; i<n; i++) {
      double magSquared = roots[0][i]*roots[0][i]+roots[1][i]*roots[1][i];
      if(roots[1][i]*roots[1][i]/magSquared<1.0e-12) { // skip small imaginary values
        temp[count] = roots[0][i];
        count++;
      }
    }
    double[] re = new double[count];
    System.arraycopy(temp, 0, re, 0, count);
    return re;
  }

  /**
   * Subtracts a constant from this polynomial.
   *
   * @param r the constant
   * @return Polynomial
   */
  public Polynomial subtract(double r) {
    return add(-r);
  }

  /**
   * Subtracts another polynomial from this polynomial.
   * @return Polynomial
   * @param p Polynomial
   */
  public Polynomial subtract(Polynomial p) {
    int n = Math.max(p.degree(), degree())+1;
    double[] coef = new double[n];
    for(int i = 0; i<n; i++) {
      coef[i] = coefficient(i)-p.coefficient(i);
    }
    return new Polynomial(coef);
  }

  /**
   * Converts this polynomial to a String.
   */
  public String toString() {
    if((coefficients==null)||(coefficients.length<1)) {
      return "Polynomial coefficients are undefined."; //$NON-NLS-1$
    }
    StringBuffer sb = new StringBuffer();
    boolean firstNonZeroCoefficientPrinted = false;
    for(int n = 0, m = coefficients.length; n<m; n++) {
      if(coefficients[n]!=0) {
        if(firstNonZeroCoefficientPrinted) {
          sb.append((coefficients[n]>0) ? " + " : " "); //$NON-NLS-1$ //$NON-NLS-2$
        } else {
          firstNonZeroCoefficientPrinted = true;
        }
        if((n==0)||(coefficients[n]!=1)) {
          sb.append(Double.toString(coefficients[n]));
        }
        if(n>0) {
          sb.append(" x^"+n);                           //$NON-NLS-1$
        }
      }
    }
    String str = sb.toString();
    if(str.equals("")) { //$NON-NLS-1$
      return "0";        //$NON-NLS-1$
    }
    return str;
  }

  /**
   * Evaluates the polynomial for the specified variable value.
   * @param x double    value at which the polynomial is evaluated
   * @return double polynomial value.
   */
  public double evaluate(double x) {
    int n = coefficients.length;
    double answer = coefficients[--n];
    while(n>0) {
      answer = answer*x+coefficients[--n];
    }
    return answer;
  }

  /**
   * Returns the value and the derivative of this polynomial
   * for the specified variable value in an array of two elements
   *
   * @param x double    value at which the polynomial is evaluated
   * @return double[0]   the value of the polynomial
   * @return double[1]   the derivative of the polynomial
   */
  public double[] valueAndDerivative(double x) {
    int n = coefficients.length;
    double[] answer = new double[2];
    answer[0] = coefficients[--n];
    answer[1] = 0;
    while(n>0) {
      answer[1] = answer[1]*x+answer[0];
      answer[0] = answer[0]*x+coefficients[--n];
    }
    return answer;
  }

}

/*
 * Open Source Physics software is free software; you can redistribute
 * it and/or modify it under the terms of the GNU General Public License (GPL) as
 * published by the Free Software Foundation; either version 2 of the License,
 * or(at your option) any later version.

 * Code that uses any portion of the code in the org.opensourcephysics package
 * or any subpackage (subdirectory) of this package must must also be be released
 * under the GNU GPL license.
 *
 * This software is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 USA
 * or view the license online at http://www.gnu.org/copyleft/gpl.html
 *
 * Copyright (c) 2007  The Open Source Physics project
 *                     http://www.opensourcephysics.org
 */
