MPolynomial.h

#pragma once

#if defined(MPolynomial_RECURSES)
#error Recursive header files inclusion detected in MPolynomial.h
#else // defined(MPolynomial_RECURSES)

#define MPolynomial_RECURSES

#if !defined MPolynomial_h

#define MPolynomial_h

// Inclusions
#include <iostream>
#include <vector>
#include "DGtal/base/Common.h"

namespace DGtal
{

  template < int n, typename TRing, 
             typename TAlloc = std::allocator<TRing> >
  class MPolynomial;
  template < int N, int n, typename TRing, 
             typename TAlloc = std::allocator<TRing> >
  class MPolynomialDerivativeComputer;
  template < int n, typename TRing, 
             typename TAlloc = std::allocator<TRing>, 
             typename TX = TRing >
  class MPolynomialEvaluator;

  template < int n,  typename TRing, typename TOwner, 
             typename TAlloc, typename TX >
  class MPolynomialEvaluatorImpl;

  template < typename TRing, typename TAlloc >
  void euclidDiv( const MPolynomial<1, TRing, TAlloc> & f, 
                  const MPolynomial<1, TRing, TAlloc> & g, 
                  MPolynomial<1, TRing, TAlloc> & q, 
                  MPolynomial<1, TRing, TAlloc> & r );

  template < typename TRing, typename TOwner, 
             typename TAlloc, typename TX >
  class MPolynomialEvaluatorImpl< 1, TRing, TOwner, TAlloc, TX >
  {
  public:
    typedef TRing Ring;
    typedef TOwner Owner;
    typedef TAlloc Alloc;
    typedef TX X;

    template <int nn, class TT, class AA, class SS>
    friend class MPolynomialEvaluator;
    
    template <int nn, class TT, class HLHL, class AA, class SS>
    friend class MPolynomialEvaluatorImpl;

private:
    const Owner & myOwner; 
    const X & myEvalPoint; 
    
    inline MPolynomialEvaluatorImpl( const Owner & owner, const X & evalpoint)
        : myOwner(owner), myEvalPoint(evalpoint)
    {}
    
    class EvalFun
    {
      const MPolynomialEvaluatorImpl<1, Ring, Owner, Alloc, X> & myOwner;
        
    public:
      inline 
      EvalFun( const MPolynomialEvaluatorImpl<1, Ring, Owner, Alloc, X> & 
               owner )
        : myOwner(owner)
      {}
      
      inline
      X operator() ( const MPolynomial<1, Ring, Alloc> & p ) const
      {
        X res = (X) 0;
        X xx = (X) 1;
        for ( int i = 0; i < (int) p.myValue.size(); ++i )
          {
            res += (X)(Ring) p.myValue[i] * xx;
            xx = xx * myOwner.myEvalPoint;
          }
        return res;
      }
    };
    
public:
    inline operator X() const
    {
        X res = (X) 0;
        myOwner.evaluate( res, EvalFun( *this ) );
        return res;
    }

    inline X operator() () const
    {
      return (X)(*this);
    }
    
  };

  template < int n, typename TRing, typename TOwner, 
             typename TAlloc, typename TX >
  class MPolynomialEvaluatorImpl
  {
  public:
    typedef TRing Ring;
    typedef TOwner Owner;
    typedef TAlloc Alloc;
    typedef TX X;
    typedef MPolynomial< n, Ring, Alloc> MPolyN;
    
    typedef MPolynomial< n - 1, X, 
                         typename Alloc::template rebind<X>::other > MPolyNM1;

    template<int nn, class TT, class AA, class SS>
    friend class MPolynomialEvaluator;
    
    template<int nn, class TT, class HLHL, class AA, class SS>
    friend class MPolynomialEvaluatorImpl;
    
  private:
    const Owner & myOwner; // The "owner"
    const X & myEvalPoint; // The evaluation point on this level
    
    inline 
    MPolynomialEvaluatorImpl(const Owner & owner, const X & evalpoint)
        : myOwner(owner), myEvalPoint(evalpoint)
    {}
    
    template < typename XX, typename Fun >
    class EvalFun
    {
      const MPolynomialEvaluatorImpl<n, Ring, Owner, Alloc, X> & myOwner;
      const Fun & myEvalFun;
      
    public:
      inline 
      EvalFun( const MPolynomialEvaluatorImpl<n, Ring, Owner, Alloc, X> & owner, 
               const Fun & evalfun)
        : myOwner(owner), myEvalFun(evalfun)
      {}
        
      inline XX operator() ( const MPolynomial<n, Ring, Alloc> & p ) const
      {
        XX res = (XX) 0;
        X xx = (X) 1;
        for ( int i = 0; i < (int) p.myValue.size(); ++i )
          {
            res += myEvalFun(p.myValue[i]) * (XX) xx;
            xx = xx * myOwner.myEvalPoint;
          }
        return res;
      }
    };

    template < typename XX, typename Fun >
    inline 
    void evaluate( XX & res, const Fun & evalfun ) const
    {
        // We have to pass evaluation on to our owner, but give a new
        // functor which now evaluates polynomials of type poly<n, T>.
      myOwner.evaluate( res, EvalFun< XX, Fun >( *this, evalfun ) );
    }
    
    class EvalFun2
    {
      const MPolynomialEvaluatorImpl<n, Ring, Owner, Alloc, X> & myOwner;
      
    public:
      inline 
      EvalFun2( const MPolynomialEvaluatorImpl<n, Ring, Owner, Alloc, X> & owner)
        : myOwner(owner)
      {
      }
      
      inline MPolyNM1 operator() (const MPolyN & p) const
        {
          MPolyNM1 res;
          X xx = (X) 1;
          for ( int i = 0; i < (int) p.myValue.size(); ++i )
            {
              res += ( (MPolyNM1) p.myValue[i] ) * xx;
              xx = xx * myOwner.myEvalPoint;
            }
          return res;
        }
    };
    
  public:
    inline
    operator MPolyNM1() const
    // operator poly<n - 1, S, typename Alloc::template rebind<S>::other>() const
    {
      MPolyNM1 res; // missing: determine allocator object
      // We need to pass evaluation on to our owner
      myOwner.evaluate( res, EvalFun2( *this ) );
      return res;
    }
    
    /*
      Continues evaluation with the next indeterminant.  Functor
      returning a "child" evaluator implementation.
    */
    template < typename XX >
    inline 
    MPolynomialEvaluatorImpl
    < n - 1, Ring, 
      MPolynomialEvaluatorImpl< n, Ring, Owner, Alloc, X >,
      Alloc, XX > operator() ( const XX & x ) const
      {
        return MPolynomialEvaluatorImpl
          < n - 1, Ring, MPolynomialEvaluatorImpl< n, Ring, Owner, Alloc, X>, 
            Alloc, XX >( *this, x );
      }
  };

  template < typename TRing, typename TAlloc, typename TX >
  class MPolynomialEvaluator< 1, TRing, TAlloc, TX>
  {
    friend class MPolynomial< 1, TRing, TAlloc >;
  public:
    typedef TRing Ring;
    typedef TAlloc Alloc;
    typedef TX X;
    typedef MPolynomial< 1, Ring, Alloc > MPoly1;

  private:
    const MPoly1 & myPoly; 
    const X & myEvalPoint; 
    
    inline 
    MPolynomialEvaluator( const MPoly1 & poly, const X & evalpoint)
      : myPoly( poly ), myEvalPoint( evalpoint )
    {}
    
  public:
    inline 
    operator X() const
    {
      X res = (X) 0;
      X xx = (X) 1;
      for ( int i = 0; i < (int) myPoly.myValue.size(); ++i )
        {
          res += ( (X) (Ring) myPoly.myValue[i] ) * xx;
          xx = xx * myEvalPoint;
        }
      return res;
    }
    
    inline
    X operator() () const
    {
      return (X) (*this);
    }
  };


  template < int n, typename TRing, typename TAlloc, typename TX >
  class MPolynomialEvaluator
  {
    friend class MPolynomial< n, TRing, TAlloc>;

    template<int nn, class TT, class HLHL, class AA, class SS>
    friend class MPolynomialEvaluatorImpl;

  public:
    typedef TRing Ring;
    typedef TAlloc Alloc;
    typedef TX X;
    typedef MPolynomial< n, Ring, Alloc > MPolyN;
    typedef MPolynomial< n - 1, X, 
                         typename Alloc::template rebind<X>::other > MPolyNM1;
  private:
    const MPolyN & myPoly; 
    const X & myEvalPoint; 
    
    inline 
    MPolynomialEvaluator( const MPolyN & poly, const X & evalpoint)
        : myPoly( poly ), myEvalPoint( evalpoint )
    {}
    
    template < typename XX, typename Fun >
    inline 
    void evaluate( XX & res, const Fun & evalfun ) const
    {
      X xx = (X) 1;
      for ( int i = 0; i < (int) myPoly.myValue.size(); ++i )
        {
          res += evalfun( myPoly.myValue[i] ) * xx;
          xx = xx * myEvalPoint;
        }
    }
    
public:
    inline 
    operator MPolyNM1() const
    {
      MPolyNM1 res( myPoly.getAllocator() );
      X xx = (X) 1;
      for ( int i = 0; i < (int) myPoly.myValue.size(); ++i )
        {
          res += MPolyNM1( myPoly.myValue[i] ) * xx;
          xx = xx * myEvalPoint;
        }
      return res;
    }
    
