#
# Minor Expansion and related routines
#
# Author: MBM 1989
#

MinorExpansion := proc(C,M,A) local d,m,n;
    m := M[Rows](A); n := M[Cols](A);
    if m <> n then ERROR(`matrix must be square`) fi;
    d := MinorExpansion1(C,M,A,[$1..n],[$1..n],n);
    MinorExpansion1 := subsop(4=NULL,op(MinorExpansion1));
    d
end:

MinorExpansion1 := proc(C,M,A,r,c,n) local i,s,t,d,m; option remember;
    if n = 1 then A[r[1]][c[1]]
    elif n = 2 then C[`-`]( C[`*`](A[r[1]][c[1]],A[r[2]][c[2]]),
	        	C[`*`](A[r[1]][c[2]],A[r[2]][c[1]]) )
    else 
        d := C[0]; 
        s := -1;
        for i to n do
	    s := -s;
	    if A[r[i]][c[1]] <> C[0] then
	        t := subsop(i=NULL,r), subsop(1=NULL,c);
	        m := MinorExpansion1(C,M,A,t,n-1);
	        t := C[`*`](A[r[i]][c[1]],m);
	        if s = 1 then d := C[`+`](d,t) else d := C[`-`](d,t) fi
	    fi
        od;
        d
    fi
end:

AdjointExpansion := proc(C,M,A)
local i,j,m,n,r,s,t,B,d,sym;
    m := M[Rows](A); n := M[Cols](A);
    if m <> n then ERROR(`matrix must be square`) fi;
    sym := M[symmetric](A);
    if sym then B := array(1..n,1..n,symmetric) else B := array(1..n,1..n) fi;
    r := [$1..n];
    for i to n do
        s := subsop(i=NULL,r);
        if sym then i else 1 fi;
	for j from " to n do
	    t := subsop(j=NULL,r);
	    B[j,i] := C[`*`]( (-1)^(i+j), MinorExpansion1(C,M,A,s,t,n-1) )
    	od
    od;
    d := MinorExpansion1(C,M,A,r,r,n);
    MinorExpansion1 := subsop(4=NULL,op(MinorExpansion1));
    [M[Matrix](B,n,n),d];
end:

CramersRule := proc(C,M,A)
local sym,a,n,d,B,i,j;
    sym := M[symmetric](A);
    a := AdjointExpansion(C,M,A);
    n := M[Rows](A);
    d := a[2];
    if d = C[0] then ERROR(`matrix is singular`) fi;
    d := C[Inv](d);
    if d = FAIL then RETURN( FAIL )  fi;
    if sym then B := array(1..n,1..n,symmetric,a[1])
	else B := array(1..n,1..n,a[1]) fi;
    for i to n do
	if sym then i else 1 fi;
	for j from " to n do B[i,j] := C[`*`](B[i,j],d) od
    od;
    M[Matrix](B,n,n);
end:

save `ME.m`;
quit
