#############################################################################
##
#W  involutive-cp.gi      GAP4 package IBNP      Gareth Evans & Chris Wensley
##

#############################################################################
##
#M  PommaretDivision( <alg> <polys> <ord> )
##
InstallMethod( PommaretDivision, "generic method for list of monomials",
    true, [ IsAlgebra, IsList, IsMonomialOrdering ], 0,
function( A, mons, ord )

    local nmons, mvars, vars, i, posi, emoni;
    nmons := Length( mons );
    mvars := ListWithIdenticalEntries( nmons, 0 );
    vars := OccuringVariableIndices( ord );
    if ( vars = true ) then
        Error( "specify the variable list in the ordering" );
    fi;
    for i in [1..nmons] do   ## for each monnomial 
        emoni := ExtRepPolynomialRatFun( mons[i] )[1];
        posi := Position( vars, emoni[1] );
        mvars[i] := [1..posi];  ## [1 .. first letter]
    od;
    return mvars;
end );

#############################################################################
##
#M  ThomasDivision( <alg> <polys> <ord> )
##
InstallMethod( ThomasDivision, "generic method for list of monomials",
    true, [ IsAlgebra, IsList, IsMonomialOrdering ], 0,
function( A, mons, ord )

    local genA, ngens, nmons, vars, e, max, j, k, p, i, m, mvars;
    genA := GeneratorsOfLeftOperatorRingWithOne( A );
    ngens := Length( genA );
    nmons := Length( mons );
    vars := OccuringVariableIndices( ord );
    if ( vars = true ) then
        Error( "specify the variable list in the ordering" );
    fi;
    ## find the maximal exponents
    max := ListWithIdenticalEntries( ngens, 0 );
    for j in [1..nmons] do
        e := ExtRepPolynomialRatFun( mons[j] )[1];
        for k in [1..Length(e)/2] do 
            p := k+k-1;
            i := Position( vars, e[p] );
            m := e[p+1];
            if ( m > max[i] ) then 
                max[i] := m;
            fi;
        od;
    od;
    Info( InfoIBNP, 2, "maximal exponents = ", max );
    ## find the multiplicative variables
    mvars := ListWithIdenticalEntries( nmons, 0 );
    for j in [1..nmons] do
        mvars[j] := [ ];
        e := ExtRepPolynomialRatFun( mons[j] )[1];
        for k in [1..Length(e)/2] do
            p := k+k-1;
            i := Position( vars, e[p] );
            m := e[p+1];
            if ( m = max[i] ) then 
                Add( mvars[j], i );
            fi;
        od;
    od;
    Info( InfoIBNP, 2, "multiplicative variables = ", mvars );
    return mvars;
end );

