supportSet=method()
supportSet(Module,Module,Symbol):=(M,N,y)->(
     --Input = Modules M and N over a quotient ring 
     --and an unassigned symbol.
     --Output = An ideal defining the support set of M.
     --Method:  Follows the algorithm of David Jorgensen.  
     --After constructing a hypersurface, using new variables 
     --and the generators of the defining ideal, we compute 
     --the appropriate Exts and take the intersection of their 
     --annihlators. Finally we lift the ideal to the ring it 
     --should live in.
     R  := ring(M);                         
     n1 := numgens(R);                 
     I := ideal(R);     
     n2 := numgens(I);	
     --Constructing the Hypersurface.  We need to make 
     --sure Hyp is a homogeneous elt.     	  
     Degs:=flatten degrees source gens I; 
     maxDeg := 1+ max Degs;            
     newDegs := apply(Degs,d->maxDeg-d);     
     S := coefficientRing(R)[gens(R),y_1..y_(n2),
	  MonomialOrder=>ProductOrder{n1,n2},
	  Degrees=>join(flatten (monoid R).Options.Degrees,newDegs)];
     --S is the new poly ring, it has the same vars of R plus the new ones.
     IS := substitute(I,S);     	       	  
     Hyp:=ideal(sum(apply(n2,i->S_(n1+i)*IS_i))); 
     T := S/Hyp;	       	    	
     --Get presentations of M and N over the Hypersurface.      	   
     MT := substitute(presentation M,T);
     NT := substitute(presentation N,T);     	       	   
     IT := substitute(gens IS,T);	            	   
     M1 := (coker MT)**coker(IT);
     N1 := (coker NT)**coker(IT);	
     L := coker(matrix{apply(n1,i->T_i)}); 
     --Computing the Exts and their annihilators.
     E1 := prune(Tor_(n1)(M1,N1) ** L);   
     A1 := newAnn E1;
     E2 := prune(Tor_(n1+1)(M1,N1) ** L);
     A2 := newAnn E2;
     A := intersect(A1,A2);
     --Lifting the support variety into the correct ring.
     AA := lift(A,S);
     J := ideal(selectInSubring(n1,gens AA));
     U := coefficientRing(R)[T_(n1)..T_(n2+n1-1),Degrees=>newDegs];
     intersect(substitute(J,U),ideal(vars U))
     )

newAnn = (M) -> (
     --Input = A module.
     --Output = the annihilator of that module.
     --Method = rather than computing (source P):(target P) all at once
     --we compute it one column at a time.  
     if M == 0 then ideal(1_(ring M))
               else (
     P:=presentation M;
     F:=target P;
     intersect(apply(numgens(F),i->(
       m1:=matrix{F_i};
	  I:=ideal(modulo(m1,P));
	  I))))
     )


TEST///
restart
load "supportVarietyModule.m2"
R=QQ[a,b,c]/ideal(a^2,b^2,c^2)
M= coker matrix{{0,0,0,0,0,0,-b,0,0,0,0,0},{0,0,0,0,0,0,a,b,0,0,0,0},{0,0,0,0,-b,0,c,0,0,0,0,0},{0,0,0,0,a,b,0,c,0,0,0,0},{0,0,-b,a,0,0,0,0,0,0,0,0},{0,0,a,0,0,0,0,0,0,0,0,0},{0,a,c,0,0,0,0,0,0,0,0,0},{a,b,0,c,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,-c,0,0,-b},{0,0,0,0,0,0,0,0,0,-c,0,a},{0,0,0,0,0,0,0,0,a,b,c,0},{0,0,0,0,0,0,0,0,0,0,0,c}}
time s=supportVariety(M,symbol y)

R=ZZ/32003[a,b,c,d]/ideal(a*b-c*d,a*c-b*d,a*d-b*c)
M=coker matrix{{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-c,0,-b,-c},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-c,a,d},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,a,b,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-c,0,-b,-b,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-c,a,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,a,b,0,d,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,-c,0,-b,-d,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,-c,a,c,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,a,b,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,0,-d,b,d,0,-b,d,0,-d,0,-a,-d},{0,0,0,0,0,0,0,0,0,0,-b,-c,0,0,0,0,0,0,0,0,d,0,0,0},{0,0,0,0,0,0,0,0,0,0,a,d,0,0,0,0,0,0,0,0,0,d,0,0},{0,-d,0,0,0,0,-b,0,0,-b,0,0,0,0,0,0,d,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,-b,-d,0,0,0,d,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,a,c,0,0,0,0,d,0,0,0,0,0,0,0,0,0,0},{0,0,-d,0,0,0,a,0,0,0,0,0,0,0,0,0,0,d,0,0,0,0,0,0},{-d,0,0,0,0,-c,0,0,0,0,c,0,0,0,0,0,0,0,0,0,0,0,d,0},{0,0,-b,0,0,0,c,0,0,0,0,0,0,0,0,0,0,0,d,0,0,0,0,0},{0,  0,  0,  -c, -d, 0,  0,  c,  0,  0,  0,  0,  0,  0,  d,  0,  0,  0,  0,  0,  0,  0,  0,  0},{0,  0,  0,  a,  b,  0,  0,  0,  c,  0,  0,  0,  0,  0,  0,  d,  0,  0,  0,  0,  0,  0,  0,  0},{0,  a,  b,  0,  0,  0,  0,  0,  0,  c,  0,  0,  0,  0,  0,  0,  0,  0,  0,  d,  0,  0,  0,  0},{a,  0,  0,  0,  0,  b,  0,  0,  0,  0,  0,  c,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  d}}  
time s=supportVariety(M,symbol y)

R=ZZ/32003[a,b,c,d]/ideal(a*b-c*d,a*c-b*d,a^3-c^3) 
M=coker matrix{{a+c}}                 
time s=supportVariety(M,symbol y)

Q=QQ[a,b,c]
I=ideal(a^2,b^2,c^2)
R=Q/I
M= coker matrix{{-c,0,0,0,0,0,0,0,0,-b,0,a},
               {0,-c,0,0,0,0,0,0,0,a,b,0},
               {0,0,-c,0,0,0,-b,0,a,0,0,0},
               {0,0,0,-c,0,0,a,0,0,0,0,b},
               {0,0,0,0,-c,0,0,-b,0,0,a,0},
               {0,0,0,0,0,-c,0,a,b,0,0,0},
               {-b,0,0,a,0,0,0,c,0,0,0,0},
               {0,0,-b,0,0,a,0,0,0,c,0,0},
               {0,0,a,b,0,0,0,0,0,0,c,0},
               {0,a,0,0,-b,0,c,0,0,0,0,0},
               {0,0,0,0,a,b,0,0,0,0,0,c},
               {a,b,0,0,0,0,0,0,c,0,0,0}}
time s=supportVariety(M,symbol y)

R=ZZ/32003[a,b,c,d]/ideal(a^2,b^2,c^2,d^2)
M=coker value get "bigguy.m2";
time s=supportVariety(M,symbol y)

R=QQ[a,b,c,Degrees=>{4,5,6}]/ideal(b^2-a*c,a^3-c^2)
isHomogeneous R
M=coker matrix{{a,b}}
time s=supportVariety(M,symbol y)
///