--author:Jessica Sidman
--6/24/09
--An implementation of the algorithm for computing the secant ideal of a monomial
--ideal as presented in "Combinatorial Secant Varieties" by Bernd Sturmfels and 
--Seth Sullivant.
--Note, I must be of class monomialIdeal.

--If I is a square-free monomial ideal, we can use
rSecantEdgeIdeal = (I,r) ->
(
  S:= ring I;
  v:= flatten entries vars S;
  n:=#v;
  m :=monomialIdeal apply(v, i -> i^2); 
  L:={};
  for i from 0 to n-1 do L= append(L, r); 
  mr:=monomialIdeal apply(v, i -> i^(r+1));
monomialIdeal( ( gens (mr: (dual I)^r) % mr) % m)
     )

secantMonomialIdeal = (I,r) ->
(
  S:= ring I;
  v:= flatten entries vars S;
  n:=#v;
  --We construct a vector a whose i-th component is max{rd_i-r+1,1}, where d_i
  --is the largest exponent of x_i appearing in a minimal generator of I.
  d:=flatten exponents lcm flatten entries gens I;
  a:=apply(d, i -> max(r*i-r+1, 1));
    --We go modulo the ideal m in the end.     
  L:={};
  for i from 0 to n-1 do (
       L = append(L, (v_i)^(a_i+1));
       )  ;
  m:=monomialIdeal L;
  monomialIdeal( gens dual( (dual(I, a))^r, r*a)  % m)
     )

--test ideal J

S = ZZ/32003[a..e]
I = monomialIdeal(a*b, a*c, a*d, b*c, b*d, c*d)
r=3
secantMonomialIdeal(I,3)