newPackage(
    "FGLM",
    Version => "1.1.0",
    Date => "September 3, 2019",
    Authors => {
        { Name => "Dylan Peifer",   Email => "djp282@cornell.edu", HomePage => "https://math.cornell.edu/~djp282" },
        { Name => "Mahrud Sayrafi", Email => "mahrud@umn.edu",     HomePage => "https://math.umn.edu/~mahrud" }
        },
    Keywords => {"Groebner Basis Algorithms"},
    Headline => "Groebner bases via the FGLM algorithm"
    )

-*
Copyright (C) 2019 Dylan Peifer and Mahrud Sayrafi

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
*-

export {"fglm"}

importFrom_Core { "RawMutableMatrix", "raw", "rawLUincremental", "rawTriangularSolve" }

-------------------------------------------------------------------------------
--- top level functions
-------------------------------------------------------------------------------

---------------------------------------------------------------------------
-- See Section 2.4.4, Algorithm 2.5 of Thibaut Verron's thesis for details:
-- https://thibautverron.github.io/doc/2016-these.pdf
---------------------------------------------------------------------------
fglm = method()
fglm(Ideal,         Ring) := GroebnerBasis => (I1, R2) -> fglm(gb I1, R2)
fglm(GroebnerBasis, Ring) := GroebnerBasis => (G1, R2) -> (
    -- G1 = a Groebner basis for I1 in R1
    -- R2 = a polynomial ring
    -- returns Groebner basis G2 of the ideal generated by G1 in the ring R2

    R1 := ring G1;
    kk := coefficientRing R1;
    I1 := ideal gens G1; -- TODO: make a github issue add gb to cache
    -- FIXME: this is a hecking hack
    gbOpt := new GroebnerBasisOptions from {HardDegreeLimit => null, Syzygies => false, SyzygyRows => 0};
    I1.generators.cache#gbOpt = G1;
    if R1 == I1 then return forceGB(sub(gens G1, R2));
    if dim I1 > 0 then error "expected zero-dimensional ideal";
    if #gens R1 != #gens R2 then error "expected the same number of generators";

    M := multiplicationMatrices G1;
    m := numcols M#0;
    n := #gens R2;

    -- elements in the groebner basis
    G2 := new MutableHashTable from {}; -- leading term => gb element
    -- elements in the staircase
    B2 := new MutableHashTable from {1_R2 => true};
    B2' := matrix {{1_R2}};

    -- initiating LU-decomposition matrices
    P := toList(0..m-1);
    v := transpose mutableMatrix {{1_kk} | toList ((m-1):0)};
    LU := mutableMatrix map(kk^m, kk^(m+1), 0);
    lambda := transpose mutableMatrix {toList (m:0_kk)};
    LUincremental(P, LU, v, 0);

    -- normal form translation table
    -- Note: we want dense mutable matrices
    -- TODO: get rid of this and use LU
    V := new MutableHashTable from {1_R2 => v};

    -- list of elements between the staircase and groebner basis generators
    -- Note: use a heap in the engine for this
    S := new MutableHashTable from for i to n - 1 list R2_i * 1_R2 => (i, 1_R2);

    -- ~2s up to this point
    s := 1;
    while #S > 0 do (
	-- ~14s in this line
	-- TODO: use O(1) min for fun and profit
	(elt, vals) := min pairs S;
	remove(S, elt);
	if any(keys G2, lt -> elt % lt == 0) then continue;
	(i, mu) := vals;
	v = M#i * V#mu;

	-- ~6.2s in this line
	-- TODO: can we not copy v?
	P = LUincremental(P, LU, submatrix(v, , {0}), s);
	if s == m or (s < m and LU_(s, s) == 0) then (
	    -- ~1.3s in this branch
	    backSub(submatrix(LU, toList(0..s-1), toList(0..s)), lambda, s);
	    -- TODO: don't remake a matrix every time
	    g := elt - B2' * matrix submatrix(lambda, toList(0..s-1),);
	    G2#elt = g;
	    ) else (
	    -- ~1.6s in this branch
	    s = s + 1;
	    -- TODO: add elt to VS and row reduce here
	    V#elt = v;
	    B2#elt = true;
	    B2' = B2' | matrix {{elt}};
	    -- Add the product of elt and generators of R2 to S
	    for j to n - 1 do if not B2#?(R2_j * elt) then S#(R2_j * elt) = (j, elt);
	    );
	);
    forceGB(matrix {values G2})
    )

