LIB "tst.lib";
tst_init();


//======================  Example 8.6 =============================
LIB "finvar.lib";        
ring R = (0,a), (x(0..4)), dp;
minpoly = a4+a3+a2+a+1;  // need fifth roots of unity
matrix Si[5][5] = 0,0,0,0,1,
                  1,0,0,0,0,
                  0,1,0,0,0,
                  0,0,1,0,0,
                  0,0,0,1,0;
matrix Ta[5][5];
Ta[1,1] = 1;  Ta[2,2] = a;  Ta[3,3] = a2;  Ta[4,4] = a3;  
Ta[5,5] = a4;
int aa = timer;          // time in seconds
matrix P,S,IS = invariant_ring(Si,Ta,intvec(0,0,1));

//->   Generating the entire matrix group. Whenever a new group element is
//->   found, the corresponding ring homomorphism of the Reynolds operator
//->   and the corresponding term of the Molien series is generated.

//->   Group element 3 has been found.
//->                        [...]
//->   Group element 125 has been found.

//->   Now we are done calculating Molien series and Reynolds operator.

//->   We can start looking for primary invariants...

//->   Computing primary invariants in degree 5:
//->   We find: x(0)*x(1)*x(2)*x(3)*x(4)
//->   We find: x(0)^3*x(2)*x(3)+x(0)*x(1)*x(3)^3+x(0)*x(2)^3*x(4)+[...]
//->   We find: x(0)*x(1)^3*x(2)+x(1)*x(2)^3*x(3)+x(0)^3*x(1)*x(4)+[...]
//->   Computing primary invariants in degree 10:
//->   We find: x(0)^10+x(1)^10+x(2)^10+x(3)^10+x(4)^10
//->   We find: -x(0)^5*x(1)^5+x(0)^5*x(2)^5-[...]

//->   We found all primary invariants.

//->   Polynomial telling us where to look for secondary invariants:
//->    x(0)^30+3*x(0)^25+24*x(0)^20+44*x(0)^15+24*x(0)^10+3*x(0)^5+1

//->   In degree 0 we have: 1

//->   Searching in degree 5, we need to find 3 invariant(s)...
//->   We find: x(0)^2*x(1)^2*x(3)+x(0)*x(2)^2*x(3)^2+[...]
//->   We find: x(0)^2*x(1)*x(2)^2+x(1)^2*x(2)*x(3)^2+[...]
//->   We find: x(0)^5+x(1)^5+x(2)^5+x(3)^5+x(4)^5 .

//->   Searching in degree 10, we need to find 24 invariant(s)...
//->                        [...]   
//->   Searching in degree 15, we need to find 44 invariant(s)...
//->                        [...]   
//->   Searching in degree 20, we need to find 24 invariant(s)...
//->                        [...]   
//->   Searching in degree 25, we need to find 3 invariant(s)...
//->                        [...]
//->   Searching in degree 30, we need to find 1 invariant(s)...
//->   We find: x(0)^15*x(1)^10*x(3)^5+x(0)^10*x(1)^15*x(3)^5+[...]

//->   We're done!

ideal HMQ = invariant_basis(5,Si,Ta);
print(HMQ);
//->   x(0)*x(1)*x(2)*x(3)*x(4),
//->   x(0)^2*x(1)*x(2)^2+x(1)^2*x(2)*x(3)^2+x(0)^2*x(3)^2*x(4)+[...]
//->   x(0)^2*x(1)^2*x(3)+x(0)*x(2)^2*x(3)^2+x(1)^2*x(2)^2*x(4)+[...]
//->   x(0)^3*x(2)*x(3)+x(0)*x(1)*x(3)^3+x(0)*x(2)^3*x(4)+[...]
//->   x(0)*x(1)^3*x(2)+x(1)*x(2)^3*x(3)+x(0)^3*x(1)*x(4)+[...]
//->   x(0)^5+x(1)^5+x(2)^5+x(3)^5+x(4)^5


kill R,aa;
//==================  Remark 8.7(new Session) =========================
if (not(defined(invariant_ring))){ LIB "finvar.lib"; }
ring R = 101, (x(0..4)), dp;
matrix Si[5][5] = 0,0,0,0,1,
                  1,0,0,0,0,
                  0,1,0,0,0,
                  0,0,1,0,0,
                  0,0,0,1,0;
number a = 36;              // primitive fifth root of unity
matrix Ta[5][5];
Ta[1,1] = 1;  Ta[2,2] = a;  Ta[3,3] = a^2;  Ta[4,4] = a^3;
Ta[5,5] = a^4;
int aa = timer;             // time in seconds
matrix P,S,IS = invariant_ring(Si,Ta,intvec(0,0,0));
size(S);
//->   100
print(S[100]);
//->   [x(0)^15*x(1)^10*x(3)^5+x(0)^10*x(1)^15*x(3)^5+[...]
ideal HMQ = invariant_basis(5,Si,Ta);