#############################################################################
##
#M  JanetDivision( <alg> <polys> <ord> )
##
InstallMethod( JanetDivision, "generic method for list of monomials",
    true, [ IsAlgebra, IsList, IsMonomialOrdering ], 0,
function( A, mons, ord )

    local  genA, ngens, nmons, vpos, vars, ordfn, expos, j, e, k, i, 
           mvars, done, ej, L, subj, ei, ok;
    genA := ShallowCopy( GeneratorsOfLeftOperatorRingWithOne( A ) );
    ngens := Length( genA );
    Info( InfoIBNP, 2, "in JanetDivision, mons = ", mons );
    nmons := Length( mons );
    vpos := List( genA, g -> ExtRepPolynomialRatFun(g)[1][1] );
    vars := OccuringVariableIndices( ord );
    if ( vars = true ) then
        Error( "specify the variable list in the ordering" );
    fi;
    ordfn := MonomialComparisonFunction( ord );
    ## create list of exponents, including zeroes
    expos := ListWithIdenticalEntries( nmons, 0 );
    for j in [1..nmons] do
        expos[j] := ListWithIdenticalEntries( ngens, 0 );
        e := ExtRepPolynomialRatFun( mons[j] )[1];
        Info( InfoIBNP, 2, "e(mons[j]) = ", e );
        for k in [1..Length(e)/2] do
            i := Position( vpos, e[k+k-1] );
            expos[j][i] := e[k+k];
        od;
    od;
    Info( InfoIBNP, 2, "expos = ", expos );
    ## now set up the lists of multiplicative variables
    mvars := ListWithIdenticalEntries( nmons, 0 );
    for j in [1..nmons] do
        mvars[j] := [ ];
    od;
    ## test for x_1 .. x_{n-1}
    for j in [1..ngens-1] do
        done := ListWithIdenticalEntries( nmons, false );
        for k in [1..nmons] do
            if not done[k] then
                done[k] := true;
                ej := expos[k][j];
                L := [k];
                subj := expos[k]{[j+1..ngens]};
                for i in [k+1..nmons] do
                    if ( subj = expos[i]{[j+1..ngens]} ) then
                        done[i] := true;
                        Add( L, i );
                        if ( expos[i][j] > ej ) then
                            ej := expos[i][j];
                        fi;
                    fi;
                od;
                for i in L do
                    if ( subj = expos[i]{[j+1..ngens]} ) 
                        and ( expos[i][j] = ej ) then
                        Add( mvars[i], j );
                    fi;
                od;
            fi;
        od;
    od;
    ## now test for x_n
    ei := expos[1][ngens];
    for k in [2..nmons] do
        if ( expos[k][ngens] > ei ) then
            ei := expos[k][ngens];
        fi;
    od;
## Error("here");
    for k in [1..nmons] do
        if ( expos[k][ngens] = ei ) then
            Add( mvars[k], ngens );
        fi;
    od;
    return mvars;
end );

#############################################################################
##
#M  DivisionRecord( <args> )
##
BindGlobal( "DivisionRecord", 
function( arg )
    local nargs, A, polys, ord;
    nargs := Length( arg );
    if not ( nargs = 3 ) then 
        Error( "expecting arguments [ A, polys, ord ]" ); 
    fi;
    A := arg[1];
    polys := arg[2];
    ord := arg[3];
    if not IsNearAdditiveMagma( A ) then 
        Error( "expecting an algebra as first parameter" ); 
    fi;
    if IsCommutative( A ) then 
        return DivisionRecordCP( A, polys, ord ); 
    else 
        return DivisionRecordNP( A, polys, ord );
    fi;
end );

############################################################################
##
#M IsDivisionRecord
## 
InstallMethod( IsDivisionRecord, "generic method for records", true, 
    [ IsRecord ], 0,
function( r )
    return ( IsBound\.( r, RNamObj( "div" ) ) ) 
       and ( IsBound\.( r, RNamObj( "mvars" ) ) )
       and ( IsBound\.( r, RNamObj( "polys" ) ) );
end );

#############################################################################
##
#M  DivisionRecordCP( <alg> <polys> <ord> )
##
InstallMethod( DivisionRecordCP, "generic method for list of monomials",
    true, [ IsAlgebra, IsList, IsMonomialOrdering ], 0,
function( A, polys, ord )

    ## Overview: returns local overlap-based multiplicative variables
    ## Detail: this function implements various algorithms
    ## described in the thesis "Noncommutative Involutive Bases"
    ## for finding multiplicative variable for a set of polynomials 

    local  pl, genA, ordfn, ngens, npolys, mons, mvars;

    pl := InfoLevel( InfoIBNP );
    if( polys = [ ] ) then 
        return [ ];
    fi;
    genA := ShallowCopy( GeneratorsOfLeftOperatorRingWithOne( A ) );
    ordfn := MonomialComparisonFunction( ord );
    Sort( genA, ordfn );
    ngens := Length( genA );
    ## set up arrays
    npolys := Length( polys );
    mons := List( polys, p -> LeadingMonomialOfPolynomial( p, ord ) );
    Info( InfoIBNP, 2, "mons = ", mons );
    if ( CommutativeDivision = "Pommaret" ) then 
        mvars := PommaretDivision( A, mons, ord );
    elif ( CommutativeDivision = "Thomas" ) then
        mvars := ThomasDivision( A, mons, ord );
    elif ( CommutativeDivision = "Janet" ) then
        mvars := JanetDivision( A, mons, ord );
    else
        Error( "invalid CommutativeDivision" );
    fi;
    return rec( div := CommutativeDivision, mvars := mvars, polys := polys );
end );

