/*******************************************************
* Unit:    matrix            release 0.10              *
* Purpose: general manipulation with matrices          *
* Licency: GPL or LGPL                                 *
* Copyright: (c) 1998-2021 Jaroslav Fojtik             *
********************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
//#include <alloc.h>

#ifndef __MATRIX__
 #include "matrix.h"
#endif


//------Basic procedures of the class matrix----------


//the constructor for the class Matrix without arguments
template <class FLOAT>
inline matrix<FLOAT>::matrix(void)
{
 size2D=size1D=0;
 data=NULL;
}


//the constructor for the class Matrix from one number
template <class FLOAT>
inline matrix<FLOAT>::matrix(FLOAT num)
{
 data=NULL;
 resize(1,1);
 if(data!=NULL) **data=num;
}


//the constructor for the class Matrix from another matrix
template <class FLOAT>
matrix<FLOAT>::matrix(const matrix<FLOAT> & m)
{
FLOAT **R1,*C1,**R2,*C2;
unsigned x,y;

 data=NULL;
 resize(m.size1D,m.size2D);
 if(data==NULL) {size1D=size2D=0;return;}

 R1=data;
 R2=m.data;
 for(y=0;y<size2D;y++)
   {
   C1=*R1++;
   C2=*R2++;
   for(x=0;x<size1D;x++)
   	{
	*C1++ = *C2++;
        }
   }
}


//the tunelling constructor for the class Matrix from temp_matrix
template <class FLOAT>
matrix<FLOAT>::matrix(temp_matrix<FLOAT> & tmp_m)
{
 data=tmp_m.data;       tmp_m.data=NULL;
 size1D=tmp_m.size1D;   tmp_m.size1D=0;
 size2D=tmp_m.size2D;   tmp_m.size2D=0;
}


//the constructor for the class Matrix with known size
template <class FLOAT>
inline matrix<FLOAT>::matrix(TSIZE x, TSIZE y)
{
 data=NULL;
 resize(x,y);
 if(data==NULL) size1D=size2D=0;
 return;
}


//the constructor of matrix from 2 dimensional array of data
template <class FLOAT>
matrix<FLOAT>::matrix(TSIZE x, TSIZE y, FLOAT **DATA)
{
FLOAT **R1,*C1,*C2;

 data=NULL;
 if(DATA==NULL) x=y=0;
 resize(x,y);
 if(data==NULL) {size1D=size2D=0;return;}

 R1=data;
 for(y=0;y<size2D;y++)
   {
   C1=*R1++;
   C2=*DATA++;
   for(x=0;x<size1D;x++)
	{
	*C1++ = *C2++;
	}
   }
}


//the constructor for the class Matrix with known size and initialization
template <class FLOAT>
inline matrix<FLOAT>::matrix(TSIZE x, TSIZE y, FLOAT contents)
{
 data=NULL;
 resize(x,y);
 fill(contents);
 return;
}


//this method allocates a new data structure with a given size
template <class FLOAT>
void matrix<FLOAT>::resize(TSIZE x, TSIZE y)
{
 if(data!=NULL)
	{
	if((x==size1D)&&(y==size2D)) return; //previous allocation is OK
	erase();
	}
 if(x==0 || y==0) return;
 if((data=(FLOAT **)calloc(y,sizeof(FLOAT *))) == NULL) //pointers to rows
	{
	size2D=size1D=0;		//allocation error
	RaiseError(MatrixId|No_Memory,this);
	return;
	}
 //memset(data,0,y*sizeof(FLOAT *)) is not needed, calloc clears memory;
 size2D=y;
 size1D=x;

 while(y-->0)
	{
	if( (data[y]=(FLOAT *)calloc(x,sizeof(FLOAT))) == NULL)
		{
		erase();		//allocation error
		RaiseError(MatrixId|No_Memory,this);
		return;
		}
	}
 return;
}

//this procedure will erase all contents of the matrix
template <class FLOAT>
void matrix<FLOAT>::erase(void)
{
TSIZE i;

 if(data!=NULL)
	{
	for(i=0;i<size2D;i++)
		{
		if(data[i]!=NULL) free(data[i]);
		}
	free(data);
	}
 size1D=size2D=0;
 data=NULL;
}


//the destructor of the class matrix
template <class FLOAT>
inline matrix<FLOAT>::~matrix(void)
{
 erase();
}


template <class FLOAT>
matrix<FLOAT> &matrix<FLOAT>::operator= (const FLOAT & F)
{
 resize(1,1);
 if(data!=NULL) **data=F;
return(*this);
}


//the move operator for Matrix object
template <class FLOAT>
matrix<FLOAT> &matrix<FLOAT>::operator= (const matrix<FLOAT> & m)
{
FLOAT **R1,*C1,**R2,*C2;
unsigned x,y;

 resize(m.size1D,m.size2D);
 if(data==NULL) return(*this);

 R1=data;
 R2=m.data;
 for(y=0;y<size2D;y++)
   {
   C1=*R1++;
   C2=*R2++;
   for(x=0;x<size1D;x++)
	{
	*C1++ = *C2++;
	}
   }
return(*this);
}


//the tunelling move operator for Matrix object
template <class FLOAT>
matrix<FLOAT> &matrix<FLOAT>::operator= (temp_matrix<FLOAT> & tmp_m)
{
 erase();
 data=tmp_m.data;       tmp_m.data=NULL;
 size1D=tmp_m.size1D;   tmp_m.size1D=0;
 size2D=tmp_m.size2D;   tmp_m.size2D=0;
return(*this);
}


//comparing two matrices
template <class FLOAT>
int matrix<FLOAT>::compare(const matrix<FLOAT> & m) const
{
FLOAT **R1,*C1,**R2,*C2;
TSIZE x,y;

 if((size1D!=m.size1D)||(size2D!=m.size2D)) return(0);	//not equal dimensions

 R1=data;
 R2=m.data;
 for(y=0;y<size2D;y++)
   {
   C1=*R1++;
   C2=*R2++;
   for(x=0;x<size1D;x++)
   	{
        if( *C1++ != *C2++) return(0); //element is not same
        }
   }

 return(1);     //yes, matrices are equal
}


//this method enables to manipulate with members of matrix
template <class FLOAT>
FLOAT & matrix<FLOAT>::member(TSIZE x, TSIZE y) const
{
static FLOAT PlaceHolder;

 if((x>=size1D)||(y>=size2D))
	{
	RaiseError(MatrixId|Bad_Index,(void *)this);	//Invalid index
	return(PlaceHolder);
	}
 return(data[y][x]);
}


template <class FLOAT>
void matrix<FLOAT>::row(TSIZE y, FLOAT *DataRow)
{
FLOAT *R;
TSIZE x;
 if(y>=size2D)
	{
	RaiseError(MatrixId|Bad_Index,(void *)this);	//Invalid index
	return;
	}
 R=data[y];
 x=size1D;
 while(x-->0)
   {
   *R++=*DataRow++;
   }
}


//Adding Operator
template <class FLOAT>
matrix<FLOAT> &matrix<FLOAT>::operator+=(const matrix<FLOAT> & m)
{
FLOAT **R1,*C1;
FLOAT **R2,*C2;
TSIZE x,y;

 if(m.size1D!=size1D || m.size2D!=size2D)
 	{
	erase();		//wrong dimensions
	RaiseError(MatrixId|Incompatible_Size,this);
        return(*this);
	}

 R1=data;
 R2=m.data;
 for(y=0;y<size2D;y++)
   {
   C1=*R1++;
   C2=*R2++;
   for(x=0;x<size1D;x++)
   	{
        *C1++ += *C2++;
        }
   }

return(*this);
}


template <class FLOAT>
matrix<FLOAT> &matrix<FLOAT>::operator+=(const FLOAT & num)
{
FLOAT **R1,*C1;
TSIZE x,y;

 R1=data;
 if(R1!=NULL)
   for(y=0;y<size2D;y++)
     {
     C1=*R1++;
     for(x=0;x<size1D;x++)
	{
	*C1++ += num;
	}
     }

return(*this);
}


template <class FLOAT>
temp_matrix<FLOAT> operator+ (const matrix<FLOAT> & m1, const matrix<FLOAT> & m2)
{
FLOAT **R,*C,**R1;
FLOAT *C1,**R2,*C2;
TSIZE x,y;

 if((m1.size1D!=m2.size1D)||(m1.size2D!=m2.size2D))
	{
	RaiseError(MatrixId|Incompatible_Size,(void *)&m1);
	return temp_matrix<FLOAT>(); 		//wrong dimensions
	}

 temp_matrix<FLOAT> tmp(m1.size1D,m1.size2D);
 if(tmp.data==NULL) return(tmp);

 R=tmp.data;
 R1=m1.data;
 R2=m2.data;
 for(y=0;y<m1.size2D;y++)
   {
   C=*R++;
   C1=*R1++;
   C2=*R2++;
   for(x=0;x<tmp.size1D;x++)
	{
	*C++ = *C1++ + *C2++;
	}
   }

 return tmp;
}


template <class FLOAT>
temp_matrix<FLOAT> operator+ (const matrix<FLOAT> & m1, const FLOAT & num)
{
FLOAT **R,*C,**R1,*C1;
TSIZE x,y;

 if(m1.data==NULL)
	{
	return temp_matrix<FLOAT>(); 	//wrong dimensions
	}

 temp_matrix<FLOAT> tmp(m1.size1D,m1.size2D);
 if(tmp.data==NULL) return(tmp);

 R=tmp.data;
 R1=m1.data;
 for(y=0;y<m1.size2D;y++)
   {
   C=*R++;
   C1=*R1++;
   for(x=0;x<tmp.size1D;x++)
   	{
        *C++ = *C1++ + num;
        }
   }

 return(tmp);
}


//Substracting Operators
template <class FLOAT>
matrix<FLOAT> &matrix<FLOAT>::operator-=(const matrix<FLOAT> & m)
{
FLOAT **R1,*C1,**R2,*C2;
TSIZE x,y;

 if(m.size1D!=size1D && m.size2D==size2D)
 	{
	erase();		//wrong dimensions
	RaiseError(MatrixId|Incompatible_Size,this);
        return(*this);
        }

 R1=data;
 R2=m.data;
 for(y=0;y<size2D;y++)
   {
   C1=*R1++;
   C2=*R2++;
   for(x=0;x<size1D;x++)
   	{
        *C1++ -= *C2++;
        }
   }

return(*this);
}


template <class FLOAT>
matrix<FLOAT> &matrix<FLOAT>::operator-=(const FLOAT & num)
{
FLOAT **R1,*C1;
TSIZE x,y;

 R1=data;
 for(y=0;y<size2D;y++)
   {
   C1=*R1++;
   for(x=0;x<size1D;x++)
   	{
        *C1++ -= num;
        }
   }

return(*this);
}


template <class FLOAT>
temp_matrix<FLOAT> operator- (const matrix<FLOAT> & m1, const matrix<FLOAT> & m2)
{
FLOAT **R,*C,**R1,*C1,**R2,*C2;
TSIZE x,y;

 if((m1.size1D!=m2.size1D)||(m1.size2D!=m2.size2D))
	{
	RaiseError(MatrixId|Incompatible_Size,(void *)&m1);
	return temp_matrix<FLOAT>();		//wrong dimensions
        }

 temp_matrix<FLOAT> tmp(m1.size1D,m1.size2D);

 if(tmp.data==NULL) return(tmp);

 R=tmp.data;
 R1=m1.data;
 R2=m2.data;
 for(y=0;y<m1.size2D;y++)
   {
   C=*R++;
   C1=*R1++;
   C2=*R2++;
   for(x=0;x<tmp.size1D;x++)
   	{
        *C++ = *C1++ - *C2++;
        }
   }

 return tmp;
}


template <class FLOAT>
temp_matrix<FLOAT> operator- (const matrix<FLOAT> & m1, const FLOAT & num)
{
FLOAT **R,*C,**R1,*C1;
TSIZE x,y;

 if(m1.data==NULL) return matrix<FLOAT>();	//empty matrix

 temp_matrix<FLOAT> tmp(m1.size1D,m1.size2D);
 if(tmp.data==NULL) return(tmp);

 R=tmp.data;
 R1=m1.data;
 for(y=0;y<m1.size2D;y++)
   {
   C=*R++;
   C1=*R1++;
   for(x=0;x<tmp.size1D;x++)
   	{
        *C++ = *C1++ - num;
        }
   }

 return tmp;
}


template <class FLOAT>
temp_matrix<FLOAT> operator- (const FLOAT & num, const matrix<FLOAT> & m1)
{
FLOAT **R,*C,**R1,*C1;
unsigned x,y;

 if(m1.data==NULL) return temp_matrix<FLOAT>();	//empty matrix

 temp_matrix<FLOAT> tmp(m1.size1D,m1.size2D);
 if(tmp.data==NULL) return(tmp);

 R=tmp.data;
 R1=m1.data;
 for(y=0;y<m1.size2D;y++)
   {
   C=*R++;
   C1=*R1++;
   for(x=0;x<tmp.size1D;x++)
	{
	*C++ = num - *C1++;
	}
   }

 return tmp;
}


template <class FLOAT>
temp_matrix<FLOAT> operator-(const matrix<FLOAT> & m)
{
FLOAT **R1,*C1,**R2,*C2;
TSIZE x,y;
temp_matrix<FLOAT> tmp(m.size1D,m.size2D);

 if(m.data==NULL || tmp.data==NULL) return(tmp);  //Matrix is empty

 R1=tmp.data;
 R2=m.data;
 for(y=0;y<m.size2D;y++)
   {
   C1=*R1++;
   C2=*R2++;
   for(x=0;x<m.size1D;x++)
   	{
        *C1++ = - *C2++;
        }
   }

return(tmp);
}


//Multiplying matrixes
template <class FLOAT>
temp_matrix<FLOAT> matrix<FLOAT>::multiply(const matrix<FLOAT> & m) const
{
if(size1D!=m.size2D)
	{
	RaiseError(MatrixId|Incompatible_Size,(void *)this);
	return temp_matrix<FLOAT>(); //wrong dimensions
	}

TSIZE i,j,k;
FLOAT **R,*C;
FLOAT **R2,**R1;
const FLOAT *C1;
temp_matrix<FLOAT> tmp(m.size1D,size2D);
 if(tmp.data==NULL) return(tmp);  //no memory - error is reported from constructor

 R=tmp.data;
 R1=data;
 R2=m.data;
 for(i=0;i<size2D;i++)
    {
    C=*R++;
    for(j=0;j<m.size1D;j++)
	 {
	 *C=0;
	 C1=*R1;

	 for(k=0;k<size1D;k++)
		    {
		    *C += *C1++ * R2[k][j];
		    }
	 C++;
	 }
    R1++;
    }

return(tmp);
}


template <class FLOAT>
temp_matrix<FLOAT> matrix<FLOAT>::multiply(const FLOAT & num) const
{
FLOAT **R,*C;
FLOAT **R1;
const FLOAT *C1;
TSIZE x,y;

 if(data==NULL)
	{
	return temp_matrix<FLOAT>(); 	//empty matrix
	}

temp_matrix<FLOAT> tmp(size1D,size2D);
 if(tmp.data==NULL) return(tmp);  //no memory - error is reported from constructor

 R=tmp.data;
 R1=data;
 for(y=0;y<size2D;y++)
   {
   C=*R++;
   C1=*R1++;
   for(x=0;x<tmp.size1D;x++)
   	{
        *C++ = *C1++ * num;
        }
   }

 return(tmp);
}


template <class FLOAT>
matrix<FLOAT> &matrix<FLOAT>::operator*=(const FLOAT & num)
{
FLOAT **R1,*C1;
TSIZE x,y;

 R1=data;
 for(y=0;y<size2D;y++)
   {
   C1=*R1++;
   for(x=0;x<size1D;x++)
	{
	*C1++ *= num;
	}
   }
return(*this);
}


template <class FLOAT>
matrix<FLOAT> &matrix<FLOAT>::operator*=(const matrix<FLOAT> & m)
{
  *this= *this * m;
return *this;
}


//dividing matrices
template <class FLOAT>
matrix<FLOAT> &matrix<FLOAT>::operator/=(const FLOAT & num)
{
FLOAT **R1,*C1;
TSIZE x,y;

 R1=data;
 for(y=0;y<size2D;y++)
   {
   C1=*R1++;
   for(x=0;x<size1D;x++)
	{
	*C1++ /= num;
	}
   }

return(*this);
}


template <class FLOAT>
matrix<FLOAT> &matrix<FLOAT>::operator/=(const matrix<FLOAT> & m)
{
matrix<FLOAT> Inv(m.inv());
  *this *= Inv;
return(*this);
}


template <class FLOAT>
temp_matrix<FLOAT> operator/ (const FLOAT &F, const matrix<FLOAT> & m)
{
temp_matrix<FLOAT> Inv(m.inv());
  Inv*=F;
return Inv;
}


template <class FLOAT>
temp_matrix<FLOAT> operator/ (const matrix<FLOAT> & m1, const matrix<FLOAT> & m2)
{
temp_matrix<FLOAT> Inv(m2.inv());
  Inv=m1 * Inv;
return Inv;
}


template <class FLOAT>
temp_matrix<FLOAT> matrix<FLOAT>::divide(const FLOAT &f) const
{
FLOAT **R,*C;
FLOAT **R1,*C1;
TSIZE x,y;

 if(data==NULL) return temp_matrix<FLOAT>(); 	//empty matrix

temp_matrix<FLOAT> tmp(size1D,size2D);
 if(tmp.data==NULL) return(tmp);

 R=tmp.data;
 R1=data;
 for(y=0;y<size2D;y++)
   {
   C=*R++;
   C1=*R1++;
   for(x=0;x<tmp.size1D;x++)
   	{
	*C++ = *C1++ / f;
	}
   }

 return(tmp);
}


//this method fill all matrix with a given value f
template <class FLOAT>
void matrix<FLOAT>::fill(FLOAT f)
{
FLOAT *C1;
TSIZE i,y;

 if(data==NULL) return;

 y=size2D;
 while(y-->0)
 	{
        C1=data[y];
        for(i=0;i<size1D;i++,C1++)
         	{
                *C1=f;
                }
        }
 return;
}


//this method produces an eigen matrix
template <class FLOAT>
void matrix<FLOAT>::ones(void)
{
FLOAT *C1;
TSIZE i,y;

 if(data==NULL) return;

 y=size2D;
 while(y-->0)
	{
	C1=data[y];
	for(i=0;i<size1D;i++,C1++)
		{
		*C1=(i==y)?1:0;
		}
	}
 return;
}


//this method returns a transposed matrix
template <class FLOAT>
temp_matrix<FLOAT> transpose( matrix<FLOAT>  & m )
{
FLOAT *C1;
TSIZE i,y;
temp_matrix<FLOAT> tmp(m.size2D,m.size1D);  //order is changed

 if(m.data==NULL || tmp.data==NULL) return tmp;

 y=m.size2D;
 while(y-->0)
	{
	for(C1=m.data[y],i=0;i<m.size1D;i++,C1++)
		{
		tmp.data[i][y]=*C1;
		}
        }
 return tmp;
}


//this procedure check consistency of the variable type matrix
template <class FLOAT>
int matrix<FLOAT>::check(void)
{
TSIZE y;

 if(data==NULL)
	{
        if(size1D!=0 || size2D!=0)
        	{
                size1D=size2D=0;
                return(MatrixId+Invalid_Size);
                }
        return(0);
        }

 if(size1D==0 || size2D==0)
 	{
        data=NULL;
        size1D=size2D=0;
	return(MatrixId+Wrong_Pointer);
        }

 for(y=0;y<size2D;y++)
 	{
        if(data[y]==NULL)
		{
                erase();
                return(MatrixId+Bad_Allocation);
		}
	}
return(0);
}


//This function calculates sum of squares all elements inside matrix
template <class FLOAT>
FLOAT matrix<FLOAT>::norm(void) const
{
FLOAT **R,*C;
FLOAT Norm=0;
TSIZE x,y;

 R=data;
 for(y=0;y<size2D;y++)
   {
   C=*R++;
   for(x=0;x<size1D;x++)
	{
	Norm+=(*C)*(*C);
	C++;
	}
   }
 return(Norm);
}


//------Advanced procedures of the class matrix----------
template<class FLOAT>
temp_matrix<FLOAT> matrix<FLOAT>::triangle(const matrix<FLOAT>& M, int *flips) const
{
FLOAT **o;
FLOAT **p, **q;
FLOAT *l, *r, *s;
FLOAT pivot_el, tmp;
TSIZE i,j, col,row;

TSIZE n = size2D;
TSIZE d = M.size1D;
TSIZE m = n+d;

temp_matrix<FLOAT> A(m,n);
 if(A.data==NULL) return A;

 o = data;
 p = A.data;
 q = M.data;
 for(i=0;i<n;i++)
   {
   l = *p++;
   r = *o++;
   for(j=0;j<n;j++) *l++ = *r++;
   if(q)
     {
     s = *q++;
     for(j=0;j<d;j++) *l++ = *s++;
     }
   }

 flips = 0;

 for(col=0, row=0; row<n; row++, col++)
   {	 // search for row j with maximal absolute entry in current col
   j = row;

   pivot_el = Abs(A.data[j][col]);
   for(i=row+1; i<n; i++)
     {
     tmp = Abs(A.data[i][col]);
     if(pivot_el < tmp) {j = i;pivot_el = tmp;}
     }

   if( n > j && j > row)
     {
     A.fliprows(j,row);
     if(flips) flips++;
     }

   q  = A.data + row;
   tmp = (*q)[col];     //A.data[row][col];

   if(Abs(tmp) < EPSILON) // matrix has not full rank
     {
     RaiseError(MatrixId|NoFullRank,(void *)this); //matrix::solve: matrix has not full rank.
     return temp_matrix<FLOAT>();
     }

   for (p = &A.data[n-1]; p != q; p--)
     {
     l = *p+col;
     s = *p+m;
     r = *q+col;

     if (*l != 0.0)
       {
       pivot_el = *l/tmp;
       while(l < s) *l++ -= *r++ * pivot_el;
       }
     }
  }
return A;
}


template<class FLOAT>
temp_matrix<FLOAT> matrix<FLOAT>::solve(const matrix<FLOAT>& M) const
{
FLOAT **p, **q;
FLOAT *l, *r, *s;

if (size1D!=size2D || size2D!=M.size2D)
	{
	RaiseError(MatrixId|Incompatible_Size,(void *)this); //Solve: wrong dimensions
	return temp_matrix<FLOAT>();
	}

 int      n = size2D;
 int      d = M.size1D;
 int      m = n+d;
 int      row, col;
 matrix<FLOAT> A(triangle(M));

 if(A.Data() == NULL)
    {
    RaiseError(MatrixId|NoFullRank,(void *)this); //matrix::solve: matrix has not full rank.
    return matrix<FLOAT>();
    }

 for(col = n-1, p = &A.data[n-1]; col>=0; p--, col--)
   {
   s = *p+m;
   FLOAT tmp = (*p)[col];

   for(l=*p+n; l < s; l++) *l /=tmp;
   for(q = A.data; q != p; q++ )
     {
     tmp = (*q)[col];
     l = *q+n;
     r = *p+n;
     while(r < s)  *l++ -= *r++ * tmp;
     }
   }

 temp_matrix<FLOAT> result(d,n);

 for(row=0; row<n; row++)
    {
    l = A[row]+n;
    r = result[row];
    for(col=0; col<d; col++)
	*r++ = *l++;
    }

return result;
}


/*This function calculates determinant of given matrix*/
template<class FLOAT>
FLOAT matrix<FLOAT>::det(void) const
{
TSIZE i;
FLOAT Det;

 if(size2D!=size1D)
     {
     RaiseError(MatrixId|NotSquare,(void *)this);
     return 0;
     }
 if(size2D==0 || size1D==0) return 0;

 matrix<FLOAT> A;
 A=triangle(matrix<FLOAT>());

 Det=A.member(0,0);
 for(i=1;i<A.size2D;i++)
	Det*=A.member(i,i);

return Det;
}


