-- file name: intersection.m2
-- intersect needs a degree limit

intersection = method(Dispatch => Thing, TypicalValue=>Ideal, Options => {DegreeLimit=>{}})

intersection Sequence := 
intersection List := o -> L -> (
     if not all(L, x -> instance(x,Ideal))
     then error "expected a list of ideals";
     B := directSum apply(L, generators);
     A := map(target B, 1, (i,j) -> 1);
     ideal syz(A|B, SyzygyRows => 1, DegreeLimit=>o.DegreeLimit)
     )

inter = (d, L) -> (
     -- L is a list of ideals
     -- d an integer
     -- compute the intersection of the ideals in degrees <= d.
     ans := L_0;
     for i from 1 to #L-1 do time ans = intersection(ans,L_i,DegreeLimit=>d);
     ans)

intersectSpaces = (d,L) -> (
     -- L is a list of ideals
     -- d an integer
     -- return the ideal generated by the intersection of the degree d parts of
     --  each ideal.
     -- Assumption: each ideal in L is generated in degree d.
     S := ring L_0;
     k := coefficientRing S;
     monoms := basis(d,S);
     time L = apply(L, I -> (
	       (mn,cf) := coefficients(gens I, Monomials=>monoms);
	       image substitute(cf,k)));
     time M := intersect L;
     ideal(monoms * (gens M))
     )
-- what about intersection of subspaces of polynomials of the same degree?
-- The above will do more work that required, since it will make
-- spairs that are not used...
end
restart
load "intersection.m2"

R = ZZ/32003[a..d]
I = (ideal random(R^1, R^{-2,-3,-5}))
J = (ideal random(R^1, R^{-2,-3,-6}))
K = ideal random(R^1, R^{-2,-3,-8})
time intersect(I,J,K);
betti oo
time intersection(I,J,K,DegreeLimit=>6);
time intersection(I,J,K,DegreeLimit=>7);
time intersection(I,J,K);