//==================  Example 8.8 (continued Session) ===================
ring P4 = 101, (x(0..4)), dp;
ideal HMQ = fetch(R,HMQ);
ideal CI = HMQ[1], HMQ[2];
list DEG = primdecGTZ(CI);
size(DEG);
//->   10
ideal I6 = DEG[6][1]; I6;
//->   I6[1]=x(1)^2*x(3)+x(2)*x(4)^2
//->   I6[2]=x(0)
degree(std(I6));
//->   // dimension (proj.)  = 2
//->   // degree (proj.)   = 3
ideal I10 = DEG[10][1]; I10;
//->   I10[1]=x(2)^2
//->   I10[2]=x(2)*x(3)^2+x(0)*x(4)^2
//->   I10[3]=x(0)*x(2)
//->   I10[4]=x(0)^2
degree(std(I10));
//->   // dimension (proj.)  = 2
//->   // degree (proj.)   = 2
ideal IX = intersect(DEG[3][1],DEG[4][1],DEG[8][1],DEG[9][1], 
                     DEG[10][1]);
degree(std(IX));
//->   // dimension (proj.)  = 2
//->   // degree (proj.)   = 10
resolution FIX = mres(IX,0);
print(betti(FIX),"betti");
//->              0     1     2     3     4
//->   ------------------------------------
//->       0:     1     -     -     -     -
//->       1:     -     -     -     -     -
//->       2:     -     -     -     -     -
//->       3:     -     -     -     -     -
//->       4:     -     3     -     -     -
//->       5:     -    15    35    20     -
//->       6:     -     -     -     -     2
//->   ------------------------------------
//->   total:     1    18    35    20     2

module N = transpose(FIX[3]);
homog(N);
//->   1
intvec deg_N = attrib(N,"isHomog");
attrib(N,"isHomog",deg_N-3);        // set degrees
resolution FN = mres(N,0);
print(betti(FN),"betti");
//->              0     1     2     3     4     5
//->   ------------------------------------------
//->      -3:    20    35    15     -     -     -
//->      -2:     -     -     4     -     -     -
//->      -1:     -     -     -     -     -     -
//->       0:     -     -     5    15    10     2
//->   ------------------------------------------
//->   total:    20    35    24    15    10     2
matrix NN = FN[2];
matrix PRESMHM[35][19] = NN[1..35,1..19];
PRESMHM = transpose(PRESMHM);
resolution FMHM = mres(PRESMHM,0);
print(betti(FMHM),"betti");
//->              0     1     2     3
//->   ------------------------------
//->       0:     4     -     -     -
//->       1:    15    35    20     -
//->       2:     -     -     -     2
//->   ------------------------------
//->   total:    19    35    20     2

matrix zero[1][15];
matrix ran = random(100,1,4);
matrix psi = transpose(concat(zero,ran));
matrix pres = PRESMHM + module(psi);
module dir = transpose(pres);
resolution fdir = mres(dir,2);
print(betti(fdir),"betti");
//->              0     1     2
//->   ------------------------
//->       0:    35    15     -
//->       1:     -     3     -
//->       2:     -     -     -
//->       3:     -     -     -
//->       4:     -     -     -
//->       5:     -     -     1
//->   ------------------------
//->   total:    35    18     1
ideal IA = groebner(flatten(fdir[2]));
int codimIA = nvars(P4) - dim(IA);
ideal sIA = minor(jacob(IA),codimIA)+IA;
nvars(P4) - dim(groebner(sIA));
//->   5

