#ifndef _INCLUDED_L3_POLY_H
#define _INCLUDED_L3_POLY_H

// *Polynominal* Copyright (C) Krzysztof Bosak 1999-03-03...1999-03-03.
// All rights reserved.
// kbosak@box43.pl
// http://www.kbosak.prv.pl

#include <string.h>
#include "l3_stack.h"
#include "l3_array.h"
#include "l3_ptr.h"

template<class real_number>
class Polynominal
{	// Coefficients are stored as follows: c0 + c1*x + c2*x^2 +c3*x^3
protected:
	autoarray<real_number> _coeff;

public:
	inline Polynominal()
		: _coeff(1)
	{
		_coeff[0]=0;
	}
	virtual ~Polynominal()
	{
	}
	inline Polynominal(const real_number * const coeff, int degree)
	{
		assert(coeff!=NULL);
		assert(degree>=0);
		assert(coeff[degree]!=0);
		_coeff.setsize(degree+1);
		for(int i=0; i<_coeff.size(); i++)
			_coeff[i]=coeff[i];
	}
	inline explicit Polynominal(const autoarray<real_number>& coeff)
		: _coeff(coeff)
	{
		assert(coeff.size()>=1);
	}
	inline Polynominal<real_number>& operator+=(const Polynominal<real_number>& poly)
	{
		const int coeff_size=(_coeff.size()<=poly._coeff.size() ? _coeff.size() : poly._coeff.size());
		SmartPtr<real_number> ptr(&_coeff[0], &_coeff[0], coeff_size);
		SmartPtrReadOnly<real_number> ptr2(&poly._coeff[0], &poly._coeff[0], coeff_size);
		for(int i=0; i<coeff_size; i++)
		{
			(*ptr)+=(*ptr2);
			++ptr;
			++ptr2;
		}
		return *this;
	}
	inline const Polynominal<real_number> operator+(const Polynominal<real_number>& poly) const
	{
		Polynominal<real_number> temp(*this);
		return temp+=poly;
	}
	inline Polynominal<real_number>& operator-=(const Polynominal<real_number>& poly)
	{
		const int coeff_size=(_coeff.size()<=poly._coeff.size() ? _coeff.size() : poly._coeff.size());
		SmartPtr<real_number> ptr(&_coeff[0], &_coeff[0], coeff_size);
		SmartPtrReadOnly<real_number> ptr2(&poly._coeff[0], &poly._coeff[0], coeff_size);
		for(int i=0; i<coeff_size; i++)
		{
			(*ptr)-=(*ptr2);
			++ptr;
			++ptr2;
		}
		return *this;
	}
	inline Polynominal<real_number> operator-() const
	{
		Polynominal<real_number> ret(*this);
		autoarray<real_number>& array_ref=ret._coeff;
		for(int i=0; i<_coeff.size(); i++)
			array_ref[i]=-array_ref[i];
		return ret;
	}
	inline Polynominal<real_number>& MultiplyByX()
	{
		_coeff.resize(_coeff.size()+1);
		for(int i=_coeff.size()-1; i>=1; i--)
		{
			_coeff[i]=_coeff[i-1];
		}
		_coeff[0]=0;
		return *this;
	}
	inline int Degree() const
	{
		return _coeff.size()-1;
	}
	inline real_number Coeff(int power_index) const
	{
		return _coeff[power_index];
	}
	friend ostream& operator<<(ostream& str, const Polynominal<real_number>& poly)
	{
		const int imax=poly._coeff.size();
		for(int i=0; i<imax; i++)
		{
			str<<poly._coeff[i];
			if(i==1)
			{
				str<<"x";
			}
			else if(i>1)
			{
				str<<"x^"<<i;
			}
			if(i<imax-1)
			{
				str<<" + ";
			}
		}
		return str;
	}
	void Out(int style=0, char * const arg=NULL) const
	{	// Output of whole polynominal
		// 0==ascii output; 1==cpp code; 2=='x' --> arg Polynominal;
		assert(style>=0);
		assert(style<=2);
		assert(arg!=NULL);
		SmartPtrReadOnly<real_number> ptr(&_coeff[0], &_coeff[0], _coeff.size());
		for(int j=0; j<_coeff.size(); j++, ptr++)
		{
			real_number c=*ptr;
			if(c!=0)
			{

				if(c>0)
				{
					cout<<" +";
				}
				else if(c<0)
				{
					cout<<" -";
					c=-c;
				}

				if(c!=1 || j==0)
				{
					cout<<c;
				}
				if(j)
				{
					if(style==1)//cpp code
					{
						if(c!=1)
							cout<<'*';
						if(j!=1)
							cout<<"pow(x,"<<static_cast<char>(j+'0')<<")";
						else
							cout<<'x';
					}
					else if(style==2)//arg Polynominal
					{
						if(c!=1)
							cout<<'*';
						if(j!=1)
							cout<<"pow("<<arg<<"),"<<static_cast<char>(j+'0')<<")";
						else
							cout<<"cos(x)";
					}//ascii output
					else
					{
						cout<<'x';
						if(j!=1)
							cout<<'^'<<static_cast<char>(j+'0');
					}
				}
			}
		}
	}
	inline real_number Value(const real_number& x) const
	{	// Value of polynominal at x
		const int degree=Degree();
		SmartPtrReadOnly<real_number> ptr(&_coeff[degree], &_coeff[0], _coeff.size());
		real_number result=*ptr;
		for(int i=degree; i>0; i--)
		{
			--ptr;
			result*=x;
			result+=(*ptr);
		}
		return result;
	}
	inline real_number ValueDerivative(const real_number& x) const
	{	// Value of derivative of polynominal at x
		const int degree=Degree();
		SmartPtrReadOnly<real_number> ptr(&_coeff[degree], &_coeff[0], _coeff.size());
		real_number resultd=(*ptr)*degree;
		for(int i=degree-1; i>0; i--)
		{
			--ptr;
			resultd*=x;
			resultd+=(*ptr)*i;
		}
		return resultd;
	}
	inline real_number Value_ValueDerivative(const real_number& x) const
	{	// Value of polynominal at x/ derivative of polynominal at x (useful for quick finding roots)
		const int degree=Degree();
		SmartPtrReadOnly<real_number> ptr(&_coeff[degree], &_coeff[0], _coeff.size());
		real_number result=*ptr;
		real_number resultd=*ptr * degree;
		--ptr;
		result*=x;
		result+=(*ptr);
		for(int i=degree-1; i; i--)
		{
			resultd*=x;
			resultd+=(*ptr)*i;
			ptr--;
			result*=x;
			result+=(*ptr);

		}
		assert(resultd!=0);
		if(resultd!=0)
			return result/resultd;
		else
			return 0;
	}
	inline Polynominal<real_number>& operator*=(const real_number& scalar)
	{
		for(int i=0; i<_coeff.size(); i++)
		{
			_coeff[i]*=scalar;
		}
		return *this;
	}
	inline Polynominal<real_number> GetDerivative() const
	{	// Returns first derivative of polynominal
		const int degree=Degree();
		autoarray<real_number> temp(degree);
		SmartPtr<real_number> ptr1(&temp[0], &temp[0], temp.size());
		SmartPtrReadOnly<real_number> ptr2(&_coeff[1], &_coeff[0], _coeff.size());
		for(int d=1, i=0; i<degree; i++, d++)
		{
			(*ptr1)=(*ptr2)*d;
			++ptr1;
			++ptr2;
		}
		return Polynominal<real_number>(&temp[0], degree-1);
	}
	void Derivate()
	{	// Derivates itself
		SmartPtr<real_number> ptr1(&_coeff[0], &_coeff[0], _coeff.size());
		SmartPtr<real_number> ptr2(&_coeff[1], &_coeff[0], _coeff.size());
		const int degree=Degree();
		for(int d=1, i=0; i<degree; i++, d++)
		{
			(*ptr1)=(*ptr2)*d;
			++ptr1;
			++ptr2;
		}
		_coeff.resize(_coeff.size()-1);
	}
};

