-- MATH 531 -- Sept 27 --  Radical of a 0-dim. ideal
reduceUP = f -> ( -- reduce univariate polynomial
     -- (we assume that the input depends on one variable only)
     df := first select(1, first entries diff f, g->g!=0);
     f // gcd(f,df)
     )
intersect1 = (I,i) -> (
     E := QQ[x_1..x_(i-1), x_(i+1)..x_n, x_i, 
	  MonomialOrder=>Eliminate(n-1)];
     eI := substitute(I,E); 
     use ring I;
     substitute(first first entries selectInSubring(1, gens gb eI), ring I)
     )
radical0 = I -> (
     for i from 1 to n do (
	  p := intersect1(I,i);
	  I = I + ideal reduceUP p;
     	  );
     I
     )

///
restart
M2dir = "/home/x/leykin/M2/" 
load (M2dir|"lab4.m2")

-- test reduceUP
R = QQ[x_1,x_2];
f = x_2^2*(x_2-1)^5;
factor reduceUP f

-- test intersect1
n = 2;
R = QQ[x_1..x_n];
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); 
dim I
intersect1(I,1)
intersect1(I,2)     

-- test radical0
rI := radical0 I
-- compare the staircases (look at the leading monomials)
first entries gens gb I / leadMonomial 
first entries gens gb rI / leadMonomial
///