#ifndef __POLYNOMIAL_HPP_INCLUDED__
#define __POLYNOMIAL_HPP_INCLUDED__

#include <vector>
#include <cmath>
#include <stdexcept>
#include <ostream>

class Polynome;

bool operator==(const Polynome& p, const Polynome& q);
inline Polynome operator*(const Polynome& p, const Polynome& q);
double puissance(double x, size_t n);

class Polynome
{

public:

// Constructeurs -------------------------------------------------------------------------------

    Polynome(double a = 0.)
            :m_termes(1,a)
    {}

    Polynome(const std::vector<double>& tab)
            :m_termes(tab)
    {
        simplifie();
    }

    template<size_t N>
    Polynome(const double (&tab)[N])
           :m_termes(&tab[0], &tab[N])
    {
        simplifie();
    }

// Accès aux termes ------------------------------------------------------------------------------

    inline const double& operator[](size_t n) const
    {
        return m_termes[n];
    }

    inline double& operator[](size_t n)
    {
        return m_termes[n];
    }

    inline const double& at(size_t n) const
    {
        return m_termes.at(n);
    }

    inline double& at(size_t n)
    {
        return m_termes.at(n);
    }

// Degre -----------------------------------------------------------------------------------

    inline size_t degre() const
    {
        return m_termes.size() - 1;
    }

    inline size_t length() const
    {
        return m_termes.size();
    }

// Evaluation -----------------------------------------------------------------------------------

    double evalue(double x) const
    {
        double resultat = 0.;

        for (size_t i(0); i<m_termes.size(); ++i)
        {
            resultat += m_termes[i] * puissance(x,i);
        }

        return resultat;
    }

    inline double operator()(double& x) const
    {
        return evalue(x);
    }

//Operateurs ----------------------------------------------------------------------------

    Polynome& operator+=(const Polynome& q)
    {
        size_t i = 0;
        for (; i<q.length() && i<m_termes.size(); ++i)
        {
            m_termes[i] += q[i];
        }

        for (; i<q.length() ; ++i)
            m_termes.push_back(q[i]);

        return *this;
    }

    Polynome& operator-=(const Polynome& q)
    {
        size_t i = 0;
        for (; i<q.length() && i<m_termes.size(); ++i)
        {
            m_termes[i] -= q[i];
        }

        for (; i<q.length() ; ++i)
            m_termes.push_back(-q[i]);

        return *this;
    }

    Polynome& operator*=(const Polynome& q)
    {
        if (q == 0. || *this == 0.)
        {
            m_termes.assign(1,0.);
            return *this;
        }

        const size_t len(q.degre() + degre() + 1);
        std::vector<double> vec(len,0.);

        for (size_t k(0); k<len; ++k)
            for (size_t i(0); i< this->length();++i)
            {
                const size_t j = k-i;
                if (j<q.length())
                    vec[k] += m_termes[i]*q[j];
            }

        *this = Polynome(vec);
        return *this;
    }

    inline Polynome& operator/=(const Polynome& q)
    {
        *this = division(*this,q).first;
        return *this;
    }

    inline Polynome& operator%=(const Polynome& q)
    {
        *this = division(*this,q).second;
        return *this;
    }

//Derivée ------------------------------------------------------------------------------

    Polynome derivee() const
    {
        if (degre() == 0)
            return 0;

        Polynome ret(*this);
        for (size_t i(0); i<ret.degre(); ++i)
            ret[i] = ret[i+1] * (i+1);

        ret.m_termes.pop_back();

        return ret;
    }

private:

//Simplification -----------------------------------------------------------------------

    void simplifie()
    {
        size_t counter = 0;
        typedef std::vector<double>::reverse_iterator Iterator;
        for (Iterator it = m_termes.rbegin(); it != m_termes.rend() && *it == 0. ; ++it)
        {
            ++counter;
        }

        m_termes.resize(m_termes.size() - counter);
    }

    typedef std::pair<Polynome, Polynome >  Paire;

    static Paire division(const Polynome& dividende, const Polynome& diviseur)
    {
        if (diviseur == 0.)
            throw std::runtime_error("Polynome: Division par 0");

        if (dividende.degre() < diviseur.degre())
            return Paire(0, dividende);

        if (dividende == diviseur)
            return Paire(1,0);

        Polynome reste(dividende);
        Polynome quotient;

        while (reste.degre() > diviseur.degre())
        {
            std::vector<double> vec(reste.degre() - diviseur.degre() + 1,0);
            vec.back() = reste.m_termes.back() / diviseur.m_termes.back();

            Polynome p(vec);

            reste -= p*diviseur;
            reste.simplifie();

            quotient += p;
        }

        return Paire(quotient,reste);
    }

//Attributs-----------------------------------------------------------------------------

    std::vector<double> m_termes;

//Amities ------------------------------------------------------------------------------

    friend bool operator==(const Polynome& p, const Polynome& q);

};

// Flux -----------------------------------------------------------------------------------

std::ostream& operator<<(std::ostream& flux, const Polynome& poly)
{
    if (poly.degre() == 0 && poly[0] == 0.)
    {
        flux << 0.;
        return flux;
    }

    if (poly[poly.length() - 1] < 0.)
        flux << "- ";

    for (size_t i(1);i<poly.length();++i)
    {
        if (poly[poly.length() - i] != 0.)
            flux << std::abs(poly[poly.length() - i]) << "x^" << poly.length() - i;

        if (poly[poly.length() - (i+1)] > 0.)
            flux << " + ";
        else if (poly[poly.length() - (i+1)] < 0.)
            flux << " - ";
    }

    if (poly[0] != 0.)
        flux << std::abs(poly[0]);

    return flux;
}

//Opérateurs de comparaison ------------------------------------------------------------------------

bool operator==(const Polynome& p, const Polynome& q)
{
    return p.m_termes == q.m_termes;
}

inline bool operator!=(const Polynome& p, const Polynome& q)
{
    return !(p == q);
}

// Opérateurs binaires ---------------------------------------------------------------------------

inline Polynome operator+(const Polynome& p, const Polynome& q)
{
    return Polynome(p) += q;
}

inline Polynome operator-(const Polynome& p, const Polynome& q)
{
    return Polynome(p) -= q;
}

inline Polynome operator*(const Polynome& p, const Polynome& q)
{
    return Polynome(p) *= q;
}

inline Polynome operator/(const Polynome& p, const Polynome& q)
{
    return Polynome(p) /= q;
}

inline Polynome operator%(const Polynome& p, const Polynome& q)
{
    return Polynome(p) %= q;
}

//Opérateurs unaires -------------------------------------------------------------------------------

Polynome operator-(const Polynome& p)
{
    Polynome result(p);
    for (size_t i(0); i<result.length(); ++i)
        result[i] = -result[i];
    return result;
}

Polynome operator+(const Polynome& p)
{
    Polynome result(p);
    for (size_t i(0); i<result.length(); ++i)
        result[i] = +result[i];
    return result;
}

//Opérateurs logiques ----------------------------------------------------------------------------

inline bool operator!(const Polynome& p)
{
    return p == 0.;
}

//Dérivée -----------------------------------------------------------------------------------

inline Polynome derivee(const Polynome& p)
{
    return p.derivee();
}

//Fonction assistante -----------------------------------------------------------------------------

double puissance(double x, size_t n)
{
    double resultat = 1.;

    for (size_t i(0); i<n ; ++i)
        resultat *= x;
    return resultat;
}

#endif //__POLYNOM_HPP_INCLUDED__

#ifndef __POLYNOMIAL_HPP_INCLUDED__
#define __POLYNOMIAL_HPP_INCLUDED__

/**
 * \file Polynom.hpp
 * \author Nanoc
 * \brief A Polynomial class
 * \date 19 April 2009
 * \version 1.1
 *
 * The operators are implemented naively. Could be improved
 * \todo Add a solve() function
 */

#include <vector>
#include <cmath>
#include <stdexcept>
#include <ostream>

/**
 * \class Polynomial<>
 *
 *  \brief A class to represent polynomials of any degree
 *
 *  All the usual operators are overloaded. One can evaluate the polynomial at any point or compute its derivative
 *  Access to all terms is also granted
 *
 */
template<typename T>
class Polynomial
{

public:

// Constructors -------------------------------------------------------------------------------

