// Filename: polynomial.cpp
// Author:   Mike Rust and Ian Weiner
// Modified: Michael P. Schubmehl
// Date(s):  07/07/2000
// Purpose:    Implements polynomial.hh, the interface to the polynomial class
//           used by dataSet (see dataset.hh) to interpolate between values in
//           the dataSet. Polynomials support integration and differentiation,
//           and can be of arbitrary degree.
//
// Notes:      Implements polynomial.hh, and requires iostream.h for stream-
//           handling routines. Uses stdlib.h for the NULL constant, and math.h
//           for the pow(x,y) routine.

#include <iostream.h>                         // For << operator
#include <stdlib.h>                           // For NULL constant
#include <math.h>                             // For powering
#include "polynomial.hh"                      // Header file for this class

// Table of Contents **********************************************************
// --public--
//   polynomial();                            // Constructors
//   polynomial(int degree_);
//   polynomial(int degree_, double* coefficients_);
//
//   double evaluate(double x);               // Returns value of polynomial,
//                                            //   evaluated at x.
//   double evaluateDerivative(double x);     // Returns the value of the
//                                            //   derivative at x.
//   double evaluateIntegral(double x0,       // Returns the integral from
//                           double x1);      //   x0 to x1.
//
//
//   polynomial& operator=(const polynomial& source); // Assignment operator
//
//   void set(int degree_,                    // Set new degree, coefficients
//             double* coefficients_);
//
// other:
// ostream& operator<<(ostream& stream,       // Operator to output polynomial
//                     const polynomial& source);
// END Table of Contents ******************************************************


// class polynomial ***********************************************************
polynomial::polynomial() {
  // Input:  None.
  // Output: Returns a polynomial with no degree or coefficients.

  degree = -1;
  coefficients = NULL;
}

polynomial::polynomial(int degree_) {
  // Input:  An integer degree.
  // Output: Returns a polynomial of specified degree, with zero coefficients.

  if (degree_ >= 0) {                         // Check for a valid degree
    degree = degree_;
    coefficients = new double[degree+1];
    for(int i = 0; i < degree + 1; i++) {     // Initialize coefficients
      coefficients[i] = 0;
    }
  }
  else {
    cerr << "Error: Invalid degree for a polynomial." << endl;
    exit(1);                                  // Terminate on error
  }
}

polynomial::polynomial(int degree_, double* coefficients_) {
  // Input:  An integer degree and an array of double coefficients that has
  //         length degree + 1.
  // Output: Returns a polynomial with specified degree and coefficients.

  if (degree_ >= 0) {                         // Check for a valid degree
    degree = degree_;
    coefficients = new double[degree+1];
    for(int i = 0; i < degree + 1; i++) {     // Set up coefficients
      coefficients[i] = coefficients_[i];
    }
  }
  else {
    cerr << "Error: Invalid degree for a polynomial." << endl;
    exit(1);                                  // Terminate on error
  }
}

polynomial::~polynomial() {
  // Input:  None.
  // Output: Deletes the coefficient array and destroys the polynomial.

  delete[] coefficients;
}

double polynomial::evaluate(double x) {
  // Input:  A value at which to evaluate the polynomial.
  // Output: The value of the polynomial at x.

  double result = 0;
  for(int i = 0; i < degree + 1; i++) {       // Add the ith term (a_i * x^i)
    result += coefficients[i] * pow(x, i);    //   to the final result.
  }
  return result;
}

double polynomial::evaluateDerivative(double x) {
  // Input:  A value at which to evaluate the derivative of the polynomial.
  // Output: Returns the derivative of the polynomial at x.

  double result = 0;
  for(int i = 1; i < degree + 1; i++) {       // Add ith term of derivative
    result += i * coefficients[i] * pow(x, i-1);
  }
  return result;
}

double polynomial::evaluateIntegral(double x0, double x1) {
  // Input:  Values x0 and x1 between which to integrate the polynomial.
  // Output: Returns the integral of the polynomial between x0 and x1.

  double resultx0 = 0;                        // Integral from 0 to x0
  double resultx1 = 0;                        // Integral from 0 to x1
  for(int i = 0; i < degree + 1; i++) {       // Add ith terms of integrals
    resultx0 += 1/(double)(i+1) * coefficients[i] * pow(x0, i+1);
    resultx1 += 1/(double)(i+1) * coefficients[i] * pow(x1, i+1);
  }
  return resultx1 - resultx0;                 // Return integral from x0 to x1
}

polynomial& polynomial::operator=(const polynomial& source) {
  // Input:  A polynomial to be copied into this one.
  // Output: Sets the values of this polynomial to match the source.

  delete[] coefficients;                      // Free up memory used by old
                                              //   coefficients.
  degree = source.degree;                     // Copy degree
  coefficients = new double[degree+1];        // Allocate coefficient space
  for(int i = 0; i < degree + 1; i++) {
    coefficients[i] = source.coefficients[i]; // Copy coefficients
  }

  return *this;                               // Return self, with new values
}

void polynomial::set(int degree_, double* coefficients_) {
  // Input:  A degree and new coefficients for this polynomial.
  // Output: Sets this polynomial to have degree and coefficients as given.

  delete[] coefficients;                      // Free old coefficients

  degree = degree_;                           // Set degree
  coefficients = new double[degree+1];        // Allocate for coefficients
  for(int i = 0; i < degree + 1; i++) {       // Set coefficients
    coefficients[i] = coefficients_[i];
  }
}
// END class polynomial *******************************************************

ostream& operator<<(ostream& stream, const polynomial& source) {
  // Input:  An ostream and a polynomial.
  // Output: Outputs that polynomial to the ostream, and returns the stream.

  for(int i = source.degree; i > 1; i--) {
    stream << source.coefficients[i] << "x^" << i << " + ";
  }

  if (source.degree > 0) {
    stream << source.coefficients[1] << "x + ";
  }

  stream << source.coefficients[0];

  return stream;
}