###################################################################################
######      Compute NullSpace over Q, and solve linear equations over Q      ######
###################################################################################
#
# Mark van Hoeij, 2020.
#
# Return the NullSpace of a matrix over Q. The optional argument SkipCheck
# is used if you want to skip the correctness check.
#
# Example:
#
# M := Matrix( [[1,2,0,3], [1,2,0,4]] );
#
# RationalDenseNullSpace(M, true);  # Skip correctness check.
#
# RationalDenseNullSpace(M);        # Proven output.
#
# Input must be over Q, do type checking before calling this.
#
# The program also returns the list of pivots and non-pivots (i.e. basic and free
# variables) because most applications will need that data.

RationalDenseNullSpace := proc(
    M::Matrix,
    {skiptest::truefalse:=false},
    {skipcheck::truefalse:=false},
$)
local A, B, N, c, d, i, inv, j, k, m, npiv, p, piv, r, rm, sm, rnk, sig;

    if not skipcheck then
        if not type(M,'Matrix'('rational')) then
            error "the input must be a Matrix with rational entries";
        end if;
    end if;
    r,c := op(1,M);
    userinfo(4,{'LREtools','SolveTools'},cat("Solving a ", r, " by ", c, " system"));
    k := kernelopts(wordsize);
    d := 'integer'[k/8];
    p := 2^(k/2 - 4); # Max prime is prevprime( 2^(k/2) )
    m := 1;

    if 0 in {rtable_dims(M)} then
        # this 1. is trivial 2. triggers wierd errors in LinearAlgebra:-Modular
        return LinearAlgebra:-NullSpace(M);
    end if;

    A := Matrix(r, c, 'datatype'=d);

    do
        p := nextprime(p);
        try
            (LinearAlgebra:-Modular:-Mod(p, M, A));
        catch "division by zero":
            userinfo(1,{'LREtools','SolveTools'},"discarding prime", p);
            next;
        end try;

        # Only do a full Row Reduction if m > 1 (then rank < c is likely)
        LinearAlgebra:-Modular:-RowReduce(p,A,r,c,c,0,0,'rnk','sig',0, evalb(m > 1));
        if m=1 then
            userinfo(4,{'LREtools','SolveTools'},cat(c-rnk," solutions"))
        end if;

        if rnk = c then
            return Matrix(c,0), [seq(i,i=1..c)], []
        elif m=1 then
            # NullSpace non-trivial so we need a full Row Reduction:
            LinearAlgebra:-Modular:-RowReduce(p,A,r,c,c,0,0,0,0,0,true)
        end if:

        if m = 1 or rnk > rm or (rnk = rm and sig < sm) then
            rm := rnk;
            sm := sig;
            if m>1 then
                userinfo(1,{'LREtools','SolveTools'},
                    "discarding previous prime(s)", m);
            end if;
            m := 1;
        elif rnk < rm or (rnk = rm and sig > sm) then
            # Discarding a prime should be highly unlikely,
            # if it ever occurs then search for a bug.
            userinfo(1,{'LREtools','SolveTools'},"discarding prime", p);
            next;
        end if;

        if m = 1 then
            piv := Vector(rnk); # pivots
            npiv := Vector(c-rnk);  # non-pivots
            inv := Vector(c);
            i := 1; k := 1;
            for j to c do
                if i>r or A[i,j]=0 then
                    npiv[k] := j;
                    inv[j] := k;
                    k := k+1;
                else
                    piv[i] := j;
                    inv[j] := i;
                    i := i+1;
                end if;
            end do;
            inv := convert(inv,list);
            piv := convert(piv,list);
            npiv := convert(npiv,list);
        end if;

        B := LinearAlgebra:-SubMatrix(A, 1..rnk, npiv);
        if m=1 then
            N := B;
        else
            N := chrem([N,B], [m, p]);
            B := 0; # hopefully this helps gc()?
        end if;

        m := m * p;
        B := -iratrecon(N, m, 'scaled');
        if B <> FAIL and length(m) > 3 + max(map(length,B)) then
            j := LinearAlgebra:-IdentityMatrix(c-rnk);
            B := <seq(LinearAlgebra:-Row(`if`(member(i,piv),B,j),inv[i]), i=1..c)>;
            if (_npassed>1 and skiptest) then
                return B, piv, npiv;
            end if;

            local res := LinearAlgebra:-MatrixMatrixMultiply(M, B);
            if rtable_scanblock(res, [rtable_dims(res)], 'NonZeros') = 0 then
                return B, piv, npiv;
            else
                userinfo(1,{'LREtools','SolveTools'},"test failed, trying an additional prime");
                #error "test failed, additional prime needed";
            end if;
        end if;

    end do;

