cmatrix& cmatrix::neg() // 矩陣自身求負
{
	isneg = !isneg;	// 取負標志求反
	return (*this);
}

cmatrix cmatrix::t()		// 矩陣轉置，產生新的矩陣
{
	cmatrix m(*this);
	m.trans();
	return m;
}

cmatrix cmatrix::operator!()		// 矩陣共軛，產生新的矩陣
{
	cmatrix m(*this);
	return m.conj();
}

cmatrix cmatrix::tc()		// 矩陣共軛轉置，產生新的矩陣
{
	cmatrix m(*this);
	return m.transconj();
};

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

cmatrix cmatrix::operator-(COMPLEX a)	// 矩陣減常數，產生新的矩陣
{
	return (*this)+(-a);
}

cmatrix operator-(COMPLEX a, cmatrix m) // 常數減矩陣，產生新的矩陣
{
	return (-m)+a;
}

void cmatrix::swapr(size_t r1, size_t r2, size_t k) // 交換矩陣r1和r2兩行
{
	COMPLEX a;
	for(size_t i=k; i<colnum; i++) {
		a = value(r1, i);
		set(r1, i, value(r2, i));
		set(r2, i, a);
	}
}

void cmatrix::swapc(size_t c1, size_t c2, size_t k) // 交換c1和c2兩列
{
	COMPLEX a;
	for(size_t i=k; i<colnum; i++) {
		a = value(i, c1);
		set(i, c1, value(i, c2));
		set(i, c2, a);
	}
}

DOUBLE cmatrix::maxabs(size_t &r, size_t &c, size_t k)
 // 求第k行和第k列后的主元及位置
{
	DOUBLE a=0.0;
	DOUBLE br,bi;
	for(size_t i=k;i<rownum;i++)
	for(size_t j=k;j<colnum;j++) {
		br = ::real(value(i,j));
		bi = ::imag(value(i,j));
		br = br*br+bi*bi;
		if(a < br) {
			r=i;c=j;a=br;
		}
	}
	return a;
}

size_t cmatrix::zgsxy(cmatrix & m, int fn) // 進行主高斯消元運算，fn為參數，缺省為0
 /* 本矩陣其實是常數陣，而矩陣m必須是方陣
	運算過程其實是對本矩陣和m同時作行初等變換，
	運算結果m的對角線相乘將是行列式，而本矩陣則變換成
	自己的原矩陣被m的逆陣左乘，m的秩被返回，如果秩等于階數
	則本矩陣中的內容已經是唯一解
 */
{
	if(rownum != m.rownum || m.rownum != m.colnum) // 本矩陣行數必須與m相等
						// 且m必須是方陣
		throw TMESSAGE("can not divide!");
	lbuffer * bb = getnewlbuffer(rownum); // 產生一維數為行數的長整數緩存區
	lbuffer & b = (*bb); // 用引用的辦法使下面的程序容易懂
	size_t is;
	COMPLEX a;
	size_t i,j,rank=0;
	for(size_t k=0; k<rownum; k++) { // 從第0行到第k行進行主高斯消元
		if(m.maxabs(is, i, k)==0) // 求m中第k級主元，主元所在的行，列存在is,i中
				// 返回的主元值放在a中
			break; // 如果主元為零，則m不可逆，運算失敗
		rank = k+1; // rank存放當前的階數
		b.retrieve(k) = i;  // 將長整數緩存區的第k個值設為i
		if(i != k)
			m.swapc(k, i); // 交換m中i,k兩列
		if(is != k) {
			m.swapr(k, is, k); // 交換m中i,k兩行,從k列以后交換
			swapr(k, is); // 交換本矩陣中i,k兩行
		}
		a = m.value(k,k); // 取出主元值
		for (j=k+1;j<rownum;j++) // 本意是將m的第k行除以主元
				// 但只需把第k行的第k+1列以上除以主元即可
				// 這樣還保留了主元作行列式運算用
			m.set(k,j,m.value(k,j)/a);
		for (j=0;j<colnum;j++) // 將本矩陣的第k行除以主元
			set(k,j,value(k,j)/a);
		// 上面兩步相當于將m和本矩陣構成的增廣矩陣第k行除以主元
		// 下面對增廣矩陣作行基本初等變換使第k行的其余列均為零
		// 但0值無必要計算，因此從第k+1列開始計算
		for(j=k+1;j<rownum;j++) // j代表列，本矩陣的行數就是m的列數
		for(i=0;i<rownum;i++) //i代表行，依次對各行計算，k行除外
			if(i!=k)
				m.set(i,j,m.value(i,j)-m.value(i,k)*m.value(k,j));
		// 對本矩陣亦作同樣的計算
		for(j=0;j<colnum;j++)
		for(i=0;i<rownum;i++)
			if(i!=k)
				set(i,j,value(i,j)-m.value(i,k)*value(k,j));
	} // 主高斯消元循環k結束
	if(fn == 1) {
		for(j=0; j<rank; j++)
		for(i=0; i<rownum; i++)
			if(i==j) m.set(i,i,1.0);
			else
				m.set(i,j,0.0);
		for(k = rank; k>0; k--)
			m.swapc(k-1,(size_t)b[k-1]);
	}
	for(k = rank; k>0; k--) // 將本矩陣中的各行按b中內容進行調節
		if(b[k-1] != k-1)
			swapr(k-1,(size_t)b[k-1]); // 行交換
	delete bb; // 釋放長整數緩存
	return rank; // 返回mm的秩
}