    template <typename XX>
    inline 
    MPolynomialEvaluatorImpl
    < n - 1, Ring, MPolynomialEvaluator< n, Ring, Alloc, X >, Alloc, XX > 
    operator() (const XX & x) const
    {
      return MPolynomialEvaluatorImpl
        < n - 1, Ring, MPolynomialEvaluator< n, Ring, Alloc, X >, Alloc, XX > 
        ( *this, x );
    }
  };

  template < typename TRing, typename TAlloc >
  class MPolynomial< 0, TRing, TAlloc >
  {
  public:
    typedef TRing Ring;
    typedef TAlloc Alloc;

  private:
    Alloc myAllocator;
    Ring myValue;

  public:
    inline 
    MPolynomial( const Ring & v = 0, const Alloc & allocator = Alloc() )
      : myAllocator(allocator), myValue(v)
    {}
    
    inline 
    MPolynomial( const Alloc & allocator )
      : myAllocator(allocator), myValue( Ring() )
    {}
    
    inline bool isZero() const
    {
      return myValue == 0;
    }
    
    inline operator const Ring & () const
    {
        return myValue;
    }
    
    inline MPolynomial & operator = (const Ring & v)
    {
      myValue = v;
      return *this;
    }

    inline Ring operator() () const
    {
      return myValue;
    }
    
    inline MPolynomial operator * (const Ring & v) const
    {
        return MPolynomial(myValue * v);
    }
    
    inline MPolynomial operator / (const Ring & v) const
    {
        return MPolynomial(myValue / v);
    }
    
    inline MPolynomial operator + (const Ring & v) const
    {
        return MPolynomial(myValue + v);
    }
    
    inline MPolynomial operator - (const Ring & v) const
    {
        return MPolynomial(myValue - v);
    }
    
    inline MPolynomial operator - () const
    {
      return MPolynomial(-myValue);
    }
    
    inline MPolynomial & operator *= (const Ring & v)
    {
      myValue *= v;
      return *this;
    }
    
    inline MPolynomial & operator /= (const Ring & v)
    {
      myValue /= v;
      return *this;
    }
    
    inline MPolynomial & operator += (const Ring & v)
    {
      myValue += v;
      return *this;
    }
    
    inline MPolynomial & operator -= (const Ring & v)
    {
      myValue -= v;
      return *this;
    }
    
    inline bool operator == (const Ring & v) const
    {
        return myValue == v;
    }
    
    inline bool operator != (const Ring & v) const
    {
        return myValue != v;
    }
    
    void selfDisplay( std::ostream & s, int /*N = 0*/ ) const
    {
      s << myValue;
    }
    
    inline void swap( MPolynomial & p )
    {
      std::swap(myValue, p.myValue);
    }
    
    Alloc getAllocator() const
    {
      return myAllocator;
    }
  };


  template < typename T, typename TAlloc = std::allocator<T>, 
             bool usePointers = false>
  class IVector
  {
  public:
    typedef TAlloc Alloc;
    typedef typename std::vector<T, Alloc>::size_type Size;
  private:
    std::vector<T, Alloc> myVec;
    
  public:
    IVector( const Alloc & allocator = Alloc() )
      : myVec(allocator)
    {}
    
    IVector( Size aSize, const Alloc & allocator = Alloc() )
      : myVec( aSize, T(), allocator )
    {}
    
    IVector( Size aSize, const T & entry, const Alloc & allocator = Alloc() )
      : myVec(aSize, entry, allocator)
    {}
    
    Size size() const
    {
      return myVec.size();
    }
    
    void resize( Size aSize, const T & entry = T() )
    {
      myVec.resize(aSize, entry);
    }
    
    const T & operator[] ( Size i ) const
    {
      return myVec[i];
    }
    
    T & operator[] ( Size i )
    {
      return myVec[i];
    }
    
    const T & back() const
    {
      return myVec.back();
    }
    
    T & back()
    {
      return myVec.back();
    }
    
    void swap(IVector & v)
    {
      myVec.swap( v.myVec );
    }
    
    Alloc get_allocator() const
    {
        return myVec.get_allocator();
    }
    
    Alloc getAllocator() const
    {
        return myVec.get_allocator();
    }
  };

  template < typename T, typename TAlloc>
  class IVector< T, TAlloc, true >
  {
  public:
    typedef TAlloc Alloc;
    // typedef typename std::vector<typename Alloc::pointer, typename Alloc::template rebind<typename Alloc::pointer>::other>::size_type Size;
    typedef unsigned long int Size;

  private:
    Alloc myAllocator;
    std::vector<typename Alloc::pointer, typename Alloc::template rebind<typename Alloc::pointer>::other> myVec;
    
