source: git/Singular/LIB/zeroset.lib @ 7de8e4

// Last change 12.02.2001 (Eric Westenberger)
///////////////////////////////////////////////////////////////////////////////
version="$Id: zeroset.lib,v 1.19 2009-04-06 09:17:01 seelisch Exp $";
category="Symbolic-numerical solving";
info="
LIBRARY:  zeroset.lib      Procedures for roots and factorization
AUTHOR:   Thomas Bayer,    email: tbayer@mathematik.uni-kl.de,@*
          http://wwwmayr.informatik.tu-muenchen.de/personen/bayert/@*
          Current address: Hochschule Ravensburg-Weingarten

OVERVIEW:
 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 algebraic number. Written in the scope of the
 diploma thesis (advisor: Prof. Gert-Martin Greuel) 'Computations of moduli
 spaces of semiquasihomogeneous singularities and an implementation in Singular'.
 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.
 
 NOTE:
 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.

PROCEDURES:
 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

AUXILIARY PROCEDURES:
 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 algebraic 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 ground field
 sqfrNormMain(f)  norm of f (f must be squarefree)
 containedQ(data, f)  f in data ?
 sameQ(a, b)    a == b (list a,b)
";

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)
// return 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); where f is a polynomial
PURPOSE: compute all roots of f in a finite extension of the ground field
         without multiplicities.
RETURN:  ring, a polynomial ring over an extension field of the ground field,
         containing a list 'roots' and polynomials 'newA' and 'f':
  @format
  - 'roots' is the list of roots of the polynomial f (no multiplicities)
  - if the ground field 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')
  @end format
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 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); where f is a polynomial
PURPOSE: compute all roots of f in a finite extension of the ground field
         without multiplicities.
RETURN:  list, all entries are polynomials
  @format
  _[1] = roots of f, each entry is a polynomial
  _[2] = 'newA' - if the ground field is Q(b) and the extension field
         is Q(a), then 'newA' is the representation of b in Q(a)
  _[3] = minpoly of the algebraic extension of the ground field
  @end format
ASSUME:  basering = Q[x,a] ideal mpoly must be defined, it might be 0!
NOTE:    might change the 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 ground field
    result[1] = roots;
    result[2] = var(2);
    result[3] = mpoly[1];
    return(result);
  }

  // process the nonlinear factors, i.e., extend the ground field
  // 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:   zeroSet(I [,opt] ); I=ideal, opt=integer
PURPOSE: compute the zero-set of the zero-dim. ideal I, in a finite extension
         of the ground field.
RETURN:  ring, a polynomial ring over an extension field of the ground field,
         containing a list 'zeroset', a polynomial 'newA', and an
         ideal 'id':
  @format
  - 'zeroset' is the list of the zeros of the ideal I, each zero is an ideal.
  - if the ground field is Q(b) and the extension field is Q(a), then
    'newA' is the representation of b 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')
  @end format
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(RZSB);
  setring 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;
");
  setring RZSB;
  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); where f is a polynomial
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 !
NOTE:    outdated, use / instead
"
{
  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
PURPOSE: 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); where f,g are polynomials;
PURPOSE: compute the quotient q and remainder r s.t. f = g*q + r, deg(r) < deg(g)
RETURN:  list of polynomials
  @format
  _[1] = quotient  q
  _[2] = remainder r
  @end format
ASSUME:  basering = Q[x] or Q(a)[x]
NOTE: outdated, use div/mod instead
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); where f,g are polynomials
PURPOSE: compute the quotient q and remainder r s.th. f = g*q + r, deg(r) < deg(g)
RETURN:  list of polynomials
  @format
  _[1] = quotient  q
  _[2] = remainder r
  @end format
ASSUME:  basering = Q[x,a] and ideal mpoly is defined (it might be 0),
         this represents the ring Q(a)[x] together with its minimal polynomial.
NOTE: outdated, use div/mod instead
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); where f,g are polynomials
PURPOSE: compute the remainder of the division of f by g, i.e. a polynomial r
         s.t. f = g*q + r, deg(r) < deg(g).
RETURN:  poly
ASSUME:  basering = Q[x] or Q(a)[x]
NOTE: outdated, use mod/reduce instead
"
{
  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); where f,g are polynomials
PURPOSE: compute the remainder r s.t. f = g*q + r, deg(r) < deg(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 its minimal polynomial.
NOTE: outdated, use mod/reduce instead
"
{
  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 egcdMain(poly f, poly g)
"USAGE:   egcdMain(f, g); where f,g are polynomials
PURPOSE: 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 its minimal polynomial.
NOTE: outdated, use gcd instead
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
PURPOSE: 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 = gcd(@f, @g);
  setring(@RGCDB);
  poly h = imap(@RGCD, @h);
  kill @RGCD;
  return(h);
}

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

proc sqfrNorm(poly f)
"USAGE:   sqfrNorm(f); where f is a polynomial
PURPOSE: compute the norm of the squarefree polynomial f in Q(a)[x].
RETURN:  list with 3 entries
  @format
  _[1] = squarefree norm of g (poly)
  _[2] = g (= f(x - s*a)) (poly)
  _[3] = s (int)
  @end format
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(SNB);  // new ring with parameter 'a'
                                // replaced by a variable
  setring 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); where f is a polynomial
PURPOSE: compute the norm of the squarefree polynomial f in Q(a)[x].
RETURN:  list with 3 entries
  @format
  _[1] = squarefree norm of g (poly)
  _[2] = g (= f(x - s*a)) (poly)
  _[3] = s (int)
  @end format
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
"
{
  def SNRMB = basering;
  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 = SNRMB, 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
  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 factorMain(poly f)
"USAGE:   factorMain(f); where f is a polynomial
PURPOSE: compute the factorization of the squarefree poly f over Q(a)[t],
         minpoly  = p(a).
RETURN:  list with 2 entries
  @format
  _[1] = factors, first is a constant
  _[2] = multiplicities (not yet implemented)
  @end format
ASSUME:  basering = Q[x,a], representing Q(a)[x]. An ideal mpoly must
         be defined, representing the minimal polynomial (it might be 0!).
NOTE: outdated, use factorize instead
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
PURPOSE: compute the zero-set of the zero-dim. ideal I, in a simple extension
         of the ground field.
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 ground field 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.

  def altring=basering;
  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 ground field 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);
PURPOSE: compute the zero-set of the zero-dim. ideal I, in a finite extension
         of the ground field (without multiplicities).
RETURN:  list, all entries are polynomials
         _[1] = zeros, each entry is an ideal
         _[2] = newA; if the ground field is Q(a') this is the rep. of a' w.r.t. a
         _[3] = minpoly of the algebraic extension of the ground field (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
"
{
  def altring=basering;
  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)
"
{
  def altring=basering;
  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:   NonLinearZeroSetMain(I, wt);
PURPOSE: solves the (nonlinear) univariate polynomials in I
         of the ground field (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;
PURPOSE: 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;
PURPOSE: 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
PURPOSE: 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;
PURPOSE: 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]); data=list; f=any type; opt=integer
PURPOSE:  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 @code{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); a,b=list/intvec
PURPOSE:  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
PURPOSE: 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 !
NOTE: outdated, use reduce instead
"
{
  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;
PURPOSE: 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 !
"
{
  def altring=basering;
  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);
PURPOSE: 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
"
{
  def altring=basering;
  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);
  setring altring;
  return(TR);
}

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

static proc NewBaseRing()
"USAGE:   NewBaseRing();
PURPOSE: 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]); }
  setring BR;
  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;

*/