#include <iostream.h>
#include <stdlib.h>
#include <math.h>
//#define NDEBUG // uncomment to remove checking of assert()
#include <assert.h>
#define DEFAULT_ALLOC 2



template <class ElType> class matrix;

template <class ElType>
class vector
{
	friend class matrix<ElType>;
	 ElType * data;
	 int len;
	 public:
	 int length()const;
	 vector();
	 vector(int n);
	 ~vector(){ delete [] data;}
	 //Copy operator
	 vector(const vector<ElType>& v) ;
	 //assignment operator
	vector<ElType>& operator =(const vector<ElType> &original);
	ElType& operator[](int i)const  ;
	void zero();
	vector<ElType> operator+(const vector<ElType>& v);
	vector<ElType> operator-(const vector<ElType>&v);
	void  rprint()const;  //print entries on a single line
	void resize(int n);
	int operator==(const vector<ElType>& v)const;
	friend   vector<ElType> operator*(ElType c,vector<ElType>& v );
	friend   vector<ElType> operator*(vector<ElType>& v,ElType c );
	friend ostream& operator<<(ostream& s,vector<ElType>& v);
};
template <class ElType>
void vector<ElType>::zero()
{
	for(int i=0;i<len;i++) data[i]=(ElType)0;
}
template <class ElType>
int vector<ElType>::length()const
{
	return len;
}
template <class ElType>
ElType& vector<ElType>::operator[](int i)const
{
	assert(i>=0 && i < len);
	return data[i];
}

template <class ElType>
vector<ElType>::vector()
{
	data=new ElType[ DEFAULT_ALLOC];
	assert(data!=0);
	len=  DEFAULT_ALLOC;
}
template <class ElType>
vector<ElType>::vector(int n)
{
	data = new ElType[len=n];
	assert(data!=0);
}
template <class ElType>
vector<ElType>::vector(const vector<ElType>& v)
{
	data=new ElType[len=v.len];
	assert(data!=0);
	for(int i=0;i<len;i++) data[i]=v.data[i];
}
template <class ElType>
vector<ElType>& vector<ElType>::operator =(const vector<ElType> &original)
{
		if(this != &original)
		{
			delete [] data;
			data= new ElType[len=original.len];
			assert(data!=0);
			for(int i=0;i<len;i++) data[i]=original.data[i];
		}
		return *this;
  }
template <class ElType>
vector<ElType> vector<ElType>::operator+(const vector<ElType>& v)
{
	vector<ElType> sum(len);
	for(int i=0;i<len;i++) sum[i] = data[i]+v.data[i];
	return sum;

}
template <class ElType>
vector<ElType> vector<ElType>::operator-(const vector<ElType>& v)
{
		vector<ElType> sum(len);
		for(int i=0;i<len;i++) sum[i] = data[i]-v.data[i];
		return sum;
}
template <class ElType>
void  vector<ElType>::rprint()const  //print entries on a single line
{
		int i;
		cout << "VECTOR: ";
		cout << "(";
		for(i=0;i<len-1;i++) cout << data[i] << ",";
		cout << data[len-1] << ")" << endl;
		return;
}
template <class ElType>
void vector<ElType>::resize(int n)
{
		delete[]data;
		data = new ElType[len=n];
		assert(data !=0);
}
template <class ElType>
int vector<ElType>::operator==(const vector<ElType>& v)const
{
		if(len != v.len) return 0;
		for(int i=0;i<len;i++) if(data[i]!=v.data[i]) return 0;
		return 1;
}
template <class ElType>
vector<ElType> operator*(ElType c,vector<ElType>& v )
{
		vector<ElType> ans(v.len);
		for(int i=0;i<v.len;i++) ans[i]=c*v[i];
		return ans;
}
template <class ElType>
vector<ElType> operator*(vector<ElType>& v,ElType c )
	{
		vector<ElType> ans(v.len);
		for(int i=0;i<v.len;i++) ans[i]=c*v[i];
		return ans;
	}
