// Rascal, the Advanced Scientific CALculator
// Copyright (C) 2001, 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 TAYLOR_HPP_INCLUDED
#define TAYLOR_HPP_INCLUDED

#include <iostream.h>

template <class T>
class mtaylor
{
   private:
      int   n;
      T     *a;
   public:
      mtaylor(void) : n(0),a(0) {}
      mtaylor(T x,int iN=1) :a(new T[iN])
      {
         if(a)
            n=iN;
         else
            n=0;
      
         int i; 
         if(n) 
            a[0]=x; 
         for(i=1;i<n;i++) 
            a[i]=0; 
      }
      mtaylor(T x,T y,int iN=2) : a(new T[iN])
      { 
         if(a)
            n=iN;
         else
            n=0;
            
         int i; 
         if(n) 
            a[0]=x; 
         if(n>1) 
            a[1]=y; 
         for(i=2;i<n;i++) 
            a[i]=0; 
      }
      mtaylor(const mtaylor<T> & b) : a(new T[b.n])
      { 
         if(a)
            n=b.n;
         else
            n=0;
         int i; 
         for(i=0;i<n;i++) 
            a[i]=b.a[i]; 
      }
      
      ~mtaylor(void) { if(a) delete [] a; }
      
      mtaylor<T> & operator =(const mtaylor<T> & b)
      {
         T * new_a=new T[b.n];

         n=new_a?b.n:0;

         int i;
         for(i=0;i<n;i++)
            new_a[i]=b.a[i];
            
         delete [] a;
            
         a=new_a;
         
         return *this;
      }
      
      void minsize(int nn)
      {
         if(nn<=n)
            return;
            
         T * new_a=new T[nn];
         
         int i;
         for(i=0;i<n;i++)
            new_a[i]=a[i];
         for(;i<nn;i++)
            new_a[i]=0;
            
         delete [] a;
            
         a=new_a;
         
         n=new_a?nn:0;
      }
      
      friend inline mtaylor<T> & operator +=(mtaylor<T> & a,const mtaylor<T> & b) { int i; a.minsize(b.n); for(i=0;i<a.n && i<b.n;i++) a.a[i]=a.a[i]+b.a[i]; return a;}
      friend inline mtaylor<T> & operator +=(mtaylor<T> & a,const T & b) { if(a.n) a.a[0]=a.a[0]+b; return a;}
      friend inline mtaylor<T> operator +(mtaylor<T> a,const mtaylor<T> & b) { return a+=b; }
      friend inline mtaylor<T> operator +(mtaylor<T> a,const T & b) { return a+=b; }
      friend inline mtaylor<T> operator +(const T & b,mtaylor<T> a) { return a+=b; }
   
      friend inline mtaylor<T> operator -(mtaylor<T> a) { int i; for(i=0;i<a.n;i++) a.a[i]=-a.a[i]; return a;}
      
      friend inline mtaylor<T> & operator -=(mtaylor<T> & a,const mtaylor<T> & b) { int i; a.minsize(b.n);for(i=0;i<a.n && i<b.n;i++) a.a[i]=a.a[i]-b.a[i]; return a;}
      friend inline mtaylor<T> & operator -=(mtaylor<T> & a,const T & b) { if(a.n) a.a[0]=a.a[0]-b; return a;}
      friend inline mtaylor<T> operator -(mtaylor<T> a,const mtaylor<T> & b) { return a-=b; }
      friend inline mtaylor<T> operator -(mtaylor<T> a,const T & b) { return a-=b; }
      friend inline mtaylor<T> operator -(const T & b,mtaylor<T> a) { return -(a-=b); }

      friend inline mtaylor<T> & operator *=(mtaylor<T> & a,const mtaylor<T> & b) 
      {  
         int i,j; 
         a.minsize(b.n);
         for(i=a.n-1;i>=0;i--)
         {
            a.a[i]=a.a[i]*b.a[0];
            for(j=1;j<=i && j<b.n;j++)
               a.a[i]=a.a[i]+a.a[i-j]*b.a[j]; 
         }   
         return a;
      }
      friend inline mtaylor<T> & operator *=(mtaylor<T> & a,const T & b) { int i;for(i=0;i<a.n;i++) a.a[i]=a.a[i]*b; return a;}
      friend inline mtaylor<T> operator *(mtaylor<T> a,const mtaylor<T> & b) { return a*=b; }
      friend inline mtaylor<T> operator *(mtaylor<T> a,const T & b) { return a*=b; }
      friend inline mtaylor<T> operator *(const T & b,mtaylor<T> a) { return a*=b; }