#############################################################################
##
#M  IPolyReduce( <alg> <poly> <drec> <ord> )
##
BindGlobal( "IPolyReduce",
function( arg )
    local nargs, A, p, drec, ord, polys, mvars;
    nargs := Length( arg );
    if not ( nargs = 4 ) then 
        Error( "expecting arguments [ A, pol, overlaprec, ordering ]" ); 
    fi;
    A := arg[1];
    p := arg[2];
    drec := arg[3];
    polys := drec.polys;
    mvars := drec.mvars;
    ord := arg[4];
    if not IsNearAdditiveMagma( A ) then 
        Error( "expecting an algebra as first parameter" ); 
    fi;
    if IsCommutative( A ) then
        ## add some tests
        return IPolyReduceCP( A, p, drec, ord ); 
    else
        ## add more tests
        return IPolyReduceNP( A, p, drec, ord );
    fi;
end );

#############################################################################
##
#M  IPolyReduceCP( <alg> <poly> <drec> <ord> )
##
InstallMethod( IPolyReduceCP, 
    "generic method for a poly, record with fields [polys, mult vars]",
    true, [ IsPolynomialRing, IsPolynomial, IsRecord, IsMonomialOrdering ], 0,
function( A, pol, drec, ord )
##
## Overview: Reduces 2nd arg w.r.t. 3rd arg (polys and vars)
##
## Detail: drec is a record containing a list 'polys' of polynomials
## and their multiplicative variables 'vars'.
## Given a polynomial _pol_, this function involutively
## reduces the polynomial with respect to the given polys
## with associated left and right multiplicative variables vars.

    return CombinedIPolyReduceCP( A, pol, drec, ord, false );
end );

#############################################################################
##
#M  LoggedIPolyReduce( <alg> <poly> <drec> <ord> )
##
BindGlobal( "LoggedIPolyReduce",
function( arg )
    local nargs, A, p, drec, ord, polys, mvars;
    nargs := Length( arg );
    if not ( nargs = 4 ) then 
        Error( "expecting arguments [ A, pol, overlaprec, ordering ]" ); 
    fi;
    A := arg[1];
    p := arg[2];
    drec := arg[3];
    polys := drec.polys;
    mvars := drec.mvars;
    ord := arg[4];
    if not IsNearAdditiveMagma( A ) then 
        Error( "expecting an algebra as first parameter" ); 
    fi;
    if IsCommutative( A ) then
        ## add some tests
        return CombinedIPolyReduceCP( A, p, drec, ord, true ); 
    else
        ## add more tests
        return CombinedIPolyReduceNP( A, p, drec, ord, true );
    fi;
end );

#############################################################################
##
#M  LoggedIPolyReduceCP( <alg> <poly> <drec> <ord> )
##
InstallMethod( LoggedIPolyReduceCP, 
    "generic method for a poly, record with fields [polys, mult vars]",
    true, [ IsPolynomialRing, IsPolynomial,
            IsRecord, IsMonomialOrdering ], 0,
function( A, pol, drec, ord )
##
## Overview: Reduces 2nd arg w.r.t. 3rd arg (polys and vars)
## just as in IPolyReduceCP, but returns a record contsaining the reduced
## polynomial and also logged information showing how the result is 
## obtained from the original polynomial.
##
    return CombinedIPolyReduceCP( A, pol, drec, ord, true );
end );