template <class ElType>
ostream& operator<<(ostream& s,vector<ElType>& v)
{
	  s << "(";
	for(int i=0;i<v.len-1;i++) s << v.data[i] << ", ";
	s << v.data[v.len-1]<<")"<<endl;
	return s;
}

template <class ElType>
class matrix
{
	vector<ElType> *m;
	int rows,cols;
	public:
	matrix();
	matrix( int r, int c);
	matrix(const matrix<ElType> &s);
	~matrix();
	matrix& operator =(const matrix<ElType>& s);
	vector<ElType>& operator[](const int i);
	vector<ElType> operator*(const vector<ElType>&);
	friend matrix<ElType>  operator*(const ElType&, const matrix<ElType>&);
	friend matrix<ElType> operator*(const matrix<ElType>&, const ElType&);
	matrix<ElType> operator*(const matrix<ElType>& a);
	matrix<ElType> operator+(const matrix<ElType>& a);
	matrix<ElType> operator-(const matrix<ElType>& a);
	matrix<ElType> transpose();
	//matrix<ElType> inverse();
	friend ostream& operator<<(ostream& s,matrix<ElType>& m);
	friend void ludcmp(matrix<ElType>& a,vector<int>& indx,double &d);
	friend void lubksb(matrix<ElType>&a,vector<int>& indx,vector<ElType>&b);
};
template <class ElType>
matrix<ElType>::matrix()
{
	m = new vector<ElType>[DEFAULT_ALLOC];
	assert(m !=0);
	rows=cols=DEFAULT_ALLOC;
	for(int i=0;i<rows;i++)
	{
		vector<ElType> v;
		m[i]= v;
	}
}

template <class ElType>
matrix<ElType>::matrix(int r, int c)
{
	m= new vector<ElType>[r];
	assert(m != 0);
	rows=r;
	cols=c;
	for(int i=0;i<r;i++)
	{
		vector<ElType> v(cols);
		m[i]=v;
	}
}
template <class ElType>
matrix<ElType>::matrix(const matrix<ElType> &s)
{
	int i;
	rows=s.rows;
	m = new vector<ElType>[rows];
	assert(m!=0);
	cols =s.cols;
	for(i=0;i<rows;i++)
	{
	  m[i]=s.m[i];
	}
}
template <class ElType>
matrix<ElType>::~matrix()
{
	delete [] m;
}

