//Author: Courtesy of  Wayne Rowcliffe
// polynomial.cpp 

#include <cmath>
#include <cstdlib>
#include <iostream>
#include <string>
#include <sstream>
#include "complexx.h"
#include "polynomial.h"

using namespace std; 


// helper declarations
	polynomial deflate(double*,int,complex);
	
	double unifRand();


// constructors & destructor (21 pts) 
	
	// (2 pts) default constructor initializes a zero polynomial
	polynomial::polynomial() : 
	degree(0), realCoeffs(true), Coeffs(0), rootsFound(false), roots(0), realRoot(0), numRealRoots(0) {
		rCoeffs = new double[1];
		rCoeffs[0] = 0;
		
	}
	
	// (4 pts) constructor over an array of coefficients in increasing degress.                 
	polynomial::polynomial(int deg, complex *coefs) : 
	degree(deg), realCoeffs(false), rCoeffs(0), rootsFound(false), roots(0), realRoot(0), numRealRoots(0) {
		Coeffs = new complex[deg+1];
		for(int i = 0; i <= deg; i++) {
			Coeffs[i] = coefs[i];
		}
	}  
	
	// (5 pts) constructor over an array of real coefficients
	polynomial::polynomial(int deg, double *coefs) :
	degree(deg), realCoeffs(true), Coeffs(0), rootsFound(false), roots(0), realRoot(0), numRealRoots(0) {
		rCoeffs = new double[deg+1];
		for(int i = 0; i <= deg; i++) {
			rCoeffs[i] = coefs[i];
		}
	}  
	
	// (5 pts) copy constructor
	polynomial::polynomial(const polynomial& p) {
		copy(p);
	}     
	
	// (5 pts)
	polynomial::~polynomial() {
		destruct();
	}                        


// utility functions (5 pts)

	// (1 pts) degree of the polynomial
	int polynomial::getDegree() const {
		return degree;
	}        
	
	// (1 pts) true if all coefficients are real              
	bool polynomial::realPoly() const {
		return realCoeffs;
	}                    
	
	 // (1 pts) store all coefficients in the array coefs. 
	void polynomial::putCoeffs(complex *coefs) const {
		if(!realCoeffs) {
			for(int i = 0; i <= degree; i++) {
				coefs[i] = Coeffs[i];
			}
		} else {
			for(int i = 0; i <= degree; i++) {
				coefs[i] = rCoeffs[i];
			}
		}
	}     
	
	// (1 pts) store in the array rCoeffs if all real coefficients;
	//         update realCoeffs
	void polynomial::putRealCoeffs(double *rCoefs) const {
		if(realCoeffs) {
			for(int i =0; i <= degree; i++) {
				rCoefs[i] = rCoeffs[i];
			}
		}
	}   
	
	
	// (1 pts) update the coefficient of a polynomial term with the given degree
	void polynomial::update(int deg, complex coef) {
		// Changin a coefficient means the roots are invalid.
		if(rootsFound) {
			rootsFound = false;
			delete [] roots;
			delete [] realRoot;
		}
		// If the coefficients are real and the updated coefficient is real, use real coeffs
		if(realCoeffs && (coef.imagPart() < eta)) {
			rCoeffs[deg] = coef.realPart();
		// Otherwise copy everything over to complex coefficients
		} else {
			if(realCoeffs) {
				Coeffs = new complex[degree+1];
				for(int i = 0; i <= degree; i++) {
					Coeffs[i] = rCoeffs[i];
				}
				realCoeffs = false;
				delete [] rCoeffs;
				rCoeffs = 0;
			}
			Coeffs[deg] = coef;
		}
	}  


