Index: Singular/LIB/all.lib
===================================================================
--- Singular/LIB/all.lib	(revision d268a947ea45dfab74cae98ecb5b6e00a67224bf)
+++ Singular/LIB/all.lib	(revision 4e461ff1df145d7679c725489c921162884348f8)
@@ -1,6 +1,6 @@
-// $Id: all.lib,v 1.27 2000-09-14 12:37:08 Singular Exp $
+// $Id: all.lib,v 1.28 2000-12-11 13:14:08 Singular Exp $
 ///////////////////////////////////////////////////////////////////////////////
 
-version="$Id: all.lib,v 1.27 2000-09-14 12:37:08 Singular Exp $";
+version="$Id: all.lib,v 1.28 2000-12-11 13:14:08 Singular Exp $";
 info="
 LIBRARY:  all.lib   Load all libraries
@@ -35,4 +35,6 @@
   spcurve.lib   PROCEDURES FOR CM CODIMENSION 2 SINGULARITIES
   algebra.lib   PROCEDURES FOR COMPUTING WITH ALGBRAS AND MAPS
+  rinvar.lib    PROCEDURES FOR INVARIANT RINGS OF REDUCTIVE GROUPS
+  zeroset.lib   PROCEDURES FOR ROOTS AND FACTORIZATION
 ";
 
@@ -72,3 +74,5 @@
 LIB "toric.lib";
 LIB "mregular.lib"
+LIB "rinvar.lib"
+LIB "zeroset.lib";
 
