macro( Rank = `Matrix/Rank` );

#
#--> Rank(M,A)
#
# Compute the rank of the matrix A in the domain M
#
# Author: MBM 1990
#

Rank := proc(M,A) local m,n,r;
	m := M[Rows](A);
	n := M[Cols](A);
	r := MatrixRankHeuristic1(M,A);
if r = FAIL then r := MatrixRankHeuristic2(M,A) fi;
	if r <> FAIL and r = min(m,n) then RETURN(r) fi;
	ERROR(`not implemented`)
end:

	
MatrixRankHeuristic1 := proc(M,A) local R,m,n,x,f,p,r,c,i,j,t,k;

	R := M[CoefficientRing]:
	if not hasOperation(R,ModularHomomorphism) then RETURN(FAIL) fi;
	m := M[Rows](A);
	n := M[Cols](A);
	x := array(1..m,1..n);
	f := R[ModularHomomorphism]();
	p := f[2]; f := f[1];
	for i to m do
	    for j to n do
		x[i,j] := f(A[i][j]);
		if t = FAIL then RETURN(FAIL) fi;
	    od
	od;

	r := 1; c := 1;
	while r <= m and c <= n do
	    for i from r to n while x[i,c] = 0 do od;
	    if i > n then c := c + 1; next fi;
	    if i <> n then
		for j from c to n do
		    t := x[i,j]; x[i,j] := x[r,j]; x[r,j] := t
		od;
	    fi;
	    i := traperror(modp(1/x[r,c],p));
	    if i = lasterror then RETURN(FAIL) fi;
	    for k from r+1 to m do
		if x[k,c] = 0 then next fi;
		t := modp(x[k,c]*i,p);
	        for j from c+1 to n do x[k,j] := modp(x[k,j]-t*x[r,j],p) od;
	    od;
	    r := r+1; c := c+1
	od;
	min(r,c)-1
end:

MatrixRankHeuristic2 := proc(M,A) local R,m,n,x,F,f,r,c,i,j,t,k;
	R := M[CoefficientRing]:
	if not hasOperation(R,ModularMapping) then RETURN(FAIL) fi;
	m := M[Rows](A);
	n := M[Cols](A);
	x := array(1..m,1..n);
	F := R[ModularMapping]();
	f := F[1];
	F := F[2];
	for i to m do
	    for j to n do
		x[i,j] := f(A[i][j]);
	    od
	od;

	r := 1; c := 1;
	while r <= m and c <= n do
	    for i from r to n while x[i,c] = F[0] do od;
	    if i > n then c := c + 1; next fi;
	    if i <> n then
		for j from c to n do
		    t := x[i,j]; x[i,j] := x[r,j]; x[r,j] := t
		od;
	    fi;
	    i := traperror(F[Inv](x[r,c]));
	    if i = lasterror then RETURN(FAIL) fi;
	    for k from r+1 to m do
		if x[k,c] = F[0] then next fi;
		t := F[`*`](x[k,c],i);
	        for j from c+1 to n do
		    x[k,j] := F[`-`](x[k,j],F[`*`](t,x[r,j]))
		od;
	    od;
	    r := r+1; c := c+1
	od;
	min(r,c)-1
end:

save `Rank.m`;
quit
