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

/////////////////////////////////////////////////////////////////////////////
version="version rinvar.lib 4.1.2.0 Feb_2019 "; // $Id$
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 Address: Institut fuer Informatik, TU Muenchen
OVERVIEW:
 Implementation based on Derksen's algorithm. Written in the scope 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)

SEE ALSO: qhmoduli_lib, zeroset_lib
";

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

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

proc EquationsOfEmbedding(ideal embedding, int nrs)
"USAGE:   EquationsOfEmbedding(embedding, s); ideal embedding; int s;
PURPOSE: 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;
PURPOSE: 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;
  string ringSTR1, ringSTR2, order;

  dbPrt = printlevel-voice+2;
  dbprint(dbPrt, "Image Group of " + string(Grp) + ", action =  "
                                   + string(Gaction));
  def RIGB = basering;
  string @mPoly = string(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

  def RIGR=addNvarsTo(basering,newVars,"z",1); setring RIGR;
  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);
  execute("minpoly = number(" + @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, actionid;
  dbprint(dbPrt+1, "
// 'ImageGroup' created a new ring.
// To see the ring, type (if the name 'R' was assigned to the return value):
     show(R);
// To access the ideal of the image of the input group and to access the new
// action of the group, type
     setring R;  groupid; actionid;
");
  setring RIGB;
  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,intvec wt)
"USAGE:   HilbertWeights(I, w); ideal I, intvec wt
PURPOSE: compute the weights of the "slack" variables 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,intvec wt)
"USAGE:   HilbertSeries(I, w); ideal I, intvec wt
PURPOSE: 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
PURPOSE: 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,ideal F, list #)
"USAGE:   ImageVariety(ideal I, F [, w]);ideal I; F is a list/ideal, intvec w.
PURPOSE: 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 ground field, 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;
  def RARB = basering;
  nrNewVars = ncols(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.

  //string @mPoly = string(minpoly);
  def RAR1=addNvarsTo(basering,nrNewVars,"Y",1); setring RAR1;
  list RAR2l=ringlist(RAR1);
  list RAR2ll=RAR2l[2];
  RAR2ll=RAR2ll[size(RAR2ll)-nrNewVars+1..size(RAR2ll)];
  RAR2l[2]=RAR2ll;
  RAR2l[3]=list(list("dp",1:nrNewVars),list("C",0));
  def RAR2=ring(RAR2l);

  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));

  setring RAR2;
  ideal imageid = imap(RAR1, J2);
  export(imageid);
     dbprint(dbPrt+1, "
// 'ImageVariety' created a new ring.
// To see the ring, type (if the name 'R' was assigned to the return value):
     show(R);
// To access the ideal of the image variety, type
     setring R;  imageid;
");
  setring RARB;
  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,def Gaction, int nrs)
