#include <stdlib.h>
#include "toolsmith.h"
#include "Vector.h"
#include "Matrix.h"
#include "Matrixlt.h"


/************************************************************
 * Matrixlt class member functions                          *
 ************************************************************/

void Matrixlt::set(int n)
{
   int i;

   if(n < 0)
      error("Matrixlt:set(int)", "Matrix dimensions must be non-negative", ERREXIT);

   N = n;
   p = new Vector [n];
   for(i = 0; i < n; i++) {
      p[i].length(i + 1);
   }
}

void Matrixlt::clear()
{
   if(N > 0) delete [] p;
}

Matrixlt::Matrixlt()
{
   N = 0;
}

Matrixlt::Matrixlt(int n)
{
   set(n);
}

Matrixlt::Matrixlt(const Matrixlt &X)
{
   int i;

   N = X.N;
   p = new Vector [N];
   for(i = 0; i < N; i++) {
      p[i] = X.p[i];
   }
}

Matrixlt::~Matrixlt()
{
   clear();
}

Vector Matrixlt::dim() const
{
   Vector w(2);
   
   w[0] = w[1] = N;
   
   return w;
}

Matrixlt Matrixlt::dim(int n)
{
   clear();
   set(n);

   return *this;
}

Matrixlt::operator Matrix() const
{
   int i, j;
   Matrix Y(N, N);

   for(i = 0; i < N; i++) {
      for(j = 0; j <= i; j++) {
         Y[i][j] = p[i][j];
      }
   }

   return Y;
}

Matrix Matrixlt::operator !() const
{
   int i, j;
   Matrix Y(N, N);

   for(i = 0; i < N; i++) {
      for(j = 0; j <= i; j++) {
         Y[i][j] = !p[i][j];
      }
      for( ; j < N; j++) {
         Y[i][j] = true;
      }
   }

   return Y;
}

Matrixlt Matrixlt::operator -() const
{
   int i;
   Matrixlt Y(N);
   Vector *ap = Y.p, *bp = p;

   for(i = 0; i < N; i++) {
      *ap++ = -(*bp++);
   }

   return Y;
}

Matrix operator %(const Vector &v, const Matrixlt &X)
{
   int i, j, n = X.N;
   double sum;
   Matrix Y(1, X.N);

   if(n != v.length())
      error("Matrix operator %", "non-conformable matrices", ERRNCM);

   for(j = 0; j < n; j++) {
      sum = 0.0;
      for(i = j; i < n; i++) {
         sum += v[i] * X.p[i][j];
      }
      Y[0][j] = sum;
   }

   return Y;
}

Matrix operator %(const Matrix &X, const Matrixlt &Y)
{
   int i, j, k, m = X.nrow(), n = Y.N;
   double sum;
   Matrix Z(m, n);

   if(n != X.ncol())
      error("Matrix operator % ", "non-conformable matrices", ERRNCM);

   for(i = 0; i < m; i++) {
      for(j = 0; j < n; j++) {
         sum = 0.0;
         for(k = j; k < n; k++) {
            sum += X[i][k] * Y.p[k][j];
         }
         Z[i][j] = sum;
      }
   }

   return Z;
}

Matrix Matrixlt::operator %(const Vector &v) const
{
   int i, j;
   double sum;
   Matrix Y(N, 1);

   if(N != v.length())
      error("Matrixlt::operator %", "non-conformable matrices", ERRNCM);

   for(i = 0; i < N; i++) {
      sum = 0.0;
      for(j = 0; j <= i; j++) {
         sum += p[i][j] * v[j];
      }
      Y[i][0] = sum;
   }

   return Y;
}

Matrix Matrixlt::operator %(const Matrix &X) const
{
   int i, j, k, m = N, n = X.ncol();
   double sum;
   Matrix Y(m, n);

   if(m != X.nrow())
      error("Matrixlt::operator %", "non-conformable matrices", ERRNCM);

   for(i = 0; i < m; i++) {
      for(j = 0; j < n; j++) {
         sum = 0.0;
         for(k = 0; k <= i; k++) {
            sum += p[i][k] * X[k][j];
         }
         Y[i][j] = sum;
      }
   }

   return Y;
}

Matrix Matrixlt::operator %(const Matrixlt &X) const
{
   int i, j, k;
   double sum;
   Matrix Y(N, N);

   if(N != X.N)
      error("Matrixlt::operator %", "non-conformable matrices", ERRNCM);

   for(i = 0; i < N; i++) {
      for(j = 0; j <= i; j++) {
         sum = 0.0;
         for(k = j; k <= i; k++) {
            sum += p[i][k] * X.p[k][j];
         }
         Y[i][j] = sum;
      }
   }

   return Y;
}


