// This is a part of the Object Oriented Mathematical Library
// Copyright (C) 1999-2008 Renishaw plc
//
// This source file is the confidential property of Renishaw plc.
// Reproduction or transmission in whole or in part, in any form or
// by any means, electronic, mechanical or otherwise, is prohibited
// without the prior written consent of the copyright owner.
//
// This is a part of the Object Oriented Mathematical Library
// Copyright (C) 1999-2008 Renishaw plc


// objalgeb.h -- algebraic objects
//=========================================================================================================================

#ifndef __OBJALGEB_H__
#define __OBJALGEB_H__

#include "OOMathsExport.h"
#include "ObjExcep.h"
//=========================================================================================================================
// CMaAlgebra - math library algebra (abstract base class)
//=========================================================================================================================

class OOMATHS_API CMaAlgebra : public CMaObject
{
public:
// destructor
   virtual ~CMaAlgebra();

protected:
// constructor
   CMaAlgebra();
};

//=========================================================================================================================
// CMaInterval - math library interval (concrete class)
//=========================================================================================================================

class OOMATHS_API CMaInterval : public CMaAlgebra
{
public:
// constructor
   CMaInterval(double d1 = 0.0, double d2 = Ma_2PI, int lopen = 0, int ropen = 0);

// destructor
   ~CMaInterval();

// attributes
   int IsLeftOpen() const;
   int IsRightOpen() const;
   void GetBounds(double* d1, double* d2) const;

// operations
   int IsIn(double value) const;
   int IsOverlapping(const CMaInterval& other) const;

private:
// implementation
   int    m_nLopen;
   int    m_nRopen;
   double m_Bound1;
   double m_Bound2;
};

//=========================================================================================================================
// CMaMatrix - math library matrix (concrete class)
//=========================================================================================================================

class OOMATHS_API CMaMatrix : public CMaAlgebra
{
public:
// constructors
   CMaMatrix();                         // default zero matrix[4][4]
   CMaMatrix(int nRow, int nColumn);    // variable-sized zero matrix[rows][columns]
   CMaMatrix(const CMaMatrix& other);   // copy

// destructor
   ~CMaMatrix();

// attributes
   inline int GetRowSize() const{return m_nRow;}					// get number of rows
   inline int GetColumnSize() const{return m_nColumn;}				// get number of columns
   const double** GetData() const{return (const double**)m_pM;}		// return two dimensional C array
   double** GetData(){ return m_pM; }

// operations
   // assignment and data element
   CMaMatrix& operator= (const CMaMatrix& other);
   inline double& Element(int i, int j){CheckRowColumnIndex(i,j);return m_pM[i][j];}
   inline double  Element(int i, int j) const{CheckRowColumnIndex(i,j);return m_pM[i][j];}
   double  GetElement(int i, int j) const{CheckRowColumnIndex(i,j);return m_pM[i][j];}
   void    SetElement(int i, int j, double value){CheckRowColumnIndex(i,j);m_pM[i][j] = value;}

   // vector element
   CMaVector GetRowVector   (int i) const;
   CMaVector GetColumnVector(int j) const;
   void SetRowVector   (int i, const CMaVector& vector);
   void SetColumnVector(int j, const CMaVector& vector);

   // access to subset
   CMaMatrix SubMatrix(int i1, int i2, int j1, int j2) const;
      // note:
      //    return a matrix with elements
      //       [i1][j1]..[i1][j2], j1 <= j2
      //       ...
      //       [i2][j1]..[i2][j2], i1 <= i2

   void SetSubMatrix(int i1, int i2, int j1, int j2, const CMaMatrix& matrix);

   // comparison
   EMaBoolean IsZero(double tolerance = MaGetZero()) const;
   EMaBoolean IsEqual(const CMaMatrix& other, double tolerance = MaGetTol()) const;
   EMaBoolean IsSymmetric(double tolerance = MaGetTol()) const;
   EMaBoolean IsOrthogonal(double tolerance = MaGetLengthTol()) const;

