jz.cpp

距阵算法C++编写的lib库

// jz.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "mjz.h"
#include <stdio.h>
#include <math.h>
//#include "stdafx.h"
#include <string.h>
#include "mjz.h"

Mjz::Mjz(int m,int n)
{
	double *M;
	if(m>JZMax)	row=JZMax;else	row=m;
	if(n>JZMax)	col=JZMax;else	col=n;
	size=row*col;
	M=new double[size];
	pHead=M;
	int i,j;
	for (i=0;i<row;i++)
		for (j=0;j<col;j++)
			*M++=0.0;
}
Mjz::Mjz(int m)
{
	double *M;
	if(m>JZMax)	row=JZMax;else	row=m;
	col=row;
	size=row*col;
	M=new double[size];
	pHead=M;
	int i,j;
	for (i=0;i<row;i++)
		for (j=0;j<col;j++)
			*M++=0.0;
}

bool Mjz::SetAt(int n,double data)
{
	//if (IsIn(n))
	//{
		int i=0;
		double *M;
		M=pHead;
		while(i++ < n) M++;
		*M=data;
		return true;
	//}
	//return false;
}
bool Mjz::SetAt(int m,int n,double data)
{
	//if (IsIn(n))
	//{
		int i=0,mn;
		double *M;
		M=pHead;
		mn=m*row+n;
		while(i++ < mn) M++;
		*M=data;
		return true;
	//}
	//return false;
}

bool Mjz::GetAt(int n,double &data)
{
	//if (IsIn(n))
	//{
		int i=0;
		double *M;
		M=pHead;
		while(i++ < n) M++;
		data=*M;
		return true;
	//}
	//return false;
}
bool Mjz::GetAt(int m,int n,double &data)
{
	//if (IsIn(n))
	//{
		int i=0,mn;
		double *M;
		M=pHead;
		mn=m*row+n;
		while(i++ < mn) M++;
		data=*M;
		return true;
	//}
	//return false;
}
double Mjz::GetAt(int n)
{
	//if (IsIn(n))
	//{
		int i=0;
		double *M;
		M=pHead;
		while(i++ < n) M++;
		return *M;
	//}
	//return RetErr;
}
double Mjz::GetAt(int m,int n)
{
	//if (IsIn(n))
	//{
		int i=0,mn;
		double *M;
		M=pHead;
		mn=m*row+n;
		while(i++ < mn) M++;
		return *M;
	//}
	//return RetErr;
}
void Mjz::I()
{
	double *M;
	int i,j;
	size=row*row;
	M=new double[size];
	pHead=M;
	for (i=0;i<row;i++)
		for (j=0;j<row;j++)
			if(i==j)
				*M++=1.0;
			else
				*M++=0.0;
}

void Mjz::I(int m)
{
	double *M;
	if(m>JZMax)	row=JZMax;else	row=m;
	col=row;
	size=row*row;
	M=new double[size];
	pHead=M;
	int i,j;
	for (i=0;i<row;i++)
		for (j=0;j<col;j++)
			if(i==j)
				*M++=1.0;
			else
				*M++=0.0;
}
void Mjz::free()
{
	//M=pHead;
	try {
		if (pHead)
		{
		delete((double*)pHead);
		row=0;
		col=0;
		size=0;
		pHead=NULL;
		}
		//delete this;
	} catch (...) {
		;
	}
	
}

bool Mjz:: ReadTextFile(char* filename)
{
	FILE *f;
	if((f=fopen(filename,"r"))==NULL) {printf("\nErr:ReadFile[%s]",filename);return false;}
	int H,L,i,j;
	double data;
	fscanf(f,"%i %i",&H,&L);
	if (H>JZMax)H=JZMax;
	if (L>JZMax)L=JZMax;
	printf("H=%i,L=%i",H,L);
//	Mjz jz(H,L);
	if(H>JZMax)	row=JZMax;else	row=H;
	if(L>JZMax)	col=JZMax;else	col=L;
	size=row*col;
	pHead=new double[size];
	////
	for(i=0;i<H;i++)
		for(j=0;j<L;j++)
		{
			fscanf(f,"%lf",&data);
			SetAt(i*H+j,data);
		//	printf("H=%.3f",data);
		};
//	*this=jz;
	fclose(f);
	return true;
}
bool Mjz:: WriteTextFile(char* filename)
{
	FILE *f;
	if((f=fopen(filename,"w+"))==NULL) return false;
	int i,j;
	double data;
	fprintf(f,"%i %i\n",row,col);
	for(i=0;i<row;i++)
	{
		fprintf(f,"%s\n","");
		for(j=0;j<col;j++)
		{
			data=GetAt(i*row+j);
			fprintf(f,"%10.3lf",data);
		//	printf("H=%.3f",data);
		}
	}
	fclose(f);
	return true;
}