    void create( Size begin, Size end, typename Alloc::const_reference entry)
    {
      for (Size i = begin; i < end; ++i)
        {
          myVec[i] = myAllocator.allocate(sizeof(T));
          myAllocator.construct(myVec[i], entry);
        }
    }
    
    void free(Size begin, Size end)
    {
      for (Size i = begin; i < end; ++i)
        {
          myAllocator.destroy(myVec[i]);
          myAllocator.deallocate(myVec[i], sizeof(T));
        }
    }
    
    template<class A>
    void copy_from(const std::vector<typename Alloc::pointer, A> & source)
    {
      for (Size i = 0; i < myVec.size(); ++i)
        {
          myVec[i] = myAllocator.allocate(sizeof(T));
          myAllocator.construct(myVec[i], *source[i]);
        }
    }
    
  public:
    IVector(const Alloc & allocator = Alloc())
      : myAllocator(allocator), myVec(allocator)
    {}
    
    IVector(Size aSize, const Alloc & allocator = Alloc())
      : myAllocator(allocator), myVec(aSize, 0, allocator)
    {
      create(0, aSize, T());
    }
    
    IVector(Size aSize, const T & entry, const Alloc & allocator = Alloc())
      : myAllocator(allocator), myVec(aSize, 0, allocator)
    {
      create(0, aSize, entry);
    }
    
    IVector(const IVector & v)
      : myVec(v.size())
    {
      copy_from(v.myVec);
    }
    
    ~IVector()
    {
      free(0, (Size)myVec.size());
    }
    
    IVector & operator = (const IVector & v)
    {
      if (&v != this)
        {
          free(0, (Size)myVec.size());
          myVec.resize(v.size());
          copy_from(v.myVec);
        }
      return *this;
    }
    
    Size size() const
    {
      return (Size)myVec.size();
    }
    
    void resize(Size aSize, const T & entry = T())
    {
      Size oldsize = (Size)myVec.size();
      if (oldsize > aSize)
        free(aSize, oldsize);
      myVec.resize(aSize);
      if (oldsize < aSize)
        create(oldsize, aSize, entry);
    }
    
    const T & operator[] (Size i) const
    {
      return *myVec[i];
    }
    
    T & operator[] (Size i)
    {
      return *myVec[i];
    }
    
    const T & back() const
    {
      return *myVec.back();
    }
    
    T & back()
    {
      return *myVec.back();
    }
    
    void swap(IVector & v)
    {
      myVec.swap(v.myVec);
    }
    
    Alloc get_allocator() const
    {
      return myVec.get_allocator();
    }

    Alloc getAllocator() const
    {
      return myVec.get_allocator();
    }
  };
  

  // template class MPolynomial
  template < int n, typename TRing, class TAlloc >
  class MPolynomial
  {
    // ----------------------- friends ---------------------------------------
    
    template<int NN, int nn, typename TT, typename AA>
    friend class MPolynomialDerivativeComputer;
    
    friend void euclidDiv<TRing, TAlloc>
    ( const MPolynomial<1, TRing, TAlloc> &, 
      const MPolynomial<1, TRing, TAlloc> &, 
      MPolynomial<1, TRing, TAlloc> &, 
      MPolynomial<1, TRing, TAlloc> & );
    
    template<int nn, typename TT, typename AA, typename SS>
    friend class MPolynomialEvaluator;
    
    template<int nn, typename TT, typename HLHL, typename AA, typename SS>
    friend class MPolynomialEvaluatorImpl;

  public:
    typedef TRing Ring;
    typedef TAlloc Alloc;
    typedef MPolynomial< n - 1, Ring, Alloc > MPolyNM1; 
    typedef IVector< MPolyNM1, 
                     typename Alloc::template rebind<MPolyNM1 >::other, (n>1) >
    Storage;
    typedef typename Storage::Size Size;

    // ----------------------- Datas  ----------------------------
private:
    static MPolyNM1 myZeroPolynomial;
    Storage myValue;

    inline 
      MPolynomial( bool, Size s, const Alloc & /*allocator*/)
      : myValue(s) 
    {}
    
public:
    inline void normalize()
    {
      Size dp1 = myValue.size();
      while ( dp1 )
        {
          if (myValue[dp1 - 1].isZero())
            --dp1;
          else
            break;
        }
      if (dp1 < myValue.size())
        myValue.resize(dp1);
    }

    // ----------------------- Standard services ----------------------------
  public:
      
    inline MPolynomial( const Alloc & allocator = Alloc() )
      : myValue( allocator )
    {}
    
    inline MPolynomial( const Ring & v, const Alloc & allocator = Alloc() )
      : myValue( 1, MPolyNM1( v ), allocator )
    {}
    