   // conversion
   operator CMaComplexMatrix() const;

   // arithmetic
   CMaMatrix operator+ (const CMaMatrix& other) const;
   CMaMatrix operator- (const CMaMatrix& other) const;
   CMaMatrix operator* (const CMaMatrix& other) const;
   CMaVector operator* (const CMaVector& other) const;
   CMaMatrix operator* (double number)          const;
   CMaMatrix operator/ (double number)          const;

   CMaMatrix operator- () const;

   CMaMatrix& operator+= (const CMaMatrix& other);
   CMaMatrix& operator-= (const CMaMatrix& other);
   CMaMatrix& operator*= (double number);
   CMaMatrix& operator/= (double number);

   CMaMatrix PreMultiplyTranspose() const;    // this.Transpose() * this
   CMaMatrix PostMultiplyTranspose() const;   // this * this.Transpose()

   OOMATHS_API friend CMaMatrix operator* (double number, const CMaMatrix& matrix);

   // decomposition
   void LUD(CMaIntArray* Index, int* Parity, int* Singular);
      // note: on lower and upper triangular decomposition
      //    method: Crout's with implicit pivoting
      //       The actual matrix decomposed is a row permutation of the original matrix.
      //    L + U -> this, where L's diagonal unity elements are not stored
      //    Index: row permutation, Parity (1/-1): number of row interchanges is even or odd
      //    Singular (0/1): singular matrix due to zero row vector or zero pivot element
   void LUD2(CMaMatrix* P, CMaMatrix* L, CMaMatrix* U) const;
      // note: explicit expression of LUD
      //    P*this = L*U
   void SVD(CMaVector* W, CMaMatrix* V);
      // note: on singular value decomposition
      //    this = U*W*V', U -> this (column-orthogonal), Wi >= 0, V orthogonal

   // other related calculation
   double Determinant()  const;
   double Trace() const;
   bool BinaryDump(char* filename) const;
   bool BinaryRead(char* filename);
   CMaMatrix Transpose() const;
   CMaVector EigenValues  (EMaStatus* status) const;
   CMaMatrix EigenVectors (EMaStatus* status) const;   // columns of the returned matrix are eigenvectors
   void EigenValuesVectors(EMaStatus* status, CMaVector* values, CMaMatrix* vectors) const;
   CMaMatrix Inverse(EMaStatus* status) const;
      // note: on square matrix
      //    This has to be square for: LUD, Determinant, Trace, EigenValues, EigenVectors and Inverse
      // note: on eigen values, eigen vectors, and inverse matrix
      //    It exists when status == Ma_SUCCESS

   // static
   static CMaMatrix Identity(int Size = 4);      // Size > 0
   static CMaMatrix DiagonalMatrix(const CMaVector& vector);

private:
// implementation
   void Init(int nRow, int nColumn);
   void Copy(const CMaMatrix& other);

   int m_nRow, m_nColumn;   // number of rows and columns (r, c)
   double** m_pM;           // m[r][c]

private:
   /// Inline private functions for speed gain
   inline void CheckRowColumnIndex(const int& iRow, const int& jCol)const
   {
	   if (CMaException::IsAssertEnabled())
		   assert(0 <= iRow && iRow < m_nRow && 0 <= jCol && jCol < m_nColumn);

	   if (iRow < 0 || m_nRow <= iRow)
		   throw CMaOperationException(CMaOperationException::badIndex,
		   "operation", "CMaMatrix::Element()", 0, m_nRow-1);
	   if (jCol < 0 || m_nColumn <= jCol)
		   throw CMaOperationException(CMaOperationException::badIndex,
		   "operation", "CMaMatrix::Element()", 0, m_nColumn-1);
   }
};

//=========================================================================================================================
// CMaComplex - math library complex number in Cartesian coordinates (concrete class)