-------------------------------------------------------------------------------
-- Fast incremental LU decomposition code
-------------------------------------------------------------------------------
-- Inputs:
--   permutation list of size n
--   n x n+1 mutable matrix LU
--   n x 1 mutable matrix v
--   integer m
-- Consequences:
--   Updates LU according to v
--   Updates v to r = v % L
-- Outputs:
--   Updated permutation list
LUincremental = method()
LUincremental(List, MutableMatrix, MutableMatrix, ZZ) := List => (P, LU, v, m) -> (
    rawLUincremental(P, raw LU, raw v, m)
    )

-------------------------------------------------------------------------------
-- Fast back substitution code
-------------------------------------------------------------------------------
-- Inputs:
--   n x n+1 mutable matrix (U|v)
--   n x 1 mutable matrix x
--   integer m
-- Consequences:
--   Updates v to r = v % U
--   Updates x such that v = Ux + r
-- Notes: assumes U is upper triangular, but the actual values above the diagonal need not be zero
-- If LUdecomposition => true, assumes the diagonal to be 1.
backSub = method(Options => {Strategy => null})
backSub(MutableMatrix, MutableMatrix, ZZ) := RawMutableMatrix => opts -> (U, x, m) -> (
    if opts.Strategy === "incremental" then
      rawTriangularSolve(raw U, raw x, m, 3)
    else
      rawTriangularSolve(raw U, raw x, m, 2)
    )

-*
backSub(MutableMatrix, MutableMatrix, ZZ) := Nothing => (U, x, n) -> opts -> (
    for i from 1 to n do (
	x_(n-i, 0) = U_(n-i, n) / U_(n-i,n-i);
	columnAdd(U, n, -x_(n-i, 0), n-i);
	);
    )
*-

-------------------------------------------------------------------------------
-- Fast forward substitution code
-------------------------------------------------------------------------------
-- Inputs:
--   n x n+1 mutable matrix (L|v)
--   n x 1 mutable matrix x
--   integer m
-- Consequences:
--   Updates v to r = v % L
--   Updates x such that v = Lx + r
-- Notes: assumes L is lower triangular, but the actual values below the diagonal need not be zero.
-- If LUdecomposition => true, assumes the diagonal to be 1. This option is mainly used for LUincremental.
forwardSub = method(Options => {Strategy => null})
forwardSub(MutableMatrix, MutableMatrix, ZZ) := RawMutableMatrix => opts -> (L, x, m) -> (
    if opts.Strategy === "incremental" then
      rawTriangularSolve(raw L, raw x, m, 1)
    else
      rawTriangularSolve(raw L, raw x, m, 0)
    )

-*
forwardSub(MutableMatrix, MutableMatrix, ZZ) := Nothing => opts -> (L, x, n) -> (
    for i to n - 1 do (
	x_(i, 0) = if opts.Strategy === "incremental" then L_(i, n) else L_(i, n) / L_(i,i);
	columnAdd(L, n, -x_(i, 0), i);
	);
    )
*-