      friend inline mtaylor<T> & operator /=(mtaylor<T> & a,const mtaylor<T> & b) 
      {  
         int i,j; 
         a.minsize(b.n);
         for(i=0;i<a.n;i++)
         {
            for(j=1;j<=i && j<b.n;j++)
               a.a[i]=a.a[i]-a.a[i-j]*b.a[j]; 
            a.a[i]=a.a[i]/b.a[0];
         }   
         return a;
      }
      friend inline mtaylor<T> & operator /=(mtaylor<T> & a,const T & b) { int i; for(i=0;i<a.n;i++) a.a[i]=a.a[i]/b; return a;}
      friend inline mtaylor<T> operator /(mtaylor<T> a,const mtaylor<T> & b) { return a/=b; }
      friend inline mtaylor<T> operator /(mtaylor<T> a,const T & b) { return a/=b; }
      friend inline mtaylor<T> operator /(const T & b,mtaylor<T> a) { return mtaylor<T>(b,a.n)/a; }
      
      friend inline ostream & operator <<(ostream & o,const mtaylor<T> & a) { int i; for(i=0;i<a.n-1;i++) o << a.a[i] << ' '; return (o << a.a[a.n-1]); }
      
      inline T   operator ()(int i) const { return a[i]; }
      inline T & operator ()(int i)       { return a[i]; }
      inline int dim(void) const { return n; }

      
      friend inline int operator==(const mtaylor<T> & a,const mtaylor<T> & b) 
      { 
         return (a.n!=0) && (b.n!=0) && (a.a[0]==b.a[0]);
      }
      friend inline int operator!=(const mtaylor<T> & a,const mtaylor<T> & b) 
      { 
         return (a.n==0) || (b.n==0) || (a.a[0]!=b.a[0]);
      }
      
      friend inline mtaylor<T> sqrt(const mtaylor<T> & a)
      {
         if(a.n<=0)
            return a; // nonsense
         mtaylor<T> result(sqrt(a.a[0]),a.n);
         int   k,j;
         T     factor=T(1)/(T(2)*result.a[0]);
         for(k=1;k<a.n;k++)
         {
            result.a[k]=a.a[k];
            if(k&1) // odd k
            {
               for(j=1;j<=(k-1)/2;j++)
                  result.a[k]=result.a[k]-T(2)*result.a[j]*result.a[k-j];

            } else  // even k
            {
               result.a[k]=result.a[k]-result.a[k/2]*result.a[k/2];
               for(j=1;j<k/2;j++)
                  result.a[k]=result.a[k]-T(2)*result.a[j]*result.a[k-j];
            }
            result.a[k]=result.a[k]*factor;
         }
         return result;
      }

      friend inline mtaylor<T> sqr(const mtaylor<T> & a)
      {
         if(a.n<=0)
            return a; // nonsense
         mtaylor<T> result(sqr(a.a[0]),a.n);
         int   k,j;
         for(k=1;k<a.n;k++)
         {
            if(k&1) // odd k
            {
               result.a[k]=a.a[0]*a.a[k];
               for(j=1;j<=(k-1)/2;j++)
                  result.a[k]=result.a[k]+a.a[j]*a.a[k-j];
               result.a[k]=T(2)*result.a[k];
            } else  // even k
            {
               result.a[k]=sqr(a.a[k/2]);
               for(j=0;j<k/2;j++)
                  result.a[k]=result.a[k]+T(2)*a.a[j]*a.a[k-j];
            }
         }
         return result;
      }

      friend inline mtaylor<T> exp(const mtaylor<T> & a)
      {
         if(a.n<=0)
            return a; // nonsense
         mtaylor<T> result(exp(a.a[0]),a.n);
         int   k,j;
         for(k=1;k<a.n;k++)
         {
            result.a[k]=result.a[0]*T(k)*a.a[k];
            for(j=1;j<k;j++)
               result.a[k]=result.a[k]+result.a[j]*T(k-j)*a.a[k-j];
               
            result.a[k]=result.a[k]/T(k);
         }
         return result;
      }
      