    /**
     *  \brief Constructor. Buils a Polynomial from a single value
     */
    Polynomial(const T& a = static_cast<T>(0))
            :m_terms(1,a)
    {}

    /**
     *  \brief Constructor. Builds a Polynomial from a vector of terms
     */
    Polynomial(const std::vector<T>& tab)
            :m_terms(tab)
    {
        simplify();
    }

    /**
     *  \brief Constructor. Buils a Polynomial from a "C-like" array
     */
    template<size_t N>
    Polynomial(const T (&tab)[N])
         :m_terms(&tab[0], &tab[N])
    {
       simplify();
    }

// Term access ------------------------------------------------------------------------------

    /**
     *  \brief Access to a term (read only)
     */
    inline const T& operator[](size_t n) const
    {
        return m_terms[n];
    }

    /**
     *  \brief Access to a term
     */
    inline T& operator[](size_t n)
    {
        return m_terms[n];
    }

    /**
     *  \brief Access to a term with rangecheck (read only)
     */
    inline const T& at(size_t n) const
    {
        return m_terms.at(n);
    }

    /**
     *  \brief Access to a term with rangecheck
     */
    inline T& at(size_t n)
    {
        return m_terms.at(n);
    }

// Degree -----------------------------------------------------------------------------------

    /**
     *  \brief Returns the degree of the polynomial
     */
    inline size_t degree() const
    {
        return m_terms.size() - 1;
    }

    /**
     *  \brief Returns the length of the internal array (should be degree()+1)
     */
    inline size_t length() const
    {
        return m_terms.size();
    }

// Evaluation -----------------------------------------------------------------------------------

    /**
     *  \brief Evaluates the polynom at a given point x
     */
    template <typename S>
    S evaluate(const S& x) const
    {
        S result = static_cast<S>(0);

        for (size_t i(0); i<m_terms.size(); ++i)
        {
            result += m_terms[i] * power(x,i);
        }

        return result;
    }

    /**
     *  \brief Evaluates the polynom at a given point x
     */
    template <typename S>
    inline S operator()(const S& x) const
    {
        return evaluate(x);
    }

//Operators ----------------------------------------------------------------------------

    /**
     *  \brief Sum term by term of two polynomials
     */
    Polynomial<T>& operator+=(const Polynomial<T>& q)
    {
        size_t i = 0;
        for (; i<q.length() && i<m_terms.size(); ++i)
        {
            m_terms[i] += q[i];
        }

        for (; i<q.length() ; ++i)
            m_terms.push_back(q[i]);

        return *this;
    }

    /**
     *  \brief Difference term by term of two polynomials
     */
    Polynomial<T>& operator-=(const Polynomial<T>& q)
    {
        size_t i = 0;
        for (; i<q.length() && i<m_terms.size(); ++i)
        {
            m_terms[i] -= q[i];
        }

        for (; i<q.length() ; ++i)
            m_terms.push_back(-q[i]);

        return *this;
    }

    /**
     *  \brief Product of two Polynomials
     */
    Polynomial<T>& operator*=(const Polynomial<T>& q)
    {
        if (q == Polynomial<T>(0) || *this == Polynomial<T>(0))
        {
            m_terms.assign(1,static_cast<T>(0));
            return *this;
        }

        const size_t len(q.degree() + degree() + 1);
        std::vector<T> vec(len,static_cast<T>(0));

        for (size_t k(0); k<len; ++k)
            for (size_t i(0); i< this->length();++i)
            {
                const size_t j = k-i;
                if (j<q.length())
                    vec[k] += m_terms[i]*q[j];
            }

        *this = Polynomial<T>(vec);
        return *this;
    }

    /**
     *  \brief Quotient of the Euclidian division
     */
    inline Polynomial<T>& operator/=(const Polynomial<T>& q)
    {
        *this = division(*this,q).first;
        return *this;
    }

    /**
     *  \brief Remainder of the Euclidian division
     */
    inline Polynomial<T>& operator%=(const Polynomial<T>& q)
    {
        *this = division(*this,q).second;
        return *this;
    }

//Derivative -------------------------------------------------------------------------

    /**
     *  \brief Returns the derivative of the polynom
     */
    Polynomial<T> derivative() const
    {
        if (degree() == 0)
            return 0;

        Polynomial<T> ret(*this);
        for (size_t i(0); i<ret.degree(); ++i)
            ret[i] = ret[i+1] * (i+1);

        ret.m_terms.pop_back();

        return ret;
    }

private:

//Simplification -----------------------------------------------------------------------

    /**
     *  \brief Simplifies the polynom by removing all 0-terms after the highest non-zero term
     */
    void simplify()
    {
        size_t counter = 0;
        typedef typename std::vector<T>::reverse_iterator Iterator;
        for (Iterator it = m_terms.rbegin(); it != m_terms.rend() && *it == static_cast<T>(0) ; ++it)
        {
            ++counter;
        }

        m_terms.resize(m_terms.size() - counter);
    }

    /**
     *  \brief Notation simplification
     */
    typedef std::pair<Polynomial<T>, Polynomial<T> >  Pair;

    /**
     *  \brief Returns a Pair containing the quotient and the remainder of a division
     */
    static Pair division(const Polynomial& dividend, const Polynomial& divisor)
    {
        if (divisor == Polynomial<T>(static_cast<T>(0)))
            throw std::runtime_error("Polnomial: Division by 0");

        if (dividend.degree() < divisor.degree())
            return Pair(0, dividend);

        if (dividend == divisor)
            return Pair(1,0);

        Polynomial<T> remainder(dividend);
        Polynomial quotient;

        while (remainder.degree() > divisor.degree())
        {
            std::vector<T> vec(remainder.degree() - divisor.degree() + 1,0);
            vec.back() = remainder.m_terms.back() / divisor.m_terms.back();

            Polynomial<T> p(vec);

            remainder -= p*divisor;
            remainder.simplify();

            quotient += p;
        }

        return Pair(quotient,remainder);
    }

//Attributes-----------------------------------------------------------------------------

    std::vector<T> m_terms; /*!< The coefficients */

//Friends--------------------------------------------------------------------------------

    friend bool operator==(const Polynomial<T>& p, const Polynomial<T>& q);

};

// Stream -----------------------------------------------------------------------------------
/**
 *  \brief Writes the polynom in a stream
 */
template <typename T>
std::ostream& operator<<(std::ostream& stream, const Polynomial<T>& poly)
{
    if (poly.degree() == 0 && poly[0] == 0)
    {
        stream << static_cast<T>(0);
        return stream;
    }

    if (poly[poly.length() - 1] < static_cast<T>(0))
        stream << "- ";

    for (size_t i(1);i<poly.length();++i)
    {
        if (poly[poly.length() - i] != static_cast<T>(0))
            stream << std::abs(poly[poly.length() - i]) << "x^" << poly.length() - i;

        if (poly[poly.length() - (i+1)] > static_cast<T>(0))
            stream << " + ";
        else if (poly[poly.length() - (i+1)] < static_cast<T>(0))
            stream << " - ";
    }

    if (poly[0] != static_cast<T>(0))
        stream << std::abs(poly[0]);

    return stream;
}

//Comparison operators ------------------------------------------------------------------------
/**
 *  \brief Compares two polynoms to see if they are equal
 */
template <typename T>
bool operator==(const Polynomial<T>& p, const Polynomial<T>& q)
{
    return p.m_terms == q.m_terms;
}

/**
 *  \brief Compares two polynoms to see if they are different
 */
template <typename T>
inline bool operator!=(const Polynomial<T>& p, const Polynomial<T>& q)
{
    return !(p == q);
}

// Binary Operators ---------------------------------------------------------------------------

template<typename T>
inline Polynomial<T> operator+(const Polynomial<T>& p, const Polynomial<T>& q)
{
    return Polynomial<T>(p) += q;
}

template<typename T>
inline Polynomial<T> operator-(const Polynomial<T>& p, const Polynomial<T>& q)
{
    return Polynomial<T>(p) -= q;
}

template<typename T>
inline Polynomial<T> operator*(const Polynomial<T>& p, const Polynomial<T>& q)
{
    return Polynomial<T>(p) *= q;
}