    inline MPolynomial( const MPolyNM1 & pdm1, 
                        const Alloc & allocator = Alloc() )
      : myValue( 1, pdm1 )
    {}

    template < typename Ring2, typename Alloc2 >
    inline MPolynomial( const MPolynomial<n, Ring2, Alloc2> & p, 
                        const Alloc & allocator = Alloc() )
      : myValue( p.degree() + 1, MPolyNM1(), allocator )
    {
      for ( Size i = 0; i < myValue.size(); ++i )
        myValue[i] = p[i];
      normalize();
    }
    
    template < typename Ring2, typename Alloc2 >
    inline 
    MPolynomial & operator= 
    (const MPolynomial< n, Ring2, Alloc2> & p )
    {
      myValue.resize(p.degree() + 1);
      for ( Size i = 0; i < myValue.size(); ++i )
        myValue[i] = p[i];
      normalize();
      return *this;
    }

    inline void swap( MPolynomial & p )
    {
        myValue.swap(p.myValue);
    }
    
    Alloc getAllocator() const
    {
        return myValue.getAllocator();
    }

    // ----------------------- MPolynomial services ----------------------------
  public:

    inline int degree() const
    {
      return (int)(myValue.size() - 1);
    }
    
    inline const MPolyNM1 & leading() const
    {
      return myValue.size() ? myValue.back() : myZeroPolynomial;
    }
    
    inline bool isZero() const
    {
      return myValue.size() == 0;
    }
    
    inline MPolyNM1 & operator[] ( Size i )
    {
      if (i >= myValue.size())
        myValue.resize(i + 1);
      return myValue[i];
    }
    
    inline const MPolyNM1 & operator[] (Size i) const
    {
      return i < myValue.size() ? myValue[i] : myZeroPolynomial;
    }

    inline 
    MPolynomialEvaluator<n, Ring, Alloc, Ring> operator() (const Ring & x) const
    {
      return MPolynomialEvaluator<n, Ring, Alloc, Ring>( *this, x );
    }

    template <typename Ring2>
    inline 
    MPolynomialEvaluator<n, Ring, Alloc, Ring2> operator() (const Ring2 & x) const
    {
      return MPolynomialEvaluator<n, Ring, Alloc, Ring2>( *this, x );
    }
    
    // ----------------------- Operators --------------------------------------
  public:

    inline MPolynomial operator * (const Ring & v) const
    {
      MPolynomial r(*this);
      for ( Size i = 0; i < myValue.size(); ++i )
        r[i] *= v;
      return r;
    }
    
    inline MPolynomial operator / (const Ring & v) const
    {
      MPolynomial r(*this);
      for ( Size i = 0; i < myValue.size(); ++i )
        r[i] /= v;
      return r;
    }

    inline MPolynomial & operator *= (const MPolynomial & p)
    {
      return *this = *this * p;
    }
    
    inline MPolynomial & operator *= (const Ring & v)
    {
      MPolynomial r( *this );
      for ( Size i = 0; i < myValue.size(); ++i )
        myValue[i] *= v;
      return *this;
    }
    
    inline MPolynomial & operator /= (const Ring & v)
    {
      for ( Size i = 0; i < myValue.size(); ++i )
        myValue[i] /= v;
      return *this;
    }
    
    friend 
    inline MPolynomial operator * (const Ring & v, const MPolynomial & p)
    {
      MPolynomial r(p);
      for ( Size i = 0; i < p.myValue.size(); ++i )
        r[i] *= v;
      return r;
    }
    
    inline MPolynomial operator - () const
    {
      MPolynomial r( true, myValue.size(), getAllocator() );
      for ( Size i = 0; i < myValue.size(); ++i )
        r[i] = -myValue[i];
      return r;
    }
    
    inline MPolynomial operator + (const MPolynomial & q) const
    {
      MPolynomial r(*this);
      if (r.myValue.size() < q.myValue.size())
        r.myValue.resize(q.myValue.size());
      for ( Size i = 0; i < q.myValue.size(); ++i )
        r[i] += q[i];
      r.normalize();
      return r;
    }
    
    inline MPolynomial operator - (const MPolynomial & q) const
    {
      MPolynomial r(*this);
      if (r.myValue.size() < q.myValue.size())
        r.myValue.resize(q.myValue.size());
      for ( Size i = 0; i < q.myValue.size(); ++i )
        r[i] -= q[i];
      r.normalize();
      return r;
    }
    
    inline MPolynomial & operator += (const MPolynomial & q)
    {
      if (myValue.size() < q.myValue.size())
        myValue.resize(q.myValue.size());
      for ( Size i = 0; i < q.myValue.size(); ++i )
        myValue[i] += q[i];
      normalize();
      return *this;
    }
    