template <class ElType>
matrix<ElType>& matrix<ElType>::operator =(const matrix<ElType> &s)
{
	if(this != &s)
	{
		delete []m;
		rows= s.rows;
		cols=s.cols;
		m = new vector<ElType>[rows];
		assert(m !=0);
		for(int i=0;i<rows;i++) m[i]=s.m[i];
	}
	return *this;
}
template <class ElType>
vector<ElType>& matrix<ElType>::operator[](const int i)
{
	assert(i>=0 && i < rows);
	return m[i];
}
template <class ElType>
vector<ElType> matrix<ElType>::operator*(const vector<ElType>& v)
{
	int i,j;
	assert(cols == v.len);
	vector<ElType> ans(rows);
	for(i=0;i<rows;i++)
	{
		ans.data[i]=0.0;
		for(j=0;j<cols;j++) ans.data[i] += m[i][j]*v.data[j];
	}
	return ans;
}
template <class ElType>
matrix<ElType> operator*(const ElType& x,const matrix<ElType>& s)
{
	matrix<ElType> ans(s.rows,s.cols);
	for(int i=0;i<ans.rows;i++)
	  {
		ans.m[i]= x*s.m[i];
	  }
	return ans;
}
template <class ElType>
matrix<ElType> matrix<ElType>::transpose()
{
  matrix<ElType> ans(cols,rows);
  for(int i=0;i<rows;i++)
	{
	  for(int j=0;j<cols;j++) ans[j][i]=m[i][j];
	  }
     return ans;
}
template <class ElType>
matrix<ElType> operator*(const matrix<ElType>& s,const ElType& x)
{
	matrix<ElType> ans(s.rows,s.cols);
	for(int i=0;i<ans.rows;i++)
	  {
		ans.m[i]= x*s.m[i];
	  }
	return ans;
}
template <class ElType>
matrix<ElType>  matrix<ElType> ::operator*(const matrix<ElType>&  a)
{

	assert(cols == a.rows);

	matrix<ElType>  ans(rows,a.cols);
	for(int i=0;i<rows;i++)
	{
		for(int j=0;j<a.cols;j++)
		{
			ans.m[i][j]=0.0;
			for(int k=0;k<cols;k++)
			{
				ans.m[i][j] += m[i][k]*a.m[k][j];
			}
		}
	}
	return ans;
}
template <class ElType>
matrix<ElType>  matrix<ElType> ::operator+(const matrix<ElType> & a)
{
	int i,j;

	assert(rows== a.rows);
	assert(cols== a.cols);

	matrix<ElType>  ans(a.rows,a.cols);
	for(i=0;i<a.rows;i++)
	{
		for(j=0;j<a.cols;j++)
		  {
			ans.m[i][j] = m[i][j] + a.m[i][j];  //faster than assigning vectors?
		}
	}
	return ans;
}
template <class ElType>
matrix<ElType> matrix<ElType>::operator-(const matrix<ElType>& a)
{
	int i,j;
	assert(rows == a.rows);
	assert(cols == a.cols);
	matrix ans(rows,cols);
	for(i=0;i<rows;i++)
	{
		for(j=0;j<cols;j++)
		ans.m[i][j] = m[i][j] - a.m[i][j];
	}
	return ans;
}
template <class ElType>
ostream& operator<<(ostream& s,matrix<ElType>& m)
{
	for(int i=0; i<m.rows;i++) s << m[i];
	return s;
}

#define TINY 1.0e-20;
//we assume fabs(ElType) is defined
//assignment of doubles to ElType is defined
template <class ElType>
void ludcmp(matrix<ElType>& a, vector<int>& indx,double& d)
{
	int i,imax,j,k;
	ElType  big,dum,sum,temp;
	int n=a.rows;
	vector<ElType> vv(n);
	assert(a.rows == a.cols);
	d=1.0;
	for (i=0;i<n;i++)
	{
		big=0.0;
		for (j=0;j<n;j++) if ((temp=fabs(a[i][j])) > big) big=temp;
		if (big == 0.0) cerr << "Singular matrix in routine LUDCMP" << endl;
		vv[i]=1.0/big;
	}
	for (j=0;j<n;j++)
	{
		for (i=0;i<j;i++)
		{
			sum=a[i][j];
			for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
			a[i][j]=sum;
		}
		big=0.0;
		for (i=j;i<n;i++)
		{
			sum=a[i][j];
			for (k=0;k<j;k++) sum -= a[i][k]*a[k][j];
			a[i][j]=sum;
			if ( (dum=vv[i]*fabs(sum)) >= big)
			{
				big=dum;
				imax=i;
			}
		}
		if (j != imax)
		{
			for (k=0;k<n;k++)
			{
				dum=a[imax][k];
				a[imax][k]=a[j][k];
				a[j][k]=dum;
			}
			d = -(d);
			vv[imax]=vv[j];
		}
		indx[j]=imax;
		if (a[j][j] == 0.0) a[j][j]=TINY;
		if (j != n-1) {
			dum=1.0/(a[j][j]);
			for (i=j+1;i<n;i++) a[i][j] *= dum;
		}
	}
}
#undef TINY
template <class ElType>
void lubksb(matrix<ElType>& a,vector<int>& indx,vector<ElType>& b)
{
	int i,ip,j;
	ElType sum;
	int n=a.rows;
	for (i=0;i<n;i++)
	{
		ip=indx[i];
		sum=b[ip];
		b[ip]=b[i];
		for (j=0;j<=i-1;j++) sum -= a[i][j]*b[j];
		b[i]=sum;
	}
	for (i=n-1;i>=0;i--)
	{
		sum=b[i];
		for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
		b[i]=sum/a[i][i];
	}
}