      friend inline mtaylor<T> sin(const mtaylor<T> & a)
      {
         if(a.n<=0)
            return a; // nonsense
         mtaylor<T> sresult(sin(a.a[0]),a.n);
         mtaylor<T> cresult(cos(a.a[0]),a.n);
         int   k,j;
         for(k=1;k<a.n;k++)
         {
            sresult.a[k]=cresult.a[0]*a.a[k]*T(k);
            for(j=1;j<k;j++)
               sresult.a[k]=sresult.a[k]+cresult.a[j]*T(k-j)*a.a[k-j];
               
            sresult.a[k]=sresult.a[k]/T(k);
            
            if(k<a.n-1)
            {
               cresult.a[k]=sresult.a[0]*a.a[k]*T(k);
               for(j=1;j<k;j++)
                  cresult.a[k]=cresult.a[k]+sresult.a[j]*T(k-j)*a.a[k-j];
               
               cresult.a[k]=cresult.a[k]/T(-k);
            }
         }
         return sresult;
      }

      friend inline mtaylor<T> cos(const mtaylor<T> & a)
      {
         if(a.n<=0)
            return a; // nonsense
         mtaylor<T> sresult(sin(a.a[0]),a.n);
         mtaylor<T> cresult(cos(a.a[0]),a.n);
         int   k,j;
         for(k=1;k<a.n;k++)
         {
            if(k<a.n-1)
            {
               sresult.a[k]=cresult.a[0]*a.a[k]*T(k);
               for(j=1;j<k;j++)
                  sresult.a[k]=sresult.a[k]+cresult.a[j]*T(k-j)*a.a[k-j];
               
               sresult.a[k]=sresult.a[k]/T(k);
            }
            cresult.a[k]=sresult.a[0]*a.a[k]*T(k);
            for(j=1;j<k;j++)
               cresult.a[k]=cresult.a[k]+sresult.a[j]*T(k-j)*a.a[k-j];
               
            cresult.a[k]=cresult.a[k]/T(-k);
         }
         return cresult;
      }
      
      friend inline mtaylor<T> tan(const mtaylor<T> & a)
      {
         if(a.n<=0)
            return a; // nonsense
         mtaylor<T> sresult(sin(a.a[0]),a.n);
         mtaylor<T> cresult(cos(a.a[0]),a.n);
         int   k,j;
         for(k=1;k<a.n;k++)
         {
            sresult.a[k]=cresult.a[0]*a.a[k]*T(k);
            for(j=1;j<k;j++)
               sresult.a[k]=sresult.a[k]+cresult.a[j]*T(k-j)*a.a[k-j];
               
            sresult.a[k]=sresult.a[k]/T(k);
          
            cresult.a[k]=sresult.a[0]*a.a[k]*T(k);
            for(j=1;j<k;j++)
               cresult.a[k]=cresult.a[k]+sresult.a[j]*T(k-j)*a.a[k-j];
               
            cresult.a[k]=cresult.a[k]/T(-k);
         }

         return sresult/cresult;
      }

      friend inline mtaylor<T> cot(const mtaylor<T> & a)
      {
         if(a.n<=0)
            return a; // nonsense
         mtaylor<T> sresult(sin(a.a[0]),a.n);
         mtaylor<T> cresult(cos(a.a[0]),a.n);
         int   k,j;
         for(k=1;k<a.n;k++)
         {
            sresult.a[k]=cresult.a[0]*a.a[k]*T(k);
            for(j=1;j<k;j++)
               sresult.a[k]=sresult.a[k]+cresult.a[j]*T(k-j)*a.a[k-j];
               
            sresult.a[k]=sresult.a[k]/T(k);
          
            cresult.a[k]=sresult.a[0]*a.a[k]*T(k);
            for(j=1;j<k;j++)
               cresult.a[k]=cresult.a[k]+sresult.a[j]*T(k-j)*a.a[k-j];
               
            cresult.a[k]=cresult.a[k]/T(-k);
         }
         
         return cresult/sresult;
      }
      