class OOMATHS_API CMaComplex : public CMaAlgebra
{
public:
// constructor
   CMaComplex(double a = 0.0, double b = 0.0);   // construct a + ib

// destructor
   ~CMaComplex();

// attributes
   double r() const;                             // get real part
   double i() const;                             // get imaginary part
   void GetData(double* a, double* b) const;

// operations
   // conversion
   CMaComplexPolar ToPolarForm() const;

   // comparison
   EMaBoolean IsZero(double tolerance = MaGetZero()) const;
   EMaBoolean IsEqual(const CMaComplex& other, double tolerance = MaGetTol()) const;
   EMaBoolean IsReal(double tolerance = MaGetTol()) const;        // is a real number
   EMaBoolean IsImaginary(double tolerance = MaGetTol()) const;   // is an imaginary number

   // arithmetic
   CMaComplex operator+ (const CMaComplex& other) const;
   CMaComplex operator- (const CMaComplex& other) const;
   CMaComplex operator* (const CMaComplex& other) const;
   CMaComplex operator/ (const CMaComplex& other) const;
   CMaComplex operator* (double number)           const;
   CMaComplex operator/ (double number)           const;

   CMaComplex operator- () const;

   CMaComplex& operator+= (const CMaComplex& other);
   CMaComplex& operator-= (const CMaComplex& other);
   CMaComplex& operator*= (double number);
   CMaComplex& operator/= (double number);

   OOMATHS_API friend CMaComplex operator* (double number, const CMaComplex& complex);

   double Modulus() const;                       // absolute value
   CMaComplex Conjugate() const;                 // construct a - bi
   CMaComplex Sqrt() const;                      // square root
   CMaComplex Power(int n) const;                // raise to the nth power (n>0)
   CMaComplexArray Roots(int n) const;           // nth roots[i], i = 0..n-1 (n>1)

private:
// implementation
   double m_Real;   // real part
   double m_Imag;   // imaginary part
};

//=========================================================================================================================
// CMaComplexPolar - math library complex number in polar coordinates (concrete class)
//=========================================================================================================================

class OOMATHS_API CMaComplexPolar : public CMaAlgebra
{
public:
// constructor
   CMaComplexPolar(double modulus = 1.0, double argument = 0.0);
      // note: on constructed complex number
      //    modulus (cos (argument) + i sin(argument)), where argument in radians

// destructor
   ~CMaComplexPolar();

// attributes
   double Modulus() const;
   double Argument() const;
   void GetData(double* modulus, double* argument) const;

// operations
   // conversion
   CMaComplex ToCartesianForm() const;

   // arithmetic
   CMaComplexPolar operator* (const CMaComplexPolar& other) const;
   CMaComplexPolar operator/ (const CMaComplexPolar& other) const;

   CMaComplexPolar Power(int n) const;           // raise to the nth power (n>0)
   CMaComplexPolarArray Roots(int n) const;      // nth roots[i], i = 0..n-1 (n>1)

private:
// implementation
   double m_M;   // modulus of complex number
   double m_A;   // argument of complex number
};

//=========================================================================================================================
// CMaComplexVector - math library complex vector (concrete class)
//=========================================================================================================================

class OOMATHS_API CMaComplexVector : public CMaAlgebra
{
public:
// constructors
   CMaComplexVector(int Size = 3);                    // size > 0
   CMaComplexVector(const CMaComplexVector& other);   // copy

// destructor
   ~CMaComplexVector();

// attributes
   int GetSize() const;                 // get dimension
   const CMaComplex* GetData() const;   // return C array
   CMaComplex* GetData();

   CMaVector Real() const;              // get real part
   CMaVector Imaginary() const;         // get imaginary part

// operations
   // assignment and data element
   CMaComplexVector& operator= (const CMaComplexVector& other);
   CMaComplex& operator[] (int Index);
   CMaComplex  operator[] (int Index) const;

   // access to subset
   CMaComplexVector SubVector(int i, int j) const;
      // note:
      //    return a vector with elements [i]..[j], i <= j