const Mjz& Mjz ::operator=(double x)
{
	int i,j;
	for (i=0;i<row;i++)
		for (j=0;j<col;j++)
			SetAt(i*row+j,x);
	return (Mjz) *this;
}

const Mjz& Mjz ::operator=(Mjz MjzSrc)
{
	int n,i,j;
	double data;
	Mjz M=Mjz(MjzSrc.row,MjzSrc.col);
	for (i=0;i<row;i++)
		for (j=0;j<col;j++)
		{
			n=i*row+j;
			data=MjzSrc.GetAt(n);
			M.SetAt(n,data);
		}
return M;
}
const Mjz& Mjz ::Set(double M[][DefCol])
{
	int n,i,j;
	double data;
	//Mjz(MjzSrc.row,MjzSrc.col);
	for (i=0;i<row;i++)
		for (j=0;j<col;j++)
		{
			n=i*row+j;
			data=M[i][j];
			SetAt(n,data);
		}
return (Mjz) *this;
}
const Mjz& Mjz ::operator+(double A)
{
	Mjz Mjzb(row,col);
	int n,i,j;
	double data;
	for (i=0;i<row;i++)
		for (j=0;j<col;j++)
		{
			n=i*row+j;
			data=GetAt(n)+A;
			Mjzb.SetAt(n,data);
		}

return Mjzb;
}
const Mjz& Mjz ::operator+(Mjz MjzSrc)
{
Mjz Mjzb(row,col);
if(RowColEq(MjzSrc))
{
	int n,i,j;
	double data;
	for (i=0;i<row;i++)
		for (j=0;j<col;j++)
		{
			n=i*row+j;
			data=GetAt(n)+MjzSrc.GetAt(n);
			Mjzb.SetAt(n,data);
		}
}
return Mjzb;
}
const Mjz& Mjz ::operator-(double A)
{
	Mjz Mjzb(row,col);
	int n,i,j;
	double data;
	for (i=0;i<row;i++)
		for (j=0;j<col;j++)
		{
			n=i*row+j;
			data=GetAt(n)-A;
			Mjzb.SetAt(n,data);
		}

return Mjzb;
}
const Mjz& Mjz ::operator-(Mjz MjzSrc)
{
Mjz Mjzb(row,col);
if(RowColEq(MjzSrc))
{
	int n,i,j;
	double data;
	for (i=0;i<row;i++)
		for (j=0;j<col;j++)
		{
			n=i*row+j;
			data=GetAt(n)-MjzSrc.GetAt(n);
			Mjzb.SetAt(n,data);
		}
}
return Mjzb;
}
const Mjz& Mjz ::operator*(double A)
{
	Mjz Mjzb(row,col);
	int n,i,j;
	double data;
	for (i=0;i<row;i++)
		for (j=0;j<col;j++)
		{
			n=i*row+j;
			data=GetAt(n)*A;
			Mjzb.SetAt(n,data);
		}
return Mjzb;
}

const Mjz& Mjz ::operator*(Mjz MjzSrc)
{
Mjz Mjzb(row,MjzSrc.col);
if(ColEq(MjzSrc))
{
	int n,al,bl,i,j,l,COL,Col;
	double A,B,data;
	COL=col;
	Col=MjzSrc.col;
    for (i=0;i<row;i++)
	  for (j=0;j<Col;j++)
	  {
		   data=0.0;
			for (l=0;l<COL;l++)
			{
				al=i*row+l;
				bl=l*row+j;
				A=GetAt(al);
				B=MjzSrc.GetAt(bl);
				data+=A*B;
			}
			n=i*row+j;
			Mjzb.SetAt(n,data);
	  }
}
//return (Mjz) *this;
return Mjzb;
}

const Mjz& Mjz ::operator/(double A)
{
	int n,i,j;
	double data;
	Mjz Mjzb(row,col);
	for (i=0;i<row;i++)
		for (j=0;j<col;j++)
		{
			n=i*row+j;
			data=GetAt(n)/A;
			Mjzb.SetAt(n,data);
		}
return Mjzb;
}

