// ************
// polynomial.h
// ************
// 
// Oct 15, 2009: removal of the second line of comment on putRealCoeffs.

#ifndef POLYNOMIAL_H
#define POLYNOMIAL_H

#include <cmath>
#include <iostream>
#include "complexx.h"

using namespace std; 

// r = a + b i is considered as a root of a polynomial p(x) if |p(r)| < epsilon 
const double epsilon = 0.000000001;  

// a root r = a + b i is considered 
//    1) a real root if |b| < eta and |p(a)| < epsilon; 
//    2) an imaginary root otherwise if |a| < eta and |p(b)| < epsilon.  
const double eta = 0.000001;   

// terminate Muller's method if root estimates change by less than this amount. 
const double TOL = 0.00001;    

class polynomial   
{
    // overloaded input & output (10 pts) 
    friend istream& operator>>(istream& istr, polynomial& p);       // (6 pts) determine degree from input
    friend ostream& operator<<(ostream& ostr, const polynomial& p); // (4 pts)


public:

    // constructors & destructor (21 pts) 
    polynomial();                        // (2 pts) default constructor initializes a zero polynomial
    polynomial(int deg, complex *coefs); // (4 pts) constructor over an array of coefficients in increasing degress. 
    polynomial(int deg, double *coefs);  // (5 pts) constructor over an array of real coefficients
    polynomial(const polynomial& p);     // (5 pts) copy constructor
    ~polynomial();                       // (5 pts) 
	

    // utility functions (5 pts)
    int getDegree() const;               // (1 pts) degree of the polynomial
    bool realPoly() const;               // (1 pts) true if all coefficients are real 
    void putCoeffs(complex *coefs) const;      // (1 pts) store all coefficients in the array coefs.
    void putRealCoeffs(double *rCoefs) const;  // (1 pts) store in the array rCoeffs if all real coefficients; 
    void update(int deg, complex coef);  // (1 pts) update the coefficient of a polynomial term with the given degree


    // arithmetic operators (30 pts)
    polynomial operator+(const polynomial& p) const;  // (4 pts) addition
    polynomial operator-(const polynomial& p) const;  // (4 pts) subtraction
    polynomial operator*(const polynomial& p) const;  // (8 pts) multiplication
    polynomial operator/(const polynomial& p) const;  // (12 pts) division;
    polynomial operator%(const polynomial& p) const;  // (2 pts) remainder;


    // assignment (5 pts) 
    polynomial& operator=(const polynomial& p);       


    // evaluation and differentiation (10 pts) 
	// if evalDeriv == false, evaluate the polynomial at x = t only. 
	// otherwise, also evaluate its derivative at x = t. 
    complex eval(complex t, bool evalDeriv = false, complex* deriv = 0) const;  


    // root display (4 pts) 
    void printRoots();                    // (1 pt) output all roots. 
    void printRealRoots();                // (1 pt) output real roots only. 
    void putRoots(complex *rts);          // (1 pt) store the roots in the array rts[]. 
    int putRealRoots(double *rts);        // (1 pt) return the number k of real roots on a return 
    // from the function call, only rts[0], ..., rts[k-1] are filled with roots.  

private: 
    int degree;         // degree of the polynomial 
    bool realCoeffs;    // true if all coefficients are real
    double *rCoeffs;    // store coefficients if realCoeffs == true
    complex *Coeffs;    // store coefficients if realCoeffs == false

    bool rootsFound;    // true if all roots have already been found
    complex *roots;     // dynamic array to store the roots
    bool *realRoot;     // realRoot[i] is true when roots[i] is a real root
    int numRealRoots;   // number of real roots 


    // root finding (65 pts) 
    void findRoots();         // (10 pts) find all roots and store them in roots[].

    // closed-form roots (degree <= 4).
    void linearRoots();       // (3 pts)
    void quadraticRoots();    // (5 pts) 
    void cubicRoots();        // (12 pts) 
    void quarticRoots();      // (10 pts) 

    // numerical roots (degree > 4) 
    complex muller() const;   // (11 pts) Muller's method to find one root 
    void newton(complex& r0); // (10 pts) Newton's method to polish a root estimate r0.
    void polishRoots();       // (4 pts) use Newton's method to polish all roots. 
    
    // helper functions
    
    polynomial deflate(complex root);
    void destruct();
    void copy(const polynomial &p);
    
    string getDegreeString(int degree) const;
};


// overload the input operator
istream& operator>>(istream& istr, polynomial& p);        

// overload the output operator
ostream& operator<<(ostream& ostr, const polynomial& p); 


#endif