source: git/Singular/LIB/rinvar.lib @ 371848

// Last change 10.12.2000 (TB)
///////////////////////////////////////////////////////////////////////////////
version="$Id: rinvar.lib,v 1.5 2000-12-22 14:43:04 greuel Exp $";
category="Invariant theory";
info="
LIBRARY:  rinvar.lib      Invariant Rings of Reductive Groups
AUTHOR:   Thomas Bayer,   tbayer@in.tum.de 
          http://wwwmayr.informatik.tu-muenchen.de/personen/bayert/
          Current Adress: Institut fuer Informatik, TU Muenchen
OVERVIEW: 
 Implementation based on Derksen's algorithm. Written in the frame of the 
 diploma thesis (advisor: Prof. Gert-Martin Greuel) 'Computations of moduli 
 spaces of semiquasihomogenous singularities and an implementation in Singular'

PROCEDURES:
 HilbertSeries(I, w);           Hilbert series of the ideal I w.r.t. weight w
 HilbertWeights(I, w);          weighted degrees of the generators of I
 ImageVariety(I, F);            ideal of the image variety F(variety(I))
 ImageGroup(G, F);              ideal of G w.r.t. the induced representation
 InvariantRing(G, Gaction);     generators of the invariant ring of G
 InvariantQ(f, G, Gaction);     decide if f is invariant w.r.t. G
 LinearizeAction(G, Gaction);   linearization of the action 'Gaction' of G
 LinearActionQ(action,s,t);     decide if action is linear in var(s..nvars)
 LinearCombinationQ(base, f);   decide if f is in the linear hull of 'base'
 MinimalDecomposition(f,s,t);   minimal decomposition of f (like coef)
 NullCone(G, act);              ideal of the nullcone of the action 'act' of G
 ReynoldsImage(RO,f);           image of f under the Reynolds operator 'RO'
 ReynoldsOperator(G, Gaction);  Reynolds operator of the group G
 SimplifyIdeal(I[,m,s]);        simplify the ideal I (try to reduce variables)
 TransferIdeal(R,name,nA);      transfer the ideal 'name' from R to basering

SEE ALSO: qhmoduli.lib, zeroset.lib
";

LIB "presolve.lib";
LIB "elim.lib";
LIB "zeroset.lib";

///////////////////////////////////////////////////////////////////////////////

proc EquationsOfEmbedding(ideal embedding, int nrs)
"USAGE:   EquationsOfEmbedding(embedding, s); ideal embedding; int s;
PUROPSE: compute the ideal of the variety parameterized by 'embedding' by
         implicitation and change the variables to the old ones.
RETURN:  ideal
ASSUME:  nvars(basering) = n, size(embedding) = r and s = n - r.
         The polynomials of embedding contain only var(s + 1 .. n).
NOTE:    the result is the Zariski closure of the parameterized variety
EXAMPLE: example  EquationsOfEmbedding; shows an example
"
{
  ideal tvars;

  for(int i = nrs + 1; i <= nvars(basering); i++) { tvars[i - nrs] = var(i); }

  def RE1 = ImageVariety(ideal(0), embedding);  // implicitation of the parameterization
  // map F = RE1, tvars;
        map F = RE1, maxideal(1);
  return(F(imageid));
}
example
{"EXAMPLE:";  echo = 2;
  ring R   = 0,(s(1..5), t(1..4)),dp;
  ideal emb = t(1), t(2), t(3), t(3)^2;
  ideal I = EquationsOfEmbedding(emb, 5);
  I;
}

///////////////////////////////////////////////////////////////////////////////

proc ImageGroup(ideal Grp, ideal Gaction)
"USAGE:   ImageGroup(G, action); ideal G, action;
PUROPSE: compute the ideal of the image of G in GL(m,K) induced by the linear
         action 'action', where G is an algebraic group and 'action' defines
         an action of G on K^m (size(action) = m).
RETURN:  ring, a polynomial ring over the same ground field as the basering,
         containing the ideals 'groupid' and 'actionid'.
         - 'groupid' is the ideal of the image of G (order <= order of G)
         - 'actionid' defines the linear action of 'groupid' on K^m.
NOTE:    'action' and 'actionid' have the same orbits
         all variables which give only rise to 0's in the m x m matrices of G
         have been omitted.
ASSUME:  basering K[s(1..r),t(1..m)] has r + m variables, G is the ideal of an
         algebraic group and F is an action of G on K^m. G contains only the
         variables s(1)...s(r). The action 'action' is given by polynomials
         f_1,...,f_m in basering, s.t. on the ring level we have
         K[t_1,...,t_m] --> K[s_1,...,s_r,t_1,...,t_m]/G
         t_i   -->  f_i(s_1,...,s_r,t_1,...,t_m)