// arithmetic operators (30 pts)

	// (4 pts) addition
	polynomial polynomial::operator+(const polynomial& p) const {
		int maxDeg = (degree > p.degree) ? degree : p.degree;
		bool real = (realCoeffs && p.realCoeffs);
		polynomial p2;
		// Decide which values to use
		if(real) {
			// Set up
			double *temp = new double[maxDeg+1];
			for(int i = 0; i <= degree; i++) {
				temp[i] = rCoeffs[i];
			}
			
			// Do the math
			for(int i = 0; i <= p.degree; i++) {
				temp[i] += p.rCoeffs[i];
			}
			
			// Reduce degree until highest degree is non-zero
			while((maxDeg > 0) && abs(temp[maxDeg]) < eta) {
				maxDeg--;
			}
			p2 = polynomial(maxDeg,temp);
			delete [] temp;
		// Same as above except complex
		} else {
			complex *lhs = 0;
			complex *rhs = 0;
			if(realCoeffs) {
				lhs = new complex[degree +1];
				for(int i = 0; i <= degree; i++) {
					lhs[i] = complex(rCoeffs[i],0);
				}
			} else {
				lhs = Coeffs;
			}
			if(p.realCoeffs) {
				rhs = new complex[p.degree +1];
				for(int i = 0; i <= p.degree; i++) {
					rhs[i] = complex(p.rCoeffs[i],0);
				}
			} else {
				rhs = p.Coeffs;
			}
			complex *temp = new complex[maxDeg +1];
			for(int i = 0; i <= degree; i++) {
				temp[i] = lhs[i];
			}
			for(int i = 0; i <= p.degree; i++) {
				temp[i] = temp[i] + rhs[i];
			}
			
			while((maxDeg > 0) && abs(temp[maxDeg].magnitude()) < eta) {
				maxDeg--;
			}
			p2 = polynomial(maxDeg, temp);
			
			delete [] temp;
			if(realCoeffs) {
				delete [] lhs;
			}
			if(p.realCoeffs) {
				delete [] rhs;
			}
		}
		return p2;
	}  
	
	// (4 pts) subtraction
	polynomial polynomial::operator-(const polynomial& p) const {
		int maxDeg = (degree > p.degree) ? degree : p.degree;
		bool real = (realCoeffs && p.realCoeffs);
		polynomial p2;
		// decide which values to use
		if(real) {
			// Set up
			double *temp = new double[maxDeg+1];
			for(int i = 0; i <= degree; i++) {
				temp[i] = rCoeffs[i];
			}
			
			// Do the math
			for(int i = 0; i <= p.degree; i++) {
				temp[i] -= p.rCoeffs[i];
			}
			
			// Reduce degree until highest degree is non-zero
			while((maxDeg > 0) && abs(temp[maxDeg]) < eta) {
				maxDeg--;
			}
			p2 = polynomial(maxDeg,temp);
			delete [] temp;
		// Same as above except complex
		} else {
			complex *lhs = 0;
			complex *rhs = 0;
			if(realCoeffs) {
				lhs = new complex[degree +1];
				for(int i = 0; i <= degree; i++) {
					lhs[i] = complex(rCoeffs[i],0);
				}
			} else {
				lhs = Coeffs;
			}
			if(p.realCoeffs) {
				rhs = new complex[p.degree +1];
				for(int i = 0; i <= p.degree; i++) {
					rhs[i] = complex(p.rCoeffs[i],0);
				}
			} else {
				rhs = p.Coeffs;
			}
			complex *temp = new complex[maxDeg +1];
			for(int i = 0; i <= degree; i++) {
				temp[i] = lhs[i];
			}
			for(int i = 0; i <= p.degree; i++) {
				temp[i] = temp[i] - rhs[i];
			}
			
			while((maxDeg > 0) && abs(temp[maxDeg].magnitude()) < eta) {
				maxDeg--;
			}
			p2 = polynomial(maxDeg, temp);
			
			delete [] temp;
			if(realCoeffs) {
				delete [] lhs;
			}
			if(p.realCoeffs) {
				delete [] rhs;
			}
		}
		return p2;
	}  
	
	// (8 pts) multiplication
	polynomial polynomial::operator*(const polynomial& p) const {
		int maxDeg = degree + p.degree;
		bool real = (realCoeffs && p.realCoeffs);
		polynomial p2;
		// Decide which method to use
		if(real) {
			// Set up the array
			double *temp = new double[maxDeg +1];
			for(int i = 0; i <= maxDeg; i++) {
				temp[i] = 0;
			}
			
			// Do the math
			for(int i = 0; i <= degree; i++) {
				for(int j = 0; j <= p.degree; j++) {
					temp[i+j] += rCoeffs[i] * p.rCoeffs[j];
				}
			}
			
			// Reduce degree until the highest degree is non-zero
			while((maxDeg > 0) && abs(temp[maxDeg]) < eta) {
				maxDeg--;
			}
			p2 = polynomial(maxDeg, temp);
			delete [] temp;
		// Same as above except complex
		} else {
			complex *lhs = 0;
			complex *rhs = 0;
			if(realCoeffs) {
				lhs = new complex[degree +1];
				for(int i = 0; i <= degree; i++) {
					lhs[i] = complex(rCoeffs[i],0);
				}
			} else {
				lhs = Coeffs;
			}
			if(p.realCoeffs) {
				rhs = new complex[p.degree +1];
				for(int i = 0; i <= p.degree; i++) {
					rhs[i] = complex(p.rCoeffs[i],0);
				}
			} else {
				rhs = p.Coeffs;
			}
			complex *temp = new complex[maxDeg +1];
			for(int i = 0; i <= maxDeg; i++) {
				temp[i] = complex(0,0);
			}
			for(int i = 0; i <= degree; i++) {
				for(int j = 0; j <= p.degree; j++) {
					temp[i + j] = temp[i +j] + (lhs[i] * rhs[j]);
				}
			}
			
			while((maxDeg > 0) && abs(temp[maxDeg].magnitude()) < eta) {
				maxDeg--;
			}
			p2 = polynomial(maxDeg, temp);
			
			delete [] temp;
			if(realCoeffs) {
				delete [] lhs;
			}
			if(p.realCoeffs) {
				delete [] rhs;
			}
		}
		return p2;
	}  
	
	// (12 pts) division
	polynomial polynomial::operator/(const polynomial& p) const {
		int maxDeg = degree - p.degree;
		bool real = (realCoeffs && p.realCoeffs);
		polynomial p2;
		if(maxDeg >= 0) {
			if(real) {
				// Create the arrays
				double *temp = new double[maxDeg +1];
				for(int i = 0; i <= maxDeg; i++) {
					temp[i] = 0;
				}
				double *numerator = new double[degree+1];
				for(int i = 0; i <= degree; i++) {
					numerator[i] = rCoeffs[i];
				}
				
				// Do the division
				for(int i = degree; i >= p.degree; i--) {
					double factor = numerator[i]/p.rCoeffs[p.degree];
					temp[i-p.degree] = factor;
					for(int j = 0; j <= p.degree; j++) {
						numerator[i+j-p.degree] -= factor * p.rCoeffs[j];
					}
				}
				
				// Count down degrees until we find one that is non-zero
				while((maxDeg > 0) && abs(temp[maxDeg]) < eta) {
					maxDeg--;
				}
				p2 = polynomial(maxDeg, temp);
				
				delete [] temp;
				delete [] numerator;
			// Same as above except complex
			} else {
				complex *numerator = 0;
				complex *rhs = 0;
				if(realCoeffs) {
					numerator = new complex[degree +1];
					for(int i = 0; i <= degree; i++) {
						numerator[i] = complex(rCoeffs[i],0);
					}
				} else {
					numerator = new complex[degree + 1];
					for(int i = 0; i <= degree; i++) {
						numerator[i] = Coeffs[i];
					}
				}
				if(p.realCoeffs) {
					rhs = new complex[p.degree +1];
					for(int i = 0; i <= p.degree; i++) {
						rhs[i] = complex(p.rCoeffs[i],0);
					}
				} else {
					rhs = p.Coeffs;
				}
				
				complex *temp = new complex[maxDeg +1];
				for(int i = 0; i <= maxDeg; i++) {
					temp[i] = complex(0,0);
				}
				
				for(int i = degree; i >= p.degree; i--) {
					complex factor = numerator[i]/rhs[p.degree];
					temp[i-p.degree] = factor;
					for(int j = 0; j <= p.degree; j++) {
						numerator[i+j-p.degree] = numerator[i+j-p.degree] - (factor * rhs[j]);
					}
				}
								
				while((maxDeg > 0) && abs(temp[maxDeg].magnitude()) < eta) {
					maxDeg--;
				}
				p2 = polynomial(maxDeg, temp);
				
				delete [] temp;
				delete [] numerator;
				if(p.realCoeffs) {
					delete [] rhs;
				}
			}
		}
		return p2;	
	}  
	
	// (2 pts) remainder
	polynomial polynomial::operator%(const polynomial& p) const {
		return *this - (*this/p) * p;
	}  