cmatrix& cmatrix::operator/=(cmatrix m) // 利用重載的除法符號/=來解方程
	// 本矩陣設為d,則方程為mx=d,考慮解寫成x=d/m的形式，
	// 而方程的解也存放在d中，則實際編程時寫d/=m
{
	if(zgsxy(m)!=rownum) // 如秩不等于m的階數，則方程無解
		throw TMESSAGE("can not divide!");
	return *this;
}

cmatrix& cmatrix::fft(int l)
{
	if(colnum != 1) throw TMESSAGE("fft need vector!");
	size_t n,i;
	int k;
	n = 1;
	for(k=0;k<10000;k++) {
		if(n>=rownum) break;
		n*=2;
	}
	if(rownum < n) {	// 將數組擴大為2的乘冪維
		cmatrix m(n);
		for(i=0; i<rownum; i++)
			m.set(i,value(i,0));
		for(i=rownum; i<n; i++)
			m.set(i,0.0);
		(*this)=m;
	}
	size_t it,m,is,j,nv;
	int l0;
	COMPLEX q,s,v,podd;
	DOUBLE p;
	cmatrix f(n);
	for (it=0; it<n; it++)
		{ m=it; is=0;
		  for (i=0; i<k; i++)
			 { j=m/2; is=2*is+(m-2*j); m=j;}
		  f.set(it,value(is,0));
		}
	set(0,1.0);
	cbuffer & pp = *buf;
	cbuffer & ff = *(f.buf);
	p = 2.0*M_PI/n;
	pp[1] = COMPLEX(cos(p),-sin(p));
	if (l!=0) pp[1] = ::conj(pp[1]);
	for (i=2; i<n; i++)
		pp[i] = pp[1]*pp[i-1];
	for (it=0; it<n-1; it+=2)
		{  v = ff[it];
			ff[it] += ff[it+1];
			ff[it+1] = v - ff[it+1];
		}
	m=n/2; nv=2;
	for (l0=k-2; l0>=0; l0--)
		{ m/=2; nv*=2;
		  for (it=0; it<=(m-1)*nv; it+=nv)
			 for (j=0; j<=(nv/2)-1; j++)
				{
					podd = pp[m*j]*ff[it+j+nv/2];
					ff[it+j+nv/2] = ff[it+j]-podd;
					ff[it+j] += podd;
				}
		}
	 if (l!=0)
		for (i=0; i<n; i++)
		  ff[i]/=n;
	(*this) = f;
	return (*this);
}

cmatrix fft(matrix& m)
{
	cmatrix mm(m);
	return mm.fft();
}

matrix convolute(matrix&a, matrix& b)	// 求矩陣a與b的離散卷積
{
	if(a.colnum != 1 || b.colnum != 1)throw TMESSAGE("must be a vector");
	size_t nn = a.rownum+b.rownum;
	size_t n;
	int k;
	n = 1;
	for(k=0;k<10000;k++) {
		if(n>=nn) break;
		n*=2;
	}
	matrix aa(n),bb(n);
	size_t i;
	aa = 0.0;
	bb = 0.0;
	for(i=0; i<a.rownum; i++)
		aa.set(i,a(i));
	for(i=0; i<b.rownum; i++)
		bb.set(i,b(i));
	cmatrix aafft(fft(aa)),bbfft(fft(bb));
	for(i=0; i<n; i++)
		aafft.set(i,aafft(i)*bbfft(i));
	nn--;
	matrix c(nn),d;
	d = (aafft.fft(1)).real();
	for(i=0;i<nn;i++)
		c.set(i,0,d(i));
	return c;
}