    inline MPolynomial & operator -= (const MPolynomial & q)
    {
      if (myValue.size() < q.myValue.size())
        myValue.resize(q.myValue.size());
      for ( Size i = 0; i < q.myValue.size(); ++i )
        myValue[i] -= q[i];
      normalize();
      return *this;
    }
    
    inline MPolynomial operator + (const MPolyNM1 & q) const
    {
      MPolynomial r(*this);
      if (r.myValue.size() < 1)
        r.myValue.resize(1);
      r[0] += q;
      r.normalize();
      return r;
    }
    
    inline MPolynomial operator - (const MPolyNM1 & q) const
    {
      MPolynomial r(*this);
      if (r.myValue.size() < 1)
        r.myValue.resize(1);
      r[0] -= q;
      r.normalize();
      return r;
    }
    
    inline MPolynomial operator + (const Ring & v) const
    {
      MPolynomial r(*this);
      if (r.myValue.size() < 1)
        r.myValue.resize(1);
      r[0] += v;
      r.normalize();
      return r;
    }

    inline MPolynomial operator - (const Ring & v) const
    {
      MPolynomial r(*this);
      if (r.myValue.size() < 1)
        r.myValue.resize(1);
      r[0] -= v;
      r.normalize();
      return r;
    }

    inline MPolynomial & operator += (const MPolyNM1 & q)
    {
      if (myValue.size() < 1)
        myValue.resize(1);
      myValue[0] += q;
      normalize();
      return *this;
    }
    
    inline MPolynomial & operator -= (const MPolyNM1 & q)
    {
      if (myValue.size() < 1)
        myValue.resize(1);
      myValue[0] -= q;
      normalize();
      return *this;
    }
    
    inline MPolynomial & operator += (const Ring & v)
    {
      if (myValue.size() < 1)
        myValue.resize(1);
      myValue[0] += v;
      normalize();
      return *this;
    }
    
    inline MPolynomial & operator -= (const Ring & v)
    {
      if (myValue.size() < 1)
        myValue.resize(1);
      myValue[0] -= v;
      normalize();
      return *this;
    }

    inline MPolynomial operator * (const MPolynomial & p) const
    {
      int d = p.degree() + degree();
      if (d < -1)
        d = -1;
      MPolynomial r( true, d + 1, getAllocator() );
      if (!isZero() && !p.isZero())
        for ( Size i = 0; i < r.myValue.size(); ++i )
          for ( Size j = 0; j < myValue.size(); ++j )
            if (i < j + p.myValue.size())
              r[i] += myValue[j] * p[i - j];
      r.normalize();
      return r;
    }

    // ----------------------- Comparison operators ---------------------------
  public:

    inline bool operator == (const MPolynomial & q) const
    {
      if (myValue.size() != q.myValue.size())
        return false;
      for (Size i = 0; i < myValue.size(); ++i)
        if (myValue[i] != q[i])
          return false;
      return true;
    }
    
    inline bool operator != (const MPolynomial & q) const
    {
      return !(*this == q);
    }
    
    inline bool operator == (const Ring & v) const
    {
      if ((v == 0) && (myValue.size() == 0))
        return true;
      if (myValue.size() != 1)
        return false;
      return myValue[0] == v;
    }

    inline bool operator != (const Ring & v) const
    {
      return !(*this == v);
    }
    

      
    // ----------------------- Interface --------------------------------------
  public:

    void selfDisplay ( std::ostream & s, int N = n ) const
    {
      if (isZero())
        s << (Ring) 0;
      else
        {
          Size nonzero = 0;
          for (Size i = 0; i < myValue.size(); ++i)
            if (!myValue[i].isZero())
              ++nonzero;
          if (nonzero > 1) s << "(";
          bool first = true;
          for (Size i = 0; i < myValue.size(); ++i)
            if (!myValue[i].isZero())
              {
                if (first) first = false;
                else       s << " + ";
                myValue[i].selfDisplay(s, N);
                if (i > 0)
                  {
                    s << " ";
                    s << "X_" << N - n;
                    if (i > 1) s << "^" << i;
                  }
              }
          if (nonzero > 1)
            s << ")";
        }
    }
    
    bool isValid() const;

     
    // ------------------------- Datas --------------------------------------
  private:
    // ------------------------- Hidden services ----------------------------
  protected:
    

  }; // end of class MPolynomial


  template <int N, typename TRing, class TAlloc>
  std::ostream&
  operator<< ( std::ostream & out, 
               const MPolynomial<N, TRing, TAlloc> & object );


  // ------------------------- monomial services ----------------------------

  template <int n, typename Ring, typename Alloc>
  class Xe_kComputer
  {
  public:
    Xe_kComputer() {}