-------------------------------------------------------------------------------
-- See Section 2.4.4, Algorithm 2.4 of Thibaut Verron's thesis for details:
-- https://thibautverron.github.io/doc/2016-these.pdf
-- Applies more "tricks"
-------------------------------------------------------------------------------
-- TODO: Move to engine
-- TODO: return MutableMatrix
multiplicationMatrices' = method()
multiplicationMatrices'(GroebnerBasis) := List => (G) -> (
    -- G = a GroebnerBasis
    -- returns the matrices giving multiplication by variables in R/I

    R := ring G;
    I := ideal gens G;
    -- FIXME: this is a hecking hack
    gbOpt := new GroebnerBasisOptions from {HardDegreeLimit => null, Syzygies => false, SyzygyRows => 0};
    I.generators.cache#gbOpt = G;
    B := first entries sub(basis (R/I), R); -- TODO: find a way to avoid recomputing GB
    N := new MutableHashTable from for b in B list b => b;
    F := flatten for x in gens R list flatten apply(B, b -> if not N#?(x * b) then x * b else {});
    F = sort F;

    for mu in F do (
	i := position(first entries leadTerm G, g -> mu == leadMonomial g);
	if i =!= null then (
	    g := (gens G)_i_0;
	    N#mu = mu - g // leadCoefficient g; -- Verron has typo in line 9
	    ) else (
	    j := position(F, mu' -> mu % mu' == 0 and any(gens R, x -> mu == x * mu'));
	    mu' := F#j;
	    (gs, cs) := coefficients N#mu';
	    N#mu = sum apply(first entries gs, first entries transpose cs, (g, c) -> N#(g * mu // mu') * c);
	    );
	);

    for x in gens R list lift(last coefficients(matrix{apply(x * B, elt -> N#elt)}, Monomials => B), coefficientRing R)
    )

-- TODO: Make this into a more general function that gets
-- (f: ring elt, S: quotient ring, B: basis) -> (M: multiplication matrix of f)
multiplicationMatrices = method()
multiplicationMatrices(GroebnerBasis) := List => (G) -> (
    -- G = a GroebnerBasis
    -- returns the matrices giving multiplication by variables in R/I

    R := ring G;
    I := ideal gens G;
    -- FIXME: this is a hecking hack
    gbOpt := new GroebnerBasisOptions from {HardDegreeLimit => null, Syzygies => false, SyzygyRows => 0};
    I.generators.cache#gbOpt = G;
    B := basis (R/I); -- TODO: find a way to avoid recomputing GB

    for x in gens R list mutableMatrix lift(last coefficients(x * B, Monomials => B), coefficientRing R)
    )

-------------------------------------------------------------------------------
--- Helper functions for tests
-------------------------------------------------------------------------------

cyclic = method(Options => {CoefficientRing => ZZ/32003, MonomialOrder => GRevLex})
cyclic(ZZ) := Ideal => opts -> (n) -> (
    R := (opts.CoefficientRing)[vars(0..n-1), MonomialOrder => opts.MonomialOrder];
    F := toList apply(1..n-1, d -> sum(0..n-1, i -> product(d, k -> R_((i+k)%n))))
         | {product gens R - 1};
    ideal F
    )

katsura = method(Options => {CoefficientRing => ZZ/32003, MonomialOrder => GRevLex})
katsura(ZZ) := Ideal => opts -> (n) -> (
    n = n-1;
    R := (opts.CoefficientRing)[vars(0..n), MonomialOrder => opts.MonomialOrder];
    L := gens R;
    u := i -> (
	 if i < 0 then i = -i;
	 if i <= n then L_i else 0_R
	 );
    f1 := -1 + sum for i from -n to n list u i;
    F := toList prepend(f1, apply(0..n-1, i -> - u i + sum(-n..n, j -> (u j) * (u (i-j)))));
    ideal F
    )

test = (I1, MO2) -> (
    R1 := ring I1;
    R2 := (coefficientRing R1)(monoid ([gens R1], MonomialOrder => MO2));
    G2 := gb(sub(I1, R2));
    elapsedTime G2' := fglm(I1, R2);
    assert(gens G2 == gens G2')
    )

-------------------------------------------------------------------------------
--- documentation
-------------------------------------------------------------------------------
beginDocumentation()

doc ///
Key
  FGLM
Headline
  Compute Groebner bases via the FGLM algorithm
Description
  Text
    FGLM is a Groebner basis conversion algorithm. This means it takes a
    Groebner basis of an ideal with respect to one monomial order and changes it
    into a Groebner basis of the same ideal over a different monomial
    order. Conversion algorithms can be useful since sometimes when a Groebner
    basis over a difficult monomial order (such as lexicographic or an
    elimination order) is desired, it can be faster to compute a Groebner basis
    directly over an easier order (such as graded reverse lexicographic) and
    then convert rather than computing directly in the original order. Other
    examples of conversion algorithms include the Groebner walk and
    Hilbert-driven Buchberger.

    FGLM performs conversion by doing linear algebra in the quotient ring R/I,
    where I is the ideal generated by the original Groebner basis in the
    polynomial ring R. This requires that I is zero-dimensional.

    In Macaulay2, monomial orders must be given as options to rings. For
    example, the following ideal has monomial order given by graded reverse
    lexicographic (which is also the default order in Macaulay2).

  Example
    R1 = QQ[x,y,z, MonomialOrder => GRevLex]
    I1 = ideal(x*y + z - x*z, x^2 - z, 2*x^3 - x^2*y*z - 1)
  Text
    If we want a Groebner basis of I1 with respect to lexicographic order
    we could substitute the ideal into
    a new ring with that order and compute directly,
  Example
    R2 = QQ[x,y,z, MonomialOrder => Lex];
    I2 = sub(I1, R2);
    gens gb I2  -- performs computation in R2
  Text
    but it may be faster to compute directly in the first order and then use
    FGLM.
  Example
    G1 = gb I1;  -- performs computation in R1
    gens fglm(G1, R2)

  Text
    @HEADER2{"References"}@
  Text
    Further background and details can be found in the following resources:
  Text
    @UL {
        "Cox, Little, O'Shea - Using Algebraic Geometry (2005)",
        "Faugere, Gianni, Lazard, Mora - Efficient Computation of Zero-dimensional
         Groebner Bases by Change of Ordering (1993)",
        "Gerdt, Yanovich - Implementation of the FGLM Algorithm and Finding Roots
         of Polynomial Involutive Systems (2003)"
     }@
  Text
    @HEADER2{"Acknowledgement"}@
  Text
    The C++ implementation of the algorithms in Version 1.2.0 was sponsored by an IMA Coding Sprint at the Cornell University.
Caveat
  The ideal generated by the Groebner basis must be zero-dimensional.
SeeAlso
  groebnerBasis
///

doc ///
Key
  fglm
  (fglm, GroebnerBasis, Ring)
  (fglm, Ideal, Ring)
Headline
  convert a Groebner basis
Usage
  H = fglm(G, R)
  H = fglm(I, R)
Inputs
  G: GroebnerBasis
     the starting Groebner basis
  I: Ideal
     the starting ideal
  R: Ring
     a ring with the target monomial order
Outputs
  H: GroebnerBasis
     the new Groebner basis in the target monomial order
Description
  Text
    FGLM takes a Groebner basis of an ideal with respect to one monomial order
    and changes it into a Groebner basis of the same ideal over a different
    monomial order. The initial order is given by the ring of G and the target
    order is the order in R. When given an ideal I as input a Groebner basis of
    I in the ring of I is initially computed directly, and then this Groebner
    basis is converted into a Groebner basis in the ring R.
  Example
    R1 = QQ[x,y,z];
    I1 = ideal(x^2 + 2*y^2 - y - 2*z, x^2 - 8*y^2 + 10*z - 1, x^2 - 7*y*z);
    R2 = QQ[x,y,z, MonomialOrder => Lex];
    fglm(I1, R2)
Caveat
  The ideal I generated by G must be zero-dimensional. The target ring R must be
  the same ring as the ring of G or I, except with different monomial order. R
  must be a polynomial ring over a field.
SeeAlso
  FGLM
  groebnerBasis
///

-------------------------------------------------------------------------------
--- tests
-------------------------------------------------------------------------------

TEST ///
  debug needsPackage "FGLM"
  R1 = ZZ/101[x,y,z]
  I1 = ideal(x*y + z - x*z, x^2 - z, 2*x^3 - x^2*y*z - 1)
  test(I1, Lex)
///

TEST ///
  debug needsPackage "FGLM"
  R1 = QQ[x,y,z]
  I1 = ideal(x^2 + 2*y^2 - y - 2*z, x^2 - 8*y^2 + 10*z - 1, x^2 - 7*y*z)
  test(I1, Lex)
///

TEST ///
  debug needsPackage "FGLM"
  R1 = QQ[x,y,z]
  I1 = ideal(x^2 + y^2 + z^2 - 2*x, x^3 - y*z - x, x - y + 2*z)
  test(I1, Lex)
///

TEST ///
  debug needsPackage "FGLM"
  R1 = QQ[x,y,z]
  I1 = ideal(x*y + z - x*z, x^2 - z, 2*x^3 - x^2*y*z - 1)
  test(I1, Lex)
///

TEST ///
  -- katsura6
  -- gb: 0.123865
  -- fglm: 0.115399
  debug needsPackage "FGLM"
  I = katsura(6, MonomialOrder=>Lex)
  G1 = elapsedTime gb I
  I = katsura(6)
  R = newRing(ring I, MonomialOrder=>Lex)
  G2 = elapsedTime fglm(I, R)
  assert(sub(gens G2, ring G1) == gens G1)
///

TEST ///
  -- cyclic6
  -- gb: 0.280165
  -- fglm: 0.885988
  debug needsPackage "FGLM"
  I = cyclic(6, MonomialOrder=>Lex)
  G1 = elapsedTime gb I
  I = cyclic(6)
  R = newRing(ring I, MonomialOrder=>Lex)
  G2 = elapsedTime fglm(I, R)
  assert(sub(gens G2, ring G1) == gens G1)
///

end--

-------------------------------------------------------------------------------
--- Development sections
-------------------------------------------------------------------------------

restart
uninstallPackage "FGLM"
restart
installPackage "FGLM"

restart
needsPackage "FGLM"
elapsedTime check FGLM -- ~3.2 seconds

viewHelp "FGLM"

-------------------------------------------------------------------------------
--- Profiling
-------------------------------------------------------------------------------
-- cyclic7
-- gb: 1354.44
-- fglm: 353.37 -> 37.76 -> 39.67 -> 25.83 (+ 0.48s gb in GRevLex w/ F4)
restart
gbTrace = 1
debug needsPackage "FGLM"
backSub = profile backSub
LUincremental  = profile LUincremental
multiplicationMatrices = profile multiplicationMatrices

--I1 = cyclic(7, MonomialOrder=>Lex)
--G1 = elapsedTime gb I1;
I2 = cyclic(7)
R2 = newRing(ring I2, MonomialOrder=>Lex)
G0 = elapsedTime forceGB(groebnerBasis(I2, Strategy => "F4"));
G2 = elapsedTime fglm(G0, R2);
--assert(gens G1 == gens G2)

profileSummary
--backSub: 35 times, used .279453 seconds
--LUincremental: 959 times, used 6.17558 seconds
--multiplicationMatrices: 1 times, used 1.99896 seconds

--T = new MutableList from toList(20:0);
--ticker := cpuTime(); ctr := 0;
--ctr = 1; T#ctr = T#ctr + cpuTime() - ticker; ticker = cpuTime();
peek T

-------------------------------------------------------------------------------
--- Longer tests
-------------------------------------------------------------------------------

-- cyclic7
-- gb: 1354.44
-- fglm: 353.37 -> 43.53 -> 31.84
restart
debug needsPackage "FGLM"
I = cyclic(7, MonomialOrder=>Lex)
G1 = elapsedTime gb I
I = cyclic(7)
R = newRing(ring I, MonomialOrder=>Lex)
G2 = elapsedTime fglm(I, R)


-- katsura7
-- gb: 6.78653
-- fglm: 0.779608 -> 0.325
restart
debug needsPackage "FGLM"
I = katsura(7, MonomialOrder=>Lex)
G1 = elapsedTime gb I
I = katsura(7)
R = newRing(ring I, MonomialOrder=>Lex)
G2 = elapsedTime fglm(I, R)


-- katsura8
-- gb: 2305.46
-- fglm: 8.23514 -> 1.571
restart
debug needsPackage "FGLM"
I = katsura(8, MonomialOrder=>Lex)
G1 = elapsedTime gb I
I = katsura(8)
R = newRing(ring I, MonomialOrder=>Lex)
G2 = elapsedTime fglm(I, R)


-- reimer5
-- gb: 8.3658
-- fglm: 3.79775 -> 1.064
restart
needsPackage "FGLM"
kk = ZZ/32003
R1 = kk[x,y,z,t,u, MonomialOrder=>Lex]
I1 = ideal(2*x^2 - 2*y^2 + 2*z^2 - 2*t^2 + 2*u^2 - 1,
           2*x^3 - 2*y^3 + 2*z^3 - 2*t^3 + 2*u^3 - 1,
           2*x^4 - 2*y^4 + 2*z^4 - 2*t^4 + 2*u^4 - 1,
           2*x^5 - 2*y^5 + 2*z^5 - 2*t^5 + 2*u^5 - 1,
           2*x^6 - 2*y^6 + 2*z^6 - 2*t^6 + 2*u^6 - 1)
G1 = elapsedTime gb I1
R2 = kk[x,y,z,t,u]
I2 = sub(I1, R2)
G2 = elapsedTime fglm(I2, R1)

-- virasoro
-- gb: 8.91079
-- fglm: 52.1752 -> 6.033 -> 4.67
restart
needsPackage "FGLM"
kk = ZZ/32003
R1 = kk[x1,x2,x3,x4,x5,x6,x7,x8, MonomialOrder=>Lex]
I1 = ideal(8*x1^2 + 8*x1*x2 + 8*x1*x3 + 2*x1*x4 + 2*x1*x5 + 2*x1*x6 + 2*x1*x7 - x1 - 8* x2*x3 - 2*x4*x7 - 2*x5*x6,
           8*x1*x2 - 8*x1*x3 + 8*x2^2 + 8*x2*x3 + 2*x2*x4 + 2*x2*x5 + 2*x2*x6 + 2*x2* x7 - x2 - 2*x4*x6 - 2*x5*x7,
	   -8*x1*x2 + 8*x1*x3 + 8*x2*x3 + 8*x3^2 + 2*x3*x4 + 2*x3*x5 + 2*x3*x6 + 2* x3*x7 - x3 - 2*x4*x5 - 2*x6*x7,
	   2*x1*x4 - 2*x1*x7 + 2*x2*x4 - 2*x2*x6 + 2*x3*x4 - 2*x3*x5 + 8*x4^2 + 8*x4* x5 + 2*x4*x6 + 2*x4*x7 + 6*x4*x8 - x4 - 6*x5*x8,
	   2*x1*x5 - 2*x1*x6 + 2*x2*x5 - 2*x2*x7 - 2*x3*x4 + 2*x3*x5 + 8*x4*x5 - 6*x4* x8 + 8*x5^2 + 2*x5*x6 + 2*x5*x7 + 6*x5*x8 - x5,
	   -2*x1*x5 + 2*x1*x6 - 2*x2*x4 + 2*x2*x6 + 2*x3*x6 - 2*x3*x7 + 2*x4*x6 + 2* x5*x6 + 8*x6^2 + 8*x6*x7 + 6*x6*x8 - x6 - 6*x7*x8,
	   -2*x1*x4 + 2*x1*x7 - 2*x2*x5 + 2*x2*x7 - 2*x3*x6 + 2*x3*x7 + 2*x4*x7 + 2* x5*x7 + 8*x6*x7 - 6*x6*x8 + 8*x7^2 + 6*x7*x8 - x7,
	   -6*x4*x5 + 6*x4*x8 + 6*x5*x8 - 6*x6*x7 + 6*x6*x8 + 6*x7*x8 + 8*x8^2 - x8)
G1 = elapsedTime gb I1;
R2 = kk[x1,x2,x3,x4,x5,x6,x7,x8];
I2 = sub(I1, R2);
G2 = elapsedTime fglm(I2, R1);

-- chemkin
-- gb: 5916.76
-- fglm: 4.46076 -> 1.611
restart
needsPackage "FGLM"
kk = ZZ/32003
R1 = kk[w,x3,x4,y2,y3,y4,y5,z2,z3,z4,z5, MonomialOrder=>Lex]
I1 = ideal(-4*w*y2 + 9*y2^2 + z2,
           x3^2 + y3^2 + z3^2 - 1,
           x4^2 + y4^2 + z4^2 - 1,
           9*y5^2 + 9*z5^2 - 8,
           -6*w*x3*y2 + 3*x3 + 3*y2*y3 + 3*z2*z3 - 1,
           3*x3*x4 + 3*y3*y4 + 3*z3*z4 - 1,
           x4 + 3*y4*y5 + 3*z4*z5 - 1,
           -6*w + 3*x3 + 3*x4 + 8,
           9*y2 + 9*y3 + 9*y4 + 9*y5 + 8,
           z2 + z3 + z4 + z5,
           w^2 - 2)
G1 = elapsedTime gb I1
R2 = kk[w,x3,x4,y2,y3,y4,y5,z2,z3,z4,z5]
I2 = sub(I1, R2)
G2 = elapsedTime fglm(I2, R1)

-- MES test
restart
needsPackage "FGLM"
kk = ZZ/32003
A = random(kk^3, kk^10)
(P,L,U) = LUdecomposition A
  Q = id_(target A) _ P
  Q*L*U == A

A = random(kk^10, kk^3)
(P,L,U) = LUdecomposition A
  Q = id_(target A) _ P
  Q*L*U == A

A = matrix"0,0,0,1;0,1,0,1;0,0,1,1" **kk
(P,L,U) = LUdecomposition A
  Q = id_(target A) _ P
  Q*L*U == A

A = transpose matrix"0,0,0,1;0,1,0,1;0,0,1,1" **kk
(P,L,U) = LUdecomposition A
  Q = id_(target A) _ P
  Q*L*U == A

A = transpose matrix"0,0,0,1;0,1,0,1;0,0,1,1" ** QQ
(P,L,U) = LUdecomposition A
  Q = id_(target A) _ P
  Q*L*U == A

A = transpose matrix"0,0,0,1;0,1,0,1;0,0,1,1" ** RR
(P,L,U) = LUdecomposition A
  Q = id_(target A) _ P
  Q*L*U == A

-*
-- add tests for LUincremental and backSub

n = 1000
kk = ZZ/32003
M = random(kk^n,kk^n)
LU = mutableMatrix map(kk^n,kk^n,0)

P = new MutableList from toList(0..n-1)
elapsedTime for i to n - 1 do { LUincr(P, LU, M_{i}, i); };

U = matrix for i to n - 1 list toList(i:0) | first entries submatrix(LU, {i}, {i..n-1});
L = matrix LU - U + id_(kk^n);
assert(0 == M - id_(kk^n)_(new List from P) * L * U)

elapsedTime (P',L',U') = LUdecomposition M;
id_(kk^n)_P' * L' * U' == M
assert(L == L')
assert(U == U')
assert(P' == new List from P)
*-

-- Singular test
LIB "poly.lib";
ring r=32003,(a,b,c,d,e,f,g),lp;
ideal I = cyclic(nvars(r));
timer=1;
int t = timer;
stdfglm(I);
t=timer-t;
t;

-- engine test
restart
debug needsPackage "FGLM"
setRandomSeed("hello")
n = 6
kk = ZZ/101
M = mutableMatrix id_(kk^n)--random(kk^n,kk^n)
LU = mutableMatrix map(kk^n,kk^(n+1),0)

P = (0..n-1)
P = new MutableList from P
elapsedTime for i to n - 1 do LUincremental(P, LU, M_{i}, i);
elapsedTime LUdecomposition M;
--rawLUincremental(P, raw LU, raw M_{0}, 0)
incrLU(P, LU, M_{0}, 0)
incrLU(P, LU, M_{1}, 1)
incrLU(P, LU, M_{2}, 2)


restart
debug needsPackage "FGLM"
R1 = ZZ/101[x,y,z]
I1 = ideal(x*y + z - x*z, x^2 - z, 2*x^3 - x^2*y*z - 1)
test(I1, Lex)

restart
debug Core
setRandomSeed("hello")

n = 3
kk= ZZ/17
U = mutableMatrix (random(kk^n,kk^(n+1), UpperTriangular => true) + (id_(kk^n) | random(kk^n,kk^1)))
x = mutableMatrix map(kk^n,kk^1,0)
v = matrix submatrix(U,{0..n-1}, {n})

rawTriangularSolve(raw U, raw x, 3, 3)

r = matrix submatrix(U,{0..n-1},{n})
assert(matrix submatrix(U,, {0..n-1}) * matrix x + r == v)