template<typename T>
inline Polynomial<T> operator/(const Polynomial<T>& p, const Polynomial<T>& q)
{
    return Polynomial<T>(p) /= q;
}

template<typename T>
inline Polynomial<T> operator%(const Polynomial<T>& p, const Polynomial<T>& q)
{
    return Polynomial<T>(p) %= q;
}

//Unary operators -------------------------------------------------------------------------------

template<typename T>
Polynomial<T> operator-(const Polynomial<T>& p)
{
    Polynomial<T> result(p);
    for (size_t i(0); i<result.length(); ++i)
        result[i] = -result[i];
    return result;
}

template<typename T>
Polynomial<T> operator+(const Polynomial<T>& p)
{
    Polynomial<T> result(p);
    for (size_t i(0); i<result.length(); ++i)
        result[i] = +result[i];
    return result;
}

//Logical operators ----------------------------------------------------------------------------

template<typename T>
inline bool operator!(const Polynomial<T>& p)
{
    return p == Polynomial<T>(0);
}

//Derivative -----------------------------------------------------------------------------------

/**
 *  \brief Returns the derivative of the Polynomial passed
 */
template<typename T>
inline Polynomial<T> derivative(const Polynomial<T>& p)
{
    return p.derivative();
}

//Useful functions -----------------------------------------------------------------------------

/**
 *  \brief Returns the n-th power of a given x
 */
template<typename T>
T power(const T& x, size_t n)
{
    T result = static_cast<T>(1);

    for (size_t i(0); i<n ; ++i)
        result *= x;
    return result;
}

//Simplification typedef -----------------------------------------------------------------------

/**
 *  \typedef Poly
 *  A simpler name for standard polynomial using double as coefficients
 */
typedef Polynomial<double> Poly;

/**
 *  \typedef Polyd
 *  A simpler name for standard polynomial using double as coefficients
 */
typedef Polynomial<double> Polyd;

/**
 *  \typedef Polyf
 *  A simpler name for standard polynomial using float as coefficients
 */
typedef Polynomial<float> Polyf;

#endif //__POLYNOM_HPP_INCLUDED__

// array_helpers.hpp 
namespace luc_lib {
    // ...
    template <typename T, int N>
	inline std::reverse_iterator<T*> rbegin(T (&array)[N])
	{ return std::reverse_iterator<T*>(array + N); }

    template <typename T, int N>
	inline std::reverse_iterator<T*> rend  (T (&array)[N])
	{ return std::reverse_iterator<T*>(array + 0); }
}

#ifndef POLYNOMIAL_HPP
#define POLYNOMIAL_HPP

#include "../array_helpers.hpp"
#include <vector>
#include <iosfwd>
#include <cmath>

namespace luc_lib
{
    typedef size_t                 degree_type;

/*===========================================================================*/
/*================================[ is_null ]================================*/
/*===========================================================================*/
    template <typename N>
    inline 
    bool is_null(N n) {
        static const N k_err = 1e5;
        const bool null =
            (std::abs(n) <= k_err * std::numeric_limits<N>::epsilon());
        return null;
    }

    template <>
    inline 
    bool is_null<int>(int n) {
        return n == 0;
    }

/*===========================================================================*/
/*===============================[ Monomial ]================================*/
/*===========================================================================*/
   /** Monomial class.
    *
    * This class serves as an optimisation hack when defining new
    * polynomials.
    * @note default copy operations, and destructor are fine.
    */
    template <typename real_type>
    class Monomial
    {
    public:
	Monomial(real_type c, degree_type d) : m_coefficient(c), m_degree(d) {}
	explicit Monomial(degree_type d) : m_degree(d) {}

	real_type   coefficient() const { return m_coefficient; }
	degree_type degree     () const { return m_degree;      }

	Monomial& operator*=(real_type c) {
	    m_coefficient *= c;
	    return *this;
	}
	Monomial& operator^=(degree_type d) {
	    m_degree *= d;
	    return *this;
	}
    private:
        real_type      m_coefficient;
        degree_type    m_degree;
    };

    template <typename real_type>
    inline
    std::ostream& operator<<(std::ostream & os, Monomial<real_type> const& m)
    {
	return os << m.coefficient() << "X^" << m.degree() ;
    }

/*===========================================================================*/
/*==============================[ Polynomial ]===============================*/
/*===========================================================================*/
   /** Polynomial class.
    *
    * @note default copy operations, and destructor, provided thanks to
    * \c std::vector<> are fine.
    */
    template <typename real_type>
    class Polynomial 
    {
    public:
        typedef std::vector<real_type> coeffs_type;

	Polynomial(real_type c = real_type())
	    : m_coefficients(1)
	    { m_coefficients[0] = c; }

#if 1
        template <typename F>
            Polynomial(F c)
            : m_coefficients(1)
            { m_coefficients[0] = c; }
#endif

	/*explicit*/ Polynomial(Monomial<real_type> const& m)
	    : m_coefficients(m.degree()+1)
	    { m_coefficients.back() = m.coefficient();}

	explicit Polynomial(coeffs_type const& c)
	    : m_coefficients(c)
	    {}

	template <degree_type N> Polynomial(real_type const (&coeffs)[N])
	    : m_coefficients(luc_lib::rbegin(coeffs), luc_lib::rend(coeffs))
	    {
	    }

        void swap(Polynomial & other) {
            m_coefficients.swap(other.m_coefficients);
        }

	degree_type const degree() const { return m_coefficients.size()-1; }

	/** Evaluates the polynomial function.
	 * Implements Horner algorithm.
	 */
	real_type const operator()(real_type x) const
	{
	    real_type res = real_type(); // 0-init
	    degree_type i = degree()+1;
	    do {
		res = res * x + m_coefficients[--i];
	    } while (i) ;
	    return res;
	}

	/** Coefficients accessor.
	 */
	real_type operator[](degree_type d) const { return m_coefficients[d]; }

	Polynomial & operator+=(Polynomial const& rhs)
	{
	    const degree_type N = rhs.degree()+1;
	    if (m_coefficients.size() < N) { m_coefficients.resize(N); }
	    for (degree_type i=0; i!=N ; ++i)
		m_coefficients[i] += rhs[i];
	    do_trim();
	    return *this;
	}

	Polynomial & operator+=(Monomial<real_type> const& rhs)
	{
	    const degree_type N = rhs.degree()+1;
	    if (m_coefficients.size() < N) { m_coefficients.resize(N); }
	    // std::cout << *this << "(" << m_coefficients.size() << ") PLUS ";
	    m_coefficients[N-1] += rhs.coefficient();
	    // std::cout << rhs << " = " << *this<<"\n";
	    do_trim();
	    return *this;
	}

	Polynomial const operator-() const
	{
	    const degree_type N = degree()+1;
	    Polynomial res(N); // not very efficient
	    for (degree_type i=0; i!=N ; ++i)
		res.m_coefficients[i] = - m_coefficients[i];
	    return res;
	}

	Polynomial & operator-=(Polynomial const& rhs)
	{
	    const degree_type N = rhs.degree()+1;
	    if (m_coefficients.size() < N) { m_coefficients.resize(N); }
	    for (degree_type i=0; i!=N ; ++i)
		m_coefficients[i] -= rhs[i];
	    do_trim();
	    return *this;
	}

	Polynomial & operator-=(Monomial<real_type> const& rhs)
	{
	    const degree_type N = rhs.degree()+1;
	    if (m_coefficients.size() < N) { m_coefficients.resize(N); }
	    m_coefficients[N-1] -= rhs.coefficient();
	    do_trim();
	    return *this;
	}

	Polynomial & operator*=(Polynomial<real_type> const& rhs)
	{
	    Polynomial tmp = (*this) * rhs;
	    tmp.m_coefficients.swap(m_coefficients);
	    return *this;
	}

	friend Polynomial const operator*(Polynomial const& lhs, Polynomial const& rhs)
	{
	    const degree_type I = lhs.degree();
	    const degree_type J = rhs.degree();
	    /*Polynomial::*/coeffs_type res(I + J + 1);
	    for (degree_type i=0; i<=I ; ++i)
		for (degree_type j=0; j<=J ; ++j)
                    res[i+j] += lhs[i] * rhs[j];
	    return Polynomial(res);
	}