template<class real_number>
inline Polynominal<real_number> LegendreP(int n, const real_number&)
{//specific definition of Legendre's polynominal (without constant coefficient)
 //(d/dx)^n ((x^2-1)^n)
	assert(n>=0);
	Polynominal<real_number> _result;
	if(n==0)
	{
		static autoarray<real_number> init(1);
		init[0]=1;
		_result=Polynominal<real_number>(init);
		return _result;
	}
	static autoarray<real_number> init(3);
	init[0]=-1;
	init[1]=0;
	init[2]=1;
	_result=Polynominal<real_number>(init);
	Polynominal<real_number> poly2(_result);
	int i;
	for(i=2; i<=n; i++)
		_result*=poly2;
	for(i=1; i<=n; i++)
		_result.Derivate();
	return _result;
}

template<class real_number>
inline Polynominal<real_number> LegendrePAsoc(int n, int m, const real_number&)
{//specific definition of associated Legendre's polynominal (without constant coefficient)
 //(-1)^m * (1-x*x)^m * (d/dx)^m LegendreP(n)
	assert(n>=0);
	assert(m>=0);
	assert(n>=m);
	Polynominal<real_number> _result(LegendreP(n));

	static autoarray<real_number> init(3);
	init[0]=-1;
	init[1]=0;
	init[2]=1;
	Polynominal<real_number> poly2(init);
	int i;
	for(i=1; i<=m; i++)
		_result.Derivate();
	for(i=1; i<=m; i++)
		_result*=poly2;
	if(!(m&1))
		return _result;
	return -_result;
}

