CCodeDir:="./"; 

minor := proc(A::Matrix,P::prime)
# Input: an m x n matrix A of polynomials over Z and (optionally) a prime P.
# Output: (r::posint,R::list(posint),C::list(posint)) where
#  r = rank(A) and R and C are lists of row and column indices, respectively,
#  such that the r x r submatrix A[R,C] has rank r.
# Do this by evaluating all variabels at random values mod P

local n,m,p,U,X,S,B,i,j,k,row,col,inv,mu,rank,c,piv,x,NZ,ij;

    n := LinearAlgebra:-RowDimension(A);
    m := LinearAlgebra:-ColumnDimension(A);
    if nargs=1 then p := 2^62-57; else p := P fi;
    U := rand(p);
    X := indets(A);
    S := {seq( x=U(), x in X )};
    B := Matrix(n,m);
    #for i to n do for j to m do B[i,j] := Eval(A[i,j],S) mod p; od od;
    for ij in indices(A) do B[op(ij)] := Eval(A[op(ij)],S) mod p; od;
printf("minor: starting elimination\n");
    for i to n do piv[i] := i; od;
    rank := 0;
    c := 0;
    k := 1;
    while k <= n and c < m do
        c := c+1; # next column
        for i from k to n while B[piv[i],c] = 0 do od; # select non-zero pivot
        if i>n then next; fi; # move to next column
        piv[k],piv[i] := piv[i],piv[k];
        rank := rank+1; # increase rank
        row[rank] := piv[k];
        col[rank] := c;
        inv := 1/B[piv[k],c] mod p;
NZ := [seq( `if`( B[piv[k],j]=0, NULL, j ), j=c+1..m )]; # improve for sparse case
        for i from k+1 to n do
            if B[piv[i],c] = 0 then next fi;
            mu := B[piv[i],c]*inv mod p;
            #for j from c+1 to m do
            for j in NZ do
                B[piv[i],j] := B[piv[i],j] - mu*B[piv[k],j] mod p;
            od;
            B[piv[i],c] := 0;
        od;
        k := k+1; # next row
     od;
     return rank, [seq( row[j], j=1..rank )], [seq( col[j], j=1..rank )];
end:


Cminor := define_external('minor64s',
  AA::ARRAY(1..nn,1..mm,integer[8],order=C_order), nn::integer[8], mm::integer[8],
  ROWS::ARRAY(1..nn,datatype=integer[4]), COLS::ARRAY(1..mm,datatype=integer[4]),
  pp::integer[8], 
  LIB=cat(CCodeDir,"minor.so"), 
  RETURN::integer[4]):

#minor64( A, n, m, R, C, p ) -> r

minor := proc(A::Matrix,P::prime)
local n,m,p,U,X,S,B,i,j,R,C,rank,ij,x,indsA;
    n := LinearAlgebra:-RowDimension(A);
    m := LinearAlgebra:-ColumnDimension(A);
printf("minor: %d x %d\n",n,m);
    if nargs=1 then p := 2^62-57; else p := P fi;
    U := rand(p);
    X := indets(A);
    S := {seq( x=U(), x in X )};
    B := Matrix(n,m,datatype=integer[8],order=C_order);
    #for i to n do for j to m do B[i,j] := Eval(A[i,j],S) mod p; od od;
    indsA := [indices(A)];
printf("minor: #nonzero=%d\n",nops(indsA));
    for ij in indsA do B[op(ij)] := Eval(A[op(ij)],S) mod p; od;
    R := Array(1..n,datatype=integer[4]);
    C := Array(1..m,datatype=integer[4]);
printf("minor: starting elimination\n");
    rank := Cminor(B,n,m,R,C,p);
#printf("rank=%d\n",rank);
    R := [seq( R[i]+1, i=1..rank )];
    C := [seq( C[i]+1, i=1..rank )];
#printf("R=%a\n",R);
#printf("C=%a\n",C);
    return rank, R, C;
end: