// 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 MATRIX_HPP_INCLUDED
#define MATRIX_HPP_INCLUDED

// #define INDCHECK

#include <iostream.h>

template <class T>
class mmatrix
{
   public:
      T  **a;
      int N,M;

   public:
      mmatrix(int aN,int aM) : N(aN),M(aM) 
      { 
         a=new T *[N];
         int i;
         for(i=0;i<N;i++)
            a[i]=new T[M];
      }
      mmatrix(int aN,int aM,const T & x,const T & y) : N(aN),M(aM) 
      { 
         a=new T *[N];
         int i;
         for(i=0;i<N;i++)
            a[i]=new T[M];
         int j;
         for(i=0;i<N;i++)
            for(j=0;j<M;j++)
               a[i][j]=(i==j)?x:y;
      }
      mmatrix(const T & b) : N(1),M(1)
      {
         a=new T *[1];
         a[0]=new T[1];
         a[0][0]=b;
      }
      mmatrix(const mmatrix<T> &b) : N(b.N),M(b.M)
      { 
         int i,j; 
         a=new T *[N];   
         for(i=0;i<N;i++) 
         {
            a[i]=new T[M];
            for(j=0;j<M;j++) 
               a[i][j]=b.a[i][j];
         } 
      }
      
     ~mmatrix(void) { int i; for(i=0;i<N;i++) delete [] a[i]; delete [] a; }
     
     
      mmatrix<T> & operator =(const mmatrix<T> & b)
      {
         int i,j;
         if(N!=b.N || M!=b.M)
         {
            for(i=0;i<N;i++) delete [] a[i]; delete [] a;  // what if self-assigment ?!??
            N=b.N;M=b.M;
            a=new T *[N];
            for(i=0;i<N;i++)
               a[i]=new T[M];
         }
         for(i=0;i<N;i++) 
            for(j=0;j<M;j++) 
               a[i][j]=b.a[i][j];
         return *this;
      } 
      
      inline mmatrix<T> & operator +=(const mmatrix<T> & b)
      {
         int i,j;
         int iN(b.N<N?b.N:N),iM(b.M<M?b.M:M);
         for(i=0;i<iN;i++) 
            for(j=0;j<iM;j++) 
               a[i][j]=a[i][j]+b.a[i][j];
         return *this;
      }

      inline mmatrix<T> & operator -=(const mmatrix<T> & b)
      {
         int i,j;
         int iN(b.N<N?b.N:N),iM(b.M<M?b.M:M);
         for(i=0;i<iN;i++) 
            for(j=0;j<iM;j++) 
               a[i][j]=a[i][j]-b.a[i][j];
         return *this;
      }
      
      inline mmatrix<T> & foreach(T (*f)(const T &))
      {
         int i,j;
         for(i=0;i<N;i++)
            for(j=0;j<M;j++)
               a[i][j]=f(a[i][j]);
         return *this;
      }
      inline mmatrix<T> foreach(T (*f)(const T &)) const
      {
         int i,j;
         mmatrix<T> temp(N,M);
         for(i=0;i<N;i++)
            for(j=0;j<M;j++)
               temp.a[i][j]=f(a[i][j]);
         return temp;
      }
      inline mmatrix<T> & foreach(T (*f)(T))
      {
         int i,j;
         for(i=0;i<N;i++)
            for(j=0;j<M;j++)
               a[i][j]=f(a[i][j]);
         return *this;
      }
      inline mmatrix<T> & foreach(T (*f)(T)) const
      {
         int i,j;
         mmatrix<T> temp(N,M);
         for(i=0;i<N;i++)
            for(j=0;j<M;j++)
               temp.a[i][j]=f(a[i][j]);
         return *this;
      }
      
      friend inline mmatrix<T> operator +(mmatrix<T> a,const mmatrix<T> & b) { return a+=b; }
      friend inline mmatrix<T> operator -(mmatrix<T> a,const mmatrix<T> & b) { return a-=b; }
      
      
      T & operator ()(int i,int j)       
      { 
#ifndef INDCHECK      
      return a[i][j]; 
#else
      if(i>=0 && i<N && j>=0 && j<M)
         return a[i][j];
      else
         throw(new int);
#endif

      }
      T &           Y(int i,int j)       { 
#ifndef INDCHECK      
      return a[i][j]; 
#else
      if(i>=0 && i<N && j>=0 && j<M)
         return a[i][j];
      else
         throw(new int);
#endif
}
      T   operator ()(int i,int j) const { 
#ifndef INDCHECK      
      return a[i][j]; 
#else
      if(i>=0 && i<N && j>=0 && j<M)
         return a[i][j];
      else
         throw(new int);
#endif
      }
      T             y(int i,int j) const { 
#ifndef INDCHECK      
      return a[i][j]; 
#else
      if(i>=0 && i<N && j>=0 && j<M)
         return a[i][j];
      else
         throw(new int);
#endif
      }
      
      mmatrix<T> row(int i) const
      {
         if(i>=0 && i<N)
         {
            mmatrix<T> r(1,M);
            for(int j=0;j<M;j++)
               r.a[0][j]=a[i][j];
            return r;
         } else
            return mmatrix<T>(0,M);
      }