   // comparison
   EMaBoolean IsZero(double tolerance = MaGetZero()) const;
   EMaBoolean IsEqual(const CMaComplexVector& other, double tolerance = MaGetTol()) const;
   EMaBoolean IsReal(double tolerance = MaGetTol()) const;        // is a real vector
   EMaBoolean IsImaginary(double tolerance = MaGetTol()) const;   // is an imaginary vector

   // arithmetic
   CMaComplexVector operator+ (const CMaComplexVector& other) const;   // addition
   CMaComplexVector operator- (const CMaComplexVector& other) const;   // subtraction
   CMaComplex       operator% (const CMaComplexVector& other) const;   // conjugate dot product
   CMaComplexMatrix operator* (const CMaComplexVector& other) const;   // column and row vector product
   CMaComplexVector operator* (const CMaComplexMatrix& other) const;   // vector matrix product

   CMaComplexVector operator* (const CMaComplex& number) const;   // scalar product
   CMaComplexVector operator/ (const CMaComplex& number) const;   // scalar division
   CMaComplexVector operator* (double number) const;              // scalar product
   CMaComplexVector operator/ (double number) const;              // scalar division

   CMaComplexVector operator- () const;                           // unary negation

   OOMATHS_API friend CMaComplexVector operator* (double number, const CMaComplexVector& vector);
   OOMATHS_API friend CMaComplexVector operator* (const CMaComplex& number, const CMaComplexVector& vector);

   // other related calculation
   void Normalise();
   double Magnitude() const;
   CMaComplexVector Conjugate() const;
   CMaComplexVector Householder() const;

private:
// implementation
   void Init(int Size);
   void Copy(const CMaComplexVector& other);

   CMaComplex* m_pV;   // array of data
   int m_nSize;        // number of elements
};

//=========================================================================================================================
// CMaComplexMatrix - math library complex matrix (concrete class)
//=========================================================================================================================

class OOMATHS_API CMaComplexMatrix : public CMaAlgebra
{
public:
// constructors
   CMaComplexMatrix();                                // default zero matrix[4][4]
   CMaComplexMatrix(int nRow, int nColumn);           // zero matrix[rows][columns]
   CMaComplexMatrix(const CMaComplexMatrix& other);   // copy

// destructor
   ~CMaComplexMatrix();

// attributes
   int GetRowSize() const;               // get number of rows
   int GetColumnSize() const;            // get number of columns
   const CMaComplex** GetData() const;   // return two dimensional C array
   CMaComplex** GetData();

   CMaMatrix Real() const;               // get real part
   CMaMatrix Imaginary() const;          // get imaginary part

// operations
   // assignment and data element
   CMaComplexMatrix& operator= (const CMaComplexMatrix& other);
   CMaComplex& Element(int i, int j);         // element [i][j]
   CMaComplex  Element(int i, int j) const;
   CMaComplex  GetElement(int i, int j) const;
   void        SetElement(int i, int j, const CMaComplex& value);

   // vector element
   CMaComplexVector GetRowVector   (int i) const;
   CMaComplexVector GetColumnVector(int j) const;
   void SetRowVector   (int i, const CMaComplexVector& vector);
   void SetColumnVector(int j, const CMaComplexVector& vector);

   // access to subset
   CMaComplexMatrix SubMatrix(int i1, int i2, int j1, int j2) const;
      // note:
      //    return a matrix with elements
      //       [i1][j1]..[i1][j2], j1 <= j2
      //       ...
      //       [i2][j1]..[i2][j2], i1 <= i2

   void SetSubMatrix(int i1, int i2, int j1, int j2, const CMaComplexMatrix& matrix);

   // comparison
   EMaBoolean IsZero (double tolerance = MaGetZero()) const;
   EMaBoolean IsEqual(const CMaComplexMatrix& other, double tolerance = MaGetTol()) const;
   EMaBoolean IsReal (double tolerance = MaGetTol()) const;       // is a real matrix
   EMaBoolean IsImaginary(double tolerance = MaGetTol()) const;   // is an imaginary matrix

