mymath.c

一个简单的矩阵类

#include "MyMath.h"

CMatrix::CMatrix()
{
	Data=NULL;
}

CMatrix::CMatrix(int s)
{
	m=s;
	n=s;
	Data=new double[m*m];
	for(int i=0;i<m*m;i++)
		Data[i]=0;

}

CMatrix::CMatrix(int s,int t)
{

	m=s;
	n=t;
	Data=new double[m*n];
	for(int i=0;i<m*n;i++)
		Data[i]=0;
}

CMatrix::CMatrix(int s,double t)
{
	m=s;
	n=s;
	Data=new double[m*m];
	for(int i=0;i<m*m;i++)
		Data[i]=0;
	for(i=0;i<m;i++)
		Data[i*m+i]=t;
}

CMatrix::CMatrix(CMatrix &m1)
{
/*	if(Data!=NULL)
	{
		delete[] Data;
	}*/
	m=m1.m;
    n=m1.n;
    Data=new double[m*n];
	for(int i=0;i<m*n;i++)
		Data[i]=m1.Data[i];
}

CMatrix::~CMatrix()
{
	if(Data!=NULL)
		delete[] Data;
}

CMatrix CMatrix :: operator +(CMatrix &m1)
{
	CMatrix m2(m,n);
	for(int i=0;i<m*n;i++)
		m2.Data[i]=Data[i]+m1.Data[i];
    return m2;
}

CMatrix CMatrix :: operator -(CMatrix &m1)
{
	CMatrix m2(m,n);
	for(int i=0;i<m*n;i++)
		m2.Data[i]=Data[i]-m1.Data[i];
    return m2;
}

CMatrix CMatrix :: operator *(CMatrix &m1)
{
	CMatrix m2(m,m1.n);
	
	int i,j,k;
	double temp=0;
    for( i=0;i<m;i++)
	{
    	 for(j=0;j<m1.n;j++)
		 {
    		 temp=0;
	    	 for(k=0;k<n;k++)
    			 temp+=Data[i*n+k]*m1.Data [k*m1.n+j];
			 m2.Data[i*m1.n+j]=temp;
		 }
	}
	return m2;
}

CMatrix CMatrix :: operator *(double f)
{
	CMatrix m2(m,n);
    for(int i=0;i<m*n;i++)
	{
    	m2.Data [i]=Data[i]*f;
	}
	return m2;
}



CMatrix CMatrix :: operator / (CMatrix &m1)
{
int i,j,k,s,t,dVec,dCol;
//    if(m==1)
  //    return true;
	double f;
	t=m1.m ;
	if(t==1)
	{
		if(m1.Data[0]==0)
			return 0;
		CMatrix m5(m,n);
		for(i=0;i<m;i++)
			for(j=0;j<n;j++)
				m5.Data[i*n+j] = Data[i*n+j]/m1.Data[0];
		return m5;
	}

    double vMod=1;
    vMod=Mod(m1);
	CMatrix m2(t-1),m3(t),m4(m,n);
    if(vMod==0)
  	   return 0;
    for(i=0;i<t;i++)
    	  for(j=0;j<t;j++)
		  {    
		  for(k=0;k<t-1;k++)
			  for(s=0;s<t-1;s++)
			  {
				  if(k<i)
					  dVec=0;
				  else 
					  dVec=1;
				  if(s<j)
					  dCol=0;
				  else
					  dCol=1;
				  m2.Data [k*(t-1)+s]=m1.Data [(k+dVec)*t+s+dCol];
			  }
		  if((i+1+j+1)%2)
		  {
			  f=Mod(m2);
	          m3.Data [j*t+i]=-f/vMod;
		  }
		  else
		  {
			  f=Mod(m2);
			  m3.Data [j*t+i]=f/vMod;
		  }
	  }
  

	double temp=0;
    for( i=0;i<m;i++)
	{
    	 for(j=0;j<n;j++)
		 {
    		 temp=0;
	    	 for(k=0;k<n;k++)
    			 temp+=Data[i*n+k]*m3.Data [k*n+j];
			 m4.Data[i*n+j]=temp;
		 }
	}
  return m4;
}


double Mod(CMatrix &m1)
{
	int i,j,k;
	int t=m1.m;
	double temp=0;
	if(t ==1)
		return double(m1.Data [0]);
	CMatrix m2(t-1);
	for(i=0;i<t;i++)
	{  
       for(j=0;j<t-1;j++)
	      if(j<i)
			  for(k=0;k<t-1;k++)
		    	   m2.Data [k*t-k+j]=m1.Data [k*t+t+j];
		  else
              for(k=0;k<t-1;k++)
		    	   m2.Data [k*t-k+j]=m1.Data [k*t+t+j+1];
       if(i%2)
		   temp-=m1.Data [i]*Mod(m2);
	   else
		   temp+=m1.Data [i]*Mod(m2);
   
	}
    return temp;
}

CMatrix CMatrix :: operator ~()
{
	CMatrix m1(n,m);
	for(int i=0;i<m;i++)
		for(int j=0;j<n;j++)
			m1.Data[j*m+i]=Data[i*n+j];
	return m1;
}


CMatrix& CMatrix :: operator = (CMatrix &m1)
{
	if(this==&m1)
		return *this;
	if(Data!=NULL)
	{
		delete[] Data;
	}
	m=m1.m;
    n=m1.n;
    Data=new double[m*n];
	
	for(int i=0;i<m*n;i++)
		Data[i]=m1.Data[i];
	return *this;
}

void CMatrix :: Set(int s,int t)
{
	m=s;
	n=t;
	Data=new double[s*t];
	for(int i=0;i<m*n;i++)
		Data[i]=0;
}


CMatrix Chol(CMatrix &m1)
{
	double temp=0;
	CMatrix m2(m1.m,m1.n);
	
	for(int j=0;j<m1.n;j++)
	{
		temp = m1.Data [j*m1.n+j];
		for(int k=0;k<j;k++)
			temp -= m2.Data [j*m1.n+k]*m2.Data [j*m1.n+k];
		m2.Data [j*m1.n+j] = sqrt(temp);
		for(int i=j+1;i<m1.m;i++)
		{
			temp = m1.Data [i*m1.n+j];
			for(k=0;k<j;k++)
			{
				temp -= m2.Data [i*m1.n+k]*m2.Data [j*m1.n+k];
			}
			m2.Data [i*m1.n+j]=temp/m2.Data [j*m1.n+j];
		}
	}
	return m2;

}


double Uniform(double a,double b, long &seed)
{
	double t;
	seed=2045*(seed)+1;
	seed=seed-(seed/1048576)*1048576;
	t=(seed)/1048576.0;
	t=a+(b-a)*t;
	return t;
}

double Gauss(double mean,double sigma, long &seed)
{
	int i;
	double x,y;
	x=0;
	for(i=0;i<12;i++)
		x+=Uniform(0.0, 1.0, seed);
	x=x-6.0;
	y=mean+x*sigma;
	return y;
}