#define arithmetic(op)                                                         \
Matrix Matrixlt::operator op (double x) const                                  \
{                                                                              \
   int i, j;                                                                   \
   Matrix Y(N, N);                                                             \
                                                                               \
   for(i = 0; i < N; i++) {                                                    \
      for(j = 0; j <= i; j++) {                                                \
         Y[i][j] = p[i][j] op x;                                               \
      }                                                                        \
      for( ; j < N; j++) {                                                     \
         Y[i][j] = 0.0 op x;                                                   \
      }                                                                        \
   }                                                                           \
                                                                               \
   return Y;                                                                   \
}                                                                              \
                                                                               \
Matrix Matrixlt::operator op (const Vector &v) const                           \
{                                                                              \
   int i, j, k = 0, l = v.length();                                            \
   Matrix Y(N, N);                                                             \
                                                                               \
   if(l == 0 || l > (N * N))                                                   \
      error("Matrixlt::operator " #op, "non-conformable matrices", ERRNCM);    \
                                                                               \
   for(j = 0; j < N; j++) {                                                    \
      for(i = 0; i < j; i++) {                                                 \
         Y[i][j] = 0.0 op v[k];                                                \
         k = (k + 1) % l;                                                      \
      }                                                                        \
      for( ; i < N; i++) {                                                     \
         Y[i][j] = p[i][j] op v[k];                                            \
         k = (k + 1) % l;                                                      \
      }                                                                        \
   }                                                                           \
                                                                               \
   return Y;                                                                   \
}                                                                              \
                                                                               \
Matrix Matrixlt::operator op (const Matrix &X) const                           \
{                                                                              \
   int i, j;                                                                   \
   Matrix Y(N, N);                                                             \
                                                                               \
   if(N != X.nrow() || N != X.ncol())                                          \
      error("Matrixlt::operator " #op, "non-conformable matrices", ERRNCM);    \
                                                                               \
   for(i = 0; i < N; i++) {                                                    \
      for(j = 0; j <= i; j++) {                                                \
         Y[i][j] = p[i][j] op X[i][j];                                         \
      }                                                                        \
      for( ; j < N; j++) {                                                     \
         Y[i][j] = 0.0 op X[i][j];                                             \
      }                                                                        \
   }                                                                           \
                                                                               \
   return Y;                                                                   \
}                                                                              \
                                                                               \
Matrix Matrixlt::operator op (const Matrixlt &X) const                         \
{                                                                              \
   int i, j;                                                                   \
   double x = 0.0;                                                             \
   Matrix Y(N, N);                                                             \
                                                                               \
   if(N != X.N)                                                                \
      error("Matrixlt::operator " #op, "non-conformable matrices", ERRNCM);    \
                                                                               \
   for(i = 0; i < N; i++) {                                                    \
      for(j = 0; j <= i; j++) {                                                \
         Y[i][j] = p[i][j] op X.p[i][j];                                       \
      }                                                                        \
      for( ; j < N; j++) {                                                     \
         Y[i][j] = 0.0 op x;                                                   \
      }                                                                        \
   }                                                                           \
                                                                               \
   return Y;                                                                   \
}                                                                              \

arithmetic(*);
arithmetic(/);

arithmetic(+);
arithmetic(-);

arithmetic(<);
arithmetic(<=);
arithmetic(>);
arithmetic(>=);

arithmetic(==);
arithmetic(!=);

arithmetic(&&);

arithmetic(||);

#undef arithmetic


Matrixlt Matrixlt::operator =(const Matrixlt &X)
{
   int i;
   Vector *ap, *bp = X.p;

   if(N != X.N) {
      clear();
      N = X.N;
      p = new Vector [N];
   }
   ap = p;
   for(i = 0; i < N; i++) {
      *ap++ = *bp++;
   }

   return *this;
}


#define assignment(op)                                                         \
Matrixlt Matrixlt::operator op (const Matrixlt &X)                             \
{                                                                              \
   int i;                                                                      \
   Vector *ap = p, *bp = X.p;                                                  \
                                                                               \
   if(N != X.N)                                                                \
      error("Matrixlt::operator " #op, "non-conformable matrices", ERRNCM);    \
                                                                               \
   for(i = 0; i < N; i++) {                                                    \
      *ap++ op *bp++;                                                          \
   }                                                                           \
                                                                               \
   return *this;                                                               \
}                                                                              \

assignment(+=);
assignment(-=);
assignment(*=);

#undef assignment


ostream &operator <<(ostream &os, const Matrixlt &X)
{
   int i, j;
   Vector *ap = X.p;

   if(X.N > 0) {
      os << *ap++;
      for(i = 1; i < X.N; i++) {
         os << endl << *ap++;
      }
   } else {
      os << "NULL";
   }

   return os;
}