    MPolynomial<n, Ring, Alloc> 
    get( unsigned int k, unsigned int e )
    {
      MPolynomial<n, Ring, Alloc> p;
      if ( k == 0 )
        p[e] = Xe_kComputer<n-1,Ring,Alloc>().get( k-1, e ); 
      else
        p[0] = Xe_kComputer<n-1,Ring,Alloc>().get( k-1, e );
      p.normalize();
      //std::cerr << "Xe_k(" << k << "," << e << ")=" << p << std::endl;
      return p;
    }

  };

  template <typename Ring, typename Alloc>
  class Xe_kComputer<0,Ring,Alloc>
  {
  public:
    Xe_kComputer() {}

    MPolynomial<0, Ring, Alloc> 
      get( unsigned int /*k*/, unsigned int /*e*/ )
    {
      MPolynomial<0, Ring, Alloc> p = 1;
      return p;
    }
  };

  template <int n, typename Ring, typename Alloc>
  inline 
  MPolynomial<n, Ring, Alloc> 
  Xe_k( unsigned int k, unsigned int e )
  {
    return Xe_kComputer<n, Ring, Alloc>().get( k, e );
  }

  template <int n, typename Ring>
  inline 
  MPolynomial<n, Ring, std::allocator<Ring> >
  Xe_k( unsigned int k, unsigned int e )
  {
    return Xe_kComputer<n, Ring, std::allocator<Ring> >().get( k, e );
  }


  template <typename Ring, typename Alloc>
  inline 
  MPolynomial<1, Ring, Alloc> 
  mmonomial(unsigned int e)
  {
    MPolynomial<1, Ring, Alloc> p;
    p[e] = 1;
    return p;
  }

  template <typename Ring, typename Alloc>
  inline 
  MPolynomial<2, Ring, Alloc> 
  mmonomial(unsigned int e, unsigned int f)
  {
    MPolynomial<2, Ring, Alloc> p;
    p[e][f] = 1;
    return p;
  }

  template <typename Ring, typename Alloc>
  inline MPolynomial<3, Ring, Alloc> 
  mmonomial(unsigned int e, unsigned int f, unsigned int g)
  {
    MPolynomial<3, Ring, Alloc> p;
    p[e][f][g] = 1;
    return p;
  }

  template <typename Ring, typename Alloc>
  inline 
  MPolynomial<4, Ring, Alloc> 
  mmonomial(unsigned int e, unsigned int f, unsigned int g, unsigned int h)
  {
    MPolynomial<4, Ring, Alloc> p;
    p[e][f][g][h] = 1;
    return p;
  }

  template <typename Ring>
  inline 
  MPolynomial<1, Ring, std::allocator<Ring> > 
  mmonomial(unsigned int e)
  {
    MPolynomial<1, Ring, std::allocator<Ring> > p;
    p[e] = 1;
    return p;
  }

  template <typename Ring>
  inline
  MPolynomial<2, Ring, std::allocator<Ring> > 
  mmonomial(unsigned int e, unsigned int f)
  {
    MPolynomial<2, Ring, std::allocator<Ring> > p;
    p[e][f] = 1;
    return p;
  }

  template <typename Ring>
  inline
  MPolynomial<3, Ring, std::allocator<Ring> > 
  mmonomial(unsigned int e, unsigned int f, unsigned int g)
  {
    MPolynomial<3, Ring, std::allocator<Ring> > p;
    p[e][f][g] = 1;
    return p;
  }

  template <typename Ring>
  inline
  MPolynomial<4, Ring, std::allocator<Ring> >
  mmonomial
  (unsigned int e, unsigned int f, unsigned int g, unsigned int h)
  {
    MPolynomial<4, Ring, std::allocator<Ring> > p;
    p[e][f][g][h] = 1;
    return p;
  }
  

  template <int n, typename Ring, class Alloc>
  class MPolynomialDerivativeComputer<0, n, Ring, Alloc>
  {
  public:
    typedef MPolynomial<n, Ring, Alloc> MPolyN;

    static inline 
    void computeDerivative( const MPolyN & src, MPolyN & dest)
    {
      dest.myValue.resize(src.degree() >= 0 ? src.degree() : 0);
      for ( int i = 1; i <= src.degree(); ++i )
        dest[i - 1] = src[i] * (Ring)i;
      dest.normalize();
    }
  };

  template <int N, int n, typename Ring, typename Alloc>
  class MPolynomialDerivativeComputer
  {
  public:
    typedef MPolynomial<n, Ring, Alloc> MPolyN;

    static inline
    void computeDerivative(const MPolyN & src, MPolyN & dest)
    {
      dest.myValue.resize(src.degree() + 1);
      for ( int i = 0; i <= src.degree(); ++i )
        MPolynomialDerivativeComputer<N - 1, n - 1, Ring, Alloc>
          ::computeDerivative( src[i], dest[i] );
      dest.normalize();
    }
  };