// assignment (5 pts) 
	polynomial& polynomial::operator=(const polynomial& p) {
		if(&p == this) {
			return *this;
		}
		destruct();
		copy(p);
		return *this;
	}       


// evaluation and differentiation (10 pts) 
// if evalDeriv == false, evaluate the polynomial at x = t only. 
// otherwise, also evaluate its derivative at x = t. 
	complex polynomial::eval(complex t, bool evalDeriv, complex* deriv) const {
		complex a;
		complex d;
		// Choose which array to use, then do horner scheme
		if(realCoeffs) {
			a = rCoeffs[degree];
			for(int i = degree-1; i >= 0; i--) {
				if(evalDeriv) {
					d = t*d + a;
				}
				a = t*a + rCoeffs[i];
			}
			
		} else {
			a = Coeffs[degree];
			for(int i = degree-1; i >= 0; i--) {
				if(evalDeriv) {
					d = t*d + a;
				}
				a = Coeffs[i] + t*a;
			}
		}
		if(evalDeriv) {
			*deriv = d;
		}
		return a;
	}  


// root display (4 pts) 

	// (1 pt) output all roots. 
	void polynomial::printRoots() {
		findRoots();
		if(rootsFound) {
			for(int i = 0; i < degree; i++) {
				cout << "(" << roots[i] << ") , ";
			}
			cout << endl;
		}
	}              
	
	// (1 pt) output real roots only. 
	void polynomial::printRealRoots() {
		findRoots();
		if(rootsFound) {
			for(int i = 0; i < degree; i++) {
				if(realRoot[i]) {
					cout << "(" << roots[i] << ") , ";
				}
			}
			cout << endl;
		}
		
	}          
	
	// (1 pt) store the roots in the array rts[]. 
	void polynomial::putRoots(complex *rts) {
		findRoots();
		if(rootsFound) {
			for(int i = 0; i < degree; i++) {
				rts[i] = roots[i];
			}
		}
	}    
	
	// (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.
	int polynomial::putRealRoots(double *rts) {
		int realCount = 0;
		findRoots();
		if(rootsFound) {
			for(int i = 0; i < degree; i++) {
				if(realRoot[i]) {
					rts[realCount++] = roots[i].realPart();
				}
			}
		}
		return realCount;
	}  
	  