#############################################################################
##
#M  CombinedIPolyReduceCP( <alg> <poly> <drec> <ord> <logging> )
##
InstallMethod( CombinedIPolyReduceCP, 
    "generic method for a poly, record with fields [polys, mult vars]",
    true, [ IsPolynomialRing, IsPolynomial,
            IsRecord, IsMonomialOrdering, IsBool ], 0,
function( A, pol, drec, ord, logging )

    local  polys, mvars, genA, vpos, a, z, front, back, eback, 
           numpolys, mons, coeffs, reduced, mon, coeff, term,
           i, poli, moni, fac, coeffi, efac, lenf, ok, logs;

    if not IsDivisionRecord( drec ) then
        Error( "third argument should be an overlap record" );
    fi;
    Info( InfoIBNP, 3, "in IPolyReduceCP: pol = ", pol, 
                       ",  logging = ", logging );
    polys := drec.polys;
    numpolys := Length( polys );
    logs := List( [1..numpolys], i -> Zero( A ) );
    mvars := drec.mvars;
    genA := GeneratorsOfLeftOperatorRingWithOne( A );
    vpos := List( genA, g -> ExtRepPolynomialRatFun(g)[1][1] );
    a := genA[1];
    z := a-a;  ## zero polynomial;
    front := z;
    back := pol;
    eback := ExtRepPolynomialRatFun( back );
    mons := List( polys, p -> LeadingMonomialOfPolynomial( p, ord ) );
    coeffs := List( polys, p -> LeadingCoefficientOfPolynomial( p, ord ) );
    ## we now recursively reduce every term in the polynomial
    ## until no more reductions are possible
    while ( eback <> [ ] ) do
        reduced := true;
        while reduced and ( eback <> [ ] ) do
            reduced := false;
            mon := LeadingMonomialOfPolynomial( back, ord );
            coeff := LeadingCoefficientOfPolynomial( back, ord );
            i := 0;
            while ( i < numpolys ) and ( not reduced ) do  
                ## for each polynomial in the list
                i := i+1;
                poli := polys[i];
                moni := mons[i];  ## pick a test monomial
                fac := mon/moni;
                if IsPolynomial( fac ) then  ## moni divides mon
                    Info( InfoIBNP, 3, "divisor found: fac = ", fac );
                    coeffi := coeffs[i];
                    efac := ExtRepPolynomialRatFun( fac )[1];
                    Info( InfoIBNP, 3, "fac = ", fac, ",  efac = ", efac );
                    lenf := Length( efac );
                    ok := ( lenf = 0 ) or ForAll( [1,3..lenf-1], 
                            j -> Position( vpos, efac[j] ) in mvars[i] );
                    ## if an involutive division was found
                    if ok then
                        ## indicate a reduction has been carried out 
                        ## to exit the loop
                        ## pick the divisor's leading coefficient
                        ## pick 'nice' cancelling coefficients
                        term := (coeff/coeffi)*fac;
                        back := back - term*poli;
                        if logging then
                            logs[i] := logs[i] + term;
                        fi;
                        Info( InfoIBNP, 2, "reduced to: ", front+back );
                        eback := ExtRepPolynomialRatFun( back );
                        reduced := true;
                    fi;
                fi;
            od;
        od;
        if ( eback <> [ ] ) then
            ## no reduction of the lead term, so add it to front
            front := front + coeff*mon;
            back := back - coeff*mon;
            eback := ExtRepPolynomialRatFun( back );
        fi;
    od;
    if ( front <> z ) then
        if ( LeadingCoefficientOfPolynomial( front, ord ) < 0 ) then
            front := (-1)*front;
            if logging then
                logs := List( logs, t -> (-1)*t );
            fi;
        fi;
    fi;
    if not logging then
        return front;  ## return the reduced and simplified polynomial
    else
        return rec( result := front, logs := logs, polys := polys );
    fi;
end );

#############################################################################
##
#M  IAutoreduce( <args> )
##
BindGlobal( "IAutoreduce",
function( arg )
    local nargs, A, polys, ord;
    nargs := Length( arg );
    if not ( nargs = 3 ) then 
        Error( "expecting arguments [ A, polys, ord ]" ); 
    fi;
    A := arg[1];
    polys := arg[2];
    ord := arg[3];
    if not IsNearAdditiveMagma( A ) then
        ## add some tgests
        Error( "expecting an algebra as first parameter" ); 
    fi;
    if IsCommutative( A ) then
        ## add more tests
        return IAutoreduceCP( A, polys, ord ); 
    else 
        return IAutoreduceNP( A, polys, ord ); 
    fi;
end );