Index: Singular/LIB/rinvar.lib
===================================================================
--- Singular/LIB/rinvar.lib	(revision 4e461ff1df145d7679c725489c921162884348f8)
+++ Singular/LIB/rinvar.lib	(revision 4e461ff1df145d7679c725489c921162884348f8)
@@ -0,0 +1,1103 @@
+// rinvar.lib
+// Invariants of reductive groups
+//
+// Implementation by : Thomas Bayer
+// Current Adress:
+// Institut fuer Informatik, Technische Universitaet Muenchen
+// www:    http://wwwmayr.informatik.tu-muenchen.de/personen/bayert/
+// email : bayert@in.tum.de
+//
+// Last change 10.12.2000
+//
+// 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"
+// Arbeitsgruppe Algebraische Geometrie, Fachbereich Mathematik,
+// Universitaet Kaiserslautern
+///////////////////////////////////////////////////////////////////////////////
+
+version="Id: rinvar.lib,v 1.0 2000/12/10 17:32:15 Singular Exp $";
+info="
+LIBRARY:  rinvar.lib      PROCEDURES FOR INVARIANT RINGS OF REDUCTIVE GROUPS
+
+AUTHOR:   Thomas Bayer,    email: tbayer@in.tum.de
+
+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
+";
+
+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);
+}
+///////////////////////////////////////////////////////////////////////////////
Index: Singular/LIB/zeroset.lib
===================================================================
--- Singular/LIB/zeroset.lib	(revision 4e461ff1df145d7679c725489c921162884348f8)
+++ Singular/LIB/zeroset.lib	(revision 4e461ff1df145d7679c725489c921162884348f8)
@@ -0,0 +1,1532 @@
+// Zeroset.lib
+// Algorithms for finding the zero-set of a zero-dim. ideal in Q(a)[x_1,..x_n],
+// Roots and Factorization of univariate polynomials over Q(a)[t]
+// where a is an algebric number
+//
+// Implementation by : Thomas Bayer
+// Current Adress:
+// Institut fuer Informatik, Technische Universitaet Muenchen
+// www:    http://wwwmayr.informatik.tu-muenchen.de/personen/bayert/
+// email : bayert@in.tum.de
+//
+// Last change 10.12.2000
+//
+// 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"
+// Arbeitsgruppe Algebraische Geometrie, Fachbereich Mathematik,
+// Universitaet Kaiserslautern
+///////////////////////////////////////////////////////////////////////////////
+
+version="Id: zeroset.lib,v 1.0 2000/05/19 12:32:15 Singular Exp $";
+info="
+LIBRARY:  zeroset.lib      PROCEDURES FOR ROOTS AND FACTORIZATION
+
+AUTHOR:   Thomas Bayer,    email: tbayer@mathematik.uni-kl.de
+
+PROCEDURES:
+ EGCD(f, g)    gcd over an algebraic extension field of Q
+ Factor(f)    factorization of f over an algebraic extension field
+ InvertNumber(c)  inverts an element of an algenraic extension field
+ LinearZeroSet(I)  find the linear (partial) solution of I
+ Quotient(f, g)    quotient q  of f w.r.t. g (in f = q*g + remainder)
+ Remainder(f,g)    remainder of the division of f by g
+ Roots(f)    computes all roots of f in an extension field of Q
+ SQFRNorm(f)    norm of f (f must be squarefree)
+ ZeroSet(I)    zero-set of the 0-dim. ideal I
+
+ SUBPROCEDURES:
+ EGCDMain(f, g)    gcd over an algebraic extension field of Q
+ FactorMain(f)    factorization of f over an algebraic extension field
+ InvertNumberMain(c)  inverts an element of an algenraic extension field
+ QuotientMain(f, g)  quotient of f w.r.t. g
+ RemainderMain(f,g)  remainder of the division of f by g
+ RootsMain(f)    computes all roots of f, might extend the groundfield
+ SQFRNormMain(f)  norm of f (f must be squarefree)
+ ContainedQ(data, f)  f in data ?
+ SameQ(a, b)    a == b (list a,b)
+ TransferRing()    create a new basering (to use ...Main())
+ NewBaseRing()    create a new basering (return from  ...Main())
+ SimplifyData(data)  reduces entries of data w.r.t. mpoly
+ SimplifyPoly(f)  reduces coefficients of f w.r.t. mpoly
+ SimplifySolution(sol)  reduces the entries of sol w.r.t. the ideal mpoly
+
+ NOTE: this library is meant as a preliminary extension of the functionality
+ of Singular for univariate factorization of polynomials over simple algebraic
+ extensions in characteristic 0.
+ Subprocedures with postfix 'Main' require that the ring contains a variable
+ 'a' and no parameters, and the ideal 'mpoly', where 'minpoly' from the
+ basering is stored.
+";
+
+LIB "primitiv.lib";
+LIB "primdec.lib";
+
+// note : return a ring : ring need not be exported !!!
+
+// Artihmetic in Q(a)[x] without built-in procedures
+// assume basering = Q[x,a] and minpoly is represented by mpoly(a).
+// the algorithms are taken from "Polynomial Algorithms in Computer Algebra",
+// F. Winkler, Springer Verlag Wien, 1996.
+
+
+// To do :
+// squarefree factorization
+// multiplicities
+
+// Improvement :
+// a main problem is the growth of the coefficients. Try Roots(x7 - 1)
+// retrurn ideal mpoly !
+// mpoly is not monic, comes from primitive_extra
+
+// IMPLEMENTATION
+//
+// In procedures with name 'proc-name'Main a polynomial ring over a simple
+// extension field is represented as Q[x...,a] together with the ideal
+// 'mpoly' (attribute "isSB"). The arithmetic in the extension field is
+// implemented in the procedures in the procedures 'MultPolys' (multiplication)
+// and 'InvertNumber' (inversion). After addition and substraction one should
+// apply 'SimplifyPoly' to the result to reduce the result w.r.t. 'mpoly'.
+// This is done by reducing each coefficient seperately, which is more
+// efficient for polynomials with many terms.
+
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc Roots(poly f)
+"USAGE:   Roots(f); poly f
+PUROPSE: compute all roots of f in a finite extension of the groundfield
+         without multiplicities.
+RETURN:  ring, a polynomial ring over an extension field of the ground field,
+         containing a list 'roots', poly 'newA', poly 'f':
+   - 'roots' is the list of roots of the polynomial f (no multiplicities)
+         - if the groundfield is Q(a') and the extension field is Q(a), then
+     'newA' is the representation of a' in Q(a). If the basering
+           contains a parameter 'a' and the minpoly remains unchanged then
+           'newA' = 'a'. If the basering does not contain a parameter then
+           'newA' = 'a' (default).
+   - 'f' is the polynomial f in  Q(a) (a' being substituted by newA)
+ASSUME:  ground field to be Q or a simple extension of Q given by a minpoly
+EXAMPLE: example  Roots; shows an example
+"
+{
+  int dbPrt = printlevel-voice+3;
+
+  // create a new ring where par(1) is replaced by the variable
+  // with the same name or, if basering does not contain a parameter,
+  // with a new variable 'a'.
+
+  def ROB = basering;
+  def ROR = TransferRing(basering);
+  setring ROR;
+  export(ROR);
+
+  // get the polynomial f and find the roots
+
+  poly f = imap(ROB, f);
+  list result = RootsMain(f);  // find roots of f
+
+  // store the roots and the the new representation of 'a' and transform
+  // the coefficients of f.
+
+  list roots = result[1];
+  poly newA = result[2];
+  map F = basering, maxideal(1);
+  F[nvars(basering)] = newA;
+  poly fn = SimplifyPoly(F(f));
+
+  // create a new ring with minploy = mpoly[1] (from ROR)
+
+  def RON = NewBaseRing();
+  setring(RON);
+  list roots = imap(ROR, roots);
+  poly newA = imap(ROR, newA);
+  poly f = imap(ROR, fn);
+  kill(ROR);
+  export(roots);
+  export(newA);
+  export(f); dbprint(dbPrt,"
+// 'Roots' created a new ring which contains the list 'roots' and
+// the polynomials 'f' and 'newA'
+// To access the roots, newA and the new representation of f, type
+   def R = Roots(f); setring R; roots; newA; f;
+");
+  return(RON);
+}
+example
+{"EXAMPLE:";  echo = 2;
+  ring R = (0,a), x, lp;
+  minpoly = a2+1;
+  poly f = x3 - a;
+  def R1 = Roots(f);
+  setring R1;
+  minpoly;
+  newA;
+  f;
+  roots;
+  map F;
+  F[1] = roots[1];
+  F(f);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc RootsMain(poly f)
+"USAGE:   RootsMain(f); poly f
+PUROPSE: compute all roots of f in a finite extension of the groundfield
+         without multiplicities.
+RETURN:  list, all entries are polynomials
+  _[1] = roots of f, each entry is a polynomial
+  _[2] = 'newA'  - if the groundfield is Q(a') and the extension field is Q(a), then
+         'newA' is the representation of a' in Q(a)
+  _[3] = minpoly of the algebraic extension of the groundfield
+ASSUME:  basering = Q[x,a]
+   ideal mpoly must be defined, it might be 0 !
+NOTE: might change teh ideal mpoly !!
+EXAMPLE: example  Roots; shows an example
+"
+{
+  int i, linFactors, nlinFactors, dbPrt;
+  intvec wt = 1,0;    // deg(a) = 0
+  list factorList, nlFactors, nlMult, roots, result;
+  poly fa, lc;
+
+  dbPrt = printlevel-voice+3;
+
+  // factor f in Q(a)[t] to obtain the roots lying in Q(a)
+  // firstly, find roots of the linear factors,
+  // nonlinear factors are processed later
+
+  dbprint(dbPrt, "Roots of " + string(f) +  ", minimal polynomial = " + string(mpoly[1]));
+  factorList = FactorMain(f);          // Factorize f
+  dbprint(dbPrt, (" prime factors of f are : " + string(factorList[1])));
+
+  linFactors = 0;
+  nlinFactors = 0;
+  for(i = 2; i <= size(factorList[1]); i = i + 1) {  // find linear and nonlinear factors
+    fa = factorList[1][i];
+    if(deg(fa, wt) == 1) {
+      linFactors++;        // get the root from the linear factor
+      lc = LeadTerm(fa, 1)[3];
+      fa = MultPolys(InvertNumberMain(lc), fa); // make factor monic
+      roots[linFactors] = var(1) - fa;  // fa is monic !!
+    }
+    else {            // ignore nonlinear factors
+      nlinFactors++;
+      nlFactors[nlinFactors] = factorList[1][i];
+      nlMult[nlinFactors] = factorList[2][i];
+    }
+  }
+  if(linFactors == size(factorList[1]) - 1) {    // all roots of f are contained in the groundfield
+    result[1] = roots;
+    result[2] = var(2);
+    result[3] = mpoly[1];
+    return(result);
+  }
+
+  // process the nonlinear factors, i.e., extend the groundfield
+  // where a nonlinear factor (irreducible) is a minimal polynomial
+  // compute the primitive element of this extension
+
+  ideal primElem, minPolys, Fid;
+  list partSol;
+  map F, Xchange;
+  poly f1, newA, mp, oldMinPoly;
+
+  Fid = mpoly;
+  F[1] = var(1);
+  Xchange[1] = var(2);      // the variables have to be exchanged
+  Xchange[2] = var(1);      // for the use of 'primitive'
+
+  if(nlinFactors == 1) {            // one nl factor
+
+    // compute the roots of the nonlinear (irreducible, monic) factor f1 of f
+    // by extending the basefield by a' with minimal polynomial f1
+    // Then call Roots(f1) to find the roots of f1 over the new base field
+
+    f1 = nlFactors[1];
+    if(mpoly[1] != 0) {
+      mp = mpoly[1];
+      minPolys = Xchange(mp), Xchange(f1);
+      primElem = primitive_extra(minPolys);  // no random coord. change
+      mpoly = std(primElem[1]);
+      F = basering, maxideal(1);
+      F[2] = primElem[2];      // transfer all to the new representation
+      newA = primElem[2];      // new representation of a
+      f1 = SimplifyPoly(F(f1));     //reduce(F(f1), mpoly);
+      if(size(roots) > 0) {roots = SimplifyData(F(roots));}
+    }
+    else {
+      mpoly = std(Xchange(f1));
+      newA = var(2);
+    }
+    result[3] = mpoly[1];
+    oldMinPoly = mpoly[1];
+    partSol = RootsMain(f1);    // find roots of f1 over extended field
+
+    if(oldMinPoly != partSol[3]) {    // minpoly has changed ?
+      // all previously computed roots must be transformed
+      // because the minpoly has changed
+      result[3] = partSol[3];    // new minpoly
+      F[2] = partSol[2];    // new representation of algebraic number
+      if(size(roots) > 0) {roots = SimplifyData(F(roots)); }
+      newA = SimplifyPoly(F(newA)); // F(newA);
+    }
+    roots = roots + partSol[1];  // add roots
+    result[2] = newA;
+    result[1] = roots;
+  }
+  else {  // more than one nonlinear (irreducible) factor (f_1,...,f_r)
+    // solve each of them by RootsMain(f_i), append their roots
+    // change the minpoly and transform all previously computed
+    // roots if necessary.
+    // Note that the for-loop is more or less book-keeping
+
+    newA = var(2);
+    result[2] = newA;
+    for(i = 1; i <= size(nlFactors); i = i + 1) {
+      oldMinPoly = mpoly[1];
+      partSol = RootsMain(nlFactors[i]);    // main work
+      nlFactors[i] = 0;        // delete factor
+      result[3] = partSol[3];        // store minpoly
+
+      // book-keeping starts here as in the case 1 nonlinear factor
+
+      if(oldMinPoly != partSol[3]) { // minpoly has changed
+        F = basering, maxideal(1);
+        F[2] = partSol[2];    // transfer all to the new representation
+        newA = SimplifyPoly(F(newA));    // F(newA); new representation of a
+        result[2] = newA;
+        if(i < size(nlFactors)) {
+          nlFactors = SimplifyData(F(nlFactors));
+        } // transform remaining factors
+        if(size(roots) > 0) {roots = SimplifyData(F(roots));}
+      }
+      roots = roots + partSol[1];    // transform roots
+      result[1] = roots;
+    }  // end more than one nl factor
+
+  }
+  return(result);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc ZeroSet(ideal I, list #)
+"USAGE:   ZeroSetMain(ideal I [, int opt]); ideal I, int opt
+PUROPSE: compute the zero-set of the zero-dim. ideal I, in a finite extension
+         of the groundfield.
+RETURN:  ring, a polynomial ring over an extension field of the ground field,
+         containing a list 'zeroset', a polynomial 'newA', and an ideal 'id'.
+   - 'zeroset'  is list of the zeros of the ideal 'I', each zero is
+           an ideal.
+         - if the groundfield is Q(a') and the extension field is Q(a), then
+     'newA' is the representation of a' in Q(a).  If the basering
+           contains a parameter 'a' and the minpoly remains unchanged then
+           'newA' = 'a'. If the basering does not contain a parameter then
+           'newA' = 'a' (default).
+   - 'id' is the ideal 'I'  in Q(a)[x_1,...] (a' substituted by 'newA')
+ASSUME:  dim(I) = 0, and ground field to be Q or a simple extension of Q given
+         by a minpoly.
+OPTIONS: opt = 0 no primary decomposition (default)
+   opt > 0 primary decomposition
+NOTE:    If I contains an algebraic number (parameter) then 'I' must be transformed
+         w.r.t. 'newA' in the new ring.
+EXAMPLE: example ZeroSet; shows an example
+"
+{
+  int primaryDecQ, dbPrt;
+  list rp;
+
+  dbPrt = printlevel-voice+2;
+
+  if(size(#) > 0) { primaryDecQ = #[1]; }
+  else { primaryDecQ = 0; }
+
+  // create a new ring 'ZSR' with one additional variable instead of the parameter
+  // if the basering does not contain a parameter then 'a' is used as the additional variable.
+
+  def RZSB = basering;
+  def ZSR = TransferRing(basering);
+  setring ZSR;
+  export(ZSR);
+
+  // get ideal I and find the zero-set
+
+  ideal id = std(imap(RZSB, I));
+        print(dim(id));
+        if(dim(id) > 1) {       // new variable adjoined to ZSR
+                ERROR(" ideal not zerodimensional ");
+        }
+
+  list result = ZeroSetMain(id, primaryDecQ);
+
+  // store the zero-set, minimal polynomial and the new representative of 'a'
+
+  list zeroset = result[1];
+  poly newA = result[2];
+  poly minPoly = result[3][1];
+
+  // transform the generators of the ideal I w.r.t. the new representation of 'a'
+
+  map F = basering, maxideal(1);
+  F[nvars(basering)] = newA;
+  id = SimplifyData(F(id));
+
+  // create a new ring with minpoly = minPoly
+
+  def RZBN = NewBaseRing();
+  setring RZBN;
+
+  list zeroset = imap(ZSR, zeroset);
+  poly newA = imap(ZSR, newA);
+  ideal id = imap(ZSR, id);
+  kill(ZSR);
+
+  export(id);
+  export(zeroset);
+  export(newA);
+    dbprint(dbPrt,"
+// 'ZeroSet' created a new ring which contains the list 'zeroset', the ideal
+// 'id' and the polynomial 'newA'. 'id' is the ideal of the input transformed
+// w.r.t. 'newA'.
+// To access the zero-set, 'newA' and the new representation of the ideal, type
+   def R = ZeroSet(I); setring R; zeroset; newA; id;
+");
+  return(RZBN);
+}
+example
+{"EXAMPLE:";  echo = 2;
+  ring R = (0,a), (x,y,z), lp;
+  minpoly = a2 + 1;
+  ideal I = x2 - 1/2, a*z - 1, y - 2;
+  def T = ZeroSet(I);
+  setring T;
+  minpoly;
+  newA;
+  id;
+  zeroset;
+  map F1 = basering, zeroset[1];
+  map F2 = basering, zeroset[2];
+  F1(id);
+  F2(id);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc InvertNumberMain(poly f)
+"USAGE:   InvertNumberMain(f); poly f
+PURPOSE: compute 1/f if f is a number in Q(a) i.e., f is represented by a
+         polynomial in Q[a].
+RETURN:  poly 1/f
+ASSUME:  basering = Q[x_1,...,x_n,a], ideal 'mpoly' must be defined and != 0 !
+"
+{
+  if(diff(f, var(1)) != 0) { ERROR("number must not contain variable !");}
+
+  int n = nvars(basering);
+  def RINB = basering;
+  string ringSTR = "ring RINR = 0, " + string(var(n)) + ", dp;";
+  execute(ringSTR);        // new ring = Q[a]
+
+  list gcdList;
+  poly f, g, inv;
+
+  f = imap(RINB, f);
+  g = imap(RINB, mpoly)[1];
+
+  if(diff(f, var(1)) != 0) { inv = extgcd(f, g)[2]; }  // f contains var(1)
+  else {  inv = 1/f;}          // f element in Q
+
+  setring(RINB);
+  return(imap(RINR, inv));
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc MultPolys(poly f, poly g)
+"USAGE:   MultPolys(f, g); poly f,g
+PURPOSE: multiply the polynomials f and g and reduce them w.r.t. mpoly
+RETURN:  poly f*g
+ASSUME:  basering = Q[x,a], ideal mpoly must be defined, it might be 0 !
+"
+{
+  return(SimplifyPoly(f * g));
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc LeadTerm(poly f, int i)
+"USAGE:   LeadTerm(f); poly f, int i
+PUROPSE: compute the leading coef and term of f w.r.t var(i), where the last
+         ring variable is treated as a parameter.
+RETURN:  list of polynomials
+         _[1] = leading term
+         _[2] = leading monomial
+         _[3] = leading coefficient
+ASSUME:  basering = Q[x_1,...,x_n,a]
+"
+{
+  list result;
+  matrix co = coef(f, var(i));
+  result[1] = co[1, 1]*co[2, 1];
+  result[2] = co[1, 1];
+  result[3] = co[2, 1];
+  return(result);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc Quotient(poly f, poly g)
+"USAGE:   Quotient(f, g); poly f, g;
+PUROPSE: compute the quotient q and remainder r s.t. f = g*q + r, deg(r) < g
+RETURN:  list of polynomials
+         _[1] = quotient  q
+         _[2] = remainder r
+ASSUME:  basering = Q[x] or Q(a)[x]
+EXAMPLE: example  Quotient; shows an example
+"
+{
+  def QUOB = basering;
+  def QUOR = TransferRing(basering);  // new ring with parameter 'a' replaced by a variable
+  setring QUOR;
+  export(QUOR);
+  poly f = imap(QUOB, f);
+  poly g = imap(QUOB, g);
+  list result = QuotientMain(f, g);
+
+  setring(QUOB);
+  list result = imap(QUOR, result);
+  kill(QUOR);
+  return(result);
+}
+example
+{"EXAMPLE:";  echo = 2;
+ ring R = (0,a), x, lp;
+ minpoly = a2+1;
+ poly f =  x4 - 2;
+ poly g = x - a;
+ list qr = Quotient(f, g);
+ qr;
+ qr[1]*g + qr[2] - f;
+}
+
+proc QuotientMain(poly f, poly g)
+"USAGE:   QuotientMain(f, g); poly f, g
+PUROPSE: compute the quotient q and remainder r s.t. f = g*q + r, deg(r) < g
+RETURN:  list of polynomials
+         _[1] = quotient  q
+         _[2] = remainder r
+ASSUME:  basering = Q[x,a] and ideal 'mpoly' is defined (it might be 0),
+         this represents the ring Q(a)[x] together with 'minpoly'.
+EXAMPLE: example  Quotient; shows an example
+"
+{
+  if(g == 0) { ERROR("Division by zero !");}
+
+  def QMB = basering;
+  def QMR = NewBaseRing();
+  setring QMR;
+  poly f, g, h;
+  h = imap(QMB, f) / imap(QMB, g);
+  setring QMB;
+  return(list(imap(QMR, h), 0));
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc Remainder(poly f, poly g)
+"USAGE:   Remainder(f, g); poly f, g
+PUROPSE: compute the remainder of the division of f by g, i.e. a polynomial r
+         s.t. f = g*q + r, deg(r) < g.
+RETURN:  poly
+ASSUME:  basering = Q[x] or Q(a)[x]
+"
+{
+  def REMB = basering;
+  def REMR = TransferRing(basering);  // new ring with parameter 'a' replaced by a variable
+  setring(REMR);
+  export(REMR);
+  poly f = imap(REMB, f);
+  poly g = imap(REMB, g);
+  poly h = RemainderMain(f, g);
+
+  setring(REMB);
+  poly r = imap(REMR, h);
+  kill(REMR);
+  return(r);
+}
+example
+{"EXAMPLE:";  echo = 2;
+ ring R = (0,a), x, lp;
+ minpoly = a2+1;
+ poly f =  x4 - 1;
+ poly g = x3 - 1;
+ Remainder(f, g);
+}
+
+proc RemainderMain(poly f, poly g)
+"USAGE:   RemainderMain(f, g); poly f, g
+PUROPSE: compute the remainder r s.t. f = g*q + r, deg(r) < g
+RETURN:  poly
+ASSUME:  basering = Q[x,a] and ideal 'mpoly' is defined (it might be 0),
+         this represents the ring Q(a)[x] together with 'minpoly'.
+"
+{
+  int dg;
+  intvec wt = 1,0;;
+  poly lc, g1, r;
+
+  if(deg(g, wt) == 0) { return(0); }
+
+  lc = LeadTerm(g, 1)[3];
+  g1 = MultPolys(InvertNumberMain(lc), g);  // make g monic
+
+  return(SimplifyPoly(reduce(f, std(g1))));
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc EGCD(poly f, poly g)
+"USAGE:   EGCDMain(f, g); poly f, g
+PUROPSE: compute the polynomial gcd of f and g over Q(a)[x]
+RETURN:  poly h s.t. h is a greatest common divisor of f and g (not nec. monic)
+ASSUME:  basering = Q(a)[t]
+EXAMPLE: example  EGCD; shows an example
+"
+{
+  def GCDB = basering;
+  def GCDR = TransferRing(basering);  // new ring with parameter 'a' replaced by a variable
+  setring GCDR;
+  export(GCDR);
+  poly f = imap(GCDB, f);
+  poly g = imap(GCDB, g);
+  poly h = EGCDMain(f, g);  // squarefree norm of f
+
+  setring(GCDB);
+  poly h = imap(GCDR, h);
+  kill(GCDR);
+  return(h);
+}
+example
+{"EXAMPLE:";  echo = 2;
+ ring R = (0,a), x, lp;
+ minpoly = a2+1;
+ poly f =  x4 - 1;
+ poly g = x2 - 2*a*x - 1;
+ EGCD(f, g);
+}
+
+proc EGCDMain(poly f, poly g)
+"USAGE:   EGCDMain(f, g); poly f, g
+PUROPSE: compute the polynomial gcd of f and g over Q(a)[x]
+RETURN:  poly
+ASSUME:  basering = Q[x,a] and ideal 'mpoly' is defined (it might be 0),
+         this represents the ring Q(a)[x] together with 'minpoly'.
+EXAMPLE: example EGCD; shows an example
+"
+// might be extended to return s1, s2 s.t. f*s1 + g*s2 = gcd
+{
+  int i = 1;
+  poly r1, r2, r;
+
+  r1 = f;
+  r2 = g;
+
+  while(r2 != 0) {
+    r  = RemainderMain(r1, r2);
+    r1 = r2;
+    r2 = r;
+  }
+  return(r1);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc MEGCD(poly f, poly g, int varIndex)
+"USAGE:   MEGCD(f, g, i); poly f, g; int i
+PUROPSE: compute  the polynomial gcd of f and g in the i'th variable
+RETURN:  poly
+ASSUME:  f, g are polynomials in var(i), last variable is the algebraic number
+EXAMPLE: example  MEGCD; shows an example
+"
+// might be extended to return s1, s2 s.t. f*s1 + g*s2 = gc
+// not used !
+{
+  string @str, @sf, @sg, @mp, @parName;
+
+  def @RGCDB = basering;
+
+  @sf = string(f);
+  @sg = string(g);
+  @mp = string(minpoly);
+
+  if(npars(basering) == 0) { @parName = "0";}
+  else { @parName = "(0, " + parstr(basering) + ")"; }
+  @str = "ring @RGCD = " + @parName + ", " + string(var(varIndex)) + ", dp;";
+  execute(@str);
+  if(@mp != "0") { execute ("minpoly = " + @mp + ";"); }
+  execute("poly @f = " + @sf + ";");
+  execute("poly @g = " + @sg + ";");
+  export(@RGCD);
+  poly @h = EGCD(@f, @g);
+  setring(@RGCDB);
+  poly h = imap(@RGCD, @h);
+  kill(@RGCD);
+  return(h);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc SQFRNorm(poly f)
+"USAGE:   SQFRNorm(f); poly f
+PUROPSE: compute the norm of the squarefree polynomial f in Q(a)[x].
+RETURN:  list
+         _[1] = squarefree norm of g (poly)
+         _[2] = g (= f(x - s*a)) (poly)
+         _[3] = s (int)
+ASSUME:  f must be squarefree, basering = Q(a)[x] and minpoly != 0.
+NOTE:    the norm is an element of Q[x]
+EXAMPLE: example  SQFRNorm; shows an example
+"
+{
+  def SNB = basering;
+  def SNR = TransferRing(basering);  // new ring with parameter 'a' replaced by a variable
+  setring SNR;
+  export(SNR);
+  poly f = imap(SNB, f);
+  list result = SQFRNormMain(f);  // squarefree norm of f
+
+  setring(SNB);
+  list result = imap(SNR, result);
+  kill(SNR);
+  return(result);
+}
+example
+{"EXAMPLE:";  echo = 2;
+   ring R = (0,a), x, lp;
+   minpoly = a2+1;
+  poly f =  x4 - 2*x + 1;
+  SQFRNorm(f);
+}
+
+proc SQFRNormMain(poly f)
+"USAGE:   SQFRNorm(f); poly f
+PUROPSE: compute the norm of the squarefree polynomial f in Q(a)[x].
+RETURN:  list
+         _[1] = squarefree norm of g (poly)
+         _[2] = g (= f(x - s*a)) (poly)
+         _[3] = s (int)
+ASSUME:  f must be squarefree, basering = Q[x,a] and ideal 'mpoly' is equal to 'minpoly',
+         this represents the ring Q(a)[x] together with 'minpoly'.
+NOTE:   the norm is an element of Q[x]
+EXAMPLE: example  SqfrNorm; shows an example
+"
+{
+  int s = 0;
+  intvec wt = 1,0;
+  ideal mapId;
+  // list result;
+  poly g, N, N1, h;
+  string ringSTR;
+
+  mapId[1] = var(1) - var(2);    // linear transformation
+  mapId[2] = var(2);
+  map Fs = basering, mapId;
+
+  N = resultant(f, mpoly[1], var(2));  // norm of f
+  N1 = diff(N, var(1));
+  g = f;
+
+  ringSTR = "ring SNRM1 = 0, " + string(var(1)) + ", dp;";  // univariate ring
+  def SNRMB = basering;
+  execute(ringSTR);
+  poly N, N1, h;        // N, N1 do not contain 'a', use built-in gcd
+  h = gcd(imap(SNRMB, N), imap(SNRMB, N1));
+  setring(SNRMB);
+  h = imap(SNRM1, h);
+  while(deg(h, wt) != 0) {    // while norm is not squarefree
+    s = s + 1;
+    g = reduce(Fs(g), mpoly);
+    N = reduce(resultant(g, mpoly[1], var(2)), mpoly);  // norm of g
+    N1 = reduce(diff(N, var(1)), mpoly);
+    setring(SNRM1);
+    h = gcd(imap(SNRMB, N), imap(SNRMB, N1));
+    setring(SNRMB);
+    h = imap(SNRM1, h);
+  }
+  return(list(N, g, s));
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc Factor(poly f)
+"USAGE:   Factor(f); poly f
+PUROPSE: compute the factorization of the squarefree poly f over Q(a)[t], minpoly  = p(a).
+RETURN:  list
+         _[1] = factors (monic), first entry is the leading coefficient
+         _[2] = multiplicities (not yet)
+ASSUME:  basering must be the univariate polynomial ring over a field which is Q
+         or a simple extension of Q given by a minpoly.
+NOTE:    if basering = Q[t] then this is the built-in 'factorize'.
+EXAMPLE: example  Factor; shows an example
+"
+{
+  def FRB = basering;
+  def FRR = TransferRing(basering);  // new ring with parameter 'a' replaced by a variable
+  setring FRR;
+  export(FRR);
+  poly f = imap(FRB, f);
+  list result = FactorMain(f);  // squarefree norm of f
+
+  setring(FRB);
+  list result = imap(FRR, result);
+  kill(FRR);
+  return(result);
+}
+example
+{"EXAMPLE:";  echo = 2;
+ ring R = (0,a), x, lp;
+ minpoly = a2+1;
+ poly f =  x4 - 1;
+ list fl = Factor(f);
+ fl;
+ fl[1][1]*fl[1][2]*fl[1][3]*fl[1][4]*fl[1][5] - f;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc FactorMain(poly f)
+"USAGE:   FactorMain(f); poly f
+PUROPSE: compute the factorization of the squarefree poly f over Q(a)[t], minpoly  = p(a).
+RETURN:  list
+         _[1] = factors, first is a constant
+         _[2] = multiplicities (not yet)
+ASSUME:  basering = Q[x,a], represents Q(a)[x] and minpoly.
+   ideal mpoly must be defined, it might be 0 !
+EXAMPLE: example  Factor; shows an example
+"
+// extend this by a squarefree factorization !!
+// multiplicities are not valid !!
+{
+  int i, s;
+  list normList, factorList, quo_rem;
+  poly f1, h, h1, H, g, leadCoef, invCoeff;
+  ideal fac1, fac2;
+  map F;
+
+  // if no minimal polynomial is defined then use 'factorize'
+  // FactorOverQ is wrapped around 'factorize'
+
+  if(mpoly[1] == 0) {
+    // print(" factorize : deg = " + string(deg(f, intvec(1,0))));
+    factorList = factorize(f); // FactorOverQ(f);
+    return(factorList);
+  }
+
+  // if mpoly != 0 and f does not contain the algebraic number, a root of
+  // f might be contained in Q(a). Hence one must not use 'factorize'.
+
+  fac1[1] = 1;
+  fac2[1] = 1;
+  normList = SQFRNormMain(f);
+  // print(" factorize : deg = " + string(deg(normList[1], intvec(1,0))));
+  factorList = factorize(normList[1]);     // factor squarefree norm of f over Q[x]
+  g = normList[2];
+  s = normList[3];
+  F[1] = var(1) + s*var(2);      // inverse transformation
+  F[2] = var(2);
+  fac1[1] = factorList[1][1];
+  fac2[1] = factorList[2][1];
+  for(i = 2; i <= size(factorList[1]); i = i + 1) {
+    H = factorList[1][i];
+    h = EGCDMain(H, g);
+    quo_rem = QuotientMain(g, h);
+    g = quo_rem[1];
+    fac1[i] = SimplifyPoly(F(h));
+    fac2[i] = 1;        // to be changed later
+  }
+  return(list(fac1, fac2));
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc ZeroSetMain(ideal I, int primDecQ)
+"USAGE:   ZeroSetMain(ideal I, int opt); ideal I, int opt
+PUROPSE: compute the zero-set of the zero-dim. ideal I, in a simple extension
+         of the groundfield.
+RETURN:  list
+   - 'f' is the polynomial f in  Q(a) (a' being substituted by newA)
+         _[1] = zero-set (list), is the list of the zero-set of the ideal I,
+                each entry is an ideal.
+         _[2] = 'newA';  if the groundfield is Q(a') and the extension field
+                is Q(a), then 'newA' is the representation of a' in Q(a).
+                If the basering contains a parameter 'a' and the minpoly
+                remains unchanged then 'newA' = 'a'. If the basering does not
+                contain a parameter then 'newA' = 'a' (default).
+         _[3] = 'mpoly' (ideal), the minimal polynomial of the simple extension
+                of the ground field.
+ASSUME: basering = K[x_1,x_2,...,x_n] where K = Q or a simple extension of Q
+        given by a minpoly.
+  dim(I) = 0.
+NOTE:    opt = 0  no primary decomposition
+  opt > 0  use a primary decomposition
+EXAMPLE: example ZeroSet; shows an example
+"
+{
+  // main work is done in ZeroSetMainWork, here the zero-set of each ideal from the
+  // primary decompostion is coputed by menas of ZeroSetMainWork, and then the
+  // minpoly and the parameter representing the algebraic extension are
+  // transformed according to 'newA', i.e., only bookeeping is done.
+
+  int i, j, n, noMP, dbPrt;
+  intvec w;
+  list currentSol, result, idealList, primDecList, zeroSet;
+  ideal J;
+  map Fa;
+  poly newA, oldMinPoly;
+
+  dbPrt = printlevel-voice+2;
+  dbprint(dbPrt, "ZeroSet of " + string(I) + ", minpoly = " + string(minpoly));
+
+  n = nvars(basering) - 1;
+  for(i = 1; i <= n; i++) { w[i] = 1;}
+  w[n + 1] = 0;
+
+  if(primDecQ == 0) { return(ZeroSetMainWork(I, w, 0)); }
+
+  newA = var(n + 1);
+  if(mpoly[1] == 0) { noMP = 1;}
+  else {noMP = 0;}
+
+  primDecList = primdecGTZ(I);      // primary decomposition
+  dbprint(dbPrt, "primary decomposition consists of " + string(size(primDecList)) + " primary ideals ");
+  // idealList = PDSort(idealList);    // high degrees first
+
+  for(i = 1; i <= size(primDecList); i = i + 1) {
+    idealList[i] = primDecList[i][2];  // use prime component
+    dbprint(dbPrt, string(i) + "  " + string(idealList[i]));
+  }
+
+  // compute the zero-set of each primary ideal and join them.
+  // If necessary, change the groundfield and transform the zero-set
+
+  dbprint(dbPrt, "
+find the zero-set of each primary ideal, form the union
+and keep track of the minimal polynomials ");
+
+  for(i = 1; i <= size(idealList); i = i + 1) {
+    J = idealList[i];
+    idealList[i] = 0;
+    oldMinPoly = mpoly[1];
+    dbprint(dbPrt, " ideal#" + string(i) + " of " + string(size(idealList)) + " = " + string(J));
+    currentSol = ZeroSetMainWork(J, w, 0);
+
+    if(oldMinPoly != currentSol[3]) {   // change minpoly and transform solutions
+      dbprint(dbPrt, " change minpoly to " + string(currentSol[3][1]));
+      dbprint(dbPrt, " new representation of algebraic number = " + string(currentSol[2]));
+      if(!noMP) {      // transform the algebraic number a
+        Fa = basering, maxideal(1);
+        Fa[n + 1] = currentSol[2];
+        newA = SimplifyPoly(Fa(newA));  // new representation of a
+        if(size(zeroSet) > 0) {zeroSet = SimplifyZeroset(Fa(zeroSet)); }
+        if(i < size(idealList)) { idealList = SimplifyZeroset(Fa(idealList)); }
+      }
+      else { noMP = 0;}
+    }
+    zeroSet = zeroSet + currentSol[1];    // add new elements
+  }
+  return(list(zeroSet, newA, mpoly));
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc ZeroSetMainWork(ideal id, intvec wt, int sVars)
+"USAGE:   ZeroSetMainWork(I, wt, sVars);
+PUROPSE: compute the zero-set of the zero-dim. ideal I, in a finite extension
+         of the groundfield (without multiplicities).
+RETURN:  list, all entries are polynomials
+         _[1] = zeros, each entry is an ideal
+         _[2] = newA; if the groundfield is Q(a') this is the rep. of a' w.r.t. a
+         _[3] = minpoly of the algebraic extension of the groundfield (ideal)
+         _[4] = name of algebraic number (default = 'a')
+ASSUME:  basering = Q[x_1,x_2,...,x_n,a]
+         ideal mpoly must be defined, it might be 0 !
+NOTE:    might change 'mpoly' !!
+EXAMPLE: example  IdealSolve; shows an example
+"
+{
+  int i, j, k, nrSols, n, noMP;
+  ideal I, generators, gens, solid, partsolid;
+  list linSol, linearSolution, nLinSol, nonlinSolutions, partSol, sol, solutions, result;
+  list linIndex, nlinIndex, index;
+  map Fa, Fsubs;
+  poly oldMinPoly, newA;
+
+  if(mpoly[1] == 0) { noMP = 1;}
+  else { noMP = 0;}
+  n = nvars(basering) - 1;
+  newA = var(n + 1);
+
+  I = std(id);
+
+  // find linear solutions of univariate generators
+
+  linSol = LinearZeroSetMain(I, wt);
+  generators = linSol[3];      // they are a standardbasis
+  linIndex = linSol[2];
+  linearSolution = linSol[1];
+  if(size(linIndex) + sVars == n) {    // all variables solved
+    solid = SubsMapIdeal(linearSolution, linIndex, 0);
+    result[1] = list(solid);
+    result[2] = var(n + 1);
+    result[3] = mpoly;
+    return(result);
+  }
+
+  // find roots of the nonlinear univariate polynomials of generators
+  // if necessary, transform linear solutions w.r.t. newA
+
+  oldMinPoly = mpoly[1];
+  nLinSol =  NonLinearZeroSetMain(generators, wt);    // find solutions of univariate generators
+  nonlinSolutions = nLinSol[1];    // store solutions
+  nlinIndex = nLinSol[4];     // and index of solved variables
+  generators = nLinSol[5];    // new generators
+
+  // change minpoly if necessary and transform the ideal and the partial solutions
+
+  if(oldMinPoly != nLinSol[3]) {
+    newA = nLinSol[2];
+    if(!noMP && size(linearSolution) > 0) {    // transform the algebraic number a
+      Fa = basering, maxideal(1);
+      Fa[n + 1] = newA;
+      linearSolution = SimplifyData(Fa(linearSolution));  // ...
+    }
+  }
+
+  // check if all variables are solved.
+
+  if(size(linIndex) + size(nlinIndex) == n - sVars) {
+    solutions = MergeSolutions(linearSolution, linIndex, nonlinSolutions, nlinIndex, list(), n);
+  }
+
+  else {
+
+  // some variables are not solved.
+  // substitute each partial solution in generators and find the
+  // zero set of the resulting ideal by recursive application
+  // of ZeroSetMainWork !
+
+  index = linIndex + nlinIndex;
+  nrSols = 0;
+  for(i = 1; i <=  size(nonlinSolutions); i = i + 1) {
+    sol = linearSolution + nonlinSolutions[i];
+    solid = SubsMapIdeal(sol, index, 1);
+    Fsubs = basering, solid;
+    gens = std(SimplifyData(Fsubs(generators)));    // substitute partial solution
+    oldMinPoly = mpoly[1];
+    partSol = ZeroSetMainWork(gens, wt, size(index) + sVars);
+
+    if(oldMinPoly != partSol[3]) {    // minpoly has changed
+      Fa = basering, maxideal(1);
+      Fa[n + 1] = partSol[2];    // a -> p(a), representation of a w.r.t. new minpoly
+      newA = reduce(Fa(newA), mpoly);
+      generators = std(SimplifyData(Fa(generators)));
+      if(size(linearSolution) > 0) { linearSolution = SimplifyData(Fa(linearSolution));}
+      if(size(nonlinSolutions) > 0) {
+        nonlinSolutions = SimplifyZeroset(Fa(nonlinSolutions));
+      }
+      sol = linearSolution + nonlinSolutions[i];
+    }
+
+    for(j = 1; j <= size(partSol[1]); j++) {   // for all partial solutions
+      partsolid = partSol[1][j];
+      for(k = 1; k <= size(index); k++) {
+        partsolid[index[k]] = sol[k];
+       }
+      nrSols++;
+      solutions[nrSols] = partsolid;
+    }
+  }
+
+  }  // end else
+  return(list(solutions, newA, mpoly));
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc LinearZeroSetMain(ideal I, intvec wt)
+"USAGE:   LinearZeroSetMain(I, wt)
+PURPOSE: solve the univariate linear polys in I
+ASSUME:  basering = Q[x_1,...,x_n,a]
+RETURN:  list
+         _[1] = partial solution of I
+         _[2] = index of solved vars
+         _[3] = new generators (standardbasis)
+"
+{
+  int i, ok, n, found, nrSols;
+  ideal generators, newGens;
+  list result, index, totalIndex, vars, sol, temp;
+  map F;
+  poly f;
+
+  result[1] = index;      // sol[1] should be the empty list
+  n = nvars(basering) - 1;
+  generators = I;        // might be wrong, use index !
+  ok = 1;
+  nrSols = 0;
+  while(ok) {
+    found = 0;
+    for(i = 1; i <= size(generators); i = i + 1) {
+      f = generators[i];
+      vars = Variables(f, n);
+      if(size(vars) == 1 && deg(f, wt) == 1) {  // univariate,linear
+        nrSols++; found++;
+        index[nrSols] = vars[1];
+        sol[nrSols] = var(vars[1]) - MultPolys(InvertNumberMain(LeadTerm(f, vars[1])[3]), f);
+      }
+    }
+    if(found > 0) {
+      F = basering, SubsMapIdeal(sol, index, 1);
+      newGens = std(SimplifyData(F(generators)));    // substitute, simplify alg. number
+      if(size(newGens) == 0) {ok = 0;}
+      generators = newGens;
+    }
+    else {
+      ok = 0;
+    }
+  }
+  if(nrSols > 0) { result[1] = sol;}
+  result[2] = index;
+  result[3] = generators;
+  return(result);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc NonLinearZeroSetMain(ideal I, intvec wt)
+"USAGE:   ZeroSetMainWork(I, wt, sVars);
+PUROPSE: solves the (nonlinear) univariate polynomials in I
+         of the groundfield (without multiplicities).
+RETURN:  list, all entries are polynomials
+         _[1] = list of solutions
+         _[2] = newA
+         _[3] = minpoly
+         _[4] - index of solved variables
+         _[5] = new representation of I
+ASSUME:  basering = Q[x_1,x_2,...,x_n,a], ideal 'mpoly' must be defined,
+         it might be 0 !
+NOTE:    might change 'mpoly' !!
+"
+{
+  int i, nrSols, ok, n;
+  ideal generators;
+  list result, sols, index, vars, partSol;
+  map F;
+  poly f, newA;
+  string ringSTR;
+
+  def NLZR = basering;
+  export(NLZR);
+
+  n = nvars(basering) - 1;
+
+  generators = I;
+  newA = var(n + 1);
+  result[2] = newA;            // default
+  nrSols = 0;
+  ok = 1;
+  i = 1;
+  while(ok) {
+
+    // test if the i-th generator of I is univariate
+
+    f = generators[i];
+    vars = Variables(f, n);
+    if(size(vars) == 1) {
+      generators[i] = 0;
+      generators = simplify(generators, 2);    // remove 0
+      nrSols++;
+      index[nrSols] = vars[1];      // store index of solved variable
+
+      // create univariate ring
+
+      ringSTR = "ring RIS1 = 0, (" + string(var(vars[1])) + ", " + string(var(n+1)) + "), lp;";
+      execute(ringSTR);
+      ideal mpoly = std(imap(NLZR, mpoly));
+      list roots;
+      poly f = imap(NLZR, f);
+      export(RIS1);
+      export(mpoly);
+      roots = RootsMain(f);
+
+      // get "old" basering with new minpoly
+
+      setring(NLZR);
+      partSol = imap(RIS1, roots);
+      kill(RIS1);
+      if(mpoly[1] != partSol[3]) {      // change minpoly
+        mpoly = std(partSol[3]);
+        F = NLZR, maxideal(1);
+        F[n + 1] = partSol[2];
+        if(size(sols) > 0) {sols = SimplifyZeroset(F(sols)); }
+        newA = reduce(F(newA), mpoly);    // normal form
+        result[2] = newA;
+        generators = SimplifyData(F(generators));  // does not remove 0's
+      }
+      sols = ExtendSolutions(sols, partSol[1]);
+    } // end univariate
+    else {
+      i = i + 1;
+    }
+    if(i > size(generators)) { ok = 0;}
+  }
+  result[1] = sols;
+  result[3] = mpoly;
+  result[4] = index;
+  result[5] = std(generators);
+
+  kill(NLZR);
+  return(result);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc ExtendSolutions(list solutions, list newSolutions)
+"USAGE:   ExtendSolutions(sols, newSols); list sols, newSols;
+PUROPSE: extend the entries of 'sols' by the entries of 'newSols',
+         each entry of 'newSols' is a number.
+RETURN:  list
+ASSUME:  basering = Q[x_1,...,x_n,a], ideal 'mpoly' must be defined,
+         it might be 0 !
+NOTE:    used by 'NonLinearZeroSetMain'
+"
+{
+  int i, j, k, n, nrSols;
+  list newSols, temp;
+
+  nrSols = size(solutions);
+  if(nrSols > 0) {n = size(solutions[1]);}
+  else {
+    n = 0;
+    nrSols = 1;
+  }
+  k = 0;
+  for(i = 1; i <= nrSols; i++) {
+    for(j = 1; j <= size(newSolutions); j++) {
+      k++;
+      if(n == 0) { temp[1] = newSolutions[j];}
+      else {
+        temp = solutions[i];
+        temp[n + 1] = newSolutions[j];
+      }
+      newSols[k] = temp;
+    }
+  }
+  return(newSols);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc MergeSolutions(list sol1, list index1, list sol2, list index2)
+"USAGE:   MergeSolutions(sol1, index1, sol2, index2); all parameters are lists
+RETURN:  list
+PURPOSE: create a list of solutions of size n, each entry of 'sol2' must
+         have size n. 'sol1' is one partial solution (from 'LinearZeroSetMain')
+         'sol2' is a list of partial solutions (from 'NonLinearZeroSetMain')
+ASSUME:  'sol2' is not empty
+NOTE:   used by 'ZeroSetMainWork'
+{
+  int i, j, k, m;
+  ideal sol;
+  list newSols;
+
+  m = 0;
+  for(i = 1; i <= size(sol2); i++) {
+    m++;
+    newSols[m] = SubsMapIdeal(sol1 + sol2[i], index1 + index2, 0);
+  }
+  return(newSols);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc SubsMapIdeal(list sol, list index, int opt)
+"USAGE:   SubsMapIdeal(sol,index,opt); list sol, index; int opt;
+PUROPSE: built an ideal I as follows.
+         if i is contained in 'index' then set I[i] = sol[i]
+         if i is not contained in 'index' then
+         - opt = 0: set I[i] = 0
+         - opt = 1: set I[i] = var(i)
+         if opt = 1 and n = nvars(basering) then set I[n] = var(n).
+RETURN:  ideal
+ASSUME:  size(sol) = size(index) <= nvars(basering)
+"
+{
+  int k = 0;
+  ideal I;
+  for(int i = 1; i <= nvars(basering) - 1; i = i + 1) {    // built subs. map
+    if(ContainedQ(index, i)) {
+      k++;
+      I[index[k]] = sol[k];
+    }
+    else {
+      if(opt) { I[i] = var(i); }
+      else { I[i] = 0; }
+    }
+  }
+  if(opt) {I[nvars(basering)] = var(nvars(basering));}
+  return(I);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc SimplifyZeroset(data)
+"USAGE:   SimplifyZeroset(data); list data
+PUROPSE: reduce the entries of the elements of 'data' w.r.t. the ideal 'mpoly'
+         'data' is a list of ideals/lists.
+RETURN:  list
+ASSUME:  basering = Q[x_1,...,x_n,a], order = lp
+         'data' is a list of ideals
+         ideal 'mpoly' must be defined, it might be 0 !
+"
+{
+  int i;
+  list result;
+
+  for(i = 1; i <= size(data); i++) {
+    result[i] = SimplifyData(data[i]);
+  }
+  return(result);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc Variables(poly f, int n)
+"USAGE:   Variables(f,n); poly f; int n;
+PUROPSE: list of variables among var(1),...,var(n) which occur in f.
+RETURN:  list
+ASSUME:  n <= nvars(basering)
+"
+{
+  int i, nrV;
+  list index;
+
+  nrV = 0;
+  for(i = 1; i <= n; i = i + 1) {
+    if(diff(f, var(i)) != 0) { nrV++; index[nrV] = i; }
+  }
+  return(index);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+proc ContainedQ(data, f, list #)
+"USAGE:    ContainedQ(data, f [, opt]); list data; f is of any type, int opt
+PUROPSE:  test if 'f' is an element of 'data'.
+RETURN:   int
+          0 if 'f' not contained in 'data'
+          1 if 'f' contained in 'data'
+OPTIONS:  opt = 0 : use '==' for comparing f with elements from data
+          opt = 1 : use 'SameQ' for comparing f with elements from data
+"
+{
+  int opt, i, found;
+  if(size(#) > 0) { opt = #[1];}
+  else { opt = 0; }
+  i = 1;
+  found = 0;
+
+  while((!found) && (i <= size(data))) {
+    if(opt == 0) {
+      if(f == data[i]) { found = 1;}
+      else {i = i + 1;}
+    }
+    else {
+      if(SameQ(f, data[i])) { found = 1;}
+      else {i = i + 1;}
+    }
+  }
+  return(found);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+
+proc SameQ(a, b)
+"USAGE:    SameQ(a, b); list/intvec a, b;
+PUROPSE:  test a == b elementwise, i.e., a[i] = b[i].
+RETURN:   int
+          0 if a != b
+          1 if a == b
+"
+{
+  if(typeof(a) == typeof(b)) {
+    if(typeof(a) == "list" || typeof(a) == "intvec") {
+      if(size(a) == size(b)) {
+        int i = 1;
+        int ok = 1;
+        while(ok && (i <= size(a))) {
+          if(a[i] == b[i]) { i = i + 1;}
+          else {ok = 0;}
+        }
+        return(ok);
+      }
+      else { return(0); }
+    }
+    else { return(a == b);}
+  }
+  else { return(0);}
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc SimplifyPoly(poly f)
+"USAGE:   SimplifyPoly(f); poly f
+PUROPSE: reduces the coefficients of f w.r.t. the ideal 'moly' if they contain
+         the algebraic number 'a'.
+RETURN:  poly
+ASSUME:  basering = Q[x_1,...,x_n,a]
+         ideal mpoly must be defined, it might be 0 !
+"
+{
+  matrix coMx;
+  poly f1, vp;
+
+  vp = 1;
+  for(int i = 1; i < nvars(basering); i++) { vp = vp * var(i);}
+
+  coMx = coef(f, vp);
+  f1 = 0;
+  for(i = 1; i <= ncols(coMx); i++) {
+    f1 = f1 + coMx[1, i] * reduce(coMx[2, i], mpoly);
+  }
+  return(f1);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc SimplifyData(data)
+"USAGE:   SimplifyData(data); ideal/list data;
+PUROPSE: reduces the entries of 'data' w.r.t. the ideal 'mpoly' if they contain
+         the algebraic number 'a'
+RETURN:  ideal/list
+ASSUME:  basering = Q[x_1,...,x_n,a]
+         ideal 'mpoly' must be defined, it might be 0 !
+"
+{
+  int n;
+  poly f;
+
+  if(typeof(data) == "ideal") { n = ncols(data); }
+  else { n = size(data);}
+
+  for(int i = 1; i <= n; i++) {
+    f = data[i];
+    data[i] = SimplifyPoly(f);
+  }
+  return(data);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc TransferRing(R)
+"USAGE:   TransferRing(R);
+PUROPSE: creates a new ring containing the same variables as R, but without
+         parameters. If R contains a parameter then this parameter is added
+         as the last variable and 'minpoly' is represented by the ideal 'mpoly'
+         If the basering does not contain a parameter then 'a' is added and
+         'mpoly' = 0.
+RETURN:  ring
+ASSUME:  R = K[x_1,...,x_n] where K = Q or K = Q(a).
+NOTE:    Creates the ring needed for all prodecures with name 'proc-name'Main
+"
+{
+  string ringSTR, parName, minPoly;
+
+  setring(R);
+
+  if(npars(basering) == 0) {
+    parName = "a";
+    minPoly = "0";
+  }
+  else {
+    parName = parstr(basering);
+    minPoly = string(minpoly);
+  }
+  ringSTR = "ring TR = 0, (" + varstr(basering) + "," + parName + "), lp;";
+
+  execute(ringSTR);
+  execute("ideal mpoly = std(" + minPoly + ");");
+  export(mpoly);
+  return(TR);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc NewBaseRing()
+"USAGE:   NewBaseRing();
+PUROPSE: creates a new ring, the last variable is added as a parameter.
+         minpoly is set to mpoly[1].
+RETURN:  ring
+ASSUME:  basering = Q[x_1,...,x_n, a], 'mpoly' must be defined
+"
+{
+  int n = nvars(basering);
+  int MP;
+  string ringSTR, parName, varString;
+
+  def BR = basering;
+  if(mpoly[1] != 0) {
+    parName = "(0, " + string(var(n)) + ")";
+    MP = 1;
+  }
+  else {
+    parName = "0";
+    MP = 0;
+  }
+
+
+  for(int i = 1; i < n - 1; i++) {
+    varString = varString + string(var(i)) + ",";
+  }
+  varString = varString + string(var(n-1));
+
+  ringSTR = "ring TR = " + parName + ", (" + varString + "), lp;";
+  execute(ringSTR);
+  if(MP) { minpoly = number(imap(BR, mpoly)[1]); }
+  return(TR);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
+/*
+                           Examples:
+
+
+// order = 20;
+ring S1 = 0, (s(1..3)), lp;
+ideal I = s(2)*s(3), s(1)^2*s(2)+s(1)^2*s(3)-1, s(1)^2*s(3)^2-s(3), s(2)^4-s(3)^4+s(1)^2, s(1)^4+s(2)^3-s(3)^3, s(3)^5-s(1)^2*s(3);
+ideal mpoly = std(0);
+
+// order = 10
+ring S2 = 0, (s(1..5)), lp;
+ideal I = s(2)+s(3)-s(5), s(4)^2-s(5), s(1)*s(5)+s(3)*s(4)-s(4)*s(5), s(1)*s(4)+s(3)-s(5), s(3)^2-2*s(3)*s(5), s(1)*s(3)-s(1)*s(5)+s(4)*s(5), s(1)^2+s(4)^2-2*s(5), -s(1)+s(5)^3, s(3)*s(5)^2+s(4)-s(5)^3, s(1)*s(5)^2-1;
+ideal mpoly = std(0);
+
+//order = 126
+ring S3 =  0, (s(1..5)), lp;
+ideal I = s(3)*s(4), s(2)*s(4), s(1)*s(3), s(1)*s(2), s(3)^3+s(4)^3-1, s(2)^3+s(4)^3-1, s(1)^3-s(4)^3, s(4)^4-s(4), s(1)*s(4)^3-s(1), s(5)^7-1;
+ideal mpoly = std(0);
+
+// order = 192
+ring S4 = 0, (s(1..4)), lp;
+ideal I = s(2)*s(3)^2*s(4)+s(1)*s(3)*s(4)^2, s(2)^2*s(3)*s(4)+s(1)*s(2)*s(4)^2, s(1)*s(3)^3+s(2)*s(4)^3, s(1)*s(2)*s(3)^2+s(1)^2*s(3)*s(4), s(1)^2*s(3)^2-s(2)^2*s(4)^2, s(1)*s(2)^2*s(3)+s(1)^2*s(2)*s(4), s(1)^3*s(3)+s(2)^3*s(4), s(2)^4-s(3)^4, s(1)*s(2)^3+s(3)*s(4)^3, s(1)^2*s(2)^2-s(3)^2*s(4)^2, s(1)^3*s(2)+s(3)^3*s(4), s(1)^4-s(4)^4, s(3)^5*s(4)-s(3)*s(4)^5, s(3)^8+14*s(3)^4*s(4)^4+s(4)^8-1, 15*s(2)*s(3)*s(4)^7-s(1)*s(4)^8+s(1), 15*s(3)^4*s(4)^5+s(4)^9-s(4), 16*s(3)*s(4)^9-s(3)*s(4), 16*s(2)*s(4)^9-s(2)*s(4), 16*s(1)*s(3)*s(4)^8-s(1)*s(3), 16*s(1)*s(2)*s(4)^8-s(1)*s(2), 16*s(1)*s(4)^10-15*s(2)*s(3)*s(4)-16*s(1)*s(4)^2, 16*s(1)^2*s(4)^9-15*s(1)*s(2)*s(3)-16*s(1)^2*s(4), 16*s(4)^13+15*s(3)^4*s(4)-16*s(4)^5;
+ideal mpoly = std(0);
+
+ring R = (0,a), (x,y,z), lp;
+minpoly = a2 + 1;
+ideal I1 = x2 - 1/2, a*z - 1, y - 2;
+ideal I2 = x3 - 1/2, a*z2 - 3, y - 2*a;
+
+*/