// root finding (65 pts) 
	
	// (10 pts) find all roots and store them in roots[].
	void polynomial::findRoots() {
		// If we can find a closed form, please do.
		if(realCoeffs && !rootsFound) {
			switch(degree) {
				case 1 : linearRoots(); break;
				case 2 : quadraticRoots(); break;
				case 3 : cubicRoots(); break;
				case 4 : quarticRoots(); break;
			}
		}
		// Otherwise use mullers to find the roots
		if(realCoeffs && !rootsFound) {
			roots = new complex[degree];
			int rootCount = 0;
			
			// Find a root
			complex root = muller();
			
			// List it					
			roots[rootCount++] = root;
			if(abs(root.imagPart()) > eta) {
				roots[rootCount++] = complex(root.realPart(),-root.imagPart());
			}
			
			// Deflate then find roots of deflated polynomial
			if(rootCount < degree) {
				polynomial p = deflate(root);
				p.findRoots();
			
				// Copy the deflated polynomials roots into the roots array
				for(int i = 0; i < p.degree; i++) {
					roots[rootCount++] = p.roots[i];
				}
			}
			
			// Polish the roots			
			polishRoots();
			
			// Decide which roots are real
			realRoot = new bool[degree];
			for(int i = 0; i < degree; i++) {
				if(abs(roots[i].imagPart()) < eta) {
					realRoot[i] = true;
				} else { 
					realRoot[i] = false;
				}
			}
			rootsFound = true;
		}
	}