end proc:
# If one chooses to omit a deterministic check (with skiptest = true),
# then one may want a probabilistic check in the form of an extra prime.
# We skip an extra prime because currently (Maple 2020), iratrecon already
# uses twice as many primes as are needed. To use fewer primes see:
#   cs.uwaterloo.ca/~astorjoh/issac11ratrecon.pdf


##PROCEDURE(help,nospec,label="RationalDense") SolveTools[LinearSolvers][RationalDense]
##TITLE(halfline="fast linear solve over the field of rational numbers") SolveTools[LinearSolvers][RationalDense]
##
##CALLINGSEQUENCE
##- RationalDense('eq', 'var')
##- RationalDense('eq', 'var', `skiptest`)
##
##PARAMETERS
##- 'eq' : list or set of equations linear in 'var' with rational coefficients
##- 'var' : list or set of variables for 'eq'
##
##DESCRIPTION
##- The ~RationalDense~ command computes the
## solution(s) of the input dense rational linear system.
## If there exist no solutions 'NULL' is returned. It follows the same
## form of input and output as "SolveTools:-Linear".
##- If the option `skiptest` is given, any correctness checks
## are skipped and the answer has a small chance of being incorrect.
## it should only be used if the user has an efficient way to check the
## answer in their result. (The test is not slow, but does generally take
## several times more computation time than the rest of the computation.)
##- For more details on the underlying algorithm, see
## "RationalDenseNullSpace" which does the work.
##
##EXAMPLES
##> with(SolveTools[LinearSolvers]):
##-(lead=indent) Basic system
##> sys := {x+y=2, x-y=0}:
##> vars := {x,y}:
##> RationalDense(sys, vars);
##< {x = 1, y = 1}
##-(lead=indent) No solution
##> sys := {x+y=2, 2*x+2*y=0}:
##> vars := {x,y}:
##> RationalDense(sys, vars);
##< NULL
##-(lead=indent) Multiple solutions
##> sys := {x+y=2, 2*x+2*y=4}:
##> vars := {x,y}:
##> RationalDense(sys, vars);
##< {x = 2-y, y = y}
##
##SEEALSO
##- "SolveTools:-LinearSolvers:-RationalDenseNullSpace"
##- "SolveTools:-Linear"
##XREFMAP
##- RationalDenseNullSpace : Help:SolveTools:-LinearSolvers:-RationalDenseNullSpace


# Same syntax as SolveTools:-Linear.
#
# This is a dense method.
#
# The input must be over Q, do type checking before calling this:
#
RationalDense := proc(
    eq::{set,list},
    var::{set,list},
    {skiptest::truefalse:=false},
    {skipcheck::truefalse:=false},
$)
local E,i,j,v,n,V,M,e,c,t;

    E := map(i -> `if`(type(i,`=`),lhs(i)-rhs(i),i), eq);
    v := [op(var),1];
    n := nops(v);

    if E = {} or n = 1 then
        return {};
    end if;

    for i to n do
        V[v[i]] := i
    end do;

    M := Matrix(nops(E),n);

    i := 0;
    for e in E do
        i := i+1;
        try
            c := [coeffs(e, v, 't')];
        catch:
            error "input system must be linear in the variables with rational coefficients";
        end try;

        if not skipcheck then
            if remove(type,c,'rational')<>[] then
                error "input system must be linear in the variables with rational coefficients";
            end if;

            if {t} minus {op(var),1}<>{} then
                error "input system must be linear in the variables with rational coefficients";
            end if;
        end if;

        t := [t];
        for j to nops(c) do
            M[ i, V[t[j]] ] := c[j];
        end do;
    end do;

    (M, E, e) := RationalDenseNullSpace(M, _passed[3..-1]);

    # E = basic variables
    # e =  free variables
    if member(n,E) then
        # n'th variable is 1, if it is not free then no solutions.
        return NULL;
    else
        n := nops(e);
        return {seq(v[i]=add(M[i,j]*v[e[j]], j=1..n),i=E),
            seq(v[i]=v[i], i=e[1..-2])};
    end if;

end proc:

