#
# Fraction-free Gaussian elimination routines
# FractionFreeElimination(C,M,A,m,n): fraction free Gaussian elimination
# of the m by n matrix A over C
#
# Author: MBM 1989
# 
FractionFreeElimination := proc(C,M,A,m,n) local d,i,j,k,s,t,size;

    if hasCategory(C,EuclideanDomain) then
	size := C[EuclideanNorm] else size := length fi;

    s := 1; d := C[1];
    for k to m-1 do

	userinfo(1,GaussianElimination,`elimination at row`,k);
	for i from k to m while A[i,k] = C[0] do od;
	if i > m then RETURN(C[0]) fi;

       	#  Select the smallest element as the pivot.
	for j from i+1 to m do
	    if A[j,k] = C[0] then next fi;
	    if size(A[j,k]) < size(A[i,k]) then i := j fi
	od;

	userinfo(2,GaussianElimination,`pivot is`, A[i,k]);
	#  Pivot is A[i,k], interchange if necessary
	if i <> k then
	    s := - s;
	    for j from k to n do t := A[k,j]; A[k,j] := A[i,j]; A[i,j] := t od
	fi;

	#  Fraction-free elimination
	for i from k+1 to m do
    	    for j from k+1 to n do 
		t := C[`-`]( C[`*`](A[i,j],A[k,k]), C[`*`](A[k,j],A[i,k]) );
		A[i,j] := C[Div](t, d)
	    od
       	od;
       	d := A[k,k]

    od;
    s

end:

FractionFreeBackElimination := proc(C,M,A,m,n) local d,i,j,k;

    d := A[m,m];
    for k from m+1 to n do
        for i from m-1 by -1 to 1 do
	    A[i,k] := C[`*`](d,A[i,k]);
	    for j from i+1 to m do
		A[i,k] := C[`-`](A[i,k],C[`*`](A[i,j]*A[j,k]))
	    od;
	    A[i,k] := C[Div](A[i,k],A[i,i])
	od
    od;

end:

IntegralDomainDet := proc(M) local n,s,A,C;

    C := M[CoefficientRing];
    if not hasCategory(C,IntegralDomain) then
	ERROR(`1st argument must be a Matrix over an IntegralDomain`) fi;

    A := M[ArrayCoeffs](args[2]); n := M[Rows](A);
    s := FractionFreeElimination(C,M,A,n,n);
    C[`*`](s,A[n,n]);

end:

IntegralDomainSolve := proc(C,M,Q,a,b) local d,i,j,k,n,A;

    if not hasCategory(C,IntegralDomain) then
	ERROR(`1st argument must be an IntegralDomain`,C) fi;

    n := M[Rows](a);
    A := array(1..n,1..n+1);
    for i to n do for j to n do A[i,j] := a[i][j] od; A[i,n+1] := b[i] od;
    d := FractionFreeElimination(C,M,A,n,n+1);
    if d=0 then RETURN(FAIL) else d := A[n,n] fi;
    FractionFreeBackSubstitution(C,M,A,n,n+1);
    [seq( Q[Slash](A[k,n+1],d), k=1..n )]

end:

IntegralDomainAdj := proc(M,a) local i,j,k,d,A,C;

    C := M[CoefficientRing];
    if not hasCategory(C,IntegralDomain) then
	ERROR(`1st argument must be a Matrix over an IntegralDomain`) fi;

    n := M[Rows](a);
    A := array(1..n,1..2*n);
    for i to n do
	for j to n do A[i,j] := a[i][j]; A[i,n+j] := C[0] od;
	A[i,n+i] := C[1]
    od;
    d := FractionFreeElimination(C,M,A,n,2*n);
    if d = 0 then RETURN(FAIL) else d := A[n,n] fi;
    FractionFreeBackElimination(C,M,A,n,2*n);
    [[seq( [seq(A[i,j], j=n+1..2*n)], i=1..n)], d ]

end:

IntegralDomainInv := proc(C,M,Q,a) local A,d,i,I,j,n;

    n := M[Dimension];
    A := IntegralDomainAdj(M,a,d);
    I := M[ArrayCoeffs](A,n,n);
    for i to n do for j to n do I[i,j] := Q[Slash](I[i,j],d) od od;
    M[Matrix](I,n,n)

end:

save `FF.m`;
quit