// closed-form roots (degree <= 4).

	// (3 pts)
	void polynomial::linearRoots() {
		// According to the assignment handout
		if(degree == 1) {
			rootsFound = true;
			roots = new complex[1];
			realRoot = new bool[1];
			roots[0] = -rCoeffs[0]/rCoeffs[1];
			realRoot[0] = true;
		}
	}       
	
	// (5 pts) 
	void polynomial::quadraticRoots() {
		// According to the assignment handout
		if(degree == 2) {
			rootsFound = true;
			roots = new complex[2];
			realRoot = new bool[2];
			double a = rCoeffs[2];
			double b = rCoeffs[1];
			double c = rCoeffs[0];
			double det = (b*b) - (4*a*c);
			if(det < 0) {
				roots[0] = (complex(-b) + complex(det).sqrt())/ (2*a);
				roots[1] = (complex(-b) - complex(det).sqrt())/ (2*a);
				realRoot[0] = false;
				realRoot[1] = false;
			} else {
				roots[0] = (-b + sqrt(det))/ (2*a);
				roots[1] = (-b - sqrt(det))/ (2*a);
				realRoot[0] = true;
				realRoot[1] = true;
			}
		}
	}    
	
	// (12 pts) 
	void polynomial::cubicRoots() {
		// according to the assignment handout
		if(degree == 3) {
			
			rootsFound = true;
			roots = new complex[3];
			realRoot = new bool[3];
			
						
			double p = rCoeffs[2]/rCoeffs[3];
			double q = rCoeffs[1]/rCoeffs[3];
			double r = rCoeffs[0]/rCoeffs[3];
			
					
			double a = ((3*q)-(p*p))/3;
			double b = (((2*p*p*p) - (9*p*q) + (27*r))/27.0);
			
			complex inner = complex(pow(b,2)/4 + pow(a,3)/27).sqrt();
			
			complex aA = (complex(-b/2)+inner);
			
			if(abs(aA.imagPart()) < epsilon) {
				if(aA.realPart() > 0) {
					aA = complex(pow(aA.realPart(),1/3.0));
				} else {
					aA = complex(-pow(abs(aA.realPart()),1/3.0));
				}
			} else {
				aA = complex(pow(aA.magnitude(),1.0/3) * cos((1.0/3) * aA.phase()), pow(aA.magnitude(),1.0/3) * sin((1.0/3)*aA.phase()));
			}
			
			complex bB = (complex(-b/2)-inner);
			
			if(abs(bB.imagPart()) < epsilon) {
				if(bB.realPart() > 0) {
					bB = complex(pow(bB.realPart(),1/3.0));
				} else {
					bB = complex(-pow(abs(bB.realPart()),1/3.0));
				}
			} else {
				bB = complex(pow(bB.magnitude(),1.0/3) * cos((1.0/3) * bB.phase()), pow(bB.magnitude(),1.0/3) * sin((1.0/3)*bB.phase()));
			}
			
			roots[0] = (aA + bB + complex(-p/3));
			roots[1] = (complex(-.5)*(aA + bB) + complex(0,sqrt(3)/2.0)*(aA - bB) + complex(-p/3));
			roots[2] = (complex(-.5)*(aA + bB) - complex(0,sqrt(3)/2.0)*(aA - bB) + complex(-p/3));
			
			for(int i = 0; i < 3; i++) {
				if(abs(roots[i].imagPart()) < epsilon) {
					realRoot[i] = true;
				} else { 
					realRoot[i] = false;
				}
			}
		}
	}        
	
	// (10 pts) 
	void polynomial::quarticRoots() {
		// according to the assignment handout
		if(degree == 4) {
			double p = rCoeffs[3]/rCoeffs[4];
			double q = rCoeffs[2]/rCoeffs[4];
			double r = rCoeffs[1]/rCoeffs[4];
			double s = rCoeffs[0]/rCoeffs[4];
			
			double rcA[] = {(4*q*s) + -(r*r) + -(p*p*s), p*r - 4*s, -q, 1};
						
			polynomial rc(3,rcA);
			
			rc.cubicRoots();
						
			complex z;
			for(int i = 0; i < 3; i++) {
				if(rc.realRoot[i]) {
					z = rc.roots[i].realPart();
					break;
				} else if(i == 2) {
					return;
				}
			}
			
			if((complex(p*p/4 - q) + z).magnitude() < epsilon) {
				return;
			}
			
			complex rR = (complex(p*p/4 - q) + z).sqrt();
					
			
			complex first = complex(.75*p*p) - rR.pow(2) - complex(2*q);
			complex second;
			if(rR.magnitude() < epsilon) {
				second = (z.pow(2) - 4*s).sqrt() * 2;
			} else {
				second = (complex(4*p*q - 8*r) - pow(p,3))/ (rR * 4);
			}
			
			complex dD = (first + second).sqrt();
			complex eE = (first - second).sqrt();
			
			rootsFound = true;
			roots = new complex[4];
			realRoot = new bool[4];
			
			roots[0] = complex(-p/4) + (rR + dD)/2;
			roots[1] = complex(-p/4) + (rR - dD)/2;
			roots[2] = complex(-p/4) - (rR - eE)/2;
			roots[3] = complex(-p/4) - (rR + eE)/2;
			
			for(int i = 0; i < 4; i++) {
				if(abs(roots[i].imagPart()) < epsilon) {
					realRoot[i] = true;
				} else { 
					realRoot[i] = false;
				}
			}
		}
	}      