	Polynomial & operator*=(Monomial<real_type> const& rhs)
	{
	    (*this) *= rhs.coefficient();
	    m_coefficients.insert(m_coefficients.begin(), rhs.degree(), real_type());
	    return *this;
	}

	static void div(Polynomial lhs, Polynomial const& rhs, Polynomial &Q, Polynomial &R)
	    ;

	Polynomial & operator/=(Monomial<real_type> const& rhs) {
            Polynomial Q, R;
            div(*this, rhs, Q, R);
            this->swap(Q);
            return *this;
        }

	Polynomial & operator%=(Monomial<real_type> const& rhs) {
            Polynomial Q, R;
            div(*this, rhs, Q, R);
            this->swap(R);
            return *this;
        }

	friend bool operator==(Polynomial const& lhs, Polynomial const& rhs) {
	    return lhs.m_coefficients == rhs.m_coefficients;
	}

    private:
	void do_trim() {
	    while (m_coefficients.size()>1 && is_null(m_coefficients.back()))
	    {
                // std::cout << "trim:" << m_coefficients.back() << "\n";
		m_coefficients.pop_back();
	    }
	}
    private:
	coeffs_type m_coefficients;

	friend std::ostream& operator<<(std::ostream & os, Polynomial const& p)
	{
	    degree_type i = p.degree()+1;
	    if (i>1)
	    {
                bool needPlus = false;
		do
		{
		    if (! is_null(p[--i]))
		    {
			if (p[i] < 0) {
			    os << " - ";
			} else if (needPlus) {
			    os << " + ";
			}
			needPlus = true;
			if (p[i] != 1 || i==0) {
			    os << std::abs(p[i]) ; 
			}
			if (i>0) {
			    os << 'X';
			    if (i>1) {
				os << '^' << (i);
			    }
			}
		    }
		}
		while (i) ;
	    } else {
		os << (i==0 ? 0 : p[0]);
	    }
	    return os;
	}
    };


/*===========================================================================*/
/*===============================[ operators ]===============================*/
/*===========================================================================*/
    template <typename F, typename real_type>
    inline
    Monomial<real_type> const operator*(F c, Monomial<real_type> const& rhs)  {
        return Monomial<real_type>(rhs) *= c;
    }
    template <typename real_type>
    inline
    Monomial<real_type> const operator^(Monomial<real_type> const& lhs, degree_type d)  {
        return Monomial<real_type>(lhs) ^= d;
    }

/*===========================================================================*/

    template <typename real_type, typename convertible_real_type>
    inline
    Polynomial<real_type> const operator+(Polynomial<real_type> const& lhs, convertible_real_type const& rhs)  {
        return Polynomial<real_type>(lhs) += rhs;
    }
    template <typename real_type, typename convertible_real_type>
    inline
    Polynomial<real_type> const operator+(convertible_real_type const& lhs, Polynomial<real_type> const& rhs)  {
        return Polynomial<real_type>(rhs) += lhs;
    }

    template <typename real_type>
    inline
    Polynomial<real_type> const operator+(Polynomial<real_type> const& lhs, Monomial<real_type> const& rhs)  {
        return Polynomial<real_type>(lhs) += rhs;
    }
    template <typename real_type>
    inline
    Polynomial<real_type> const operator+(Monomial<real_type> const& lhs, Polynomial<real_type> const& rhs)  {
        return Polynomial<real_type>(rhs) += lhs;
    }

    template <typename real_type>
    inline
    Polynomial<real_type> const operator+(Polynomial<real_type> const& lhs, Polynomial<real_type> const& rhs)  {
        return Polynomial<real_type>(lhs) += rhs;
    }

    template <typename real_type>
    inline
    Polynomial<real_type> const operator+(Monomial<real_type> const& lhs, Monomial<real_type> const& rhs)  {
        return Polynomial<real_type>(lhs) += rhs;
    }
    template <typename real_type, typename convertible_real_type>
    inline
    Polynomial<real_type> const operator+(Monomial<real_type> const& lhs, convertible_real_type const& rhs)  {
        return Polynomial<real_type>(lhs) += rhs;
    }
    template <typename real_type, typename convertible_real_type>
    inline
    Polynomial<real_type> const operator+(convertible_real_type const& lhs, Monomial<real_type> const& rhs)  {
        return Polynomial<real_type>(rhs) += lhs;
    }

/*===========================================================================*/
    template <typename real_type>
    inline
    Polynomial<real_type> const operator-(Polynomial<real_type> const& lhs, Monomial<real_type> const& rhs)  {
        return Polynomial<real_type>(lhs) -= rhs;
    }
    template <typename real_type>
    inline
    Polynomial<real_type> const operator-(Monomial<real_type> const& lhs, Polynomial<real_type> const& rhs)  {
        return Polynomial<real_type>(-rhs) += lhs;
    }

    template <typename real_type>
    inline
    Polynomial<real_type> const operator-(Polynomial<real_type> const& lhs, Polynomial<real_type> const& rhs)  {
        return Polynomial<real_type>(lhs) -= rhs;
    }

    template <typename real_type, typename convertible_real_type>
    inline
    Polynomial<real_type> const operator-(Monomial<real_type> const& lhs, convertible_real_type const& rhs)  {
        return Polynomial<real_type>(lhs) -= rhs;
    }

    template <typename real_type>
    inline
    Polynomial<real_type> const operator-(Monomial<real_type> const& lhs, Monomial<real_type> const& rhs)  {
        return Polynomial<real_type>(lhs) -= rhs;
    }

    template <typename real_type, typename convertible_real_type>
    inline
    Polynomial<real_type> const operator-(convertible_real_type const& lhs, Monomial<real_type> const& rhs)  {
        return Polynomial<real_type>(rhs) -= lhs;
    }

/*===========================================================================*/
    template <typename real_type>
    inline
    Polynomial<real_type> const operator/(Polynomial<real_type> const& lhs, Polynomial<real_type> const& rhs)  {
        Polynomial<real_type> Q,R;
        Polynomial<real_type>::div(lhs, rhs, Q, R);
        return Q;
    }

    template <typename real_type>
    inline
    Polynomial<real_type> const operator%(Polynomial<real_type> const& lhs, Polynomial<real_type> const& rhs)  {
        Polynomial<real_type> Q,R;
        Polynomial<real_type>::div(lhs, rhs, Q, R);
        return R;
    }

#if  0
    inline
        std::ostream& operator<<(std::ostream & os, Polynomial const& p)
        {
        }
#endif

    template <typename real_type>
    inline
    void Polynomial<real_type>::div(
            Polynomial<real_type> lhs, Polynomial<real_type> const& rhs,
            Polynomial<real_type> &Q, Polynomial<real_type> &R)
    {
        Q.m_coefficients.resize(0); Q.m_coefficients.reserve(lhs.degree()-rhs.degree()+1);
        // R.m_coefficients.resize(0); R.m_coefficients.reserve(rhs.degree());

        while (lhs.degree() >= rhs.degree()) {
            const real_type   c = lhs[lhs.degree()] / rhs[rhs.degree()];
            const degree_type d = lhs.degree() - rhs.degree(); 
            const Monomial<real_type> m(c,d);
            lhs -= m * rhs;
            Q += m;
        }
        R.swap(lhs);
    }

    // const Monomial X(1,1);
#define X Monomial(1,1)    
} // namespace 


#endif // POLYNOMIAL_HPP
// Vim: let g:BTW_project='test_polynomial'

#ifndef ESSAIS_POLYNOMIAL_TEST_POLYNOMIAL_HPP__
#define ESSAIS_POLYNOMIAL_TEST_POLYNOMIAL_HPP__
#include <cxxtest/TestSuite.h>

#include <sstream>

// #define DOUBLE
#define TEMPLATE
// #define EXPR_TEMPLATE

#if defined(DOUBLE)
#  include "polynomial.hpp"
#elif defined(TEMPLATE)
#  include "polynomial_T.hpp"
#  define Polynomial Polynomial<double>
#  define Monomial   Monomial<double>
#elif defined(EXPR_TEMPLATE)
#  include "polynomial_ET.hpp"
#  define Polynomial Polynomial<double>
#  define Monomial   Monomial<double>
#endif
using namespace luc_lib;