   EMaBoolean IsDiagonal (double tolerance = MaGetTol()) const;
   EMaBoolean IsHermitian(double tolerance = MaGetTol()) const;
   EMaBoolean IsSymmetric(double tolerance = MaGetTol()) const;
   EMaBoolean IsUnitary  (double tolerance = MaGetTol()) const;
   EMaBoolean IsUpperHessenberg(double tolerance = MaGetTol()) const;
   EMaBoolean IsUpperTriangular(double tolerance = MaGetTol()) const;

   // arithmetic
   CMaComplexMatrix operator+ (const CMaComplexMatrix& other) const;
   CMaComplexMatrix operator- (const CMaComplexMatrix& other) const;
   CMaComplexMatrix operator* (const CMaComplexMatrix& other) const;
   CMaComplexVector operator* (const CMaComplexVector& other) const;
   CMaComplexMatrix operator* (const CMaComplex& number)      const;
   CMaComplexMatrix operator/ (const CMaComplex& number)      const;
   CMaComplexMatrix operator* (double number)                 const;
   CMaComplexMatrix operator/ (double number)                 const;

   CMaComplexMatrix operator- () const;

   OOMATHS_API friend CMaComplexMatrix operator* (double number, const CMaComplexMatrix& matrix);
   OOMATHS_API friend CMaComplexMatrix operator* (const CMaComplex& number, const CMaComplexMatrix& matrix);

   CMaComplexMatrix& operator+= (const CMaComplexMatrix& other);
   CMaComplexMatrix& operator-= (const CMaComplexMatrix& other);
   CMaComplexMatrix& operator*= (const CMaComplexMatrix& other);

   // decomposition
   void QRD(CMaComplexMatrix& Q, CMaComplexMatrix& R) const;
      // note: on unitary triangular decomposition
      //    this = Q*R
      //    Q: a unitary matrix, R: an upper triangular matrix

   // other related calculation
   CMaComplex Determinant() const;
   CMaComplexVector EigenValues () const;
   CMaComplexMatrix EigenVectors() const;   // columns of the returned matrix are eigenvectors
   CMaComplexMatrix ConjTranspose() const;

   // static
   static CMaComplexMatrix Identity(int Size = 4);   // Size > 0
   static CMaComplexMatrix DiagonalMatrix(const CMaComplexVector& vector);

private:
// implementation
   void Init(int nRow, int nColumn);
   void Copy(const CMaComplexMatrix& other);

   CMaComplexMatrix& Row_operation(int i, int j);                  // interchange row(i) and row(j)
   CMaComplexMatrix& Row_operation(int i, int j, CMaComplex a);    // add a*row(i) to row(j)
   CMaComplexMatrix& Row_operation(int i, CMaComplex a);           // multiply row(i) by a
   CMaComplexMatrix& Column_operation(int i, int j);               // interchange col(i) and col(j)
   CMaComplexMatrix& Column_operation(int i, int j, CMaComplex a); // add a*col(i) to col(j)
   CMaComplexMatrix& Column_operation(int i, CMaComplex a);        // multiply col(i) by a

   void static Calculate_Givens_parameters(CMaComplex *c,CMaComplex *s,CMaComplex a,CMaComplex b);
   void Givens_pre_multiplication(int i,int k,CMaComplex c,CMaComplex s);
   void Givens_post_multiplication(int i,int k,CMaComplex c,CMaComplex s);

   CMaComplexVector Calculate_2_by_2_matrix_eigenvector();
   void Hessenberg_reduction(CMaComplexMatrix& Q);
   void Schur_decomposition(CMaComplexMatrix & Q, int *conv);

   int m_nRow, m_nColumn;   // number of rows and columns (r, c)
   CMaComplex** m_pM;       // m[r][c]
};

