-- What is the following code doing?
-- You can load it into Macaulay 2 by saying 
-- load "hw.m2"
-- and run the subroutines yourself to figure out what they do.

syze = 2;
-- If you're willing to wait a while (overnight?), try changing this to syze=3

R = QQ[a_(1,1)..a_(syze,syze),b_(1,1)..b_(syze,syze)];

A = transpose genericMatrix(R,a_(1,1),syze,syze);
B = transpose genericMatrix(R,b_(1,1),syze,syze);

dp = M -> matrix apply(syze,i->apply(syze,j->(
	       if (i==j) then M_(i,j) else 0)));
lp = M -> matrix apply(syze,i->apply(syze,j->(
	       if (i>j) then M_(i,j) else 0)));

dec = I -> (
     print "old:";
     scan(flatten entries gens trim I, print); print "";
     cs = decompose I;
     scan(#cs, i->(print ("new in #" | toString(i+1) | ":");
	       scan(select(flatten entries gens cs_i, g->(g%I != 0)), print);
	       print ""; ))
     )

I1 = ideal {lp(A*B), lp(B*A)}; -- what are these equations, in words?
dec(I1); -- what does this do?

C = A*B-B*A;
C = C - dp(C);
I2 = ideal C;  -- what are these equations, in words?
dec(I2);