namespace CxxTest 
{
    CXXTEST_TEMPLATE_INSTANTIATION
        class ValueTraits<Polynomial>
        {
            std::string value;
        public:
            ValueTraits(Polynomial const& p) {
                std::ostringstream oss;
                oss << p;
                value = oss.str();
            }
            const char *asString( void ) const { return value.c_str(); }
        };
}

/*===========================================================================*/
/*==============================[ ValueTraits ]==============================*/
/*===========================================================================*/

class test_polunomial : public CxxTest::TestSuite
{
public:
    void test_epsilon()
    {
	TS_TRACE(is_null(0));
	TS_TRACE(is_null(1e-9));
	TS_TRACE(is_null(1e-10));
	TS_TRACE(is_null(1e-11));

    }

/*===========================================================================*/
/*===============================[ Monomials ]===============================*/
/*===========================================================================*/
    void test_X()
    {
        Monomial m = X;
        TS_ASSERT_EQUALS(m.degree(), 1);
        TS_ASSERT_EQUALS(m.coefficient(), 1);
    }
    void test_X3()
    {
        Monomial m = X^3;
        TS_ASSERT_EQUALS(m.degree(), 3);
        TS_ASSERT_EQUALS(m.coefficient(), 1);
    }
    void test_2X()
    {
        Monomial m = 2 * X;
        TS_ASSERT_EQUALS(m.degree(), 1);
        TS_ASSERT_EQUALS(m.coefficient(), 2);
    }
    void test_2X3()
    {
        Monomial m = 2 * X^3;
        TS_ASSERT_EQUALS(m.degree(), 3);
        TS_ASSERT_EQUALS(m.coefficient(), 2);
    }

/*===========================================================================*/
/*==============================[ Polynomials ]==============================*/
/*===========================================================================*/
    void test_XpX()
    {
        Polynomial p = X + X;
        TS_ASSERT_EQUALS(p.degree(), 1);
        TS_ASSERT_EQUALS(p[1], 2);
        TS_ASSERT_EQUALS(p[0], 0);

        TS_ASSERT_EQUALS(p(0), 0);
        TS_ASSERT_EQUALS(p(1), 2);
        TS_ASSERT_EQUALS(p(2), 4);
        TS_ASSERT_EQUALS(p(3), 6);
    }
    void test_Xp5()
    {
        Polynomial p = X + 5;
        TS_ASSERT_EQUALS(p.degree(), 1);
        TS_ASSERT_EQUALS(p[1], 1);
        TS_ASSERT_EQUALS(p[0], 5);

        TS_ASSERT_EQUALS(p(0), 5);
        TS_ASSERT_EQUALS(p(1), 6);
        TS_ASSERT_EQUALS(p(2), 7);
        TS_ASSERT_EQUALS(p(3), 8);
    }
    void test_22X2p11Xp5()
    {
        Polynomial p = (22*X^2) + (11*X) + 5;
        p = (22*X^2) + (11*X) + 5;
        TS_ASSERT_EQUALS(p.degree(), 2);
        TS_ASSERT_EQUALS(p[2], 22);
        TS_ASSERT_EQUALS(p[1], 11);
        TS_ASSERT_EQUALS(p[0], 5);

        TS_ASSERT_EQUALS(p(0), 5);
        TS_ASSERT_EQUALS(p(1), 22*1*1 + 11*1 + 5);
        TS_ASSERT_EQUALS(p(2), 22*2*2 + 11*2 + 5);
        TS_ASSERT_EQUALS(p(3), 22*3*3 + 11*3 + 5);
    }
    void test_33X322X2p11Xp5()
    {
        Polynomial p = (33*X^3) + (22*X^2) + (11*X) + 5;
        p = (33*X^3) + (22*X^2) + (11*X) + 5;
        TS_ASSERT_EQUALS(p.degree(), 3);
        TS_ASSERT_EQUALS(p[3], 33);
        TS_ASSERT_EQUALS(p[2], 22);
        TS_ASSERT_EQUALS(p[1], 11);
        TS_ASSERT_EQUALS(p[0], 5);

        TS_ASSERT_EQUALS(p(0), 5);
        TS_ASSERT_EQUALS(p(1), 33*1*1*1 + 22*1*1 + 11*1 + 5);
        TS_ASSERT_EQUALS(p(2), 33*2*2*2 + 22*2*2 + 11*2 + 5);
        TS_ASSERT_EQUALS(p(3), 33*3*3*3 + 22*3*3 + 11*3 + 5);
    }

    void test_mult_Xp1_sq()
    {
        const Polynomial p1 = (X+1);
        const Polynomial m = p1 * p1;

        TS_ASSERT_EQUALS(m.degree(), 2);
        TS_ASSERT_EQUALS(m[2], 1);
        TS_ASSERT_EQUALS(m[1], 2);
        TS_ASSERT_EQUALS(m[0], 1);
    }
    void test_mult_Xp1_dot_Xm1()
    {
        const Polynomial p1 = (X+1);
        const Polynomial p2 = (X-1);
        const Polynomial m = p1 * p2;

        TS_ASSERT_EQUALS(m.degree(), 2);
        TS_ASSERT_EQUALS(m[2], 1);
        TS_ASSERT_EQUALS(m[1], 0);
        TS_ASSERT_EQUALS(m[0], -1);
    }

    void test_div()
    {
	// const Polynomial A = (-7*X^3) + (2*X^2) + (3*X) + 1;
	// const Polynomial B = (-2*X^2) + (1*X) + 1;
	// const Polynomial eQ = (-5*X^3) + (2*X^2) + (2*X) + 1;
	// const Polynomial eR = -10*X + 9;

	const Polynomial A = (X^5) - (X^4) - (X^3) + (3*X^2) - (2*X);
	const Polynomial B = (X^2) - (X) + 1;
	const Polynomial eQ = (X^3) - (2*X) + 1;
	const Polynomial eR = X - 1;

	// const Polynomial A = (X^5) + (2*X^3) - (3*X) - 2;
	// const Polynomial B = (X^3) + (X) + 1;
	// const Polynomial eQ = (X^2) + 1;
	// const Polynomial eR = (-1*X^2) - (4*X) - 3;

        Polynomial Q, R;
	Polynomial::div(A, B, Q, R);
	TS_ASSERT_EQUALS(Q, eQ);
	TS_ASSERT_EQUALS(R, eR);
	TS_ASSERT_EQUALS(A, B*Q+R);
    }
};


#endif // ESSAIS_POLYNOMIAL_TEST_POLYNOMIAL_HPP__

# ---------------------------------------------------------------------------
# Development tools directories definition
# ---------------------------------------------------------------------------
# CC_OPT_XARCH    = -xarch=v9
# LINK_OPT_XARCH  = -xarch=v9
# CC_DIR          = /opt/SUNWspro/bin
# CXX             = ${CC_DIR}/CC
# CXXFLAGS        = -g ${CC_OPT_XARCH} -DUNIX -DSUN 
# INCLUDE_PATH    = -I${HOME}/include

INCLUDE_PATH    = -I/cygdrive/e/Dev/include
CXXFLAGS        = -g -pedantic -Wall ${INCLUDE_PATH}
# ---------------------------------------------------------------------------
# Définition de la construction de la cible principale
# ---------------------------------------------------------------------------
all: test_polynomial

# ---------------------------------------------------------------------------
test_polynomial.cpp: test_polynomial.hpp polynomial.hpp polynomial_T.hpp polynomial_ET.hpp
        cxxtestgen.pl -error-printer -o $@ $^

# test_polynomial: test_polynomial.cpp 
        # $(CPP_COMPILER) ${CC_FLAGS} ${CC_OPT_FLAGS} ${INCLUDE_PATH} -o $@ $^ ${LIBRARIES}  
        # ./$@


# ---------------------------------------------------------------------------
# cibles de nettoyage
# ---------------------------------------------------------------------------
clean: 
        rm test_polynomial test_polynomial.o test_polynomial.cpp

#ifndef POLYNOME_H
#define POLYNOME_H

#include <vector>
#include <iostream>
#include <stdarg.h>
#include <stdexcept>

#define POLYNOME_NB_ITERATIONS_RECHERCHE 10
#define POLYNOME_NB_ITERATIONS_NEWTON 1000
#define POLYNOME_NB_ITERATIONS_DICHO 10
#define POLYNOME_BIG_EPSILON 0.00000001
#define POLYNOME_EPSILON 0.000000000000000001

