------------------------------------------------------------
  --
  -- ResZ(L)  -- Resolution of double points In P1 x P1 with ACM support
  --
  -- L is a partition describing ACM support
  -- These calculations are purely combinatorial And Do Not depend
  -- upon Grobner Basis Calculations
  --
  -- L:=[L1,L2,...,Ln] is a list with L1 >= L2 >= ...
  --
  -- The algorithm is based upon the paper:
  -- "The Minimal Resolutions of Double Points In P^1 x P^1
  --      with ACM support"
  -- by Elena Guardo And Adam Van Tuyl
  --
  -----------------------------------------------------------

Define ResZ(L)
  
  -- Step 0: Verify For Valid Input

  If (Check(L) = False) Then
    PrintLn "Input must be a list of nonincreasing integers";
    PrintLn "For example, L:=[5,4,4,2,1]";
    Return(0);
  End;

  -- Values used throughout program
  
  Lambda1:=L[1];
  R:=Len(L);
  DegMat:=NewMat(R,Lambda1,1);
  For I:=1 To R Do
    For J:=1 To Lambda1 Do
      If J <= L[I] Then DegMat[I,J]:=2; End;
    End;
  End;
  
  --Step 1: Find the bigraded Betti numbers of the completion Y

  SY0:=[[2*R,0],[R,Lambda1],[0,2*Lambda1]];
  For J:=2 To R Do
    If (L[J]-L[J-1] < 0) Then 
SY0:=Concat(SY0,[[J-1,Lambda1+L[J]],[J+R-1,L[J]]]); End;
  End;
  SY1:=[[2*R,L[R]],[R,Lambda1+L[R]]];
  For J:=2 To R Do
    If (L[J]-L[J-1] < 0) Then 
SY1:=Concat(SY1,[[J-1,Lambda1+L[J-1]],[J+R-1,L[J-1]]]); End;
  End;
  
  --Step 2: Locate the corners
  
  BC:=[];
  For I:=2 To R Do
    If (L[I]-L[I-1] <0) Then Append(BC,[I,L[I]+1]); End;
  End;
  SortBy(BC,Function("LexCompare"));
  C:=[];
  For I:=1 To (Len(BC)-1) Do
    For J:=(I+1) To Len(BC) Do
      Append(C,[BC[I][1],BC[J][2]]);
    End;
  End;
  C:=Concat(C,BC);
  SortBy(C,Function("LexCompare"));
  
  -- Step 3: Compute bigraded Betti numbers
  
  SZ0:=[];   -- Degrees of the generators
  SZ1:=[];   -- Degrees of the first syzygies
  SZ2:=[];   -- Degrees of the second syzygies
  
  For L:= 1 To Len(C) Do
    I:=C[L][1];
    J:=C[L][2];
    UIJ:=Sum([DegMat[K,J] | K In 1..(I-1)]);
    VIJ:=Sum([DegMat[I,K] | K In 1..(J-1)]);
    AIJ:=Sum([DegMat[K,J] | K In I..R]);
    BIJ:=Sum([DegMat[I,K] | K In J..Lambda1]);
    Append(SZ0,[UIJ,VIJ]);
    Append(SZ1,[UIJ+AIJ,VIJ]);
    Append(SZ1,[UIJ,VIJ+BIJ]);
    Append(SZ2,[UIJ+AIJ,VIJ+BIJ]);
    DegMat:=RewriteDegMat(DegMat,I,J);
  End;

  -- Step 4: Output the bigraded Betti numbers

  PrintLn "--Degrees of 0th Syzygies --";
  SZ0:=Concat(SZ0,SY0);
  SortBy(SZ0,Function("LexCompare"));
  PrintLn SZ0;
  PrintLn "--Degrees of 1st Syzygies --";
  SZ1:=Concat(SZ1,SY1);
  SortBy(SZ1,Function("LexCompare"));
  PrintLn SZ1;
  PrintLn "--Degrees of 2nd Syzygies --";
  SortBy(SZ2,Function("LexCompare"));
  PrintLn SZ2;

End; -- ResZ

-- Check(L) - Check verifies that the input is valid -------

Define Check(L)
  If (Type(L) <> LIST) Then Return(False); End;
  C:=True;
  For I:=1 To Len(L)-1 Do
    If L[I] < L[I+1] Then C:=False; End;
  End;
  Return(C);
End;
	     
-- RewriteDegMat - RewriteDegMat changes the DegMatrix at each step

Define RewriteDegMat(DegMat,I,J)
  A:=Len(DegMat);
  B:=Len(DegMat[1]);
  For M:=1 To A Do
    For N:=1 To B Do
      If (M >= I) And (N >= J) Then
	DegMat[M,N]:=0;
      End;
    End;
  End;
  Return(DegMat);