"USAGE:   LinearizeAction(G,action,r); ideal G, action; int r
PURPOSE: 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 containing 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
"
{
  def altring = basering;
  int i, j, k, ok, loop, nrt, sizeOfDecomp, dbPrt;
  intvec wt;
  ideal action, basis, G, reduceIdeal;
  matrix decompMx;
  poly actCoeff;
  string str, order;

  dbPrt = printlevel-voice+2;
  dbprint(dbPrt, "LinearizeAction " + string(Gaction));
  def RLAR = basering;
  string @mPoly = 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("minpoly = number(" + @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, embedid, groupid;
dbprint(dbPrt+1, "
// 'LinearizeAction' created a new ring.
// To see the ring, type (if the name 'R' was assigned to the return value):
     show(R);
// To access the new action and the equivariant embedding, type
     setring R; actionid; embedid; groupid
");
  setring altring;
  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(def Gaction, int nrs)
"USAGE:   LinearActionQ(action,nrs); ideal action, int nrs
PURPOSE: check whether 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);
  LinearActionQ(Gaction, 8);
}

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

proc LinearCombinationQ(ideal I, poly f)
"USAGE:   LinearCombination(I, f); ideal I, poly f
PURPOSE: test whether 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
PURPOSE: 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:  ring R; this ring comes with the ideals 'invars' and 'groupid' and
         with the poly 'newA':
         - 'invars' contains the algebra generators of the invariant ring
         - 'groupid' is the ideal of G in the new ring
         - 'newA' is the new representation of the primitive root of the
           minimal polynomial of the ring which was active when calling the
           procedure (if the minpoly did not change, 'newA' is set to 'a').
NOTE:    the minimal polynomial of the output ring depends on some random
         choices
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 nullcone of 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
  // ground field.

  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 !
        noReynolds = 0;
        def RORN = ReynoldsOperator(groupid, action, primaryDec);
        setring RORN;
        ideal groupid = std(id);
        attrib(groupid, "isSB", 1);
        ideal action = actionid;
        setring RIRR;
        string parName, minPoly;
        if(npars(basering) == 0) {
          parName = "a";
          minPoly = "0";
        }
        else {
          parName = parstr(basering);
          minPoly = string(minpoly);
        }
        ring RA1 = create_ring(0, "(" + varstr(basering) + "," + parName + ")", "lp");
        if (minPoly!="0") { execute("ideal mpoly = std(" + minPoly + ");"); }
        ideal I = imap(RIRR,invarsGens);
        setring RORN;
        map Phi = RA1, maxideal(1);
        Phi[nvars(RORN) + 1] = newA;
        ideal invarsGens = Phi(I);
        kill Phi,RA1,RIRR;
// end of ersetzt durch

      }
      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 ");}
  }
  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;
  if (defined(ROelements)) { kill ROelements,actionid,theZeroset,id; }
  export invars;
  dbprint(dbPrt+1, "
// 'InvariantRing' created a new ring.
// To see the ring, type (if the name 'R' was assigned to the return value):
     show(R);
// To access the generators of the invariant ring type
     setring R; invars;
// Note that the input group G is stored 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'
");
  setring IRB;
  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,def action)
"USAGE:   InvariantQ(f, G, action); poly f; ideal G, action
PURPOSE: check whether 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
"
{
  def altring=basering;
  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.
PURPOSE: 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,def action)
"USAGE:   NullCone(G, action); ideal G, action
PURPOSE: 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 in general not invariant
EXAMPLE: example NullCone; shows an example
"
{
  int i, nt, dbPrt, offset, groupVars;
  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

  string @minPoly = string(minpoly);
  offset =  size(G) + nt;
  execute(ringSTR);
  execute("poly aaa = number(" + @minPoly + ");");
  if (aaa!=0) { minpoly = number(aaa); }
  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
PURPOSE: compute the Reynolds operator of the group G which acts via 'action'
RETURN:  polynomial ring R over a simple extension of the ground field 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 ideals, 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
"
{
  def ROBR = basering;
  int i, j, n, ns, primaryDec;
  ideal G1 = Grp;
  list solution, saction;
  string str;

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

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

  def RORR = zeroSet(G1, primaryDec);
  setring ROBR;
  string parName, minPoly;
  if(npars(basering) == 0) {
    parName = "a";
    minPoly = "0";
  }
  else {
    parName = parstr(basering);
    minPoly = string(minpoly);
  }
  ring RA1 = create_ring(0, "(" + varstr(basering) + "," + parName + ")", "lp");
  if (minPoly!="0") { execute("ideal mpoly = std(" + minPoly + ");"); }
  ideal Grp = imap(ROBR,Grp);
  ideal Gaction = imap(ROBR,Gaction);
  setring RORR;
  map Phi = RA1, maxideal(1);
  Phi[nvars(RORR) + 1] = newA;
  id = Phi(Grp); // id already defined by zeroSet of level 0
  ideal actionid = Phi(Gaction);
  kill parName,minPoly,Phi,RA1;
// end of ersetzt durch
  list ROelements;
  ideal Rf;
  map groupElem;
  poly h1, h2;

  for(i = 1; i <= size(theZeroset); i++) {
    groupElem = theZeroset[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, ROelements;
  setring ROBR;
  return(RORR);
}

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

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

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

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

static proc SimplifyCoefficientMatrix(matrix coefMatrix)
"USAGE:   SimplifyCoefficientMatrix(M); M matrix coming from coef(...)
PURPOSE: 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.th.
       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);
}

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

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