class Polynome{
    public:
    // Constructeurs
        Polynome();
        Polynome(const Polynome &);
        Polynome(const std::vector<double> &);
        Polynome(int, const double, ...);
        Polynome(const int, const double[]);
    // Destructeur
        ~Polynome();

    // Opérateurs
        Polynome operator=(const Polynome &);

        Polynome operator+=(const Polynome &P);
        Polynome operator-=(const Polynome &P);
        Polynome operator*=(const Polynome &P);
        Polynome operator/=(const Polynome &P) throw(std::domain_error);
        Polynome operator%=(const Polynome &P) throw(std::domain_error);

        Polynome operator+(const Polynome &) const;
        Polynome operator-(const Polynome &) const;
        Polynome operator-() const;
        Polynome operator*(const Polynome &) const;
        Polynome operator/(const Polynome &) const throw(std::domain_error);
        Polynome operator%(const Polynome &) const throw(std::domain_error);
        bool operator==(const Polynome &) const;
        bool operator!=(const Polynome &) const;
        double operator()(const double) const;

    // Arithmétique
        static void div_eucli(const Polynome &, const Polynome &, Polynome &, Polynome &) throw(std::domain_error);
        static Polynome pgcd(const Polynome &, const Polynome &);
        Polynome derive(unsigned int d=1) const;
        Polynome integre(unsigned int d=1) const;
        void unitaire();

    // Racines
        int nb_roots();
        int nb_roots(double, double);
        std::vector<double> roots();
        void sturm_compute();
        // Accesseurs
            bool get_sturmcomputed();
        double newton(double) throw(std::domain_error, std::runtime_error);
        void dicho(double &, double &, int);

    // Autres
        double eval(double) const;
        int degre() const;
        void affiche(std::ostream &) const;

    protected:
        int nb_chg_signe(double);

    // Attributs
        std::vector<double> m_coeffs;   // vecteur des coefficients
        std::vector<Polynome> Sturm;    // Suite de Sturm du polynôme
        bool sturm_computed;            // Si la suite de Sturm a été calculée
        double M;                       // Majorant de la valeur absolue de toutes les racines

};

std::ostream& operator<<(std::ostream &, const Polynome &);

#endif

#include "polynome.h"
#include <cmath>

#include <iomanip>

using namespace std;


/********************************************************
                    Constructeurs
 *******************************************************/

// Constructeur de copie
Polynome::Polynome(const Polynome &P){
    m_coeffs = P.m_coeffs;
    sturm_computed = P.sturm_computed;
    Sturm = P.Sturm;
}

// Par défault : Créer le polynome nul
Polynome::Polynome(){
    m_coeffs.push_back(0.);
    sturm_computed = false;
}

// Via un vecteur : les coefficients du polynôme sont recopié du vecteur
Polynome::Polynome(const std::vector<double> &coeffs){
    m_coeffs = coeffs;
    sturm_computed = false;
}

// Via une liste variable : Polynome(nb de coeffs, a0, a1, a2, ...)
// ex : Polynome(3, 1., -2., 4.) crée le polynôme 4x^2-2x+1
// La liste variable est constituée de double !!! Polynome(2, 1, 2) ne renverra pas 2x+1 !!!
Polynome::Polynome(int n, const double x, ...){
    sturm_computed = false;
    if (n <= 0)
        return;

    va_list varg; // Variable stockant la liste
    va_start(varg, x); // On initialise la variable
    m_coeffs.push_back(x);
    n--;
    do {
        m_coeffs.push_back(va_arg(varg, double));
    }while(--n!=0);
    va_end(varg); // On détruit la liste variable
}

// Via un tableau de double, ceci permet de profiter de l'auto-cast int -> double
// et donc de pouvoir écrire Polynome(3, (double[]) {1, 0, 1.}) pour le polynôme x^2+1
// ou bien double coef[] = {1, 0, 1.}; Polynome(3, coef)
Polynome::Polynome(const int n, const double *liste){
    sturm_computed = false;
    if (n <= 0)
        return;
    for(int i = 0 ; i < n ; i++)
        m_coeffs.push_back(liste[i]);
}



/********************************************************
 *******************  Destructeur  **********************
 *******************************************************/
Polynome::~Polynome(){
    // Rien de bien palpitant
}




/********************************************************
 ******************  Opérateurs  ************************
 *******************************************************/

// Opérateur d'affectation
Polynome Polynome::operator=(const Polynome &P){
    m_coeffs = P.m_coeffs;
    sturm_computed = P.sturm_computed;
    Sturm = P.Sturm;
    return *this;
}

Polynome Polynome::operator+=(const Polynome &P)
{
    int n = this->degre();
    int m = P.degre();
    int i;

    if (m <= n)
    {
        // On rajoute un polynome de degré plus petit ou égal
        for(i = 0 ; i <= m ; i++)
            m_coeffs[i] += P.m_coeffs[i];
    }
    else
    {
        // On rajoute un polynôme de degré plus grand
        m_coeffs.reserve(m+1);
        for(i = 0 ; i <= n ; i++)
            m_coeffs[i] += P.m_coeffs[i];
        for(i = n + 1 ; i<= m ; i++)
            m_coeffs[i] = P.m_coeffs[i];
    }
    return *this;
}

Polynome Polynome::operator-=(const Polynome &P)
{
    return (*this)+=-P;
}

Polynome Polynome::operator*=(const Polynome &P)
{
    *this = *this * P;
    return *this;
}

Polynome Polynome::operator/=(const Polynome &P) throw(std::domain_error)
{
    if(P.degre() == -1)
        throw std::domain_error("Division par zéro");
    *this = *this / P;
    return *this;
}

Polynome Polynome::operator%=(const Polynome &P) throw(std::domain_error)
{
    if(P.degre() == -1)
        throw std::domain_error("Division par zéro");
    *this = *this % P;
    return *this;
}


// Ajoute 2 polynômes entre eux
Polynome Polynome::operator+(const Polynome &P) const{
    return Polynome(*this)+=P;
}

// Retranche 2 polynômes
Polynome Polynome::operator-(const Polynome &P) const{
    return Polynome(*this)+=-P;
}

// Renvoie l'opposée du polynôme
Polynome Polynome::operator-() const{
    int n = m_coeffs.size();

    Polynome res;
    res.m_coeffs.resize(n);

    for(int i = 0 ; i < n ; i++)
        res.m_coeffs[i] = - m_coeffs[i];
    return res;
}

// Multiplie 2 polynômes entre eux : X^i * X^j = X^(i+j)
Polynome Polynome::operator*(const Polynome &P) const{
    int n = this->degre();
    int m = P.degre();
    if (n < 0 || m < 0)
        return Polynome();
    vector<double> M(n+m+1,0);
    for(int i = 0 ; i<= n ; i++)
        for(int j = 0 ; j <= m ; j++)
            M[i+j] += m_coeffs[i]*P.m_coeffs[j];
    return Polynome(M);
}

// Renvoie le quotient de la division euclidienne
Polynome Polynome::operator/(const Polynome &P) const throw(std::domain_error){
    if(P.degre() == -1)
        throw std::domain_error("Division par zéro");
    Polynome Q,R;
    div_eucli(*this, P, Q, R);
    return Q;
}

// Renvoie le reste de la division euclidienne
Polynome Polynome::operator%(const Polynome &P) const  throw(std::domain_error){
    if(P.degre() == -1)
        throw std::domain_error("Division par zéro");
    Polynome Q;
    Polynome R;
    div_eucli(*this, P, Q, R);
    return R;
}

// Test l'égalité de 2 polynômes
bool Polynome::operator==(const Polynome &P) const{
    int n = this->degre();
    int m = P.degre();
    if (n != m)
        return false;
    for(int i = 0 ; i <=n ; i++)
        if(m_coeffs[i] != P.m_coeffs[i])
            return false;
    return true;
}

// Test la différence de 2 polynômes
bool Polynome::operator!=(const Polynome &P) const{
    return !(*this == P);
}

// Confie à Poly.affiche la tache cout << P
ostream &operator<<(ostream &out, const Polynome &Poly){
    Poly.affiche(out);
    return out;
}