bool Mjz::RowEq (int n,Mjz MjzSrc)
{
	if(!ColEq(MjzSrc)) return false;
	int i,j;
	for(i=0;i<col;i++)
	{
		j=n*row+i;
		if (GetAt(j)!=MjzSrc.GetAt(j)) return false;
	}
	return true;
}
bool Mjz::RowEq (int n,double x)
{
	int i;
	double d;
	for(i=0;i<col;i++)
	{
		d=GetAt(n*row+i);
		printf("\nGetAt(%i,%i)=%.3f",n,i,d);
		if (d!=x) return false;
	}
	printf("\nGetAt==1");
	return true;
}
bool Mjz::ColEq (int n,double x)
{
	int i;
	for(i=0;i<row;i++)
		if (GetAt(i*row+n)!=x) return false;
	return true;
}
bool Mjz::ColEq (int n,Mjz MjzSrc)
{
	if(!RowEq(MjzSrc)) return false;
	int i,j;
	for(i=0;i<row;i++)
	{
		j=i*row+n;
		if (GetAt(j)!=MjzSrc.GetAt(j)) return false;
	}
	return true;
}
bool Mjz:: operator== (Mjz MjzSrc)
{
	if(!RowColEq(MjzSrc)) return false;
	int i,j,k;
	for(i=0;i<row;i++)
		for(j=0;j<col;j++)
		{
			k=i*row+j;
			if(GetAt(k)!=MjzSrc.GetAt(k)) return false;
		}

	return true;
}
bool Mjz:: operator!= (Mjz MjzSrc)
{
	if(!RowColEq(MjzSrc)) return true;
	int i,j,k;
	for(i=0;i<row;i++)
		for(j=0;j<col;j++)
		{
			k=i*row+j;
			if(GetAt(k)!=MjzSrc.GetAt(k)) return true;
		}

	return false;
}

void Mjz:: print(FILE *f)
{
	int i,j;
	fprintf(f,"\n***M[][]=");
	for (i=0;i<row;i++)
	{
		fprintf(f,"\nRow[%3i]:",i);
		for (j=0;j<col;j++)
		fprintf(f," %.3f",GetAt(i*row+j));
	}
	fprintf(f,"\n_______________");
}
//////////////////////Basic BH
void Mjz:: print(char *msg)
{
	int i,j;
	#ifdef _ads_h
		ads_printf("\n%s",msg);
	#else 
		printf("\n%s",msg);
	#endif
	#ifdef _ads_h
		ads_printf("\n***M[][]=");
	#else 
		printf("\n***M[][]=");
	#endif
	for (i=0;i<row;i++)
	{
		#ifdef _ads_h
			ads_printf("\nRow[%3i]:",i);
		#else 
			printf("\nRow[%3i]:",i);
		#endif
		for (j=0;j<col;j++)
			#ifdef _ads_h
				ads_printf(" %.3f",GetAt(i*row+j));
			#else 
				printf(" %.3f",GetAt(i*row+j));
			#endif
	}
	#ifdef _ads_h
			ads_printf("\n_______________");
	#else
			printf("\n_______________");
	#endif
}
int Mjz ::RowJH(int nn,int mm)
{
	if(nn==mm)return 1;
	int n,m;
	if(nn<0) n=0;else n=nn;
	if(mm<0) m=0;else m=mm;
	if(nn>=col) n=col-1;else n=nn;
	if(mm>=col) m=col-1;else m=mm;

	double tmp;
	int l,i,j;
	for(l=0;l<col;l++)
	{
		i=n*row+l;
		j=m*row+l;
		tmp=GetAt(j);
		SetAt(j,GetAt(i));
		SetAt(i,tmp);

	}
	return -1;
}
int Mjz ::ColJH(int nn,int mm)
{
	if(nn==mm)return 1;
	int n,m;
	if(nn<0) n=0;else n=nn;
	if(mm<0) m=0;else m=mm;
	if(nn>=col) n=col-1;else n=nn;
	if(mm>=col) m=col-1;else m=mm;

	double tmp;
	int h,i,j;
	for(h=0;h<row;h++)
	{
		i=h*row+n;
		j=h*row+m;
		tmp=GetAt(j);
		SetAt(j,GetAt(i));
		SetAt(i,tmp);

	}
	return -1;
}