class EXW
{//represents a*exp[-x]*x^w ; a, w are integers
	int _a, _w;
public:
	inline EXW(int a=0, int w=0)
	{
		assert(w>=0);
		_a=a;
		_w=w;
	}
	inline EXW d_de() const
	{// -a*exp[-x]*x^w, first term of derivative
		return EXW(-_a, _w);
	}
	inline EXW d_dx() const
	{// w*a*exp[-x]*x^(w-1), second term of derivative
		if(_w)
			return EXW(_w*_a, _w-1);
		else
			return EXW(0, 0);
	}
	inline int A() const
	{
		return _a;
	}
	inline int& A()
	{
		return _a;
	}
	inline int W() const
	{
		return _w;
	}

};

template<class real_number>
inline Polynominal<real_number> LaguerreP(int n, const real_number&)
{//specific definition of Laguerre's polynominal:
 //e^x*[(d/dx)^n (exp[-x]*x^n)]
	assert(n>=0);

	growing_arraystack<EXW> lagu, lagu2;
	lagu.push(EXW(1, n));

	//differentiate
	for(int i=0; i<n; i++)
	{
		while(!lagu.empty())
		{
			EXW exw=lagu.top();
			lagu.pop();
			if(exw.A())
			{
				lagu2.push(exw.d_de());
				if(exw.W()>0)
					lagu2.push(exw.d_dx());
			}
		}

		lagu=lagu2;
		lagu2.free();
	}

	autoarray<real_number> tmp(n+1);
	//copy coefficients to autoarray while simplyfying factors (look at +=)
	while(!lagu.empty())
	{
		EXW exw=lagu.top();
		lagu.pop();
		const real_number a=exw.A();
		if(a!=0)
			tmp[n-exw.W()]+=a;
	}
	return Polynominal<real_number>(tmp);
}

template<class real_number>
inline const Polynominal<real_number> operator*(
	const Polynominal<real_number>& poly,
	const real_number& scalar
	)
{
	Polynominal<real_number> result=poly;
	result.operator*=(scalar);
	return result;
}


#endif // _INCLUDED_L3_POLY_H