// numerical roots (degree > 4) 

	// (11 pts) Muller's method to find one root 
	complex polynomial::muller() const {
		double rad = 0;
		// Find the radius to search within
		if(abs(rCoeffs[1]) > epsilon) {
			rad = min(degree*abs(rCoeffs[0]/rCoeffs[1]),pow(abs(rCoeffs[0]/rCoeffs[degree]),1.0/degree));
		} else {
			rad = pow(abs(rCoeffs[0]/rCoeffs[degree]),1.0/degree);
		}
		// Set them up in an array so we can just cycle through the values.
		complex xk[] = {complex(unifRand()*rad,unifRand()*rad), 
						complex(unifRand()*rad,unifRand()*rad), 
						complex(unifRand()*rad,unifRand()*rad)};
		int curpos = 0;
		while(eval(xk[curpos]).magnitude() > epsilon) {
			complex &x2 = (curpos < 2) ? xk[curpos+1] : xk[0]; 
			complex &x1 = (curpos < 1) ? xk[curpos+2] : xk[curpos-1];
			complex &x0 = xk[curpos];
			if((x0 - x1).magnitude() < TOL) { break;}
			// Mullers method as given by the assignment handout
			complex p1 = (eval(x2) - eval(x1))/(x2-x1);
			complex p2 = (eval(x1) - eval(x0))/(x1-x0);
			complex p3 = (p1 - p2)/(x2- x0);
			complex &a = p3;
			complex b = p2 + p3*(x0-x1);
			complex pk = eval(x0);
			complex det = (b.pow(2) - (a * pk * 4)).sqrt();
			(curpos == 2) ? curpos = 0 : curpos++;
			xk[curpos] = x0 - (pk * 2)/(((b+det).magnitude() > (b-det).magnitude()) ? (b + det) : (b - det));
		}
		
		return xk[curpos];
	}   
	
	// (10 pts) Newton's method to polish a root estimate r0.
	void polynomial::newton(complex& r0) {
		while(eval(r0).magnitude() > epsilon) {
			complex pp;
			complex p = eval(r0,true, &pp);
			r0 = r0 - (p/ pp);
			if((p/pp).magnitude() < epsilon) {break;} 
		}
	}
	
	// (4 pts) use Newton's method to polish all roots. 
	void polynomial::polishRoots() {
		for(int i = 0; i < degree; i++) {
			newton(roots[i]);
		}
	}       



