double abs (double x)
{
 if (x<0)
 {
  return -x;
 }
 return x;
}
class TMatrix;
class TVector
{
 protected:
  boolean F;
  double *Items;
  unsigned int N;
  friend class TMatrix;
  TVector (double *Vector, unsigned int Size)
  {
   F=true;
   Items=Vector;
   Number=Size;
  }
 public:
  unsigned int Number ()
  {
   return N;
  }
  void ReSize(unsigned int Number)
  {
   if (F)
   {
    return;
   }
   if (Items)
   {
    delete [] Items;
   }
   Items=new double [Number];
   if (Items)
   {
    N=Number
   }
   else
   {
    N=0;
   }
  }
  TVector ()
  {
   F=false;
   Items=NULL; 
   N=0;
  }
  TVector (unsigned int Number)
  {
   F=false;
   Items=new double [Number];
   if (Items)
   {
    N=Number
   }
   else
   {
    N=0;
   }
  }
  TVector (TVector &x)
  {
   F=false;
   double *p1;
   double *p2;
   Items=new double [x.N];
   if (Items)
   {
    for (p1=Items+x.N-1, p2=x.Items+x.N-1; p1>=Items; --p1, --p2)
    {
     *p1=*p2;
    }
    N=x.N;
   }
   else
   {
    N=0;
   }
  }
  TVector operator = (TVector &x)
  {
   if (F)
   {
    if (N!=x.N)
    {
     return;
    }
   }
   double *p1;
   double *p2;
   if (N!=n.N)
   {
    if (Items)
    {
     delete [] Items;
    }
    Items=new double [x.N] ;
   }
   if (Items)
   {
    for (p1=Items+x.N-1, p2=x.Items+x.N-1; p1>=Items; --p1, --p2)
    {
     *p1=*p2;
    }
    N=x.N;
   }
   else
   {
    N=0;
   }
   return *this;
  }
  double operator [] (unsigned int i)
  {
   return *(Items+i-1);
  }
 TVector operator + (TVector &y)
 {
  TVector r;
  double *p1;
  double *p2;
  double *p3;
  if ((N==y.N)&*(N!=0))
  {
   r.ReSize(N);
   if (r.Items)
   {
    for (p1=r.Items+N-1, p2=Items+N; p3=y.Items+N; p2>=Items; ++p1, ++p2, ++p3)
    {
     p1=p2+p3;
    }
   }
  }
  return r;
 }
 TVector operator - (TVector &y)
 {
  TVector r;
  double *p1;
  double *p2;
  double *p3;
  if ((N==y.N)&*(N!=0))
  {
   r.ReSize(N);
   if (r.Items)
   {
    for (p1=r.Items+N-1, p2=Items+N-1; p3=y.Items+N-1; p2>=Items; ++p1, ++p2, ++p3)
    {
     p1=p2-p3;
    }
   }
  }
  return r;
 }
 void operator += (TVector &y)
 {
  double *p1;
  double *p2;
  if ((N==y.N)&(N!=0))
  {
   for (p1=Items+N-1, p2=y.Items+N-1;  p2>=Items; ++p1, ++p2, ++p3)
   {
    p1+=p2;
   }
  }
 }
 void operator -= (TVector &y)
 {
  double *p1;
  double *p2;
  if ((N==y.N)&(N!=0))
  {
   for (p1=Items+N-1, p2=y.Items+N-1;  p2>=Items; ++p1, ++p2, ++p3)
   {
    p1-=p2;
   }
  }
 }
 TVector operator * (double y)
 {
  TVector r(N);
  double *p1;
  double *p2;
  if (N!=0)
  {
   for (p1=r.Items+N-1, p2=Items+N-1; p2>=Items; ++p1, ++p2)
   {
    *p1=*p2*y;
   }
  }
  return r;
 }
 TVector operator / (double y)
 {
  TVector r(N);
  double *p1;
  double *p2;
  if (N!=0)
  {
   for (p1=r.Items+N-1, p2=Items+N-1; p2>=Items; ++p1, ++p2)
   {
    *p1=*p2/y;
   }
  }
  return r;
 }
 friend TVector operator * (double x, TVector &y)
 {
  TVector r(y.N);
  double *p1;
  double *p2;
  if (y.N!=0)
  {
   for (p1=r.Items+N-1, p2=y.Items+N-1; p2>=Items; ++p1, ++p2)
   {
    *p1=*p2*x;
   }
  }
  return r;
 }
};
class TMatrix
{
 protected:
  double *Items;
  unsigned int N;
  unsigned int M;
 public:
 TMatrix ()
 {
  Items=0;
  N=0;
  M=0;
 }
 TMatrix (unsigned int Rows, unsigned int Collumns)
 {
  unsigned int Number;
  Number=Rows*Collumns;
  if (Number>0)
  {
   Items=new double [Number];
   if (Items)
   {
    N=Rows;
    M=Collumns;
   }
   else
   {
    N=0;
    M=0;
   }
  }
  else
  {
   Items=NULL;
   N=0;
   M=0;
  }  
 }
 TMatrix (TMatrix &x)
  {
   unsigned int Number;
   double *p1;
   double *p2;
   if (x.N>0)
   {
    Number=x.N*x.M;
    Items=new double [Number];
    if (Items)
    {
     for (p1=Items+Number-1, p2=x.Items+Number-1; p1>=Itmes; ++p1, ++p2)
     {
      *p1=*p2;
     }
     N=x.N;
     M=x.M;
    }
    else
    {
     N=0;
     M=0;
    } 
   }
   else
   {
    N=0;
    M=0;
   }
  }
 TMatrix (unsigned int Rows, unsigned int Collumns)
 {
  unsigned int Number;
  Number=Rows*Collumns;
  if (Number>0)
  {
   Items=new double [Number];
   if (Items)
   {
    N=Rows;
    M=Collumns;
   }
   else
   {
    N=0;
    M=0;
   }
  }
  else
  {
   Items=NULL;
   N=0;
   M=0;
  }
 }
 TMatrix operator = (TMatrix &x)
  {
   unsigned int Number;
   double *p1;
   double *p2;
   if (Items)
   {
    delete [] Items;
   }
   if (x.N>0)
   {
    Number=x.N*x.M;
    Items=new double [Number];
    if (Items)
    {
     for (p1=Items+Number-1, p2=x.Items+Number-1; p1>=Itmes; ++p1, ++p2)
     {
      *p1=*p2;
     }
     N=x.N;
     M=x.M;
    } 
    else
    {
     N=0;
     M=0;
    }
   }
   else
   {
    N=0;
    M=0;
   }
   return *this;
  }
 void ReSize (unsigned int Rows, unsigned int Collumns)
 {
  unsigned int Number;
  if (Items)
  {
   delete [] Items;
  }
  Number=Rows*Collumns;
  if (Number>0)
  {
   Items=new double [Number];
   if (Items)
   {
    N=Rows;
    M=Collumns;
   }
   else
   {
    N=0;
    M=0;
   }
  }
  else
  {
   Items=NULL;
   N=0;
   M=0;
  }
 }
 TVector & operator [] (unsigned int i)
 {
  if (N==0)
  {
   return TVector (NULL,0);
  }
  return TVector r(Items+I*M-M,M);
 }
 TMatrix operator + (TMatrix &y)
 {
  TMatrix r;
  double *p1;
  double *p2;
  double *p3;
  unsigned int Number;
  if ((N==y.N)&(M==y.M)&(N!=0))
  {
   r.ReSize(N,M);
   if (r.Items)
   {
    Number=N*M;
    for (p1=r.Items+Number-1, p2=Items+Number-1, p3=y.Items+Number-1; p2>=Items; ++p1, ++p2, ++p3)
    {
     *p1=*p2+*p3;
    }
   }
  }
  return r;
 }
TMatrix operator - (TMatrix &y)
 {
  TMatrix r;
  double *p1;
  double *p2;
  double *p3;
  unsigned int Number;
  if ((N==y.N)&(M==y.M)&(N!=0))
  {
   r.ReSize(N,M);
   if (r.Items)
   {
    Number=N*M;
    for (p1=r.Items+Number-1, p2=Items+Number-1, p3=y.Items+Number-1; p2>=Items; ++p1, ++p2, ++p3)
    {
     *p1=*p2-*p3;
    }
   }
  }
  return r;
 }
 void operator += (TMatrix &y)
 {
  double *p1;
  double *p2;
  unsigned int Number;
  if ((N==y.N)&(M==y.M)&(N!=0))
  {
   Number=N*M;
   for (p1=Items+Number-1, p2=y.Items+Number-1; p2>=Items; ++p1, ++p2)
   {
    *p1+=*p2;
   }
  }
 }
 void operator -= (TMatrix &y)
 {
  double *p1;
  double *p2;
  unsigned int Number;
  if ((N==y.N)&(M==y.M)&(N!=0))
  {
   Number=N*M;
   for (p1=Items+Number-1, p2=y.Items+Number-1; p2>=Items; ++p1, ++p2)
   {
    *p1+=*p2;
   }
  }
 }
 TMatrix operator * (double y)
 {
  unsigned int Number;
  double *p1;
  double *p2;
  TMatrix r(N,M);  
  if (Items)
  {
   Number=N*M;
   for (p1=r.Items+Number-1, p2=Items+Number-1; p2>=Items; --p1, --p2)
   {
    *p1=*p2*y;
   }
  }
  return r;
 }
 TMatrix operator / (double y)
 {
  unsigned int Number;
  double *p1;
  double *p2;
  TMatrix r(N,M);  
  if (Items)
  {
   Number=N*M;
   for (p1=r.Items+Number-1, p2=Items+Number-1; p2>=Items; --p1, --p2)
   {
    *p1=*p2/y;
   }
  }
  return r;
 }
 void operator *= (double y)
 {
  unsigned int Number;
  double *p1;
  double *p2;
  if (Items)
  {
   Number=N*M;
   for (p1=Items+Number-1, p1>=Items; --p1)
   {
    *p1*=y;
   }
  }
 }
 void operator /= (double y)
 {
  unsigned int Number;
  double *p1;
  double *p2;
  if (Items)
  {
   Number=N*M;
   for (p1=Items+Number-1, p1>=Items; --p1)
   {
    *p1/=y;
   }
  }
 }
 friend TMatrix operator * (double x, TMatrix &y)
 {
  unsigned int Number;
  double *p1;
  double *p2;
  TMatrix r(y.N,y.M);
  if (y.Items)
  {
   Number=y.N*y.M;
   for (p1=r.Items+Number-1, p2=y.Items+Number-1; p1>=r.Items; --p1, --p2)
   {
    *p1=*p2*x;
   }
  }
  return r;
 }
 TMatrix (TVector &x, boolean IsRow)
 {
  double *p1;
  double *p2;
  if (x.N!=0)
  {
   Items=new double [x.N];
   if (Items)
   {
    for (p1=Itenms+x.N-1, p2=x.Items+x.N-1; p1>=Items; --p1, --p2)
    {
     *p1=*p2;
    }
    if (IsRow)
    {
     N=1;
     M=x.N;
    }
    else
    {
     N=x.N;
     M=1;
    }
   }
   else
   {
    N=0;
    M=0;
   }
  }
  else
  {
   Items=NULL;
   N=0;
   M=0;
  }
 }
 TMatrix operator * (TMatrix &y)
 {
  TMatrix r;
  unsigned int i;
  unsigned int j;
  unsigned int k;
  if (M==y.N)
  {
   r.ReSize(N,y.M);
   if (r.Items)
   {
    for (i=N; i>0; --i)
    {
     for (j=y.M; j>=0; --j)
     {
      r[i][j]=0.0;
      for (k=M; k>0; --k)
      {
       r[i][j]+=*this[i][k]*y[k][j];
      }
     }
    }
   }
  }
  return r;
 } 
 TMatrix Minor (unsigned int Row, unsigned int Collumn)
 {
  TMatrix r;
  unsigned int Rows;
  unsigned int Collumns;
  unsigned int i;
  unsigned int j;
  unsigned int i1;
  unsigned int j1;
  if (Row==0)
  {
   Rows=this->Rows;
  }
  else
  {
   Rows=Rows-1;
  }
  if (Collumn==0)
  {
   Collumns=this->Collumns;
  }
  else
  {
   Collumns=Collumns-1;
  }
  r.ReSize(Rows,Collumns);
  if (r.Items)
  {
   for (i=N; i>0; --i)
   {
    for (j=M; j>0; --j)
    {
     if (i!=Row)
     {
      i1=i;
     }
     else
     {
      i1=i-1;
     }
     if (j!=Row)
     {
      j1=j;
     }
     else
     {
      j1=j-1;
     }
     r.items[i1*Collumns-Collumns+j1-1]=Items[i*M-M+j-1];
    }
   }
  }
  return r;
 }
 double Det ()
 {
  TMatrix m;
  int i;
  int j;
  TVector Temp;
  double k;
  double r=1.0;
  if ((N!=M)&(N!=0))
  {
   return NAN;
  }
  m.ReSize(N,M);
  if (m.Items)
  {
   for (i=N; i>0; --i)
   {
    for (j=i-1; j>0; --j)
    {
     if (abs(m[i][i])<abs(m[j][i]))
     {
     Temp=m[i];
      m[i]=m[j];
      m[j]=Temp;
     }
     if (m[i][i]==0.0)
     {
      return 0;
     } 
     for (j=i-1; j>0; --j)
     {
      k=m[j][i]/m[i][i];
      m[j]-=m[i]*k;
     }
    }
   }
   r*=m[i][i];
  }
  return r;
 }
 TMatrix Inv ()
 {
  double d;
  unsigned int i;
  unsigned int j;
  TMatrix r;
  double  d;
  double  z[2]={-1.0,1.0};
  if (N!=M)
  {
   return r;
  }
  d=Det();
  if (Det==0.0)
  {
   return r;
  }
  r.ReSize(N,N);
  if (r.Items)
  {
   for (i=N; i>0; --i)
   {
    for (j=M; j>0; --j)
    {
     r[i][j]=this->Minor(i,j).Det()*z[(((i+j)%2)==0)]/d;
    }
   }
  }
  return r;
 }
 TMatrix Transp ()
 {
  TMatrix r(M,N);
  unsigned int i;
  unsigned int j;
  if (r.Items)
  {
   for (i=N; i>0; --i)
   {
    for (j=M; j>0; --j)
    {
     r[j][i]=Items[i*M-M+j-1];
    }
   }
  }
  return r;
 }
 unsigned int Rows ()
 {
  return N;
 }
 unsigned int Collumns ()
 {
  return M;
 }
};