End; -- RewriteDegMat

-- LexCompare is a sorting function -----------------

Define LexCompare(S,T)
  If S[1] > T[1] Then Return(True);
   Else If (S[1]=T[1]) And (S[2] > T[2]) Then Return(True);
    Else Return(False);
   End;
  End;
End; -- LexCompare

ResZ([6,5,3,1]);

--Degrees of 0th Sygygies --
[[10, 0], [8, 1], [7, 3], [6, 5], [6, 2], [5, 6], [5, 4], [4, 6], [4, 6], 
[3, 8], [3, 7], [2, 10], [2, 9], [1, 11], [0, 12]]
--Degrees of 1st Sygygies --
[[10, 1], [8, 3], [8, 2], [7, 5], [7, 4], [6, 6], [6, 6], [6, 4], [5, 7], 
[5, 6], [5, 6], [4, 8], [4, 8], [4, 7], [3, 10], [3, 9], [3, 9], [2, 11], 
[2, 11], [1, 12]]
--Degrees of 2nd Sygygies --
[[8, 4], [7, 6], [6, 7], [5, 8], [4, 9], [3, 11]]

-- doublePts is a M2 implementation of the algorthim
-- described in the paper of Guardo-Van Tuyl (2006)

doublePts = l -> (
     sy0 = {(2*#l,0),(#l,first l),(0,2*first l)};
     sy1 = {(2*#l,last l),(#l,first l + last l)};
     C0 ={};
     for k from 1 to #l-1 do (
	  if (l_k - l_(k-1) < 0) then (
       	       sy0 = join (sy0, {(k,first l + l_k), (k+#l,l_k)});
       	       sy1 = join (sy1, {(k,first l + l_(k-1)),(k+#l,l_(k-1))});
	       C0 = join (C0, {(k+1,l_k+1)});
	       );	 
	  );	 
     C0 = rsort C0;
     C1 ={};
     for k1 from 0 to #C0-2 do(
      	   for k2 from k1+1 to #C0-1 do(
		C1 = join(C1, {((C0_k1)_0,(C0_k2)_1)});
		);     
	   );	
     C = rsort(join(C0,C1));
     M = {};
     for i from 0 to #l-1 do (
	  v = {};
     	  for j from 0 to first l -1  do (
	       if j < l_i then v = join(v,{2}) else v = join(v,{1});
	       );
	  M = join(M,{v});
	  );
     M = matrix M;
     sz0 = sy0;
     sz1 = sy1;
     sz2 = {};
     for t from 0 to #C -1 do (
	  i = (C_t)_0;
	  j = (C_t)_1;
	  uij =  sum toList apply(0..i-2, k -> M_(k,j-1));
	  vij =  sum toList apply(0..j-2, k -> M_(i-1,k));
	  aij =  sum toList apply((i-1)..(#l -1), k-> M_(k,j-1));
	  bij =  sum toList apply((j-1)..(first l-1), k-> M_(i-1,k));
	  sz0 = join(sz0, {(uij,vij)});
	  sz1 = join(sz1, {(uij+aij,vij),(uij,vij+bij)});
	  sz2 = join(sz2, {(uij+aij,vij+bij)});
	  
          MP={};
     	  for a from 0 to (#l-1) do (
	       v = {};
     	       for b from 0 to (first l -1)  do (
	       	    if ((a >= i-1) and (b >=j-1)) then v = join(v,{0}) 
                    else v= join(v,{M_(a,b)});
		    );
	       MP = join(MP,{v});	  
	       );	 
	  M = matrix MP;
	  );
     print "--Degrees of the 0th syzygies:--";
     print sort sz0;
     print "--Degrees of the 1st syzygies--";
     print sort sz1;
     print "--Degrees of the 2nd syzygies--";
     print sort sz2;
     );     

i2 : doublePts(6,5,3,1,1);

--Degrees of the 0th syzygies:--
{(0, 12), (1, 11), (2, 9), (2, 10), (3, 7), (3, 8), (4, 6), (4, 6), (5, 4), 
(5, 6), (6, 2), (6, 5), (7, 3), (8, 1), (10, 0)}
--Degrees of the 1st syzygies--
{(1, 12), (2, 11), (2, 11), (3, 9), (3, 9), (3, 10), (4, 7), (4, 8), (4, 8), 
(5, 6), (5, 6), (5, 7), (6, 4), (6, 6), (6, 6), (7, 4), (7, 5), (8, 2), 
(8, 3), (10, 1)}
--Degrees of the 2nd syzygies--
{(3, 11), (4, 9), (5, 8), (6, 7), (7, 6), (8, 4)}