double Mjz::RowMul(int nn,double x)
{
	if(x==0.0)return 1.0;
	int n;
	if(nn<0) n=0;else n=nn;
	if(nn>=row) n=row-1;else n=nn;
	
	double tmp;
	int h,i;
	for(h=0;h<col;h++)
	{
		i=n*row+h;
		tmp=GetAt(i)*x;
		SetAt(i,tmp);

	}
	return x;
}
double Mjz::ColMul(int nn,double x)
{
	if(x==0.0)return 1.0;
	int n;
	if(nn<0) n=0;else n=nn;
	if(nn>=col) n=col-1;else n=nn;
	
	double tmp;
	int h,i;
	for(h=0;h<row;h++)
	{
		i=h*row+n;
		tmp=GetAt(i)*x;
		SetAt(i,tmp);

	}
	return x;
}
void Mjz::RowAdd(int on,int sn,double x)///Row(o)=Row(o)+Row(s)*x;
{
	if((x==0.0) || (on==sn))return;
	int o,s;
	if(on<0) o=0;else o=on;
	if(sn<0) s=0;else s=sn;
	if(on>=row) o=row-1;else o=on;
	if(sn>=row) s=row-1;else s=sn;
	
	double tmp;
	int h,i,j;
	for(h=0;h<col;h++)
	{
		i=o*row+h;
		j=s*row+h;
		tmp=GetAt(i)+GetAt(j)*x;
		SetAt(i,tmp);
	}
}
void Mjz::ColAdd(int on,int sn,double x)///Row(o)=Row(o)+Row(s)*x;
{
	if((x==0.0) || (on==sn))return;
	int o,s;
	if(on<0) o=0;else o=on;
	if(sn<0) s=0;else s=sn;
	if(on>=row) o=col-1;else o=on;
	if(sn>=row) s=col-1;else s=sn;
	
	double tmp;
	int h,i,j;
	for(h=0;h<row;h++)
	{
		i=h*row+o;
		j=h*row+s;
		tmp=GetAt(i)+GetAt(j)*x;
		SetAt(i,tmp);
	}
}
/////////////////////
double Mjz::HLS()
{
	double x=1.0,A,B;
	Mjz M;
	M=*this;
	int i,j,k,a,b;
	for(i=0;i<row;i++)
	  {
		a=i*row+i;//M[i,i]
		A=M.GetAt(a);
		if(A==0.0)
		{
			k=i;
			while((A==0.0) && (k++<row))
			{
				a=k*row+i;//M[i,i]
				A=M.GetAt(a);
			}
			if((A==0.0) && (k>=row)) return 0.0;
			M.RowJH(i,k);
			x*=-1.0;
		}
		for(j=i+1;j<row;j++)
		{
			b=j*row+i;
			B=M.GetAt(b);
			if(B!=0.0)
			{
			B/=A;
			M.RowSub(j,i,B);
			}
		}
	x*=A;
	}
	return x;
}
void Mjz::HLS(FILE *f)
{
	char s[100];
	double x=1.0,A,B;
	Mjz M(row,col);
	M=*this;
	int i,j,k,a,b;
	strcpy(s,"");
	fprintf(f,"\n***Sore=");
	M.print(f);
	fprintf(f,"\n***HLS Program***");
	for(i=0;i<row;i++)
	  {
		a=i*row+i;//M[i,i]
		A=M.GetAt(a);
		if(A==0.0)
		{
			k=i;
			while((A==0.0) && (k++<row))
			{
				a=k*row+i;//M[i,i]
				A=M.GetAt(a);
			}
			if((A==0.0) && (k>=row)) return;
			fprintf(f,"\nRow[%i]<-->Row[%i]:",i,k);
			M.RowJH(i,k);
			M.print(f);
			x*=-1.0;
			if (!strcmp(s,""))
				sprintf(s,"%.1f",x);
			else
				sprintf(s,"%s*-1.0",s);

		}
		for(j=i+1;j<row;j++)
		{
			b=j*row+i;
			B=M.GetAt(b);
			if(B!=0.0)
			{
			B/=A;
			printf("\nM[%i][%i]=M[%i][%i](%.3f) - M[%i][%i] (%.3f) * %.3f ",j,i,j,i,B,i,i,A,B);
			M.RowSub(j,i,B);
			M.print(f);
			}
		}
	x*=A;
	if (!strcmp(s,""))
		sprintf(s,"%.1f*%.1f",x,A);
	else
		sprintf(s,"%s*%.1f",s,A);

	}
	fprintf(f,"\nHLS=%s=%.3lf",s,x);	
}
int Mjz::R()
{
	int r;
	double A,B;
	Mjz M(row,col);
	M=*this;
	int i,j,k,a,b;
	for(i=0;i<row;i++)
	  {
		a=i*row+i;//M[i,i]
		A=M.GetAt(a);
		if(A==0.0)
		{
			k=i;
			while((A==0.0) && (k++<row))
			{
				a=k*row+i;//M[i,i]
				A=M.GetAt(a);
			}
			//if((A==0.0) && (k>=row)) return 0;
			M.RowJH(i,k);
			//x*=-1.0;
		}
		for(j=i+1;j<row;j++)
		{
			b=j*row+i;
			B=M.GetAt(b);
			if(B!=0.0)
			{
			B/=A;
			M.RowSub(j,i,B);
			}
		}
	//x*=A;
	}
	M.print();
	int rH=0,rL=0;
	for(i=0;i<row;i++)
	{
		printf("\nH[%i]==%i",i,M.RowEq(i,0));
		if(!M.RowEq(i,0.0)) rH+=1;
	}
	for(i=0;i<col;i++)
	{
		printf("\nL[%i]==%i",i,M.ColEq(i,0));
		if(!M.ColEq(i,0)) rL+=1;
	}

	printf("\nH=%i",row);
	printf("\nL=%i",col);
	printf("\nrH=%i",rH);
	printf("\nrL=%i",rL);
	r=min(rH,rL);
	rH=min(row,col);
	if(r>rH) r=rH;
	return r;
}
double Mjz::YZS(int H,int L)
{
	int i,j,ii,jj,ki,kj,Row,Col;
	Row=row-1;
	Col=col-1;
	Mjz M(Row,Col);
	for(i=0;i<row;i++)
	{
		ki=i;
		ii=i;
		if(i>H) ii=i-1;
		if(i>0 && i==H) ki=i-1;
		for(j=0;j<col;j++)
		{
			if(i!=H || j!=L)
			{
				kj=j;
				if(j>0 && j==L) kj=j-1;
				if(j<L) jj=j;
				if(j>L) jj=j-1;
				M.SetAt(ii*Row+jj,GetAt(ki*row+kj));
			}
		}
	}
		//printf("\nM[%i][%i]=",H,L);
		//M.print();
	//return 0;
	return M.HLS();
}
const Mjz& Mjz::T()
{
	int i,j;
	Mjz M(col,row);

	for(i=0;i<col;i++)
		for(j=0;j<row;j++)
		{
			M.SetAt(i*col+j,GetAt(j*row+i));
		}
	return M;
}
const Mjz& Mjz::N()
{
	
	int i,j;
	//Mjz MN(row,col);
	Mjz MN(row);
	double D=HLS();
	if(row!=col || D!=0) return MN=0.0;
	for(i=0;i<row;i++)
		for(j=0;j<col;j++)
		{
			MN.SetAt(i*row+j,Aij(j,i));
		}
		//MN.print();
	return MN/D;
}
const Mjz& Mjz::NZ()
{
	double x=1.0,A,B;
	int i,j,k,a,b;
	Mjz M(row);
	Mjz E(row);
	if(row!=col) return M=0.0;
	E.I(row);
	E.print();
	
	M=*this;
	/////Up sjx
	for(i=0;i<row;i++)
	  {
		a=i*row+i;//M[i,i]
		A=M.GetAt(a);
		if(A==0.0)
		{
			k=i;
			while((A==0.0) && (k++<row))
			{
				a=k*row+i;//M[i,i]
				A=M.GetAt(a);
			}
			if((A==0.0) && (k>=row)) return E=0.0;
			M.RowJH(i,k);
			E.RowJH(i,k);
			//x*=-1.0;
		}
		///M[i][i]=1
		M.RowDiv(i,A);
		E.RowDiv(i,A);
		for(j=i+1;j<row;j++)
		{
			b=j*row+i;
			B=M.GetAt(b);
			if(B!=0.0)
			{
			//B/=A;
			M.RowSub(j,i,B);
			E.RowSub(j,i,B);
			}
		}
	//x*=A;
	}
	/////////M[i][i]=1
	for(i=row-1;i>=0;i--)//col
	  {
		a=i*row+i;//M[i,i]
		A=M.GetAt(a);
		///M[i][i]=1
		for(j=i-1;j>=0;j--)//row
		{
			b=j*row+i;
			B=M.GetAt(b);
			if(B!=0.0)
			{
			//B/=A;
			M.RowSub(j,i,B);
			E.RowSub(j,i,B);
			}
		}
	//x*=A;
	}
	//////////////////
	return E;
}
///////////////////////////////Gol


/////////////////////////////////
int main(int argc, char* argv[])
{
	printf("Hello World!\n");
	Mjz jz;
	jz.ReadTextFile("c:\\jz.txt");
	jz.print("I");
	printf("HLS=%lf",jz.HLS());
	jz.print("I");
	return 0;
}