#
# SquareMatrix(n:posint,R:Domain)
#
# A Square Matrix domain -- Create the ring of n by n square matrices over R
#
# Author: MBM 1990 - 1992
#

macro(
	IntegralDomainDet = readlib(IntegralDomainDet, ``.gauss.`/SM/FF.m`),
	IntegralDomainAdj = readlib(IntegralDomainAdj, ``.gauss.`/SM/FF.m`),
	MatrixInverse = readlib(MatrixInverse, ``.gauss.`/SM/GE.m`),
	MatrixSolve = readlib(MatrixSolve, ``.gauss.`/SM/GE.m`),
	MatrixDeterminant = readlib(MatrixDeterminant, ``.gauss.`/SM/GE.m`),
	MinorExpansion = readlib(MinorExpansion, ``.gauss.`/SM/ME.m`),
	Cramers = readlib(CramersRule, ``.gauss.`/SM/ME.m`),
	AdjointExpansion = readlib(AdjointExpansion, ``.gauss.`/SM/ME.m`) );

SquareMatrix := proc(n) local R,S,env,Rep; option remember;

	if not type(n,integer) or n < 1 then
		ERROR(`dimension must be a positive integer`) fi;
	R := args[2];
	if not hasCategory(R,Ring) then
		ERROR(`2nd argument must be a Ring`) fi;

	S := Ring();
	S[CoefficientRing] := R;
	S[Characteristic] := R[Characteristic];
	S[DomainName] := SquareMatrix;
	defOperation( Symmetric, S &-> Boolean, S );
	defOperation( Matrix, List(List(R)) &-> S, S );
	defOperation( ArrayCoeffs, S &-> Array(R), S );
	defOperations( {Rows,Cols}, S &-> Integer, S );
	defOperations( {Row,Col}, [S,Integer] &-> List(R), S );
	defOperation( List, S &-> List(List(R)), S );
	defOperation( `*V`, [S,Vector(R)] &-> Vector(R), S ); 
	defOperation( `V*`, [Vector(R),S] &-> Vector(R), S ); 
	defOperation( Transpose, S &-> S, S );

	Rep := Matrix(R);
	env := ['dim' = n, 'M' = S, 'C' = R, 'T' = Rep];

	S[Random] := subs(env, proc() T[RandMatrix](dim,dim) end);
	S[Output] := Rep[Output];
	S[Symmetric] := Rep[Symmetric];
	S[Input] := subs(env, proc(a) local m;
		m := T['Input'](a);
		if m = FAIL or T[Rows](m) <> dim or T[Cols](m) <> dim
			then FAIL else m fi;
	    end);

	# Arithmetic operations
	S[0] := Rep[New](n,n,R[0]);
	S[1] := Rep[Diagonal]([R[1] $ n]);
	S[`=`] := subs(env, proc(x,y) local i,j;
		for i to dim do
		    for j to dim do
			if not C[`=`](x[i][j],y[i][j]) then RETURN(false) fi
		    od
		od;
		true
	    end);
	S[`+`] := subs(env, proc(x) local r,y;
		if nargs = 0 then RETURN(M[0]) else r := x fi;
		for y in [args[2..nargs]] do r := T[`+`](r,y) od;
		r
	    end);
	S[`*`] := subs(env, proc(x) local r,y;
		if nargs = 0 then RETURN(M[1]) fi;
		if type(x,integer) then
		    r := C[Coerce](x);
		    r := T[Diagonal]([r$dim]);
		else r := x
		fi;
		for y in [args[2..nargs]] do r := T[`*`](r,y) od;
		r
	    end);

	S[`.`] := Rep[`.`];
	S[`*V`] := Rep[`*V`];
	S[`V*`] := Rep[`V*`];
	S[Matrix] := Rep[Matrix];
	S[ArrayCoeffs] := Rep[ArrayCoeffs];
	S[Row] := Rep[Row];
	S[Rows] := subs(env, proc() dim end);
	S[Col] := Rep[Col];
	S[Cols] := subs(env, proc() dim end);
	S[Transpose] := Rep[Transpose];
	
	# Linear Algebra Operations
	if hasCategory(R,CommutativeRing) then
		defOperation(Det, S &-> R, S );
		defOperation(Adj, S &-> S, S );
		defOperation(Rank, S &-> Integer, S );
		defOperation(Singular, S &-> Boolean, S );
		S[Adj] := subs(env, proc(a) op(1,AdjointExpansion(C,M,a)) end);
		S[Det] := subs(env, proc(a) MinorExpansion(C,M,a) end);
		S[Rank] := Rep[Rank];
		S[Singular] := subs(env, proc(a) evalb(M[Rank](a) <> dim) end);
	fi;
	if hasCategory(R,IntegralDomain) then
		#S[Adj] := subs(env, proc(a) op(1,IntegralDomainAdj(M,a)) end);
		S[Det] := subs(env, proc(a) IntegralDomainDet(M,a) end);
	fi;
	if hasCategory(R,EuclideanDomain) then
		defOperations( {Hermite,Smith}, S &-> S, S );
		S[Hermite] := Rep[Hermite];
	fi;
	if hasCategory(R,Field) then
		defOperation( CramersRule, S &-> S, S );
		defOperation( Solve, [S,List(R)] &-> Union(List(R),FAIL), S );
		defOperation( NullSpace, S &-> FiniteSet(List(R)), S );
		S[Inv] := subs(env, proc(a) MatrixInverse(C,M,a) end);
		S[Det] := subs(env, proc(a) MatrixDeterminant(C,M,a) end);
		S[Solve] := subs(env, proc(a,b) MatrixSolve(C,M,a,b) end);
		S[CramersRule] := subs(env, proc(a) Cramers(C,M,a) end);
	fi;
	op(S)
end:

save `SM.m`:
quit
