#include "stdio.h"
#include "string.h"
#include <complex.h>
#include "d:\book\計算機\用c++編寫數學常用算法\用c++編寫數學常用算法\cmatrix.h"
#include "fstream.h"

// 缺省的構造函數，產生0行0列空復矩陣
cmatrix::cmatrix(cbuffer * b): rownum(0),colnum(0),
	isneg(0),istrans(0),isconj(0)
{
		if(b){
			buf=b;
			buf->alloc(0);
			}
		else
			buf = getnewcbuffer(0); // 如果缺省的b沒有給出，則產生一
				// 長度為行乘列個復數和緩存
}
// 構造一r行c列的復矩陣，可選的復緩存指向b
cmatrix::cmatrix(size_t r, size_t c, cbuffer * b):
		rownum(r),colnum(c),istrans(0),isneg(0),isconj(0)
{
		if(b){
			buf=b;
			buf->alloc(r*c);
			}
		else
			buf = getnewcbuffer(r*c); // 如果缺省的b沒有給出，則產生一
				// 長度為行乘列個復數和緩存
}

cmatrix::cmatrix(size_t n, cbuffer * b):
		rownum(n),colnum(1),istrans(0),isneg(0),isconj(0)
{
		if(b){
			buf=b;
			buf->alloc(n);
			}
		else
			buf = getnewcbuffer(n); // 如果缺省的b沒有給出，則產生一
				// 長度為行乘列個復數和緩存
}

cmatrix cunit(size_t n) // 產生n階單位矩陣
{
	cmatrix m(n,n);
	for(size_t i=0; i<n; i++)
	for(size_t j=0; j<n; j++)
			if(i==j) m.set(i,i,1.0);
			else m.set(i,j,0.0);
	return m;
}

// 拷貝構造函數
cmatrix::cmatrix(cmatrix& m)
{
	rownum = m.rownum; // 復制行數和列數
	colnum = m.colnum;
	istrans = m.istrans; // 復制標志
	isneg = m.isneg;
	isconj = m.isconj;
	buf = m.buf; // 指向相同的緩存
	buf->refnum++;	// 緩存的引用數增1
}

// 用數據文件構造復矩陣
cmatrix::cmatrix(const char * filename, cbuffer * b) :
	istrans(0),isneg(0),isconj(0)
{
	char label[10];
	ifstream in(filename); // 打開文件輸入
	in >> label; // 文件以cmatrix關鍵詞打頭
	if(strcmp(label, "cmatrix")!=0) throw TMESSAGE("format error!");
	in >> rownum >> colnum; // 接著是行數和列數
	if(!in.good()) throw TMESSAGE("read file error!");
	if(b) { buf=b;	// 產生所需長度的緩存
		buf->alloc(rownum*colnum);
	}
	else buf = getnewcbuffer(rownum*colnum);
	size_t line;
	for(size_t i=0; i<rownum; i++) { // 按行循環輸入
		in >> line;	// 每行以行號開始
		if(line != i+1) throw TMESSAGE("format error!");
		in.width(sizeof(label));
		in >> label; // 緊接著是冒號
		if(label[0] != ':') throw TMESSAGE("format error!");
		COMPLEX a;
		for(size_t j=0; j<colnum; j++) { // 隨后是一行的內容
			in >> a;
			set(i,j,a);
		}
		if(!in.good()) throw TMESSAGE("read file error!");
	}
}

// 從內存數據塊構造復矩陣。data指向內存數據塊，構造完畢后data指向的內存可以釋放
cmatrix::cmatrix(void * data, size_t r, size_t c, cbuffer * b):
		rownum(r),colnum(c),istrans(0),isneg(0),isconj(0)
{
	if(b){
		buf=b;
		buf->alloc(r*c);
		}
	else
		buf = getnewcbuffer(r*c);
	COMPLEX * d = (COMPLEX *)data;
	for(size_t i=0; i<r*c; i++)
		buf->set(i,d[i]);
}