// overloaded input & output (10 pts) 

	// (6 pts) determine degree from input
	istream& operator>>(istream& istr, polynomial& p) {
		// Setup
		//double r(0), i(0);
		stringbuf line;
		istr.get(line);
		stringstream poly;
		poly.str(line.str());
		int deg = 0;
		bool comp = false;
		string temp;
		stringstream tempStream;
		int tempDeg = 0;
		while(!poly.eof()) {
			poly >> temp;
			// If i appears, the polynomial is complex
			if(!(string::npos == temp.find('i'))) {
				comp = true;
			}
			// Find the maximum degree
			if((temp.length() > 0) && (temp[0] == 'x')) {
				if(temp.length() > 2) {
					tempStream.clear();
					tempStream.str(temp.substr(2));
					tempStream >> tempDeg;
					if(tempDeg > deg) {
						deg = tempDeg;
					}
				} else if(deg == 0) {
					deg = 1;
				}
			}
		}
		// Reset our stringstream
		poly.clear();
		poly.seekg(0);
	
		int maxDeg = deg;
		bool neg = false;	
		if(comp) {
			complex degVal = 1;
			complex* coefs = new complex[deg +1];
			while(!poly.eof()) {
				poly >> temp;
				// Finding a -, evaluate the previous degree, indicate that the next will be negative
				if((temp[0] == '-') && (temp.length() == 1)) {
					coefs[deg] = (neg) ? -degVal : degVal;
					deg = 0;
					degVal = 1;
					neg = true;
				// Finding a +, evaluate the previous degree, indicate that the next will be positive
				} else if((temp[0] == '+') && (temp.length() == 1)) {
					coefs[deg] = (neg) ? -degVal : degVal;
					deg = 0;
					degVal = 1;
					neg = false;
				// Find x, determine the next degree.
				} else if(temp[0] == 'x') {
					if(temp.length() > 1) {
						tempStream.clear();
						tempStream.str(temp.substr(2));
						tempStream >> deg;
					} else {
						deg = 1;
					}
				// Finding a complex value
				} else if(temp[0] == '(') {
					tempStream.clear();
					tempStream.str(temp.substr(1,temp.length() -1));
					tempStream >> degVal;
				// Finding a real or imaginary value
				} else {
					tempStream.clear();
					tempStream.str(temp);
					tempStream >> degVal;
				}
			}
			// Evaluate the final value.
			coefs[deg] = (neg) ? -degVal : degVal;
			p = polynomial(maxDeg,coefs);
			delete [] coefs;
		} else {
			double degVal = 1;
			double* coefs = new double[deg +1];
			while(!poly.eof()) {
				poly >> temp;
				// Finding a -, evaluate the previous degree, indicate that the next will be negative
				if((temp[0] == '-') && (temp.length() == 1)) {
					coefs[deg] = (neg) ? -degVal : degVal;
					deg = 0;
					degVal = 1;
					neg = true;
				// Finding a +, evaluate the previous degree, indicate that the next will be positive
				} else if((temp[0] == '+') && (temp.length() == 1)) {
					coefs[deg] = (neg) ? -degVal : degVal;
					deg = 0;
					degVal = 1;
					neg = false;
				// Find x, determine the next degree.
				} else if(temp[0] == 'x') {
					if(temp.length() > 1) {
						tempStream.clear();
						tempStream.str(temp.substr(2));
						tempStream >> deg;
					} else {
						deg = 1;
					}
				// In case there are parenthesis
				} else if(temp[0] == '(') {
					tempStream.clear();
					tempStream.str(temp.substr(1,temp.length() -1));
					tempStream >> degVal;
				// Evaluate the degree coefficient
				} else {
					tempStream.clear();
					tempStream.str(temp);
					tempStream >> degVal;
				}
			}
			coefs[deg] = (neg) ? -degVal : degVal;
			p = polynomial(maxDeg,coefs);
			delete [] coefs;
		}
		
		return istr;
	}       
	
	// (4 pts)
	ostream& operator<<(ostream& ostr, const polynomial& p) {
		// Choose the proper array to read from, then write out the 
		// Value given by getDegreeString() for each degree that is non-zero
		if(p.realCoeffs) {
			for(int i = p.degree; i >= 0 ; i--) {
				if(abs(p.rCoeffs[i]) > epsilon) {
					ostr << p.getDegreeString(i) << ' ';
				}
			}
			if(p.degree == 0 && abs(p.rCoeffs[0]) < epsilon) {
				ostr  << 0;
			}
		} else {
			for(int i = p.degree; i >= 0; i--) {
				if(abs(p.Coeffs[i].magnitude()) > epsilon) {
					ostr << p.getDegreeString(i) << ' ';
				}
			}
			if(p.degree == 0 && abs(p.Coeffs[0].magnitude()) < epsilon) {
				ostr  << 0;
			}
		}
		return ostr;
	} 
	