#############################################################################
##
#M  IAutoreduceCP( <alg> <polys> <ord> )
##
InstallMethod( IAutoreduceCP, "generic method for list of polys",
    true, [ IsAlgebra, IsList, IsMonomialOrdering ], 0,
function( A, polys, ord )

    ## Overview: Autoreduces a list of polynomials recursively 
    ## until no more reductions are possible 
    ## Detail: This function involutively reduces each member of a
    ## list of polynomials w.r.t. all the other members of the list, 
    ## removing the polynomial from the list if it is involutively
    ## reduced to 0. Iterate until no more such reductions are possible. 

    local pl, div, npolys, genA, a, z, pos, old, reduced, ordfn,
          oldPoly, nrec, newvars, vars, drec, oldvars, new, newPoly;

    pl := InfoLevel( InfoIBNP );
    div := CommutativeDivision;
    npolys := Length( polys );
    if( npolys < 2 ) then 
        ## the input basis is empty or consists of a single polynomial
        return polys;
    fi;
    genA := GeneratorsOfLeftOperatorRingWithOne( A );
    a := genA[1];
    z := a-a;  ## zero polynomial;
    ## start by reducing the final element (working backwards means
    ## that less work has to be done calculating multiplicative variables)
    ## if we are using a local division and the basis is sorted by DegRevLex,
    ## the last polynomial is irreducible so we do not have to consider it.
    ## make a copy of the input basis for traversal
    old := StructuralCopy( polys );
    ordfn := MonomialComparisonFunction( ord );
    reduced := true;
    while reduced do
        reduced := false;
        Sort( old, ordfn );
        npolys := Length( old );
        pos := npolys;
        drec := DivisionRecordCP( A, old, ord );
        oldvars := drec.mvars;
        if ( pl > 2 ) then 
            Print( "old = ", old, "\n" );
            Print( "oldvars = ", oldvars, "\n" );
        fi;
        while ( pos > 0 ) and ( not reduced ) do   ## for each poly in old
            ## extract pos-th element of the basis and reduce using the rest
            oldPoly := old[pos];
            if ( pl > 2 ) then 
                Print( "oldPoly = ", oldPoly, "\n" );
            fi;
            ## construct basis without oldPoly
            ## calculate multiplicative Variables if using a local division
            newvars := Concatenation( oldvars{[1..pos-1]},
                                      oldvars{[pos+1..npolys]} ); 
            new := Concatenation( old{[1..pos-1]}, old{[pos+1..npolys]} );
            nrec := rec( div := div, mvars := newvars, polys := new );
            ## involutively reduce old polynomial w.r.t. the truncated list
            newPoly := IPolyReduceCP( A, oldPoly, nrec, ord );
            Info( InfoIBNP, 2, "newPoly = ", newPoly );
            ## if the polynomial did not reduce to 0
            if ( newPoly = z ) then 
                ## the polynomial reduced to zero
                ## remove the polynomial from the list
                new := Concatenation( new{[1..pos-1]},
                                      new{[pos+1..npolys]} );
                pos := pos - 1;
                npolys := npolys - 1;
            elif ( oldPoly = newPoly ) then 
                ## we may proceed to look at the next polynomial 
                pos := pos - 1;
            else
                ## otherwise some reduction took place so start again
                ## add the new polynomial onto the list
                reduced := true;
                old := Concatenation( new, [ newPoly ]  );
                Info( InfoIBNP, 2, "reduction!  now old = ", old );
            fi;
        od;
    od; 
    Sort( old, ordfn );
    ##  return the fully autoreduced basis
    if ( old = polys ) then ## no change
        return true;
    else
        return old;
    fi;
end );