matrix dummy[1][3] = IA[1..3];     // the 3 quintics in IA
ideal CI2 = dummy*random(100,3,2);
ideal IA' = sat(CI2,IA)[1];
resolution FIA' = mres(IA',0);
print(betti(FIA'),"betti");
//->              0     1     2     3     4
//->   ------------------------------------
//->       0:     1     -     -     -     -
//->       1:     -     -     -     -     -
//->       2:     -     -     -     -     -
//->       3:     -     -     -     -     -
//->       4:     -     3     -     -     -
//->       5:     -     -     -     -     -
//->       6:     -     5    15    10     2
//->   ------------------------------------
//->   total:     1     8    15    10     2
int codimIA' = nvars(P4) - dim(IA');
ideal sIA' = minor(jacob(IA'),codimIA')+IA';
nvars(P4) - dim(groebner(sIA'));
//->   5

ideal QA = IA[1..3];
ideal HMlines = sat(QA,IA)[1];     // result is a Groebner basis
degree(HMlines);
//->   // dimension (proj.)  = 1
//->   // degree (proj.)   = 25


kill R,aa,P4,codimIA,codimIA',deg_N;
//================== Example 8.9 (new Session) =========================
ring R = 2, (x(1..4)), dp;
matrix A[4][4];
A[1,4] = 1;  A[2,1] = 1;  A[3,2] = 1;  A[4,3] = 1;
print(A);
//->   0,0,0,1,
//->   1,0,0,0,
//->   0,1,0,0,
//->   0,0,1,0 
if (not(defined(invariant_ring))){ LIB "finvar.lib"; }
matrix P,S = invariant_ring(A);
P;
//->   P[1,1]=x(1)+x(2)+x(3)+x(4)
//->   P[1,2]=x(1)*x(3)+x(2)*x(4)
//->   P[1,3]=x(1)*x(2)+x(2)*x(3)+x(1)*x(4)+x(3)*x(4)
//->   P[1,4]=x(1)*x(2)*x(3)*x(4)
size(S);
//->   5


kill R;
//================== Example 8.10 (new Session) =========================
if (not(defined(invariant_ring))){ LIB "finvar.lib"; }
ring R = (0,a), (x(0..3)), dp;
minpoly = a4+a3+a2+a+1;
matrix A[4][4];
A[1,1] = a;  A[2,2] = a2;  A[3,3] = a3;  A[4,4] = a4;
matrix P,S,IS = invariant_ring(A,intvec(0,0,0));
size(P);
//->   4
size(S);
//->   12

proc min_generating_set (matrix P,S)
"USAGE:  min_generating_set(P,S);   P,S matrix
ASSUME: The entries of P,S are homogeneous and ordered by ascending 
        degrees. The first entry of S equals 1. (As satisfied by 
        the first two output matrices of invariant_ring(G).)
RETURN: ideal
NOTE:   The given generators for the output ideal form a minimal 
        generating set for the ring generated by the entries of 
        P,S. The generators are homogeneous and ordered by 
        descending degrees.
"
{
  if (defined(flatten)==0) { LIB "matrix.lib"; }
  ideal I1,I2 = flatten(P),flatten(S); 
  int i1,i2 = size(I1),size(I2);
  // We order the generators by descending degrees
  // (the first generator 1 of I2 is omitted):
  int i,j,s = i1,i2,i1+i2-1;
  ideal I;
  for (int k=1; k<=s; k++)
  { 
    if (i==0) { I[k]=I2[j]; j--; } 
    else
    { 
      if (j==0) { I[k]=I1[i]; i--; } 
      else 
      {
        if (deg(I1[i])>deg(I2[j])) { I[k]=I1[i]; i--; }
        else { I[k]=I2[j]; j--; }
      }
    }
  } 
  intvec deg_I = deg(I[1..s]);
  int n = nvars(basering);
  def BR = basering;

  // Create a new ring with elimination order:
  //---------------------------------------------------------------
  // ****    this part uses the command ringlist which is      ****
  // ****     only available in SINGULAR-3-0-0 or newer        ****
  //---------------------------------------------------------------
  list rData = ringlist(BR);
  intvec wDp;
  for (k=1; k<=n; k++) { 
    rData[2][k] ="x("+string(k)+ ")"; 
    wDp[k]=1;
  }
  for (k=1; k<=s; k++) { rData[2][n+k] ="y("+string(k)+ ")"; } 
  rData[3][1] = list("dp",wDp);
  rData[3][2] = list("wp",deg_I);
  def R_aux = ring(rData);
  setring R_aux;
  //---------------------------------------------------------------

  ideal J;
  map phi = BR, x(1..n);
  ideal I = phi(I);
  for (k=1; k<=s; k++) { J[k] = y(k)-I[k]; } 
  option(redSB);
  J = std(J);

  // Remove all generators that are depending on some x(i) from J: 
  int s_J = size(J);
  for (k=1; k<=s_J; k++) { if (J[k]>=x(n)) {J[k]=0;} }

  // The monomial order on K[y] is chosen such that linear leading 
  // terms in J are in 1-1 correspondence to superfluous generators
  // in I :
  ideal J_1jet = std(jet(lead(J),1));
  intvec to_remove; 
  i=1;
  for (k=1; k<=s; k++)
  { 
    if (reduce(y(k),J_1jet)==0){ to_remove[i]=k; i++; }
  }
  setring BR;
  if (to_remove == 0) { return(ideal(I)); }
  for (i=1; i<=size(to_remove); i++)
  { 
    I[to_remove[i]] = 0;
  } 
  I = simplify(I,2);
  return(I);
}

ideal FSI = min_generating_set(P,S);
size(FSI);
//->   14
ring Rnew = 0, (x(0..3)), dp;  // coefficient field is now Q
ideal FSI = fetch(R,FSI);
ideal ZERO;
ring R1 = 0, (y(0..13)), wp(5,5,5,5,4,4,4,4,3,3,3,3,2,2);
ideal REL = preimage(Rnew,FSI,ZERO);
homog(REL);                    // check that REL is homogeneous
//->   1
size(REL);
//->   54
setring Rnew;
FSI[4];
//->   x(0)^5+x(1)^5+x(2)^5+x(3)^5
ideal F = FSI[4];
setring R1;
ideal GODEAUX = preimage(Rnew,FSI,F);
size(GODEAUX);
//->   55
GODEAUX[1];
//->   y(3)
dim(std(GODEAUX));
//->   3

tst_status(1);$