//=========================================================================================================================
// CMaPolynomial - math library real coefficient polynomial (concrete class)
//=========================================================================================================================

// formula: a0 + a1 * x + a2 * x**2 + ... + an * x**n

class OOMATHS_API CMaPolynomial : public CMaAlgebra
{
public:
// constructors
   CMaPolynomial(int degree = 0);               // degree >=0
   CMaPolynomial(const CMaPolynomial& other);   // copy

// destructor
   ~CMaPolynomial();

// attributes
   int GetDegree() const;
   const double* GetData() const;               // return C array (coefficients)
   double* GetData();

// operations
   // assignment and data element
   CMaPolynomial& operator= (const CMaPolynomial& other);
   double& operator[] (int Index);              // access coefficients[i]
   double  operator[] (int Index) const;

   // comparison
   EMaBoolean IsZero(double tolerance = MaGetZero()) const;
   EMaBoolean IsEqual(const CMaPolynomial& other, double tolerance = MaGetTol()) const;

   // arithmetic
   CMaPolynomial operator+ (const CMaPolynomial& other) const;
   CMaPolynomial operator- (const CMaPolynomial& other) const;
   CMaPolynomial operator* (const CMaPolynomial& other) const;
   CMaPolynomial operator+ (double number) const;
   CMaPolynomial operator- (double number) const;
   CMaPolynomial operator* (double number) const;
   CMaPolynomial operator/ (double number) const;

   CMaPolynomial operator- () const;

   CMaPolynomial& operator+= (double number);
   CMaPolynomial& operator-= (double number);
   CMaPolynomial& operator*= (double number);
   CMaPolynomial& operator/= (double number);

   OOMATHS_API friend CMaPolynomial operator* (double number, const CMaPolynomial& poly);

   // least-squares fitting: closest fit
   void Fit(const CMaDoubleArray& X, const CMaDoubleArray& Y);
   void Fit(const CMaDoubleArray& X, const CMaDoubleArray& Y, double* sqrErrors);

   // Chebyshev fitting: minimized maximum deviation fit
   void FitMinMax(double a, double b, double (*f)(double x));   // y = f(x), a <= x <= b

   // other related calculation
   double ValueAt(double x) const;
   double DerivativeAt(double x) const;
   void Simplify();                             // make sure coefficient a[n] != 0

private:
// implementation
   void Copy(const CMaPolynomial& other);

   int m_nN;       // degree of polynomial n
   double* m_pA;   // coefficients a[0..n]
   CMaChebyshevPolynomialExpression* m_pExp;
};

//=========================================================================================================================
// CMaChebyshevPolynomial - math library Chebyshev polynomial (concrete class)
//=========================================================================================================================

// formula: Tn(x) = cos (n arccos x), -1 <= x <= 1
// recurrence:
//    T0(x) = 1
//    T1(x) = x
//    Ti(x) = 2 * x * Ti-1(x) - Ti-2(x), i >= 2

class OOMATHS_API CMaChebyshevPolynomial : public CMaAlgebra
{
public:
// constructors
   CMaChebyshevPolynomial(int degree = 0);      // degree >= 0

// operations
   // conversions
   CMaPolynomial ToPolynomial() const;

private:
// implementation
   CMaPolynomial m_pPoly;
};

//=========================================================================================================================
// CMaChebyshevPolynomialExpression - math library Chebyshev polynomial expression (concrete class)
//=========================================================================================================================

// formula: c0 * T0(x) + c1 * T1(x) +.... + cn * Tn(x) - c0/2
//    linear combination of Chebyshev polynomials

class OOMATHS_API CMaChebyshevPolynomialExpression : public CMaAlgebra
{
public:
// constructors
   CMaChebyshevPolynomialExpression(int degree = 0);   // degree >= 0

// operations
   // Chebyshev fitting: minimized maximum deviation fit
   void FitMinMax(double a, double b, double (*f)(double x));   // y = f(x), a <= x <= b