// Evalue un polynôme en 1 point
double Polynome::operator()(const double x) const{
    return this->eval(x);
}






/********************************************************
 ******************  Arithmétique  **********************
 *******************************************************/

// Effectue la division euclidienne de A par B, le quotient va dans Q et le reste dans R
void Polynome::div_eucli(const Polynome &A, const Polynome &B, Polynome &Q, Polynome &R) throw(std::domain_error){
    if(B.degre() == -1)
        throw std::domain_error("Division par zéro");

    int n = A.degre();
    int m = B.degre();
    if (m < 0)
        return;
    if (n-m < 0){
        Q = Polynome();
        R = A;
        return;
    }

    R.m_coeffs = A.m_coeffs;
    Q.m_coeffs.resize(n-m+1);
    double pivot = 1/B.m_coeffs[m], r;

    for(int i = n-m ; i >= 0 ; i--){
        r = R.m_coeffs[i+m] * pivot;
        Q.m_coeffs[i] = r;
        if(r != 0)
            for(int j = i ; j<=i+m ; j++){
                R.m_coeffs[j] -= r * B.m_coeffs[j-i];
                if ( abs(R.m_coeffs[j]) < POLYNOME_BIG_EPSILON)
                    R.m_coeffs[j] = 0;
            }
        R.m_coeffs[i+m] = 0;
    }
}

// Renvoie le pgcd de 2 polynômes
Polynome Polynome::pgcd(const Polynome &A, const Polynome &B){
    Polynome A1 = A, B1 = B, Q, R;
    do{
        R = A1 % B1;
        A1 = B1;
        B1 = R;
    } while(R.degre() != -1);

    return A1;
}

// Renvoie la dérivée du polynôme
Polynome::Polynome Polynome::derive(unsigned int d) const{
    int n = this->degre();
    if (n+1-(int)d <= 0)
        return Polynome();
    vector<double> res(n+1-d);
    int p = 1;
    for(unsigned int i = 1 ; i <= d ; i++)
        p *= i;
    for(int i = 0; i <= n-(int)d ; i++){
        res[i] = p*m_coeffs[i+d];
        p = p * (d+i+1)/(i+1);
    }
    return Polynome(res);
}

// Renvoie la primitive du polynôme nulle en 0
Polynome::Polynome Polynome::integre(unsigned int d) const{
    int n = this->degre();
    if (n+1+d <= 0)
        return Polynome();
    vector<double> res(n+1+d);
    int p = 1;
    for(unsigned int i = 1 ; i <= d ; i++)
        p *= i;
    for(unsigned int i = 0 ; i < d ; i++)
        res[i] = 0;
    for(int i = 0 ; i <= n ; i++){
        res[i+d] = m_coeffs[i]/p;
        p = p * (d+i+1)/(i+1);
    }
    return Polynome(res);
}

// Rend le polynôme unitaire (son coefficient de plus haut degré vaudra 1)
void Polynome::unitaire(){
    int n = this->degre();
    if (n == -1)
        return;
    double div = 1 / m_coeffs[n];
    for(int i = 0 ; i <= n ; i++)
        m_coeffs[i] *= div;
    return;
}





/********************************************************
 ********************* Racines  *************************
 *******************************************************/

// Renvoie le nombre de racines réelles du polynômes (Théorème de Sturm)
int Polynome::nb_roots(){

    if( !sturm_computed )
        sturm_compute();

    return abs(nb_chg_signe(-M) - nb_chg_signe(M));
}

// Renvoie le nombre de racines réelles du polynômes dans l'intervalle [a,b] (Théorème de Sturm)
int Polynome::nb_roots(double a, double b){

    if( !sturm_computed )
        sturm_compute();

    return abs(nb_chg_signe(a) - nb_chg_signe(b));
}

// Calcul les racines réelles d'un polynôme
vector<double> Polynome::roots(){
    int nb_racines = nb_roots();

    if (nb_racines == 0)
        return vector<double>();

// Recherches les intervalles sur lesquels les racines sont isolés
// Cette recherche se fait via une dichotomie utilisée sur le théorème de Sturm
    double *bornes = new double[nb_racines + 1]; // Pour n racines, il y a n+1 intervalles
    bool *borne_fixe = new bool[nb_racines + 1];
    for(int i=0;i<nb_racines+1;i++)
        borne_fixe[i] = false;
    vector<double> racines;

    bornes[0] = -M; bornes[nb_racines] = M;        // On fixe les bornes limites [-M,M]
    borne_fixe[0] = true; borne_fixe[nb_racines] = true;

    for(int i = 1 ; i < nb_racines ; i++){

        // On recherche le premier intervalle qui contient plusieurs racines pour le scindé en 2
        // Derniere borne fixe
        int gauche = 0, droite, rac_gauche, rac_droite, rac_milieu, nb_rac;
        double x_g, x_d, milieu;
        bool continuer;

        while(borne_fixe[gauche])
            gauche++;
        gauche--;
        // Cette erreur ne devrait jamais se produire
        if (gauche > nb_racines){
            std::cerr << "Erreur lors du scindage de l'intervalle (recherche de la borne gauche\n";
            std::cerr << left << setw(20) << "Valeur des bornes" << "Fixé\n";
            for(int i = 0 ; i <= nb_racines ; i++)
                std::cerr << left << setw(20) << bornes[i] << borne_fixe[i] << "\n";
            std::cerr << flush;
            return vector<double>(0);
        }

        // Prochaine borne fixe
        droite = gauche + 1;
        while(!borne_fixe[droite])
            droite++;
        // Cette erreur ne devrait jamais se produire
        if (droite > nb_racines){
            std::cerr << "Erreur lors du scindage de l'intervalle (recherche de la borne droite\n";
            std::cerr << left << setw(20) << "Valeur des bornes" << "Fixé\n";
            for(int i = 0 ; i <= nb_racines ; i++)
                std::cerr << left << setw(20) << bornes[i] << borne_fixe[i] << "\n";
            std::cerr << flush;
            return vector<double>(0);
        }

        x_g = bornes[gauche];
        x_d = bornes[droite];
        rac_gauche = nb_chg_signe(bornes[gauche]);
        rac_droite = nb_chg_signe(bornes[droite]);
        // On cherche à couper l'intervalle [x_g, x_d] en 2 intervalles [x_g, milieu] et [x_d, milieu]
        // où chacun des 2 contient au moins une racine
        continuer = true;
        while(continuer){
            milieu = (x_g + x_d)/2;
            rac_milieu = nb_chg_signe(milieu);
            if( rac_milieu == rac_gauche)       // Il n'y a pas de racines entre x_g et milieu
                x_g = milieu;                       // donc on décale la borne gauche
            else if( rac_milieu == rac_droite)  // Il n'y a pas de racines entre milieu et x_d
                x_d = milieu;                       // donc on décale la borne droite
            else
                continuer = false;              // On a réussi à couper l'intervalle
        }

        // On scinde l'intervalle
        nb_rac = abs(rac_gauche - rac_milieu);
        bornes[gauche + nb_rac] = milieu;
        borne_fixe[gauche + nb_rac] = true;

        // Réduit la taille de l'intervalle total [-M, M] pour accélerer la convergence de la méthode de newton
        if (gauche == 0)
            bornes[gauche] = x_g;
        if (droite == nb_racines)
            bornes[droite] = x_d;
    }

    Polynome calc = *this / Polynome::pgcd(*this, this->derive()); // Permet de transformer toutes les racines de multiplicité > 1 en racines de multiplicité 1
    for(int i = 0 ; i < nb_racines ; i++){               // Ce qui permet une convergence quadratique de la méthode de newton
        int iter = 0;       // Nombre d'itérations pour le calcul d'une racine
        double rac;
        double borne_gauche = bornes[i], borne_droite = bornes[i+1];    // Bornes de recherche de la racine
        bool calcul = true; // On effectue le calcul
        while(calcul){
            try{
                iter++;
                rac = calc.newton((borne_gauche + borne_droite)/2);
                if (rac >= bornes[i] && rac <= bornes[i+1]){ // Si la racine se situe dans l'intervalle de recherche
                    racines.push_back(rac); // On l'empile sur le vecteur de solutions
                    calcul = false;         // Le calcul s'est effectué jusqu'au bout
                    calc = calc / Polynome(2, -racines[i], 1.); // Divise le polynôme de recherche par (X-rac) pour supprimer la racine (enlève un bassin d'attraction pour la méthode de newton)
                }
            }
            catch(const std::domain_error &e){
                // La dérivée s'est annulée, laisse calcul à true
            }
            catch(const std::runtime_error &e){
                // L'algorithme tourne en boucle, laisse calcul à true
            }
            if (calcul){ // Si l'on doit continuer à chercher la racine
                calc.dicho(borne_gauche, borne_droite, POLYNOME_NB_ITERATIONS_DICHO); // On divise l'intervalle de recherche par 2^POLYNOME_NB_ITERATIONS_DICHO
            }
            if (iter > POLYNOME_NB_ITERATIONS_RECHERCHE) // Si l'on a effectué trop d'itérations, on abandonne la recherche
                calcul = false;
        }
    }

    delete[] bornes;
    delete[] borne_fixe;

    return racines;
}