      friend inline mtaylor<T> sinh(const mtaylor<T> & a)
      {
         if(a.n<=0)
            return a; // nonsense
         mtaylor<T> sresult(sinh(a.a[0]),a.n);
         mtaylor<T> cresult(cosh(a.a[0]),a.n);
         int   k,j;
         for(k=1;k<a.n;k++)
         {
            sresult.a[k]=cresult.a[0]*a.a[k]*T(k);
            for(j=1;j<k;j++)
               sresult.a[k]=sresult.a[k]+cresult.a[j]*T(k-j)*a.a[k-j];
               
            sresult.a[k]=sresult.a[k]/T(k);
            
            if(k<a.n-1)
            {
               cresult.a[k]=sresult.a[0]*a.a[k]*T(k);
               for(j=1;j<k;j++)
                  cresult.a[k]=cresult.a[k]+sresult.a[j]*T(k-j)*a.a[k-j];
               
               cresult.a[k]=cresult.a[k]/T(k);
            }
         }
         return sresult;
      }

      friend inline mtaylor<T> cosh(const mtaylor<T> & a)
      {
         if(a.n<=0)
            return a; // nonsense
         mtaylor<T> sresult(sinh(a.a[0]),a.n);
         mtaylor<T> cresult(cosh(a.a[0]),a.n);
         int   k,j;
         for(k=1;k<a.n;k++)
         {
            if(k<a.n-1)
            {
               sresult.a[k]=cresult.a[0]*a.a[k]*T(k);
               for(j=1;j<k;j++)
                  sresult.a[k]=sresult.a[k]+cresult.a[j]*T(k-j)*a.a[k-j];
               
               sresult.a[k]=sresult.a[k]/T(k);
            }
            cresult.a[k]=sresult.a[0]*a.a[k]*T(k);
            for(j=1;j<k;j++)
               cresult.a[k]=cresult.a[k]+sresult.a[j]*T(k-j)*a.a[k-j];
               
            cresult.a[k]=cresult.a[k]/T(k);
         }
         return cresult;
      }

      friend inline mtaylor<T> tanh(const mtaylor<T> & a)
      {
         if(a.n<=0)
            return a; // nonsense
         mtaylor<T> sresult(sinh(a.a[0]),a.n);
         mtaylor<T> cresult(cosh(a.a[0]),a.n);
         int   k,j;
         for(k=1;k<a.n;k++)
         {
            sresult.a[k]=cresult.a[0]*a.a[k]*T(k);
            for(j=1;j<k;j++)
               sresult.a[k]=sresult.a[k]+cresult.a[j]*T(k-j)*a.a[k-j];
               
            sresult.a[k]=sresult.a[k]/T(k);

            cresult.a[k]=sresult.a[0]*a.a[k]*T(k);
            for(j=1;j<k;j++)
               cresult.a[k]=cresult.a[k]+sresult.a[j]*T(k-j)*a.a[k-j];
               
            cresult.a[k]=cresult.a[k]/T(k);
         }

         return sresult/cresult;
      }
      
      friend inline mtaylor<T> coth(const mtaylor<T> & a)
      {
         if(a.n<=0)
            return a; // nonsense
         mtaylor<T> sresult(sinh(a.a[0]),a.n);
         mtaylor<T> cresult(cosh(a.a[0]),a.n);
         int   k,j;
         for(k=1;k<a.n;k++)
         {
            sresult.a[k]=cresult.a[0]*a.a[k]*T(k);
            for(j=1;j<k;j++)
               sresult.a[k]=sresult.a[k]+cresult.a[j]*T(k-j)*a.a[k-j];
               
            sresult.a[k]=sresult.a[k]/T(k);

            cresult.a[k]=sresult.a[0]*a.a[k]*T(k);
            for(j=1;j<k;j++)
               cresult.a[k]=cresult.a[k]+sresult.a[j]*T(k-j)*a.a[k-j];
               
            cresult.a[k]=cresult.a[k]/T(k);
         }

         return cresult/sresult;
      }
      
      friend inline mtaylor<T> log(const mtaylor<T> & a)
      {
         if(a.n<=0)
            return a; // nonsense
         mtaylor<T> result(log(a.a[0]),a.n);

         int k,j;
         for(k=1;k<a.n;k++)
         {
            result.a[k]=a.a[k];
            for(j=1;j<k;j++)
            {
               result.a[k]=result.a[k]-T(j)*result.a[j]*a.a[k-j]/T(k);
            }
            result.a[k]=result.a[k]/a.a[0];
         }
         return result;         
      }

