//=============================================================
//  Numerics Library
//=============================================================
//  Polynomial.h
//
//  Template for polynomials
//
//  Copyright 1995-1996 by Scott Robert Ladd.
//=============================================================

#ifndef COYOTE_POLYNOMIAL_H
#define COYOTE_POLYNOMIAL_H

#include "stddef.h"
#include "math.h"
#include "float.h"

#include "xcomplex.h"
#include "diagnose.h"

//using namespace Coyote;

//namespace Coyote {

//---------------------------
// Polynomial exception types
//---------------------------

enum PolyError
    {
    PX_ALLOC,
    PX_TOOSMALL,
    PX_RANGE,
    PX_INCOMPAT,
    PX_OVERFLOW
    };

class PolyEx : public ExceptionBase
	{
	public:
        PolyEx
            (
            PolyError err
            )
            {
            Error = err;
            }

		PolyError WhatsWrong()
            {
            return Error;
            }

		virtual void Explain
            (
            DiagOutput & out
            );

	private:
		PolyError Error;
	};

//---------------------------
// Polynomial class
//---------------------------

template <class D>
    class Polynomial
        {
        public:
            // constructor
            Polynomial
                (
                size_t terms,
                const D * array = NULL
                );

            // copy constructor
            Polynomial
                (
                const Polynomial<D> & poly
                );

            // destructor
            ~Polynomial();

            // assignment operator
            void operator =
                (
                const Polynomial<D> & poly
                );

            void operator =
                (
                const D * array
                );

            // conversion to array
            operator D * () const;

            // increase degree
            Polynomial<D> Stretch
                (
                size_t newN
                ) const;

            // interrogation
            size_t Degree() const;

            // get coefficients
            D Get
                (
                size_t term
                ) const;

            D & operator []
                (
                size_t term
                );

            // evaluate a polynomial for given value
            D operator ()
                (
                D x
                ) const;

            // operations
            Polynomial<D> operator - () const;
            Polynomial<D> operator + () const;

            Polynomial<D> operator +
                (
                const Polynomial<D> & poly
                ) const;

            void operator +=
                (
                const Polynomial<D> & poly
                );

            Polynomial<D> operator -
                (
                const Polynomial<D> & poly
                ) const;

            void operator -=
                (
                const Polynomial<D> & poly
                );

            Polynomial<D> operator *
                (
                const Polynomial<D> & poly
                ) const;

        protected:
            D * Coeff;
            size_t N;
        };

// constructor
template <class D>
    Polynomial<D>::Polynomial
        (
        size_t terms,
        const D * array = NULL
        )
        {
        if (terms < 2)
            throw PolyEx(PX_TOOSMALL);

        N = terms;

        Coeff = new D [N];

        if (Coeff == NULL)
            throw PolyEx(PX_ALLOC);

        if (array == NULL)
            {
            for (size_t i = 0; i < N; ++i)
                Coeff[i] = (D)0;
            }
        else
            {
            for (size_t i = 0; i < N; ++i)
                Coeff[i] = array[i];
            }
        }

// copy constructor
template <class D>
    Polynomial<D>::Polynomial
        (
        const Polynomial<D> & poly
        )
        {
        N     = poly.N;
        Coeff = new D [N];

        if (Coeff == NULL)
            throw PolyEx(PX_ALLOC);

        for (size_t i = 0; i < N; ++i)
            Coeff[i] = poly.Coeff[i];
        }

// destructor
template <class D>
    Polynomial<D>::~Polynomial()
        {
        delete [] Coeff;
        }

// assignment operator
template <class D>
    void Polynomial<D>::operator =
        (
        const Polynomial<D> & poly
        )
        {
        if (N != poly.N)
            {
            delete [] Coeff;

            N = poly.N;

            Coeff = new D [N];

            if (Coeff == NULL)
                throw PolyEx(PX_ALLOC);
            }

        for (size_t i = 0; i < N; ++i)
            Coeff[i] = poly.Coeff[i];
        }

template <class D>
    void Polynomial<D>::operator =
        (
        const D * array
        )
        {
        if (array == NULL)
            return;

        for (size_t i = 0; i < N; ++i)
            Coeff[i] = array[i];
        }

// conversion to array
template <class D>
    Polynomial<D>::operator D * () const
        {
        D * result = new D [N];

        for (size_t i = 0; i < N; ++i)
            result[i] = Coeff[i];

        return result;
        }

// increase degree
template <class D>
    Polynomial<D> Polynomial<D>::Stretch
        (
        size_t newN
        )
        const
        {
        if (newN <= N)
            return (*this);

        Polynomial<D> result(newN);

        for (size_t i = 0; i < N; ++i)
            result.Coeff[i] = Coeff[i];

        return result;
        }

// interrogation
template <class D>
    inline size_t Polynomial<D>::Degree() const
        {
        return N;
        }

// get coefficients
template <class D>
    D Polynomial<D>::Get
        (
        size_t term
        )
        const
        {
        if (term >= N)
            throw PolyEx(PX_RANGE);

        return Coeff[term];
        }

template <class D>
    D & Polynomial<D>::operator []
        (
        size_t term
        )
        {
        if (term >= N)
            throw PolyEx(PX_RANGE);

        return Coeff[term];
        }

// evaluate a polynomial for given value
template <class D>
    D Polynomial<D>::operator ()
        (
        D x
        )
        const
        {
        D y = Coeff[N - 1];
        size_t i = N - 2;

        while(1)
            {
            y = x * y + Coeff[i];

            if (i == 0)
                break;

            --i;
            }

        return y;
        }

// operations
template <class D>
    Polynomial<D> Polynomial<D>::operator - () const
        {
        Polynomial<D> result(N);

        for (size_t i = 0; i < N; ++ i)
            result.Coeff[i] = -Coeff[i];

        return result;
        }

template <class D>
    inline Polynomial<D> Polynomial<D>::operator + () const
        {
        return *this;
        }

template <class D>
    Polynomial<D> Polynomial<D>::operator +
        (
        const Polynomial<D> & poly
        )
        const
        {
        if (N > poly.N)
            {
            Polynomial<D> res1(*this);

            for (size_t i = 0; i < poly.N; ++i)
                res1.Coeff[i] += poly.Coeff[i];

            return res1;
            }
        else
            {
            Polynomial<D> res2(poly);

            for (size_t i = 0; i < N; ++i)
                res2.Coeff[i] += Coeff[i];

            return res2;
            }
        }

template <class D>
    inline void Polynomial<D>::operator +=
        (
        const Polynomial<D> & poly
        )
        {
        *this = (*this) + poly;
        }

template <class D>
    Polynomial<D> Polynomial<D>::operator -
        (
        const Polynomial<D> & poly
        )
        const
        {
        if (N > poly.N)
            {
            Polynomial<D> res1(*this);

            for (size_t i = 0; i < poly.N; ++i)
                res1.Coeff[i] -= poly.Coeff[i];

            return res1;
            }
        else
            {
            Polynomial<D> res2(-poly);

            for (size_t i = 0; i < N; ++i)
                res2.Coeff[i] += Coeff[i];

            return res2;
            }
        }

template <class D>
    inline void Polynomial<D>::operator -=
        (
        const Polynomial<D> & poly
        )
        {
        *this = (*this) - poly;
        }

template <class D>
    Polynomial<D> Polynomial<D>::operator *
        (
        const Polynomial<D> & poly
        )
        const
        {
        if (N != poly.N)
            throw PolyEx(PX_INCOMPAT);

        Polynomial<D> result(2 * N - 1);

        for (size_t i = 0; i < N; ++i)
            {
            for (size_t j = 0; j < N; ++j)
                result[i + j] += Coeff[i] * poly.Coeff[j];
            }

        return result;
        }

//-------------
// declarations
//-------------

typedef Polynomial<float>       FPoly;
typedef Polynomial<double>      DPoly;
typedef Polynomial<long double> LPoly;
typedef Polynomial<FComplex>    CFPoly;
typedef Polynomial<DComplex>    CDPoly;
typedef Polynomial<LComplex>    CLPoly;

DPoly FFTMultiply
    (
    const DPoly & p1,
    const DPoly & p2
    );

// } // end namespace Coyote

#endif