// Construit la suite de Sturm
void Polynome::sturm_compute(){
    Polynome A = *this;
    Polynome B = A.derive();
    Polynome R;
    A.unitaire();
    Sturm.clear();
    Sturm.push_back(A);
    int r = 1;

    do{
        Sturm.push_back(B);
        r++;
        R = - (A % B);
        A = B;
        B = R;
    } while(R.degre() != -1);

 // Si le polynôme a des racines multiples on réduit leur multiplicité à 1
    if (Sturm[r-1].degre() != 0)
        for(int i = 0 ; i < r ; i++)
            Sturm[i] = Sturm[i] / Sturm[r-1];

// Toutes les racines seront dans [-M, M]
    int degre = this->degre();
    M = 0;
    double coeff = m_coeffs[degre];
    for(int i = 0 ; i <= degre ; i++)
        M += abs(m_coeffs[i] / coeff);
    if (M < 1)
        M = 1;

    sturm_computed = true;
}

// Indique si la suite de Sturm est construite ou non
bool Polynome::get_sturmcomputed(){
    return sturm_computed;
}

// Calcul le nombre de changement de signe utilisé dans le théorème de Sturm
int Polynome::nb_chg_signe(double x){
    int nb = 0;
    int n = Sturm.size();

    // On recherche le premier terme non nul
    int i = 0;
    while(i < n && Sturm[i].eval(x) == 0)
        i++;
    if (i == n)
        return 0; // Tout les termes sont nuls

    double terme_precedent = Sturm[i].eval(x), calc;
    for(i = i+1 ; i < n ; i++){
        calc = Sturm[i].eval(x);
        if (terme_precedent * calc < 0){
            nb++;
            terme_precedent = calc;
        }
    }

    return nb;
}

// Approxime la racine via l'algorithme de Newton
double Polynome::newton(double x) throw(std::domain_error, std::runtime_error){
    Polynome DP = this->derive();
    int iter = 0;
    double f_x = eval(x), fp_x;
    while(abs(f_x) > POLYNOME_EPSILON){
        fp_x = DP.eval(x);

        // Exceptions
        if (abs(fp_x) < POLYNOME_EPSILON)    // La dérivée est nulle
            throw std::domain_error("Division par zéro");
        if (iter > POLYNOME_NB_ITERATIONS_NEWTON){                   // On a effectué trop d'itérations
            if (abs(f_x) < POLYNOME_BIG_EPSILON)            // Mais on est presque sur un zéro
                return x;
            throw std::runtime_error("Trop d'itérations");
        }

        x = x - f_x/fp_x;
        f_x = eval(x);
        iter++;
    }
    return x;
}

// Effectue un dichotomie sur un nombre d'itérations fixés
// Elle réduit l'intervalle [b_g, b_d] autour de la racine
void Polynome::dicho(double &b_g, double &b_d, int iter){
    double val_g = eval(b_g), val_d = eval(b_d);
    double milieu, val_m;

    if (val_g * val_d > 0)      // On ne peut pas lancer la dichotomie
        return;

    while (iter > 0){
        milieu = (b_g + b_d) / 2;
        val_m = eval(milieu);
        if (val_g * val_m <= 0){ // La racine est située dans [b_g, milieu]
            b_d = milieu; // On se place sur le segment [b_g, milieu]
        }
        else{                   // La racine est située dans [milieu, b_d]
            b_g = milieu; // On se place sur le segment [milieu, b_d]
            val_g = val_m;
        }
        iter--;
    }
}


/********************************************************
 ********************** Autres  *************************
 *******************************************************/

// Evalue le polynôme en un point, méthode de Horner
double Polynome::eval(double x) const{
    double res = 0;
    int n = this->degre();
    for(int i = n ; i >= 0 ; i--){
        res*= x;
        res+= m_coeffs[i];
    }
    return res;
}

// Renvoie le degré du polynôme. Par convention le degré du polynôme nul sera -1 et non pas - l'infini
int Polynome::degre() const{
    int n = m_coeffs.size() - 1;
    while(n >= 0 && m_coeffs[n] == 0)
        n--;
    return n;
}

// Envoie l'expression litérale du polynôme sur out (utiliser par cout << P)
void Polynome::affiche(ostream &out) const{
    int i;
    double c;
    bool poly_nul = true;
    for(i=m_coeffs.size()-1 ; i!=-1 ; i--){
        c = m_coeffs[i];
        if (i == 0){
            if(poly_nul)
                out << noshowpos << c;
            else
                if(c != 0)
                    out << showpos << c;
        }
        else{
            if(poly_nul && c != 0){
                if(c == -1)
                    out << '-';
                else if(c != 1)
                    out << noshowpos << c;
                poly_nul = false;
            }
            else if(c != 0){
                if( c == 1)
                    out << '+';
                else if(c == -1)
                    out << '-';
                else
                    out << showpos << c;
            }
        }
        if(c != 0){
            if(i==1)
                out << 'X';
            else if(i > 1)
                out << "X^(" << noshowpos << i <<')';
        }
    }
    out << noshowpos;
}

#include <iostream>
#include "polynome.h"
#include <windows.h>
using namespace std;

// Affiche un vector<double>
std::ostream &operator<<(std::ostream &out, const std::vector<double> &vec){
    int n = vec.size();
    for(int i = 0 ; i < n ; i++)
        out << "[" << i << "] = " << vec[i] << endl;
    return out;
}

int main(int argc, char *argv[])
{
//    SetConsoleOutputCP( 1252 );

    Polynome A(6, (double[]) {1, 0, 0, -1, 4., 1.});
    cout << "A(X) = " << A << endl;

    Polynome B(3, 1., 0., 1.); // Création via une liste variable de double
    cout << "B(X) = " << B << endl;

    cout << endl << "Division euclidienne : A = B * Q + R" << endl;
    Polynome Q, R;
    Polynome::div_eucli(A,B,Q,R);

    cout << "Q(X) = " << Q << endl;
    cout << "R(X) = " << R << endl;

    cout << "Vérification (on teste en même temps l'opérateur d'égalité) : ";
    if (A == B * Q + R)
        cout << "ok" << endl;
    else
        cout << "erreur" << endl;

    cout << endl << "A(X) / B(X) = " << A / B << endl;
    cout << "A(X) % B(X) = " << A % B << endl;

    cout << endl << "Calcul de PGCD" << endl;
    double coeffs_A[] = {7, -7, 0, 8, -7, 0, 1};
    A = Polynome(7, coeffs_A);
    B = Polynome(6, (double[]) {-7, 0, 3, -7, 0, 3});
    cout << "A(X) = " << A << endl << "B(X) = " << B << endl;

    R = Polynome::pgcd(B,A);
    cout << "PGCD(A,B) : " << R << endl;
    R.unitaire();
    cout << "PGCD(A,B) unitarisé : " << R << endl;

    A = Polynome(6, -12., 20., -1., -9., 1., 1.);
    //A = Polynome(3, -1., 0., 1.);
    A = Polynome(11, 30., 30., -11., 35., -13., 36., 36., 20., 42., 15., -34.);
    cout << endl << "A(X) = " << A << endl;
    cout << "Nombre de racines réelles : " << A.nb_roots() << endl;

    cout << "Racines : " << endl << A.roots();

    return 0;
}