-- MATH 531 -- Oct 4 --  FGLM, Grobner Fan
find0lin = X -> ( -- attempts to find a linear combination of polynomials in X 
     S := syz(matrix{X}); 
     if S != 0 then (
	  r := select(1, entries transpose S, 
	       s->all(0..#X-1, k->s#k==(s#k)_{0}) and s#(#X-1)!=0);
	  if #r != 0 then return first r;
	  );
     null
     ) 
intersect1 = (I,i) -> ( -- a version of intersect1 that utilizes linear algebra 
     R := ring I;
     xi := R_(i-1);
     X := {}; -- store the remainders of x_i^j 
     j := 0;
     while true do (
     	  X = X | {xi^j % I}; -- the remainder of x_i^j is stored
	  L := find0lin X;
	  -- if a linear combination of the remainders is found 
	  -- then the corresponding linear combination of x_i^j is in I
	  if L =!= null 
	  then return sum(0..#L-1, k->L#k * xi^k); 
	  j = j + 1;
	  )
     )

FGLM = (G,T) -> (
     -- redefine the comparison operator in order to go around a bug in M2
     -- (it is used by sort)
     T ? T := (a,b)->( 
     if leadTerm(a+b) == a 
     then symbol >
     else if leadTerm(a+b) == b
     then symbol <
     else symbol ==
     );  
     R := ring G; n := numgens R;
     -- m is a vector of the dimensions of the "bounding box"
     m := apply(1..n, i->first degree intersect1(ideal G, i));
     -- B is the sorted list of all monomials in the box
     B := sort toList apply((n:0)..m, i->T_(toList i));
     << "sorted list of monomials in the box:" << endl << B << endl;
     H := {}; -- the new G.b. 
     S := {}; -- the new standard monomials
     remS := {}; -- the remainders of the new std.mon's w.r.t. G in the old ordering
     scan(B, m->(
	       -- we don't need to consider monomials 
	       -- that are above the (partially known) staircase
	       if any(H, h-> m % leadMonomial h == 0) then return;
	       -- compute the remainder of m w.r.t. G 
	       rem := sub(m,R)%G; 
	       L := find0lin(remS|{rem});
	       -- if the remainder of m does not depend linearly 
	       -- on the remainders of the (already known) std.monomials 
	       -- then it is standard   
	       if L === null  then ( 
	       	    S = S | {m};
		    remS = remS | {rem};
		    ) 
	       -- otherwise, we found a corner of the staircase 
	       -- (with the corresponding element of H)
	       else (
		    H = H | {sum(0..#L-1, k -> substitute(L#k,T) * (S|{m})#k)};
	       	    )
	       ));
     H
     ) 
///
restart
M2dir = "/home/x/leykin/M2/" 
load (M2dir|"lab5.m2")

R = QQ[x_1,x_2];
I = ideal(x_2^4*x_1+3*x_1^3-x_2^4-3*x_1^2, 
     x_1^2*x_2-2*x_1^2,
     2*x_2^4*x_1-x_1^3-2*x_2^4+x_1^2); 

--test intersect
intersect1(I,1)
intersect1(I,2)     

G = gens gb I
T = QQ[x_1,x_2, MonomialOrder=>Lex]
time FGLM(G,T) 
transpose gens gb sub(I,T)

-- explore Grobner bases for various weights  
R = QQ[x_1,x_2, Weights=>{100,1}]
transpose gens gb sub(I,R)
R = QQ[x_1,x_2, Weights=>{1,100}]
transpose gens gb sub(I,R)

R = QQ[x_1..x_3]
I = ideal {x_3^2-x_1+x_2-1, x_1^2-x_2*x_3+x_1, x_2^3-x_1*x_3+2} 
transpose (G = gens gb I)
transpose gens gb substitute(I,QQ[x_1..x_3,Weights=>{2,1,5}])
transpose gens gb substitute(I,QQ[x_1..x_3,MonomialOrder=>Lex])
///