#############################################################################
##
#M  InvolutiveBasis( <args> )
##
BindGlobal( "InvolutiveBasis",
function( arg )
    local nargs, A, polys, ord;
    nargs := Length( arg );
    if not ( nargs = 3 ) then 
        Error( "expecting arguments [ A, polys, ord ]" ); 
    fi;
    A := arg[1];
    polys := arg[2];
    ord := arg[3];
    if not IsNearAdditiveMagma( A ) then
        ## add some tests
        Error( "expecting an algebra as first parameter" ); 
    fi;
    if IsCommutative( A ) then
        ## add more tests
        return InvolutiveBasisCP( A, polys, ord ); 
    else 
        return InvolutiveBasisNP( A, polys, ord ); 
    fi;
end );

#############################################################################
##
#M  InvolutiveBasisCP( <alg> <polys> <ord> )
##
InstallMethod( InvolutiveBasisCP, 
    "generic method for commutative algebra + list of polys + ordering",
    true, [ IsAlgebra, IsList, IsMonomialOrdering ], 0,
function( A, polys, ord )

    ##  Implements Algorithm 9 for computing commutative involutive bases. 

    local  pl, npolys, genA, a, z, ngens, all, basis, vars, done, 
           ordfn, basis0, nbasis, drec, mvars, nvars, prolong, 
           i, j, npro, found, u, v;

    pl := InfoLevel( InfoIBNP );
    Info( InfoIBNP, 2, "Computing an Involutive Basis...");
    npolys := Length( polys );
    if( npolys < 2 ) then 
        ## the input basis is empty or consists of a single polynomial
        return polys;
    fi;
    genA := GeneratorsOfLeftOperatorRingWithOne( A );
    a := genA[1];
    z := a-a;  ## zero polynomial;
    ngens := Length( genA );
    all := [ 1..ngens ];
    basis := IAutoreduceCP( A, polys, ord );
    if ( basis = true ) then ## no reduction
        basis := polys;
    fi;
    Info( InfoIBNP, 2, "initial basis = ", basis ); 
    done := false;
    ordfn := MonomialComparisonFunction( ord );
    while ( not done ) do
        Sort( basis, ordfn );
        Info( InfoIBNP, 1, "restarting with basis:\n", basis );
        basis0 := ShallowCopy( basis );
        basis := IAutoreduceCP( A, basis0, ord );
        if ( basis = true ) then 
            basis := basis0;
        else
            Info( InfoIBNP, 1, "after autoreduction basis = \n", basis );
        fi;
        nbasis := Length( basis );
        drec := DivisionRecordCP( A, basis, ord );
        Info( InfoIBNP, 1, "division record for basis: ", drec );
        mvars := drec.mvars;
        nvars := List( mvars, v -> Difference( all, v ) );
        Info( InfoIBNP, 2, "mvars = ", mvars );
        Info( InfoIBNP, 2, "nvars = ", nvars );
        ## construct the prolongations
        prolong := [ ];
        for i in [ 1..nbasis ] do
            for j in nvars[i] do
                Add( prolong, genA[j]*basis[i] );
            od;
        od;
        Sort( prolong, ordfn );
        Info( InfoIBNP, 1, "prolongations = ", prolong );
        npro := Length( prolong ); 
        found := false;
        for i in [1..npro] do
            u := prolong[i];
            v := IPolyReduceCP( A, u, drec, ord );
            Info( InfoIBNP, 2, "u = ", u, ",  v = ", v );
            if ( v <> z ) then 
                found := true; 
                Add( basis, v );
                Info( InfoIBNP, 2, "adding ", v );
                if ( pl > 0 ) and 
                   ( Length(basis) > InvolutiveAbortLimit ) then
                      return fail;
                fi;
                drec := DivisionRecordCP( A, basis, ord );
            fi;
        od;
        if not found then
            done := true;
        fi;
    od;
    return DivisionRecordCP( A, basis, ord);
end );

#############################################################################
##
#E  involutive-cp.gi . . . . . . . . . . . . . . . . . . . . . . .  ends here
##