cmatrix::cmatrix(matrix& m, cbuffer* b) : // 由實矩陣產生一個內容相同的復矩陣
	istrans(0),isneg(0),isconj(0)
{
	rownum = m.rownum;
	colnum = m.colnum;
	if(b) { buf=b;	// 產生所需長度的緩存
		buf->alloc(rownum*colnum);
	}
	else buf = getnewcbuffer(rownum*colnum);
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		set(i,j,m.value(i,j));
}

cmatrix::cmatrix(matrix& mr, matrix& mi, cbuffer * b): // 兩個實矩陣成為實部和虛部
	istrans(0),isneg(0),isconj(0)
{
	if(mr.rownum != mi.rownum || mr.colnum != mi.colnum)
		throw TMESSAGE("dimension not same!");
	rownum = mi.rownum;
	colnum = mr.colnum;
	if(b) { buf=b;	// 產生所需長度的緩存
		buf->alloc(rownum*colnum);
	}
	else buf = getnewcbuffer(rownum*colnum);
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		set(i,j,COMPLEX(mr.value(i,j),mi.value(i,j)));
}

matrix cmatrix::real() // 由實部生成實矩陣
{
	matrix m(rownum, colnum);
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++) {
		m.set(i,j,::real(value(i,j)));
	}
	return m;
}

matrix cmatrix::image() // 由虛部生成實矩陣
{
	matrix m(rownum, colnum);
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		m.set(i,j,imag(value(i,j)));
	return m;
}

// 矩陣賦值
cmatrix& cmatrix::operator=(cmatrix& m)
{
	rownum = m.rownum; // 復制行數和列數
	colnum = m.colnum;
	if(buf == m.buf) // 如果本來就指向相同的緩存，則返回
		return (*this);
	buf->refnum--;	// 否則放棄原來的緩存，先將原緩存引用數減1
						// 如果不為0說明還有其它的矩陣引用數據
	if(!buf->refnum)delete buf; // 否則刪除緩存
	istrans = m.istrans;	// 復制標志
	isneg = m.isneg;
	isconj = m.isconj;
	buf = m.buf; // 復制緩存指針
	buf->refnum++; // 緩存引用數增1
	return (*this);
}

cmatrix& cmatrix::operator=(matrix& m)
{
	rownum = m.rownum; // 復制行數和列數
	colnum = m.colnum;
	buf->refnum--;
	if(!buf->refnum) delete buf; // 放棄原來緩存
	istrans = isneg = isconj = 0; // 標志清零
	buf = getnewcbuffer(rownum*colnum); // 申請新的緩存
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		set(i,j,m.value(i,j));
	return (*this);
}

cmatrix& cmatrix::operator=(COMPLEX a) // 將矩陣所有元素設為常數a
{
	if(rownum == 0 || colnum==0 ) return (*this);
	istrans = isneg = isconj = 0; // 標志清零
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		set(i,j,a);
	return (*this);
}

COMPLEX cmatrix::operator()(size_t r, size_t c)
 // 重載函數運算符，返回第r行c列的值
{
	if(r>= rownum || c >= colnum)
		throw TMESSAGE("range overflow");
	return value(r,c);
}

COMPLEX cmatrix::operator()(size_t r)
	// 重載函數運算符，返回第r行0列的值，用于向量
{
	return (*this)(r,0);
}


// 流輸出
ostream& operator<<(ostream& o, cmatrix& m)
{	// 先輸出關鍵詞cmatrix,后跟行數和列數
	o << "cmatrix " << m.rownum << '\t' << m.colnum << endl;
	// 每行由行號，冒號，后跟一行的復數構成
	for(size_t i=0; i<m.rownum; i++) {
		o<< (i+1) <<':';
		size_t k=8;
		for(size_t j=0; j<m.colnum; j++) {
			o<<'\t'<<m.value(i,j);
			if(--k==0) {
				k=8;
				o<<endl;
			}
		}
		o<<endl;
	}
	return o;
}