// helper methods because this code is unruly.

	polynomial polynomial::deflate(complex root) {
		polynomial p2;
		// If its a real root, deflate with horner scheme.
		if(abs(root.imagPart()) < epsilon) {
			double *d = new double[degree];
			double r = root.realPart();
			d[degree-1] = rCoeffs[degree];
			for(int i = degree-2; i >= 0; i--) {
				d[i] = r*d[i+1] + rCoeffs[i+1];
			}
			p2 = polynomial(degree-1,d);
			delete [] d;
		// Otherwise combine with the complex conjugate to form a real polynomial and divide it out.
		} else {
			double d[] = {(pow(root.realPart(),2) + pow(root.imagPart(),2)),-2*root.realPart(),1};
			
			polynomial p3(2,d);
			
			p2 = *this/p3;
		}
		return p2;
	}
	
	void polynomial::destruct() {
		// Delete things if they exist
		if(realCoeffs) {
			delete [] rCoeffs;
			rCoeffs = 0;
		} else {
			delete [] Coeffs;
			Coeffs = 0;
		}
		if(rootsFound) {
			delete [] roots;
			roots = 0;
			delete [] realRoot;
			realRoot = 0;
		}
	}
	
	void polynomial::copy(const polynomial &p) {
		// constants
		degree = p.degree;
		realCoeffs = p.realCoeffs;
		rootsFound = p.rootsFound;
		numRealRoots = p.numRealRoots;
		// create array to hold the coefficients
		if(realCoeffs) {
			Coeffs = 0;
			rCoeffs = new double[degree +1];
			for(int i = 0; i <= degree; i++) {
				rCoeffs[i] = p.rCoeffs[i];
			}
		} else {
			rCoeffs = 0;
			Coeffs = new complex[degree +1];
			for(int i = 0; i <= degree; i++) {
				Coeffs[i] = p.Coeffs[i];
			}
		}
		// Store the roots if we have them
		if(rootsFound) {
			roots = new complex[degree];
			realRoot = new bool[degree];
			for(int i = 0; i < degree; i++) {
				roots[i] = p.roots[i];
				realRoot[i] = p.realRoot[i];
			}
		} else {
			roots = 0;
			realRoot = 0;
		}
	}
	
	// Determine how the value of a given degree should be displayed
	string polynomial::getDegreeString(int deg) const {
		stringstream s;
		complex val;
		if(realCoeffs) {
			val = rCoeffs[deg];
		} else {
			val = Coeffs[deg];
		}
		// If the complex number is 0 in one of its parts only print the non-zero part
		bool comp = (abs(val.realPart()) > eta) && (abs(val.imagPart()) > eta);
		if(!comp) {
			double rVal = (abs(val.realPart()) > eta) ? val.realPart() : val.imagPart();
						
			if((deg < degree) && (rVal >= 0)) {
				s << '+';
			} else if(rVal < 0) {
				s << '-';
			}
			
			if(deg < degree) {
				s << ' ';
			}
			
			if(!(abs(rVal) == 1) || (deg == 0)) {
				s << abs(rVal);
				if(abs(val.imagPart()) > eta) {
					s << "i ";
				} else {
					s << " ";
				}
			}
		// If both a and b in a + bi are non-zero, print in parenthesis
		} else {
			if(val.realPart() >= 0) {
				if(deg < degree) {
					s << "+ ";
				}				
				s << "(" << val << ") ";
			} else {
				s << '-';
				if(deg < degree) {
					s << ' ';
				}
				s << "(" << -val << ") ";
			}
		}
		
		s << ((deg == 0) ? "" : ((deg == 1) ? "x" : "x^"));
		if(deg > 1) {
				s << deg;
		}
		
		return s.str();
	}
	
	// Random number between -1 and 1
	double unifRand() {
    	return rand() / double(RAND_MAX) -.5;
	}