template<class FLOAT>
temp_matrix<FLOAT> matrix<FLOAT>::inv() const
{
 if(size2D!=size1D)
     {
     RaiseError(MatrixId|NotSquare,(void *)this);	//matrix::inv: matrix not quadratic.
     return temp_matrix<FLOAT>();
     }

 matrix<FLOAT> I(size1D,size2D);
 I.ones();
return solve(I);
}


template<class FLOAT>
matrix<FLOAT>  matrix<FLOAT>::exp(int MaxIterations) const
{
matrix<FLOAT> result(size1D,size2D),temp(*this);
int i;
//FLOAT prev,next;

 if(size2D!=size1D)
     {
     RaiseError(MatrixId|NotSquare,(void *)this);	//matrix::exp: matrix not quadratic.
     return matrix<FLOAT>();
     }
 if(size2D==0 || size1D==0 || data==NULL) return matrix<FLOAT>();

 result.ones();
 i=2;
 do {
    //prev=result(0,0);
    result+=temp;
    temp=temp*(*this)*((FLOAT)1/i);
    //if(prev==result(0,0)) break;
    } while(++i<=MaxIterations);
return(result);
}



#ifdef Streams
template<class FLOAT> ostream &operator<<(ostream &xout, const matrix<FLOAT> &M)
{
TSIZE x,y;
 xout<<"[";
 for(y=0;y<M.rows();y++)
   {
   if(y>0) xout<<"\n";
   for(x=0;x<M.cols();x++)
	{
	if(x>0) xout<<',';
	xout<<M.member(x,y);
	}
   }
xout<<"]";
return xout;
}

/*istream &operator>>(istream &, string &)
{
} */
#endif

