// Rascal, the Advanced Scientific CALculator
// Copyright (C) 2001, 2002, Sebastian Ritterbusch (Rascal@Ritterbusch.de)
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 2
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more detauls.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
//

#ifndef POLY_HPP_INCLUDED
#define POLY_HPP_INCLUDED

#include <iostream.h>

template <class T>
class mpoly
{
   public:
      T  *a;
      int N;
      
      mpoly(void) { a=0;N=0; }
      mpoly(const T & c) { a=new T[N=1];a[0]=c; }
      mpoly(const T & c, const T & x) { a=new T[N=2];a[0]=c;a[1]=x; }
      
      mpoly(const mpoly<T> & b)
      {
         a=new T[N=b.N];
         int i;
         for(i=0;i<N;i++)
            a[i]=b.a[i];
         for(;N>0 && a[N-1]==T(0);N--);
      }
     ~mpoly(void) { delete [] a; }
      
      mpoly<T> & operator =(const mpoly<T> & c)
      {
         T *n;
         if(c.N==N)
            n=a;
         else
            n=new T[N=c.N];
         int i;
         for(i=0;i<N;i++)
            n[i]=c.a[i];
         if(n!=a)
         {
            delete [] a;
            a=n;
         }
         for(;N>0 && a[N-1]==T(0);N--);
         return *this;
      }
      int   deg(void) const
      {
         int n=N;
         for(;n>0 && a[n-1]==T(0);n--);
         return n>0?n-1:0;
      }
      int   deg(void)
      {
         for(;N>0 && a[N-1]==T(0);N--);
         return N>0?N-1:0;
      }
      mpoly<T> & foreach(T (*f)(const T &))
      {
         int i;
         for(i=0;i<N;i++)
            a[i]=f(a[i]);
         return *this;
      }
      mpoly<T> & foreach(T (*f)(T))
      {
         int i;
         for(i=0;i<N;i++)
            a[i]=f(a[i]);
         return *this;
      }
      mpoly<T> foreach(T (*f)(const T &)) const
      {
         int i;
         mpoly<T> temp(*this);
         for(i=0;i<N;i++)
            temp.a[i]=f(a[i]);
         return temp;
      }
      mpoly<T> foreach(T (*f)(T)) const
      {
         int i;
         mpoly<T> temp(*this);
         for(i=0;i<N;i++)
            temp.a[i]=f(a[i]);
         return temp;
      }
   private:
      void ensure(int n)
      {
         if(N<n)
         {
            T *u=new T[n];
            int i;
            for(i=0;i<N;i++)
               u[i]=a[i];
            for(;i<n;i++)
               u[i]=T(0);
            delete [] a;
            a=u;
            N=n;
         }
      }
   public:
      mpoly<T> & operator +=(const mpoly<T> & b)
      {
         ensure(b.N);
         int i;
         for(i=0;i<N && i<b.N;i++)
            a[i]+=b.a[i];
         for(;N>0 && a[N-1]==T(0);N--);
         return *this;
      }
      mpoly<T> & operator -=(const mpoly<T> & b)
      {
         ensure(b.N);
         int i;
         for(i=0;i<N && i<b.N;i++)
            a[i]-=b.a[i];
         for(;N>0 && a[N-1]==T(0);N--);
         return *this;
      }
      mpoly<T> & operator *=(const mpoly<T> & b)
      {
         int n;
         T *re=new T[n=b.N+N-1];
         int i,j;
         
         for(i=0;i<n;i++) re[i]=0;
         
         for(i=N-1;i>=0;i--)
            for(j=b.N-1;j>=0;j--)
               re[i+j]=re[i+j]+a[i]*b.a[j];
         delete [] a;
         a=re;
         N=n;
         for(;N>0 && a[N-1]==T(0);N--);
         return *this;
      }
      
      T operator()(const value & x) const
      {
         T res=a[N-1];
         int i;
         for(i=N-2;i>=0;i--)
            res=res*x+a[i];

         return res;
      }
      
      friend mpoly<T> operator +(mpoly<T> a,const mpoly<T> &b) { return a+=b; }
      friend mpoly<T> operator -(mpoly<T> a,const mpoly<T> &b) { return a-=b; }
      friend mpoly<T> operator *(mpoly<T> a,const mpoly<T> &b) { return a*=b; }
      friend bool     operator ==(const mpoly<T> & a,const mpoly<T> & b)
      {
         int i;
         for(i=0;i<a.N && i<b.N;i++)
            if(a.a[i]!=b.a[i])
               return false;
         for(;i<a.N;i++)
            if(a.a[i]!=T(0))
               return false;
         for(;i<b.N;i++)
            if(b.a[i]!=T(0))
               return false;
         return true;
      }
      friend bool     operator !=(const mpoly<T> & a,const mpoly<T> & b)
      {
         int i;
         for(i=0;i<a.N && i<b.N;i++)
            if(a.a[i]!=b.a[i])
               return true;
         for(;i<a.N;i++)
            if(a.a[i]!=T(0))
               return true;
         for(;i<b.N;i++)
            if(b.a[i]!=T(0))
               return true;
         return false;
      }      
      friend mpoly<T> power(mpoly<T> a,int e)
      {
         mpoly<T> res(T(1));
         while(e)
         {
            if(e&1)
               res*=a;
            e>>=1;
            if(e)
               a*=a;
         }  
         return res;
      }
      friend mpoly<T> deriv(mpoly<T> a)
      {
         int i;
         for(i=1;i<a.N;i++)
            a.a[i-1]=i*a.a[i];
         if(a.N>0)
            a.N--;
         return a;
      }      
};

#endif