// 流輸入
istream& operator>>(istream& in, cmatrix& m)
{
	char label[10];
	in.width(sizeof(label));
	in >> label;
	if(strcmp(label, "cmatrix")!=0) // 以關鍵詞cmatrix打到
		throw TMESSAGE("format error!");
	in >> m.rownum >> m.colnum; // 后跟行數和列數
	if(!in.good()) throw TMESSAGE("read file error!");
	m.buf->refnum--; // 放棄原來的緩存
	if(!m.buf->refnum) delete m.buf;
	m.buf = getnewcbuffer(m.rownum*m.colnum); // 申請新的緩存
	m.istrans = m.isneg = m.isconj = 0; // 標志清零
	size_t line;
	for(size_t i=0; i<m.rownum; i++) {
		in >> line;
		if(line != i+1) throw TMESSAGE("format error!"); // 每行以行號開始
		in.width(sizeof(label));
		in >> label;
		if(label[0] != ':') throw TMESSAGE("format error!"); // 后跟冒號
		COMPLEX a;
		for(size_t j=0; j<m.colnum; j++) { // 然后是這一行的數據
			in >> a;
			m.set(i,j,a);
		}
		if(!in.good()) throw TMESSAGE("read file error!");
	}
	return in;
}

// 數乘矩陣，原矩陣內容改變為結果
cmatrix& cmatrix::operator*=(COMPLEX a)
{
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		set(i,j,a*value(i,j));
	return (*this);
}

// 數乘矩陣，原矩陣內容不變，返回結果矩陣
cmatrix cmatrix::operator*(COMPLEX a)
{
	cmatrix m(rownum, colnum);
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		m.set(i,j,a*value(i,j));
	return m;
}

cmatrix cmatrix::operator+(COMPLEX a) // 矩陣加常數，產生新矩陣
{
	cmatrix m(rownum, colnum);
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		m.set(i,j,a+value(i,j));
	return m;

}

cmatrix& cmatrix::operator+=(COMPLEX a) // 矩陣加常數，更動原矩陣
{
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		set(i,j,a+value(i,j));
	return (*this);
}

// 矩陣相加
cmatrix cmatrix::operator+(cmatrix& m)
{
	if(rownum != m.rownum || colnum != m.colnum) // 行列必須對應相等
		throw TMESSAGE("can not do add of cmatrix");
	cmatrix mm(rownum, colnum);
	COMPLEX a;
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
	{
		a = value(i,j)+m.value(i,j);
		mm.set(i,j,a);
	}
	return mm;
}

// 矩陣相加，不產生新的矩陣，修改原矩陣的內容
cmatrix& cmatrix::operator+=(cmatrix &m)
{
	COMPLEX a;
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
	{
		a = value(i,j)+m.value(i,j);
		set(i,j,a);
	}
	return (*this);
}

cmatrix cmatrix::operator-(cmatrix& a)	// 矩陣相減
{
	return (*this)+(-a);
}

cmatrix& cmatrix::operator-=(cmatrix& a) // 矩陣自身減矩陣
{
	(*this)+=(-a);
	return (*this);
}

// 矩陣相乘，返回新產生的結果矩陣
cmatrix cmatrix::operator*(cmatrix& m)
{
	if(colnum != m.rownum) // 前矩陣的列數必須等于后矩陣的行數
		throw TMESSAGE("can not multiply!");
	cmatrix mm(rownum,m.colnum);
	COMPLEX a;
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<m.colnum; j++){
		a = 0.0;
		for(size_t k=0; k<colnum; k++)
			a += value(i,k)*m.value(k,j);
		mm.set(i,j,a);
	}
	return mm;
}

// 矩陣相乘，改變原矩陣的內容
cmatrix& cmatrix::operator*=(cmatrix& m)
{
	(*this) = (*this)*m;
	return (*this);
}

cmatrix cmatrix::operator-() // 矩陣求負，產生新的矩陣
{
	cmatrix m(*this);
	m.neg();
	return m;
}