      friend ostream & operator <<(ostream & o,const mmatrix<T> & x)
      {
         int i,j;
         for(i=0;i<x.N;i++)
         {
            for(j=0;j<x.M-1;j++)
               o << x.a[i][j] << '\t';
            o << x.a[i][j] << endl;
         }
         return o;
      }
};
      template <class S>
      inline mmatrix<S> operator *(const mmatrix<S> & x,const mmatrix<S> & y)
      {
         int i,j,k;
         int iB(x.M<y.N?x.M:y.N);
         mmatrix<S> z(x.N,y.M);
         
         for(i=0;i<x.N;i++)
            for(j=0;j<y.M;j++)
            {
               if(0<iB)
                  z(i,j)=x(i,0)*y(0,j);
               for(k=1;k<iB;k++)
                  z(i,j)=z(i,j)+x(i,k)*y(k,j);      
            }         
         return z;
      } 

      template <class S>
      inline int operator ==(const mmatrix<S> & x,const mmatrix<S> & y)
      {
         int i,j,k;
         if(x.M!=y.M || x.N!=y.N)
            return 0;
         
         for(i=0;i<x.N;i++)
            for(j=0;j<x.M;j++)
               if(x(i,j)!=y(i,j))
                  return 0;
         return 1;
      } 

      template <class S>
      inline int operator !=(const mmatrix<S> & x,const mmatrix<S> & y)
      {
         int i,j,k;
         if(x.M!=y.M || x.N!=y.N)
            return 1;
         
         for(i=0;i<x.N;i++)
            for(j=0;j<x.M;j++)
               if(x(i,j)!=y(i,j))
                  return 1;
         return 0;
      } 

      template <class S>
      inline mmatrix<S> operator *(const S & x,const mmatrix<S> & y)
      {
         int i,j,k;
         mmatrix<S> z(y.N,y.M);
         
         for(i=0;i<y.N;i++)
            for(j=0;j<y.M;j++)
               z(i,j)=x*y(i,j);
         return z;
      } 

      template <class S>
      inline mmatrix<S> operator -(const mmatrix<S> & y)
      {
         int i,j,k;
         mmatrix<S> z(y.N,y.M);
         
         for(i=0;i<y.N;i++)
            for(j=0;j<y.M;j++)
               z(i,j)=-y(i,j);
         return z;
      } 



template <class T>
mmatrix<T> trans(const mmatrix<T> & A)
{
   mmatrix<T> B(A.M,A.N);
   int i,j;
   for(i=0;i<B.N;i++)
      for(j=0;j<B.M;j++)
         B(i,j)=A(j,i);

   return B;
}

template <class T>
mmatrix<T> appendcol(const mmatrix<T> & A,const mmatrix<T> & c) // append col
{
   int maxN(A.N>c.N?A.N:c.N);
   mmatrix<T> B(maxN,A.M+c.M);
   int i,j;
   for(i=0;i<A.N;i++)
      for(j=0;j<A.M;j++)
         B(i,j)=A(i,j);
   for(i=0;i<c.N;i++)  
      for(j=0;j<c.M;j++)
         B(i,j+A.M)=c(i,j);

   return B;
}

template <class T>
mmatrix<T> appendline(const mmatrix<T> & A,const mmatrix<T> & l) // append line
{
   int maxM(A.M>l.M?A.M:l.M);
   mmatrix<T> B(A.N+l.N,maxM);
   int i,j;
   for(i=0;i<A.N;i++)
      for(j=0;j<A.M;j++)
         B(i,j)=A(i,j);

   for(i=0;i<l.N;i++)
      for(j=0;j<l.M;j++)
         B(i+A.N,j)=l(i,j);

   return B;
}

template <class T>
mmatrix<T> inv(mmatrix<T> A)
{
   int N(A.N<A.M?A.N:A.M);
   mmatrix<T> B(N,N);
   int i,j,k;
   for(i=0;i<N;i++)
      for(j=0;j<N;j++)
         B(i,j)=(i==j)?1:0;

   for(i=0;i<N;i++)
   {
      if(A(i,i)==T(0)) // trouble...
      {
         for(k=i+1;k<N && (k==i || A(k,i)==T(0));k++);
         if(k<N)       // okay, found a better row
         {
            T swap;
            for(j=0;j<N;j++)
            {
               swap=A(k,j);A(k,j)=A(i,j);A(i,j)=swap;
            }
         } else
         {
            return mmatrix<T>(0,0); // Sorry, inversion not possible.
         }  
      }
      for(j=0;j<N;j++)
         if(i!=j) 
         {
            T quot=A(j,i)/A(i,i);
            for(k=0;k<i+1;k++)
            {
               B(j,k)=B(j,k)-quot*B(i,k);
            }

            for(k=i+1;k<N;k++)
               A(j,k)=A(j,k)-quot*A(i,k);
         }
   }

   
   
   for(i=0;i<N;i++)
      for(j=0;j<N;j++)
         B(i,j)=B(i,j)/A(i,i);


   return B;
}

template <class T>
T det(mmatrix<T> A)
{
   if((A.N==1 && A.M>0) || (A.M==1 && A.N>0))
      return A(0,0);
      
   int N(A.N<A.M?A.N:A.M);
   mmatrix<T> B(N-1,N-1);
   T res;
   
   int i,j,k;
   for(i=0;i<N;i++)
   {
      for(j=0;j<i;j++)
         for(k=1;k<N;k++)  
            B(j,k-1)=A(j,k);
      j++;
      for(;j<N;j++)
         for(k=1;k<N;k++)  
            B(j-1,k-1)=A(j,k);
      
      if(i==0)      
         res=A(i,0)*det(B);
      else if(i&1)
         res=res-A(i,0)*det(B);
      else
         res=res+A(i,0)*det(B);
   }
   return res;
}

template<class T>
mmatrix<T> pow(mmatrix<T> A,int p)
{
   mmatrix<T> b(A.N,A.M,1,0);
   if(p<0)
   {
      A=inv(A);
      p=-p;
   }
   while(p!=0)
   {
      if(p&1)
         b=b*A;
      p=p/2;
      if(p)
         A=A*A;
   }
   return b;
}

#endif