   // calculation
   double ValueAt(double x) const;
   double DerivativeAt(double x) const;

   // conversion
   CMaPolynomial ToPolynomial() const;

private:
// implementation
   int m_nN;                // truncated degree
   double m_A, m_B;         // x range for function approximation
   CMaPolynomial m_pPoly;   // coefficients c[0..n], maximum degree n
};

//=========================================================================================================================
// CMaAlgEqu - math library algebraic equation (concrete class)
//=========================================================================================================================

class OOMATHS_API CMaAlgEqu : public CMaAlgebra
{
public:
// constructor
   CMaAlgEqu();

// destructor
   ~CMaAlgEqu();

// attributes
   int GetNumberRealRoots() const;
   int GetNumberComplexRoots() const;
   CMaDoubleArray GetRealRoots() const;
   CMaDoubleArray GetDistinctRealRoots() const;
   CMaComplexArray GetComplexRoots() const;

// operation
   void Solve(const CMaPolynomial& p);     // real coefficient equation
   void Solve(const CMaComplexArray& p);   // complex coefficient equation
      // note: on Solve
      //    equation form: p[0] + p[1] * x + p[2] * x**2 + ... + p[n] * x**n
      //    degree n > 0

private:
// implementation
   // default copy constructor and assignment disabled
   CMaAlgEqu(const CMaAlgEqu&);
   void operator=(const CMaAlgEqu&);

   int m_nReal;             // number of real roots
   int m_nComp;             // number of complex roots
   double*     m_pRoots;    // real roots[nReal]
   CMaComplex* m_pCroots;   // complex roots[nComp]
};

//=========================================================================================================================
// CMaLinearEquSys - math library linear equation system (concrete class)
//=========================================================================================================================

// formula: A[m][n] * x[n] = b[m], matrix A and vector b given, vector x unknown

class OOMATHS_API CMaLinearEquSys : public CMaAlgebra
{
public:
// constructor
   CMaLinearEquSys(int fit = 1);   // seek least-squares or exact solution (fit = 0)

// destructor
   ~CMaLinearEquSys();

// attributes
   // a unique solution
   EMaBoolean SolutionExists() const;
   CMaVector  GetSolution() const;

   // multiple solutions: the first vector + any linear combination of the rest
   //CMaVectorArray GetSolutionSet() const;   // for the case: fit = 0

// operations
   // set linear equality and inequality constraints for least-squares problems
   void SetEquality  (const CMaMatrix& C, const CMaVector& d);   // C * x = d
   void SetInequality(const CMaMatrix& E, const CMaVector& f);   // E * x >= f
   void SetNonnegative();          // x >= 0
   void SetLeastDistance();        // minimize |x| subject to E * x >= f
   void SetSingularSolvable();     // replace zero pivot element in LUD

   // solve
   void Solve();                                         // least distance programming
   void Solve(const CMaMatrix& A, const CMaVector& b);   // A * x = b
   void SetRightHandSide(const CMaVector& b);            // solve for another b

private:
// implementation
   // default copy constructor and assignment disabled
   CMaLinearEquSys(const CMaLinearEquSys&);
   void operator=(const CMaLinearEquSys&);

   void LUDbksb();
   void SVDbksb(const CMaVector& b);

   int m_nFit;
   int m_nType;   // constraints type 0..6:
       // null, equality, inequality, equality & inequality, nonnegative,
       // least distance, solution for singular case

   CMaVector*     m_pX;    // unique solution
   //CMaVectorArray m_pXi;   // solution set

   CMaMatrix*   m_cA;      // LUD of A[n][n] = L + U
   CMaIntArray* m_nIndx;   // LUD row permutation
   CMaVector*   m_cW;      // SVD of A = U.W.V', m_cU = m_cA
   CMaMatrix*   m_cV;
};


#endif // __OBJALGEB_H__

//=========================================================================================================================