EXAMPLE: example ImageGroup; shows an example
"
{
  int i, j, k, newVars, nrt, imageSize, dbPrt;
  ideal matrixEntries;
  matrix coMx;
  poly tVars, mPoly;
  string ringSTR1, ringSTR2, order;

  dbPrt = printlevel-voice+2;
  dbprint(dbPrt, "Image Group of " + string(Grp) + ", action =  " + string(Gaction));
  def RIGB = basering;
  mPoly = minpoly;
  tVars = 1;
  k = 0;

  // compute the representation of G induced by Gaction, i.e., a matrix
  // of size(Gaction) x size(Gaction) and polynomials in s(1),...,s(r) as entries
  // the matrix is represented as the list 'matrixEntries' where
  // the entries which are always 0 are omittet.

  for(i = 1; i <= ncols(Gaction); i++) {
    tVars = tVars * var(i + nvars(basering) - ncols(Gaction));
  }
  for(i = 1; i <= ncols(Gaction); i++){
    coMx = coef(Gaction[i], tVars);
    for(j = 1; j <= ncols(coMx); j++){
      k++;
      matrixEntries[k] = coMx[2, j];
    }
  }
  newVars = size(matrixEntries);
  nrt = ncols(Gaction);

  // this matrix defines an embedding of G into GL(m, K).
  // in the next step the ideal of this image is computed
  // note that we have omitted all variables which give give rise
  // only to 0's. Note that z(1..newVars) are slack variables

  order = "(dp(" + string(nvars(basering)) + "), dp);";
  ringSTR1 = "ring RIGR = (" + charstr(basering) + "), (" + varstr(basering) + ", z(1.." + string(newVars) + "))," + order;
  execute(ringSTR1);
  minpoly = number(imap(RIGB, mPoly));
  ideal I1, I2, Gn, G, F, mEntries, newGaction;
  G = imap(RIGB, Grp);
  F = imap(RIGB, Gaction);
  mEntries = imap(RIGB, matrixEntries);

  // prepare the ideals needed to compute the image
  // and compute the new action of the image on K^m

  for(i = 1; i <= size(mEntries); i++) { I1[i] = var(i + nvars(RIGB)) - mEntries[i]; }
  I1 = std(I1);

  for(i = 1; i <= ncols(F); i++) { newGaction[i] = reduce(F[i], I1); }
  I2 = G, I1;
  I2 = std(I2);
  Gn = nselect(I2, 1, nvars(RIGB));
  imageSize = ncols(Gn);

  // create a new basering which might contain more variables
  // s(1..newVars) as the original basering and map the ideal
  // Gn (contians only z(1..newVars)) to this ring

  ringSTR2 = "ring RIGS = (" + charstr(basering) + "), (s(1.." + string(newVars) + "), t(1.." + string(nrt) + ")), lp;";
  execute(ringSTR2);
  minpoly = number(imap(RIGB, mPoly));
  ideal mapIdeal, groupid, actionid;
  int offset;

  // construct the map F : RIGB -> RIGS

  for(i = 1; i <= nvars(RIGB) - nrt; i++) { mapIdeal[i] = 0;}      // s(i) -> 0
  offset = nvars(RIGB) - nrt;
  for(i = 1; i <= nrt; i++) { mapIdeal[i + offset] = var(newVars + i);}    // t(i) -> t(i)
  offset = offset + nrt;
  for(i = 1; i <= newVars; i++) { mapIdeal[i + offset] = var(i);}      // z(i) -> s(i)

  // map Gn and newGaction to RIGS

  map F = RIGR, mapIdeal;
  groupid = F(Gn);
  actionid = F(newGaction);
  export(groupid);
  export(actionid);
  dbprint(dbPrt, "
// 'ImageGroup' created a new ring.
// To see the ring, type (if the name of the ring is R):
     show(R);
// To access the ideal of the image-group Gn of the group G w.r.t. 'action'
// and the new action of Gn, type
     def R = ImageGroup(G, action); setring R;  groupid; actionid;
// 'groupid' is the ideal defining the image of the group G w.r.t. 'action'
// and 'actionid' is the new action of 'groupid'.
");
  return(RIGS);
}
example
{"EXAMPLE:";  echo = 2;
  ring B   = 0,(s(1..2), t(1..2)),dp;
  ideal G = s(1)^3-1, s(2)^10-1;
  ideal action = s(1)*s(2)^8*t(1), s(1)*s(2)^7*t(2);
  def R = ImageGroup(G, action);
  setring R;
  groupid;
  actionid;
}

///////////////////////////////////////////////////////////////////////////////

proc HilbertWeights(ideal I, wt)
"USAGE:   HilbertWeights(I, w); ideal I, intvec wt
PUROPSE: compute the weights of the "slack" varaibles needed for the
         computation of the algebraic relations of the generators of 'I' s.t.
         the Hilbert driven 'std' can be used.
RETURN: intvec
ASSUME: basering = K[t_1,...,t_m,...], 'I' is quasihomogenous w.r.t. 'w' and
        contains only polynomials in t_1,...,t_m
"
{
  int offset = size(wt);
  intvec wtn = wt;

  for(int i = 1; i <= size(I); i++) { wtn[offset + i] = deg(I[i], wt); }
  return(wtn);
}

///////////////////////////////////////////////////////////////////////////////

proc HilbertSeries(ideal I, wt)
"USAGE:   HilbertSeries(I, w); ideal I, intvec wt
PUROPSE: compute the polynomial p of the Hilbert Series,represented by p/q, of
         the ring K[t_1,...,t_m,y_1,...,y_r]/I1 where 'w' are the weights of
         the variables, computed, e.g., by 'HilbertWeights', 'I1' is of the
         form I[1] - y_1,...,I[r] - y_r and is quasihomogenous w.r.t. 'w'
RETURN:  intvec
NOTE:    the leading 0 of the result does not belong to p, but is needed in
         the hilbert-driven 'std'.
"
{
  int i;
  intvec hs1;
  matrix coMx;
  poly f = 1;

  for(i = 1; i <= ncols(I); i++) { f = f * (1 - var(1)^deg(I[i], wt));}
  coMx = coeffs(f, var(1));
  for(i = 1; i <= deg(f) + 1; i++) {
    hs1[i] = int(coMx[i, 1]);
  }
  hs1[size(hs1) + 1] = 0;
  return(hs1);
}
///////////////////////////////////////////////////////////////////////////////

proc HilbertSeries1(wt)
"USAGE:   HilbertSeries1(wt); ideal I, intvec wt
PUROPSE: compute the polynomial p of the Hilbert Series represented by p/q of
         the ring K[t_1,...,t_m,y_1,...,y_r]/I where I is a complete inter-
         section and the generator I[i] has degree wt[i]
RETURN:  poly
"
{
  int i, j;
  intvec hs1;
  matrix ma;
  poly f = 1;

  for(i = 1; i <= size(wt); i++) { f = f * (1 - var(1)^wt[i]);}
  ma = coef(f, var(1));
  j = ncols(ma);
  for(i = 0; i <= deg(f); i++) {
    if(var(1)^i == ma[1, j]) {
      hs1[i + 1] = int(ma[2, j]);
      j--;
    }
    else { hs1[i + 1] = 0; }
  }
  hs1[size(hs1) + 1] = 0;
  return(hs1);
}

///////////////////////////////////////////////////////////////////////////////

proc ImageVariety(ideal I, F, list #)
"USAGE:   ImageVariety(ideal I, F [, w]);ideal I; F is a list/ideal, intvec w.
PUROPSE: compute the Zariski closure of the image of the variety of I under
         the morphism F.
NOTE:    if 'I' and 'F' are quasihomogenous w.r.t. 'w' then the Hilbert-driven
         'std' is used.
RETURN:  polynomial ring over the same groundfield, containing the ideal
         'imageid'. The variables are Y(1),...,Y(k) where k = size(F)
         - 'imageid' is the ideal of the Zariski closure of F(X) where
           X is the variety of I.
EXAMPLE: example ImageVariety; shows an example
"
{
  int i, dbPrt, nrNewVars;
  intvec wt, wth, hs1;
  string ringSTR1, ringSTR2, order;

  def RARB = basering;
  nrNewVars = size(F);

  dbPrt = printlevel-voice+2;
  dbprint(dbPrt, "ImageVariety of " + string(I) + " under the map " + string(F));

  if(size(#) > 0) { wt = #[1]; }

  // create new ring for elimination, Y(1),...,Y(m) are slack variables.

  poly mPoly = minpoly;
  order = "(dp(" + string(nvars(basering)) + "), dp);";
  ringSTR1 = "ring RAR1 = (" + charstr(basering) + "), (" + varstr(basering) + ", Y(1.." + string(nrNewVars) + ")), " + order;
  ringSTR2 = "ring RAR2 = (" + charstr(basering) + "), Y(1.." + string(nrNewVars) + "), dp;";
  execute(ringSTR1);
  minpoly = number(imap(RARB, mPoly));

  ideal I, J1, J2, Fm;

  I = imap(RARB, I);
  Fm = imap(RARB, F);

  if(size(wt) > 1) {
    wth = HilbertWeights(Fm, wt);
    hs1 = HilbertSeries(Fm, wt);
  }

  // get the ideal of the graph of F : X -> Y and compute a standard basis

  for(i = 1; i <= nrNewVars; i++) { J1[i] = var(i + nvars(RARB)) - Fm[i];}
  J1 = J1, I;
  if(size(wt) > 1) {
    J1 = std(J1, hs1, wth);  // Hilbert-driven algorithm
  }
  else {
    J1 = std(J1);
  }

  // forget all elements which contain other than the slack variables

  J2 = nselect(J1, 1, nvars(RARB));

  execute(ringSTR2);
  minpoly = number(imap(RARB, mPoly));
  ideal imageid = imap(RAR1, J2);
  export(imageid);
     dbprint(dbPrt, "
// 'ImageVariety' created a new ring.
// To see the ring, type (if the name of the ring is R):
     show(R);
// To access the ideal of the image F(X), where F is a map and X is a variety
// with ideal I,type
     def R = ImageVariety(I, F); setring R;  imageid;
// 'imageid' is the ideal of the Zariski closure of F(X).
");
  return(RAR2);
}
example
{"EXAMPLE:";  echo = 2;
  ring B   = 0,(x,y),dp;
  ideal I  = x4 - y4;
  ideal F  = x2, y2, x*y;
  def R = ImageVariety(I, F);
  setring R;
  imageid;
}

///////////////////////////////////////////////////////////////////////////////

proc LinearizeAction(ideal Grp, Gaction, int nrs)
"USAGE:   LinearizeAction(G,action,r); ideal G, action; int r
PUROPSE: linearize the group action 'action' and find an equivariant embedding
         of K^m where m = size(action).
ASSUME:  G contains only variables var(1..r) (r = nrs)
basering = K[s(1..r),t(1..m)], K = Q or K = Q(a) and minpoly != 0.
RETURN:  polynomial ring contianing the ideals 'actionid', 'embedid', 'groupid'
         - 'actionid' is the ideal defining the linearized action of G
         - 'embedid' is a parameterization of an equivariant embedding (closed)
         - 'groupid' is the ideal of G in the new ring
NOTE:    set printlevel > 0 to see a trace
EXAMPLE: example LinearizeAction; shows an example
"
{
  int i, j, k, ok, loop, nrt, sizeOfDecomp, dbPrt;
  intvec wt;
  ideal action, basis, G, reduceIdeal;
  matrix decompMx;
  poly actCoeff;
  string str, order, mPoly;

  dbPrt = printlevel-voice+2;
  dbprint(dbPrt, "LinearizeAction " + string(Gaction));
  def RLAR = basering;
  mPoly = "minpoly = " + string(minpoly) + ";";
  order = ordstr(basering);
  nrt = ncols(Gaction);
  for(i = 1; i <= nrs; i++) { wt[i] = 0;}
  for(i = nrs + 1; i <= nrs + nrt; i++) { basis[i - nrs] = var(i); wt[i] = 1;}
  dbprint(dbPrt, "  basis = " + string(basis));
  if(attrib(Grp, "isSB")) { G = Grp; }
  else { G = std(Grp); }
  reduceIdeal = G;
  action = Gaction;
  loop = 1;
  i = 1;

  // check if each component of 'action' is linear in t(1),...,t(nrt).

  while(loop){
    if(deg(action[i], wt) <= 1) {
      sizeOfDecomp = 0;
      dbprint(dbPrt, "  " + string(action[i]) + " is linear");
    }
    else {   // action[i] is not linear

    // compute the minimal decomposition of action[i]
    // action[i] = decompMx[1,1]*decompMx[2,1] + ... + decompMx[1,k]*decompMx[2,k]
    // decompMx[1,j] contains variables var(1)...var(nrs)
    // decompMx[2,j] contains variables var(nrs + 1)...var(nvars(basering))

    dbprint(dbPrt, "  " + string(action[i]) + " is not linear, a minimal decomposition is :");
    decompMx = MinimalDecomposition(action[i], nrs, nrt);
    sizeOfDecomp = ncols(decompMx);
    dbprint(dbPrt, decompMx);

    for(j = 1; j <= sizeOfDecomp; j++) {   // check if decompMx[2,j] is a linear combination of basis elements
      actCoeff = decompMx[2, j];
      ok = LinearCombinationQ(basis, actCoeff, nrt + nrs);
      if(ok == 0) {

        // nonlinear element, compute new component of the action

        dbprint(dbPrt, "  the polynomial " + string(actCoeff) + " is not a linear combination of the elements of basis");
        nrt++;
        str = charstr(basering) + ", (" + varstr(basering) + ",t(" + string(nrt) + ")),";
        if(defined(RLAB)) { kill(RLAB);}
        def RLAB = basering;
        if(defined(RLAR)) { kill(RLAR);}
        execute("ring RLAR = " + str + "(" + order + ");");
        execute(mPoly);

        ideal basis, action, G, reduceIdeal;
        matrix decompMx;
        map F;
        poly actCoeff;

        wt[nrs + nrt] = 1;
        basis = imap(RLAB, basis), imap(RLAB, actCoeff);
        action = imap(RLAB, action);
        decompMx = imap(RLAB, decompMx);
        actCoeff = imap(RLAB, actCoeff);
        G = imap(RLAB, G);
        attrib(G, "isSB", 1);
        reduceIdeal = imap(RLAB, reduceIdeal), actCoeff - var(nrs + nrt);

        // compute action on the new basis element

        for(k = 1; k <= nrs; k++) { F[k] = 0;}
        for(k = nrs + 1; k < nrs + nrt; k++) { F[k] = action[k - nrs];}
        actCoeff = reduce(F(actCoeff), G);
        action[ncols(action) + 1] = actCoeff;
        dbprint(dbPrt, "  extend basering by " + string(var(nrs + nrt)));
        dbprint(dbPrt, "  new basis = " + string(basis));
        dbprint(dbPrt, "  action of G on new basis element = " + string(actCoeff));
        dbprint(dbPrt, "  decomp : " + string(decompMx[2, j]) + " -> " + string(var(nrs + nrt)));
      } // end if
      else {
        dbprint(dbPrt, "  the polynomial " + string(actCoeff) + " is a linear combination of the elements of basis");
      }
    } // end for
    reduceIdeal = std(reduceIdeal);
    action[i] = reduce(action[i], reduceIdeal);
    } // end else
    if(i < ncols(action)) { i++;}
    else {loop = 0;}
  } // end while
  if(defined(actionid)) { kill(actionid); }
  ideal actionid, embedid, groupid;
  actionid = action;
  embedid = basis;
  groupid = G;
  export(actionid);
  export(embedid);
  export(groupid);
dbprint(dbPrt, "
// 'LinearizeAction' created a new ring.
// To see the ring, type (if the name of the ring is R):
     show(R);
// To access the new action and the equivariant embedding, where G and 'action'
// are the original group and group-action contained in K[s(1..ns)] and
// K[s(1..ns),t(1..nt)] respectively, type
     def R = LinearizeAction(G, action, ns, nt); setring R; actionid; embedid; groupid
// 'actionid' is the ideal of the linearized action, 'embedid' is the ideal
// defining the equivariant embedding and 'grouid' is the ideal G.
");
  return(RLAR);
}
example
{"EXAMPLE:";  echo = 2;
  ring B   = 0,(s(1..5), t(1..3)),dp;
  ideal G =  s(3)-s(4), s(2)-s(5), s(4)*s(5), s(1)^2*s(4)+s(1)^2*s(5)-1, s(1)^2*s(5)^2-s(5), s(4)^4-s(5)^4+s(1)^2, s(1)^4+s(4)^3-s(5)^3, s(5)^5-s(1)^2*s(5);
  ideal action = -s(4)*t(1)+s(5)*t(1), -s(4)^2*t(2)+2*s(4)^2*t(3)^2+s(5)^2*t(2), s(4)*t(3)+s(5)*t(3);
  LinearActionQ(action, 5);
  def R = LinearizeAction(G, action, 5);
  setring R;
  R;
  actionid;
  embedid;
  groupid;
  LinearActionQ(actionid, 5);
}

///////////////////////////////////////////////////////////////////////////////

proc LinearActionQ(Gaction, int nrs)
"USAGE:   LinearActionQ(action,nrs,nrt); ideal action, int nrs
PUROPSE: check if the action defined by 'action' is linear w.r.t. the variables
         var(nrs + 1...nvars(basering)).
RETURN:  0 action not linear
         1 action is linear
EXAMPLE: example LinearActionQ; shows an example
"
{
  int i, nrt, loop;
  intvec wt;

  nrt = ncols(Gaction);
  for(i = 1; i <= nrs; i++) { wt[i] = 0;}
  for(i = nrs + 1; i <= nrs + nrt; i++) { wt[i] = 1;}
  loop = 1;
  i = 1;
  while(loop)
  {
    if(deg(Gaction[i], wt) > 1) { loop = 0; }
    else
    {
      i++;
      if(i > ncols(Gaction)) { loop = 0;}
    }
  }
  return(i > ncols(Gaction));
}
example
{"EXAMPLE:";  echo = 2;
  ring R   = 0,(s(1..5), t(1..3)),dp;
  ideal G =  s(3)-s(4), s(2)-s(5), s(4)*s(5), s(1)^2*s(4)+s(1)^2*s(5)-1, s(1)^2*s(5)^2-s(5), s(4)^4-s(5)^4+s(1)^2, s(1)^4+s(4)^3-s(5)^3, s(5)^5-s(1)^2*s(5);
  ideal Gaction = -s(4)*t(1)+s(5)*t(1), -s(4)^2*t(2)+2*s(4)^2*t(3)^2+s(5)^2*t(2), s(4)*t(3)+s(5)*t(3);
  LinearActionQ(Gaction, 5, 3);
}

///////////////////////////////////////////////////////////////////////////////

proc LinearCombinationQ(ideal I, poly f)
"USAGE:   LinearCombination(I, f); ideal I, poly f
PUROPSE: test if f can be written as a linear combination of the generators of I.
RETURN:  0 f is not a linear combination
         1 f is a linear combination
"
{
  int i, loop, sizeJ;
  ideal J;

  J = I, f;
  sizeJ = size(J);

  def RLC = ImageVariety(ideal(0), J);   // compute algebraic relations
  setring RLC;
  matrix coMx;
  poly relation = 0;

  loop = 1;
  i = 1;
  while(loop)
  { // look for a linear relation containing Y(nr)
    if(deg(imageid[i]) == 1)
    {
      coMx = coef(imageid[i], var(sizeJ));
      if(coMx[1,1] == var(sizeJ))
      {
        relation = imageid[i];
        loop = 0;
      }
    }
    else
    {
      i++;
      if(i > ncols(imageid)) { loop = 0;}
    }
  }
  return(i <= ncols(imageid));
}

///////////////////////////////////////////////////////////////////////////////

proc InvariantRing(ideal G, ideal action, list #)
"USAGE:   InvariantRing(G, Gact [, opt]); ideal G, Gact; int opt
PUROPSE: compute generators of the invariant ring of G w.r.t. the action 'Gact'
ASSUME:  G is a finite group and 'Gact' is a linear action.
RETURN:  polynomial ring over a simple extension of the groundfield of the
         basering (the extension might be trivial), containing the ideals
         'invars' and 'groupid' and the poly 'newA'
         - 'invars' contains the algebra-generators of the invariant ring
         - 'groupid' is the ideal of G in the new ring
         - 'newA' if the minpoly changes this is the new representation of the
           algebraic number, otherwise it is set to 'a'.
NOTE:    the delivered ring might have a different minimal polynomial
EXAMPLE: example InvariantRing; shows an example
"
{
  int i, ok, dbPrt, noReynolds, primaryDec;
  ideal invarsGens, groupid;

  dbPrt = printlevel-voice+2;

  if(size(#) > 0) { primaryDec = #[1]; }
  else { primaryDec = 0; }

  dbprint(dbPrt, "InvariantRing of " + string(G));
  dbprint(dbPrt, " action  = " + string(action));

  if(!attrib(G, "isSB")) { groupid = std(G);}
  else { groupid = G; }

  // compute the nullcone of G by means of Derksen's algorithm

  invarsGens = NullCone(groupid, action);    // compute the nullcone of the linear action
  dbprint(dbPrt, " generators of zero-fibre ideal are " + string(invarsGens));

  // make all generators of the nullcone invariant
  // if necessary, compute the Reynolds Operator, i.e., find all elements
  // of the variety defined by G. It might be necessary to extend the groundfield.

  def IRB = basering;
  if(defined(RIRR)) { kill(RIRR);}
  def RIRR = basering;
  setring RIRR;
  export(RIRR);
  export(invarsGens);
  noReynolds = 1;
  dbprint(dbPrt, " nullcone is generated by " + string(size(invarsGens)));
   dbprint(dbPrt, " degrees = " + string(maxdeg(invarsGens)));
  for(i = 1; i <= ncols(invarsGens); i++){
    ok =  InvariantQ(invarsGens[i], groupid, action);
    if(ok) { dbprint(dbPrt, string(i) + ": poly " + string(invarsGens[i]) + " is invariant");}
    else {
      if(noReynolds) {

        // compute the Reynolds operator and change the ring !

        def RORN = ReynoldsOperator(groupid, action, primaryDec);
        noReynolds = 0;
        setring RORN;
        export(RORN);
        ideal groupid = std(id);
        attrib(groupid, "isSB", 1);
        ideal action = actionid;
        ideal invarsGens = TransferIdeal(RIRR, "invarsGens", newA);
        export(invarsGens);
        kill(RIRR);

      }
      dbprint(dbPrt, string(i) + ": poly " + string(invarsGens[i]) + " is NOT invariant");
      invarsGens[i] = ReynoldsImage(ROelements, invarsGens[i]);
      dbprint(dbPrt, " --> " + string(invarsGens[i]));
    }
  }
  for(i = 1; i <= ncols(invarsGens); i++){
    ok =  InvariantQ(invarsGens[i], groupid, action);
    if(ok) { dbprint(dbPrt, string(i) + ": poly " + string(invarsGens[i]) + " is invariant"); }
    else { print(string(i) + ": Fatal Error with Reynolds ");}
  }
  kill(IRB);
  if(noReynolds == 0) {
    def RIRS = RORN;
    setring(RIRS);
    kill(RORN);export(groupid);
  }
  else {
    def RIRS = RIRR;
    kill(RIRR);
    setring(RIRS);
                export(groupid);
  }
  ideal invars = invarsGens;
  kill(invarsGens);
  export(invars);
  // export(groupid);
dbprint(dbPrt, "
// 'InvariantRing' created a new ring.
// To see the ring, type (if the name of the ring is R):
     show(R);
// To access the generators of the invariant ring of G w.r.t. the linear
// group-action 'action' of G, where G is contained in K[s(1..ns)] and
// 'action' in K[s(1..ns),t(1..nt)], type
     def R = InvariantRing(G, action); setring R; invars;
// 'invars' contains generator of the invariant ring.
// Note that G is containd in R as the ideal 'groupid', to see it, type
    groupid;
// Note that 'InvariantRing' might change the minimal polynomial
// The representation of the algebraic number is given by 'newA'
");
  return(RIRS);
}
example
{"EXAMPLE:";  echo = 2;
  ring B = 0, (s(1..2), t(1..2)), dp;
  ideal G = -s(1)+s(2)^3, s(1)^4-1;
  ideal action = s(1)*t(1), s(2)*t(2);

  def R = InvariantRing(std(G), action);
  setring R;
  invars;
}

///////////////////////////////////////////////////////////////////////////////

proc InvariantQ(poly f, ideal G, action)
"USAGE:   InvariantQ(f, G, action); poly f; ideal G, action
PUROPSE: check if the polynomial f is invariant w.r.t. G where G acts via
         'action' on K^m.
ASSUME:  basering = K[s_1,...,s_m,t_1,...,t_m] where K = Q of K = Q(a) and
         minpoly != 0, f contains only t_1,...,t_m, G is the ideal of an
         algebraic group and a standardbasis.
RETURN:  int;
         0 if f is not invariant,
         1 if f is invariant
NOTE:    G need not be finite
EXAMPLE: example InvariantQ; shows an example
"
{
  map F;

  if(deg(f) == 0) { return(1); }
  for(int i = 1; i <= size(action); i++) {
    F[nvars(basering) - size(action) + i] = action[i];
  }
  return(reduce(f - F(f), G) == 0);
}

///////////////////////////////////////////////////////////////////////////////

proc MinimalDecomposition(poly f, int nrs, int nrt)
"USAGE:   MinimalDecomposition(f,a,b); poly f; int a, b.
PUROPSE: decompose f as a sum M[1,1]*M[2,1] + ... + M[1,r]*M[2,r] where M[1,i]
         contains only s(1..a), M[2,i] contains only t(1...b) s.t. r is minimal
ASSUME:  f polynomial in K[s(1..a),t(1..b)], K = Q or K = Q(a) and minpoly != 0
RETURN:  2 x r matrix M s.t.  f = M[1,1]*M[2,1] + ... + M[1,r]*M[2,r]
EXAMPLE: example MinimalDecomposition;
"
{
  int i, sizeOfMx, changed, loop;
  list initialTerms;
  matrix coM1, coM2, coM, decompMx, auxM;
  matrix m[2][2] = 0,1,1,0;
  poly vars1, vars2;

  if(f == 0) { return(decompMx); }

  // first decompose f w.r.t. t(1..nrt)
  // then  decompose f w.r.t. s(1..nrs)

  vars1 = RingVarProduct(nrs+1..nrt+nrs);
  vars2 = RingVarProduct(1..nrs);
  coM1 = SimplifyCoefficientMatrix(m*coef(f, vars1));    // exchange rows of decomposition
  coM2 = SimplifyCoefficientMatrix(coef(f, vars2));
  if(ncols(coM2) < ncols(coM1)) {
    auxM = coM1;
    coM1 = coM2;
    coM2 = auxM;
  }
  decompMx = coM1;        // decompMx is the smaller decomposition
  if(ncols(decompMx) == 1) { return(decompMx);}      // n = 1 is minimal
  changed = 0;
  loop = 1;
  i = 1;

  // first loop, try coM1

  while(loop) {
    coM = MinimalDecomposition(f - coM1[1, i]*coM1[2, i], nrs, nrt);
    if(size(coM) == 1) { sizeOfMx = 0; }    // coM = 0
    else {sizeOfMx = ncols(coM); }      // number of columns
    if(sizeOfMx + 1 < ncols(decompMx)) {    // shorter decomposition
      changed = 1;
      decompMx = coM;
      initialTerms[1] =  coM1[1, i];
      initialTerms[2] =  coM1[2, i];
    }
    if(sizeOfMx == 1) { loop = 0;}      // n = 2 is minimal
    if(i < ncols(coM1)) {i++;}
    else {loop = 0;}
  }
  if(sizeOfMx > 1) {            // n > 2
    loop = 1;          // coM2 might yield
    i = 1;            // a smaller decomposition
  }

  // first loop, try coM2

  while(loop) {
    coM = MinimalDecomposition(f - coM2[1, i]*coM2[2, i], nrs, nrt);
    if(size(coM) == 1) { sizeOfMx = 0; }
    else {sizeOfMx = ncols(coM); }
    if(sizeOfMx + 1 < ncols(decompMx)) {
      changed = 1;
      decompMx = coM;
      initialTerms[1] =  coM2[1, i];
      initialTerms[2] =  coM2[2, i];
    }
    if(sizeOfMx == 1) { loop = 0;}
    if(i < ncols(coM2)) {i++;}
    else {loop = 0;}
  }
  if(!changed) { return(decompMx); }
  if(size(decompMx) == 1) { matrix decompositionM[2][1];}
  else  { matrix decompositionM[2][ncols(decompMx) + 1];}
  decompositionM[1, 1] = initialTerms[1];
  decompositionM[2, 1] = initialTerms[2];
  if(size(decompMx) > 1) {
    for(i = 1; i <= ncols(decompMx); i++) {
      decompositionM[1, i + 1] = decompMx[1, i];
      decompositionM[2, i + 1] = decompMx[2, i];
    }
  }
  return(decompositionM);
}
example
{"EXAMPLE:"; echo = 2;
  ring R = 0, (s(1..2), t(1..2)), dp;
  poly h = s(1)*(t(1) + t(1)^2) +  (t(2) + t(2)^2)*(s(1)^2 + s(2));
  matrix M = MinimalDecomposition(h, 2, 2);
  M;
  M[1,1]*M[2,1] + M[1,2]*M[2,2] - h;
}

///////////////////////////////////////////////////////////////////////////////

proc NullCone(ideal G, action)
"USAGE:   NullCone(G, action); ideal G, action
PUROPSE: compute the ideal of the nullcone of the linear action of G on K^n,
         given by 'action', by means of Deksen's algorithm
ASSUME:  basering = K[s(1..r),t(1..n)], K = Q or K = Q(a) and minpoly != 0,
         G is an ideal of a reductive algebraic group in K[s(1..r)],
         'action' is a linear group action of G on K^n (n = ncols(action))
RETURN:  ideal of the nullcone of G.
NOTE:    the generators of the nullcone are homogenous, but i.g. not invariant
EXAMPLE: example NullCone; shows an example
"
{
  int i, nt, dbPrt, offset, groupVars;
  poly minPoly;
  string ringSTR, vars, order;
  def RNCB = basering;

  // prepare the ring needed for the computation
  // s(1...) variables of the group
  // t(1...) variables of the affine space
  // y(1...) additional 'slack' variables

  nt = size(action);
  order = "(dp(" + string(nvars(basering) - nt) + "), dp);";
  vars = "(s(1.." + string(nvars(basering) - nt);
  vars = vars + "), t(1.." + string(nt) + "), Y(1.." + string(nt) + "))," + order;
  ringSTR = "ring RNCR = (" + charstr(basering) +  "),"  + vars; // ring for the computation

  minPoly = minpoly;
  offset =  size(G) + nt;
  execute(ringSTR);
  minpoly = number(imap(RNCB, minPoly));
  ideal action, G, I, J, N, generators;
  map F;
  poly f;

  // built the ideal of the graph of GxV -> V, (s,v) -> s(v), i.e.
  // of the image of the map GxV -> GxVxV, (s,v) -> (s,v,s(v))

  G = fetch(RNCB, G);
  action = fetch(RNCB, action);
  groupVars = nvars(basering) - 2*nt;
  offset =  groupVars + nt;
  I = G;
  for(i = 1; i <= nt; i = i + 1) {
    I = I, var(offset + i) - action[i];
  }

  J = std(I);  // takes long, try to improve

  // eliminate

  N = nselect(J, 1, groupVars);

  // substitute

  for(i = 1; i <= nvars(basering); i = i + 1) { F[i] = 0; }
  for(i = groupVars + 1; i <= offset; i = i + 1) { F[i] = var(i); }

  generators = mstd(F(N))[2];
        setring(RNCB);
  return(fetch(RNCR, generators));
}
example
{"EXAMPLE:";  echo = 2;
  ring R = 0, (s(1..2), x, y), dp;
  ideal G = -s(1)+s(2)^3, s(1)^4-1;
  ideal action = s(1)*x, s(2)*y;

  ideal inv = NullCone(G, action);
  inv;
}

///////////////////////////////////////////////////////////////////////////////

proc ReynoldsOperator(ideal Grp, ideal Gaction, list #)
"USAGE:   ReynoldsOperator(G, action [, opt); ideal G, action; int opt
PUROPSE: compute the Reynolds operator of the group G which act via 'action'
RETURN:  polynomial ring R over a simple extension of the groundfield of the
         basering (the extension might be trivial), containing a list
         'ROelements', the ideals 'id', 'actionid' and the polynomial 'newA'.
         R = K(a)[s(1..r),t(1..n)].
         - 'ROelements'  is a list of ideal, each ideal represents a
           substitution map F : R -> R according to the zero-set of G
         - 'id' is the ideal of G in the new ring
         - 'newA' is the new representation of a' in terms of a. If the
           basering does not contain a parameter then 'newA' = 'a'.
ASSUME:  basering = K[s(1..r),t(1..n)], K = Q or K = Q(a') and minpoly != 0,
         G is the ideal of a finite group in K[s(1..r)], 'action' is a linear
         group action of G
EXAMPLE: example ReynoldsOperator; shows an example
"
{
  int i, j, n, ns, primaryDec;
  ideal G1 = Grp;
  list solution, saction;
  string str;

  if(size(#) > 0) { primaryDec = #[1]; }
  else { primaryDec = 0; }

  n = nvars(basering);
  ns = n - size(Gaction);
  for(i = ns + 1; i <= n; i++) { G1 = G1, var(i);}

  def ROBR = basering;
  export(Grp);
  export(Gaction);
  def RORN = ZeroSet(G1, primaryDec);
  setring RORN;
  id = TransferIdeal(ROBR, "Grp", newA);  // defined in ZeroSet ...
  ideal actionid = TransferIdeal(ROBR, "Gaction", newA);
  list ROelements;
  ideal Rf;
  map groupElem;
  poly h1, h2;

  for(i = 1; i <= size(zeroset); i++) {
    groupElem = zeroset[i];    // element of G
    for(j = ns + 1; j <= n; j++) { groupElem[j] = var(j); } // do not change t's
    for(j = 1; j <= n - ns; j++) {
      h1 = actionid[j];
      h2 = groupElem(h1);
      Rf[ns + j] = h2;
    }
    ROelements[i] = Rf;
  }
  export(actionid);
  export(ROelements);
  return(RORN);
}

///////////////////////////////////////////////////////////////////////////////

proc ReynoldsImage(list reynoldsOp, poly f)
"USAGE:   ReynoldsImage(RO, f); list RO, poly f
PUROPSE: compute the Reynolds image of the polynomial f where RO represents
         the Reynolds operator
RETURN:  poly
"
{
  map F;
  poly h = 0;

  for(int i = 1; i <= size(reynoldsOp); i++) {
    F = basering, reynoldsOp[i];
    h = h + F(f);
  }
  return(h/size(reynoldsOp));
}

///////////////////////////////////////////////////////////////////////////////

static proc SimplifyCoefficientMatrix(matrix coefMatrix)
"USAGE:   SimplifyCoefficientMatrix(M); M matrix coming from coef(...)
PUROPSE: simplify the matrix, i.e. find linear dependencies among the columns
RETURN:  matrix M, f = M[1,1]*M[2,1] + ... + M[1,n]*M[2,n]
"
{
  int i, j , loop;
  intvec columnList;
  matrix decompMx = coefMatrix;

  loop = 1;
  i = 1;
  while(loop) {
    columnList = 1..i;            // current column
    for(j = i + 1; j <= ncols(decompMx); j++) {
      // test if decompMx[2, j] equals const * decompMx[2, i]
      if(LinearCombinationQ(ideal(decompMx[2, i]), decompMx[2, j])) {    // column not needed
        decompMx[1, i] = decompMx[1, i] +  decompMx[2, j] / decompMx[2, i] * decompMx[1, j];
      }
      else { columnList[size(columnList) + 1] = j; }
    }
    if(defined(auxM)) { kill(auxM);}
    matrix auxM[2][size(columnList)];      // built new matrix and omit
    for(j = 1; j <= size(columnList); j++) {    // the linear dependent colums
      auxM[1, j] = decompMx[1, columnList[j]];    // found above
      auxM[2, j] = decompMx[2, columnList[j]];
    }
    decompMx = auxM;
    if(i < ncols(decompMx) - 1) { i++;}
    else { loop = 0;}
  }
  return(decompMx);
}

///////////////////////////////////////////////////////////////////////////////

proc SimplifyIdeal(ideal I, list #)
"USAGE:   SimplifyIdeal(I [,m, name]); ideal I; int m, string name"
PURPOSE: simplify ideal I to the ideal I', do not change the names of the
         first m variables, new ideal I' might contain less variables.
         I' contains variables var(1..m)
RETURN: list
  _[1] ideal I'
  _[2] ideal representing a map phi to a ring with probably less vars. s.t.
       phi(I) = I'
  _[3] list of variables
  _[4] list from 'elimpart'
"
{
  int i, k, m;
  string nameCMD;
  ideal mId, In, mapId;  // ideal for the map
  list sList, result;

  sList = elimpart(I);
  In = sList[1];
  mapId = sList[5];

  if(size(#) > 0)
  {
    m = #[1];
    nameCMD = #[2];
  }
  else { m = 0;} // nvars(basering);
  k = 0;
  for(i = 1; i <= nvars(basering); i++)
  {
    if(sList[4][i] != 0)
    {
      k++;
      if(k <= m) { mId[i] = sList[4][i]; }
      else { execute("mId["+string(i) +"] = "+nameCMD+"("+string(k-m)+");");}
    }
    else { mId[i] = 0;}
  }
  map phi = basering, mId;
  result[1] = phi(In);
  result[2] = phi(mapId);
  result[3] = simplify(sList[4], 2);
  result[4] = sList;
  return(result);
}

///////////////////////////////////////////////////////////////////////////////

static proc TransferIdeal(R, string name, poly newA)
" USAGE:  TransferIdeal(R, name, newA); ring R, string name, poly newA
PUROPSE: Maps an ideal with name 'name' in R to the basering, s.t. all
         variables are fixed but par(1) is replaced by 'newA'.
RETURN:  ideal
NOTE:    this is used to transfor an ideal if the minimal polynomial has changed
"
{
  def RAB = basering;
  def RA1 = TransferRing(R);

  setring RA1;
  execute("ideal I = imap(R, " + name + ");");
  setring RAB;
  map F = RA1, maxideal(1);
  F[nvars(RAB) + 1] = newA;
  return(F(I));
}

///////////////////////////////////////////////////////////////////////////////

static proc RingVarProduct(index)
// list of indices
{
  poly f = 1;
  for(int i = 1; i <= size(index); i++)
  {
    f = f * var(index[i]);
  }
  return(f);
}
///////////////////////////////////////////////////////////////////////////////