cmatrix& cmatrix::inv()		// 用全選主元高斯-約當法求逆矩陣
{
	if(rownum != colnum) throw TMESSAGE("must be a square matrix");
	size_t n = rownum;
	lbuffer * isp = getnewlbuffer(n);
	lbuffer * jsp = getnewlbuffer(n);
	lbuffer& is = (*isp);
	lbuffer& js = (*jsp);
	size_t i,j,k,u,v;
	DOUBLE d;
	for (k=0; k<n; k++)
		{
			if((d=maxabs(u,v,k))== 0.0) {
				delete isp;
				delete jsp;
				throw TMESSAGE("can not inverse");
			}
			is[k]=u;
			js[k]=v;
/*
		 d=0.0;
		  for (i=k; i<=n-1; i++)
		  for (j=k; j<=n-1; j++)
			 { u=i*n+j;
				p=ar[u]*ar[u]+ai[u]*ai[u];
				if (p>d) { d=p; is[k]=i; js[k]=j;}
			 }
		  if (d+1.0==1.0)
			 { free(is); free(js); printf("err**not inv\n");
				return(0);
			 }
			 */
		  if (is[k]!=k) swapr(k,(size_t)is[k]);
/*
			 for (j=0; j<=n-1; j++)
				{ u=k*n+j; v=is[k]*n+j;
				  t=ar[u]; ar[u]=ar[v]; ar[v]=t;
				  t=ai[u]; ai[u]=ai[v]; ai[v]=t;
				} */
		  if (js[k]!=k) swapc(k,(size_t)js[k]);
/*			 for (i=0; i<=n-1; i++)
				{ u=i*n+k; v=i*n+js[k];
				  t=ar[u]; ar[u]=ar[v]; ar[v]=t;
				  t=ai[u]; ai[u]=ai[v]; ai[v]=t;
				} */
//		  l=k*n+k;
		  set(k,k,::conj(value(k,k))/d);
//		  ar[l]=ar[l]/d; ai[l]=-ai[l]/d;
		  for (j=0; j<n; j++)
			 if (j!=k)
				set(k,j,value(k,j)*value(k,k));
/*
				{ u=k*n+j;
				  p=ar[u]*ar[l]; q=ai[u]*ai[l];
				  s=(ar[u]+ai[u])*(ar[l]+ai[l]);
				  ar[u]=p-q; ai[u]=s-p-q;
				} */
		  for (i=0; i<n; i++)
			 if (i!=k)
				{ // v=i*n+k;
				  for (j=0; j<n; j++)
					 if (j!=k)
						set(i,j,value(i,j)-value(k,j)*value(i,k));
/*
						{ u=k*n+j;  w=i*n+j;
						  p=ar[u]*ar[v]; q=ai[u]*ai[v];
						  s=(ar[u]+ai[u])*(ar[v]+ai[v]);
						  t=p-q; b=s-p-q;
						  ar[w]=ar[w]-t;
						  ai[w]=ai[w]-b;
						} */
				}
		  for (i=0; i<n; i++)
			 if (i!=k)
				set(i,k,-value(i,k)*value(k,k));
/*				{ u=i*n+k;
				  p=ar[u]*ar[l]; q=ai[u]*ai[l];
				  s=(ar[u]+ai[u])*(ar[l]+ai[l]);
				  ar[u]=q-p; ai[u]=p+q-s;
				} */
		}
	 for (k=n; k>0; k--)
		{ if (js[k-1]!=k-1) swapr((size_t)js[k-1],k-1);
/*			 for (j=0; j<n; j++) swapr((size_t)js[k-1],k-1);
				{ u=k*n+j; v=js[k]*n+j;
				  t=ar[u]; ar[u]=ar[v]; ar[v]=t;
				  t=ai[u]; ai[u]=ai[v]; ai[v]=t;
				} */
		  if (is[k-1]!=k-1) swapc((size_t)is[k-1],k-1);
/*			 for (i=0; i<=n-1; i++)
				{ u=i*n+k; v=i*n+is[k];
				  t=ar[u]; ar[u]=ar[v]; ar[v]=t;
				  t=ai[u]; ai[u]=ai[v]; ai[v]=t;
				} */
		}
	delete isp;
	delete jsp;
	return (*this);
}

cmatrix cmatrix::operator~()	// 求逆矩陣，但產生新矩陣
{
	cmatrix m(*this);
	m.inv();
	return m;
}

cmatrix operator/(COMPLEX a, cmatrix m) // 求逆矩陣再乘常數
{
	return (~m)*a;
}

cmatrix& cmatrix::operator/=(COMPLEX a) // 所有元素乘a的倒數，自身改變
{
	(*this)*=1.0/a;
	return (*this);
}

cmatrix cmatrix::operator/(COMPLEX a) // 所有元素乘a的倒數，產生新的矩陣
{
	cmatrix m(*this);
	m/=a;
	return m;
}