      friend inline mtaylor<T> asin(const mtaylor<T> & a)
      {
         if(a.n<=0)
            return a; // nonsense
         mtaylor<T> result(asin(a.a[0]),a.n);
         mtaylor<T> invdiff(sqrt(T(1)-a*a));

         int k,j;
         for(k=1;k<a.n;k++)
         {
            result.a[k]=a.a[k];
            for(j=1;j<k;j++)
            {
               result.a[k]=result.a[k]-T(j)*result.a[j]*invdiff.a[k-j]/T(k);
            }
            result.a[k]=result.a[k]/invdiff.a[0];
         }
         return result;         
      }

      friend inline mtaylor<T> acos(const mtaylor<T> & a)
      {
         if(a.n<=0)
            return a; // nonsense
         mtaylor<T> result(acos(a.a[0]),a.n);
         mtaylor<T> invdiff(-sqrt(T(1)-a*a));

         int k,j;
         for(k=1;k<a.n;k++)
         {
            result.a[k]=a.a[k];
            for(j=1;j<k;j++)
            {
               result.a[k]=result.a[k]-T(j)*result.a[j]*invdiff.a[k-j]/T(k);
            }
            result.a[k]=result.a[k]/invdiff.a[0];
         }
         return result;         
      }

      friend inline mtaylor<T> atan(const mtaylor<T> & a)
      {
         if(a.n<=0)
            return a; // nonsense
         mtaylor<T> result(atan(a.a[0]),a.n);
         mtaylor<T> invdiff(T(1)+a*a);

         int k,j;
         for(k=1;k<a.n;k++)
         {
            result.a[k]=a.a[k];
            for(j=1;j<k;j++)
            {
               result.a[k]=result.a[k]-T(j)*result.a[j]*invdiff.a[k-j]/T(k);
            }
            result.a[k]=result.a[k]/invdiff.a[0];
         }
         return result;         
      }

      friend inline mtaylor<T> acot(const mtaylor<T> & a)
      {
         if(a.n<=0)
            return a; // nonsense
         mtaylor<T> result(acot(a.a[0]),a.n);
         mtaylor<T> invdiff(T(-1)-a*a);

         int k,j;
         for(k=1;k<a.n;k++)
         {
            result.a[k]=a.a[k];
            for(j=1;j<k;j++)
            {
               result.a[k]=result.a[k]-T(j)*result.a[j]*invdiff.a[k-j]/T(k);
            }
            result.a[k]=result.a[k]/invdiff.a[0];
         }
         return result;         
      }

      friend inline mtaylor<T> asinh(const mtaylor<T> & a)
      {
         return log(a+sqrt(a*a+T(1)));
      }
      friend inline mtaylor<T> acosh(const mtaylor<T> & a)
      {
         return log(a+sqrt(a*a-T(1)));
      }
      friend inline mtaylor<T> atanh(const mtaylor<T> & a)
      {
         return log((T(1)+a)/(T(1)-a))/T(2);
      }
      friend inline mtaylor<T> acoth(const mtaylor<T> & a)
      {
         return log((a+T(1))/(a-T(1)))/T(2);
      }
      
      friend inline mtaylor<T> pow(const mtaylor<T> & a,const T &b)
      {
         // if(a.n>1)
         //   a.a[1]=b*a.a[1]*pow(a.a[0],b-1);
         // a.a[0]=pow(a.a[0],b);
         // return a;
         
         return exp(b*log(a));
      } 

      friend inline mtaylor<T> power(const mtaylor<T> & a,int b)
      {
         // if(a.n>1)
         //   a.a[1]=b*a.a[1]*pow(a.a[0],b-1);
         // a.a[0]=pow(a.a[0],b);
         // return a;
         if(b<0)
            return power(1/a,-b);
         if(b==1)
            return a;
         if(b==0)
            return mtaylor<T>(1);
         if(b&1)
            return a*power(a,b-1);
         else
            return sqr(power(a,b/2));
      } 

};



#endif