*** Gröbner bases conversion***
----------------------------------------------
-- 
----------------------------------------------
-- For a 0-dimensional ideal an FGLM-style bases conversion algorithm relies on linear algebra. 
--
-- * Subroutines
--
-- * Main algorithm
----------------------------------------------
-- Subroutines
----------------------------------------------
-- For our running example $k=\QQ$,
R = QQ[x,y];
I = ideal(y^4*x+3*x^3-y^4-3*x^2, 
     x^2*y-2*x^2,
     2*y^4*x-x^3-2*y^4+x^2); 
-- compute a GB using relative to the default (GRevLex) monomial order
getGB = flatten@@entries@@gens@@gb
G1 = getGB I
-- A function that finds the border of a give monomial set:
border = B -> (
    R := ring matrix {B};
    B' := set {}; 
    for b in B do scan(numgens R, i ->(
	    b' := b*R_i;
	    if not member(b',B) then B'=B'+set{b'};
	    ));
    toList B'
    )
-- One of the main subroutines needs to find the next unprocessed monomial:
nextMonomial = (B,G) -> (
    C := select(border B, b'->all(G, g-> b' % leadMonomial g != 0));
    if #C>0 then min C else null
    )
-- Another one computes a polynomial $f$ as a $k$-linear combination
-- of the given list of polynomials:
linearCombination = (f,Bbar) ->  (
    M := last coefficients matrix {{f}|Bbar};
    Mf := M_{0};
    M'Bbar := submatrix'(M,{0});
    c := Mf //  M'Bbar;
    if M'Bbar*c == Mf 
    then flatten entries sub(c,coefficientRing ring f)
    else null
    )
-- Lets check the routine:
linearCombination(G1#0+3*G1#2,G1)
linearCombination(x^3+y, G1) -- returns null

----------------------------------------------
-- Main algorithm
----------------------------------------------
FGLM = (G1,mo) -> (
    I1 := ideal G1; R1 := ring I1;
    if dim I1 != 0 then error "expected a 0-dimensional ideal";
    R2 := newRing(R1,MonomialOrder=>mo); 
    m := 1_R2;
    G2 := {}; B2 := {m}; B2bar := {1_R1};
    while (m = nextMonomial(B2,G2)) =!= null do ( 
	f := sub(m,R1) % I;
	if (c:=linearCombination(f,B2bar)) === null 
	then (B2 = B2 | {m}; B2bar = B2bar | {f})
	else G2 = G2 | {m - sum(#B2, i->c#i*B2#i)}
	);
    G2
    )
-- For the running example,
G1 = getGB I
-- we compute 
FGLM(G1,Lex) 
-- and compare with the built the results of the built-in method
Rlex = newRing(R,MonomialOrder=>Lex)
getGB sub(I,Rlex) 
-- Similarly, 
for i to 5 do (
    mo = {Weights=>{i,1}};
    << "w=(" << i << ",1)    FGLM: " << FGLM(G1,mo) << endl;
    Rw = newRing(R,MonomialOrder=>mo);
    << "w=(" << i << ",1)    M2gb: " << getGB sub(I,Rw) << endl;
    print "---------------------------------------------------";
    )
-- shows the "evolution" of the bases as the weight monomial order is varied.