  template<typename Ring, class Alloc>
  class MPolynomialDerivativeComputer<0, 0, Ring, Alloc>
  {
  public:
    class ERROR_N_must_be_strictly_less_than_n_in_derivative_template;
    typedef MPolynomial<0, Ring, Alloc> MPoly0;

    static inline 
    void computeDerivative(const MPoly0 &, MPoly0 &)
    {
      ERROR_N_must_be_strictly_less_than_n_in_derivative_template();
    }
  };
  
  template<int N, typename Ring, class Alloc>
  class MPolynomialDerivativeComputer<N, 0, Ring, Alloc>
  {
  public:
    class ERROR_N_must_be_strictly_less_than_n_in_derivative_template;
    typedef MPolynomial<0, Ring, Alloc> MPoly0;

    static inline
    void computeDerivative(const MPoly0 &, MPoly0 &)
    {
      ERROR_N_must_be_strictly_less_than_n_in_derivative_template();
    }
  };

  template <int N, int n, typename Ring, typename Alloc>
  inline 
  MPolynomial<n, Ring, Alloc> 
  derivative( const MPolynomial<n, Ring, Alloc> & p )
  {
    MPolynomial<n, Ring, Alloc> res( p.getAllocator() );
    MPolynomialDerivativeComputer<N, n, Ring, Alloc>::computeDerivative(p, res);
    return res;
  }

  template<int N, int n, typename Ring>
  inline
  MPolynomial<n, Ring, std::allocator<Ring> > 
  derivative
  (const MPolynomial<n,Ring,std::allocator<Ring> > & p)
   {
     MPolynomial<n, Ring, std::allocator<Ring> > res( p.getAllocator() );
     MPolynomialDerivativeComputer<N, n, Ring, std::allocator<Ring> >
       ::computeDerivative( p, res );
     return res;
   }
   
   template<typename Ring, typename Alloc>
   void 
   euclidDiv( const MPolynomial<1, Ring, Alloc> & f, 
              const MPolynomial<1, Ring, Alloc> & g,
              MPolynomial<1, Ring, Alloc> & q, 
              MPolynomial<1, Ring, Alloc> & r )
   {
     if (f.degree() < g.degree())
       {
         // Ignore the trivial case
         q = MPolynomial<1, Ring, Alloc>(f.getAllocator());
         r = f;
         return;
       }
    q = MPolynomial<1, Ring, Alloc>( true, f.degree() - g.degree() + 1, 
                                     f.getAllocator() );
    r = f;
    for (int i = q.degree(); i >= 0; --i)
      {
        q[i] = r[i + g.degree()] / g.leading();
        for (int j = g.degree(); j >= 0; --j)
          r[i + j] -= q[i] * g[j];
      }
    r.normalize();
    // Note that the degree of q is already correct.
   }

   template <typename Ring>
   void 
   euclidDiv( const MPolynomial<1, Ring, std::allocator<Ring> > & f,
              const MPolynomial<1, Ring, std::allocator<Ring> > & g, 
              MPolynomial<1, Ring, std::allocator<Ring> > & q, 
              MPolynomial<1, Ring, std::allocator<Ring> > & r )
   {
     euclidDiv<Ring, std::allocator<Ring> >(f, g, q, r);
   }

   template<typename Ring, typename Alloc>
   MPolynomial<1, Ring, Alloc> 
   gcd( const MPolynomial<1, Ring, Alloc> & f, 
        const MPolynomial<1, Ring, Alloc> & g )
   {
     if (f.isZero())
       {
         if (g.isZero()) return f; // both are zero
         else            return g / g.leading(); // make g monic
       }
     MPolynomial<1, Ring, Alloc> 
       d1(f / f.leading()), 
       d2(g / g.leading()), 
       q(f.getAllocator()), 
       r(f.getAllocator());
     while (!d2.isZero())
       {
         euclidDiv(d1, d2, q, r);
         d1.swap(d2);
         d2 = r;
         d2 /= r.leading(); // make r monic
       }
     return d1;
   }
  
  template<typename Ring>
  MPolynomial<1, Ring, std::allocator<Ring> > 
  gcd( const MPolynomial<1, Ring, std::allocator<Ring> > & f, 
       const MPolynomial<1, Ring, std::allocator<Ring> > & g )
  {
    return gcd<Ring, std::allocator<Ring> >(f, g);
  }

} // namespace DGtal


// Includes inline functions.
#include "DGtal/math/MPolynomial.ih"

//                                                                           //

#endif // !defined MPolynomial_h

#undef MPolynomial_RECURSES
#endif // else defined(MPolynomial_RECURSES)