LIB "tst.lib";
tst_init();


//======================  Exercise 3.1 =============================
proc is_reg_sequence (ideal I)
"USAGE:   is_reg_sequence(I);  I ideal, 
RETURN:  1 if the given (ordered) list of generators for I is a 
           regular sequence;
         0 otherwise.
"
{
  int i;
  ideal J;
  while(i<size(I))
  {
    i++;
    if (size(reduce(quotient(J,I[i]),J))!=0) 
    { 
      return(0); 
    }
    J = groebner(J+I[i]);
  }
  if (size(reduce(1,J))==0) { return(0); }
  return(1);
}

ring R = 0, (x,y,z), dp;
ideal I = (x-1)*z, (x-1)*y, x;
is_reg_sequence (I); 
//-> 0
I = (x-1)*z, x, (x-1)*y;
is_reg_sequence (I);
//-> 1


kill R;
//======================  Exercise 3.2 =============================
if (not(defined(isCM))){ LIB "homolog.lib"; }
ring R1 = 0, (x,y,z), dp;
ideal I =  xy, yz, xz;
ring R1_loc = 0, (x,y,z), ds;
ideal I = imap(R1,I);  // ideal generated by I in localized ring
isCM(I); 

ring R2 = 0, (s,t,x,y,z,w), dp;
ideal I =  x-s4, y-s3t, z-st3, w-t4;
ideal IC = eliminate(I,st);
ring R2_loc = 0, (x,y,z,w), ds;
ideal IC = imap(R2,IC);
isCM(IC); 


kill R1,R1_loc,R2,R2_loc;
//======================  Exercise 3.3 =============================
proc truncate(module phi, int d)
"USAGE:   truncate(phi,d);  phi module, d int 
ASSUME:  phi comes assigned with an admissible degree vector as an
         attribute
RETURN:  module
NOTE:    Output is a presentation matrix for the truncation of 
         coker(phi) at d.
"
{
  if ( typeof(attrib(phi,"isHomog"))=="string" )
  {  
    ERROR("No admissible degree vector assigned");
  }
  else 
  {
    intvec v=attrib(phi,"isHomog");
  }
  int s = nrows(phi);
  int i,m,dummy;
  module L;
  for (i=1; i<=s; i++)
  {
    if (d>v[i]) 
    { 
      L = L+maxideal(d-v[i])*gen(i); 
    }
    else 
    { 
      L = L+gen(i); 
    }
  }
  L = modulo(L,phi);
  L = prune(L);
  if (size(L)==0) {return(L);}

  // it only remains to set the degrees for L:
  // ------------------------------------------
  m = v[1];
  for(i=2; i<=size(v); i++)
  {
    if(v[i]<m)
    {
      m = v[i];
    }
  }
  dummy = homog(L);
  intvec vv = attrib(L,"isHomog");
  if (d>m) 
  { 
    vv = vv+d-m; 
  }
  attrib(L,"isHomog",vv);
  return(L);
}

proc CM_regularity (module phi)
"USAGE:   CM_regularity(phi);    phi module
ASSUME:  phi comes assigned with an admissible degree vector as an
         attribute
RETURN:  integer 
NOTE:    Output is the Castelnuovo-Mumford regularity of coker(phi).
"
{
  if ( typeof(attrib(phi,"isHomog"))=="string" )
  {  
    ERROR("No admissible degree vector assigned");
  }
  def L = mres(phi,0);
  intmat BeL = betti(L);
  int r = nrows(module(matrix(BeL)));
  int shift = attrib(BeL,"rowShift");  // See Section 3.4
  return(r+shift-1);
}

ring R = 0, (w,x,y,z), dp;
module I = [xz,0,-w,-1,0], [-yz2,y2, 0,-w,0], [y2z,0,-z2,0,-x],
           [y3,0,-yz,-x,0], [-z3,yz,0,0,-w], [-yz2,y2,0,-w,0],
           [0,0,-wy2+xz2,-y2,x2];  
homog(I);
CM_regularity(I);
//->  3 
def T2 = mres(truncate(I,2),0);
print (betti(T2),"betti");
//->             0     1     2     3
//->  ------------------------------
//->      2:    19    36    23     6
//->      3:     -     -     1     -
//->  ------------------------------
//->  total:    19    36    24     6

def T3 = mres(truncate(I,3),0);
print (betti(T3),"betti");
//->             0     1     2     3
//->  ------------------------------
//->      3:    40    91    71    19
//->  ------------------------------
//->  total:    40    91    71    19


kill R;
//======================  Exercise 3.4 =============================
//===== Procedures are stored in the library file sol_3_4.lib ======
//==================================================================
LIB "sol_3_4.lib";
example ker_Mod;
example Ext_Mod;


//======================  Exercise 3.5 =============================
ring S = 32003, x(0..4), dp;
resolution kos = nres(maxideal(1),0);
print(betti(kos),"betti");
matrix alpha0 = random(32002,10,5);
matrix pres = module(alpha0)+kos[4];
matrix dir = transpose(pres);
resolution fdir = mres(dir,2);
print(betti(fdir),"betti");
if (not(defined(flatten))) { LIB "matrix.lib"; }
ideal I = flatten(fdir[2]);
resolution FI = mres(I,0);
print(betti(FI),"betti");

ring S1 = 32003, x(0..4), dp;
resolution kos = nres(maxideal(1),0);
betti(kos);
matrix gammatilde = random(32002,20,19);
matrix kos1 = matrix(kos[1]);
matrix kos2 = kos[2];
if (not(defined(dsum))){ LIB"matrix.lib"; }
matrix kos2pluskos1pluskos1 = dsum(kos2,kos1,kos1);
module delta = kos2pluskos1pluskos1*gammatilde;
attrib(delta,"isHomog",intvec(-1,-1,-1,-1,-1,-1,-1));
resolution fdelta = mres(delta,0);
print(betti(fdelta),"betti"); 
//->             0     1     2     3     4     5
//->  ------------------------------------------
//->     -1:     7    19    25    15     3     -
//->      0:     -     -     -     2     3     1
//->  ------------------------------------------
//->  total:     7    19    25    17     6     1

matrix psi = matrix(fdelta[3]);
matrix talpha1 = random(32002,3,15);
matrix zero[3][2];
talpha1 = concat(talpha1,zero);
matrix kos5 = kos[5];
matrix tphi = transpose(dsum(kos5,kos5,kos5));
matrix talpha1tilde = talpha1*transpose(psi);
matrix talpha0 = lift(tphi,talpha1tilde);

matrix dir = transpose(concat(psi,transpose(talpha0)));
resolution fdir = mres(dir,2);
print(betti(fdir),"betti");
ideal I = groebner(flatten(fdir[2]));
resolution FI = mres(I,0);
print(betti(FI),"betti");
// ---------- Check Smoothness ------------
int codimI = nvars(S1)-dim(I);
nvars(S1) - dim(groebner(minor(jacob(I),codimI) + I));
//-> 5
// ----------------------------------------

tst_status(1);$