source: git/Singular/LIB/finvar.lib @ 400884

// $Header: /exports/cvsroot-2/cvsroot/Singular/LIB/finvar.lib,v 1.2 1997-07-03 15:49:54 Singular Exp $
////////////////////////////////////////////////////////////////////////////////
// send bugs and comments to agnes@math.uni-sb.de
LIBRARY:  finvar.lib             LIBRARY TO CALCULATE INVARIANT RINGS & MORE
                                   by Agnes Eileen Heydtmann, send bugs and
                                   comments to agnes@math.uni-sb.de

  rey_mol(G1,G2,...[,int]);      Reynolds operator and Molien series of the
                                 finite matrix group generated by G1,G2,...
  part_mol(M,n[,p]);             n terms of partial expansion of Molien series M
  eval_rey(RO,p);                evaluate poly p under Reynolds operator RO
  rad_con(p,I);                  check radical containment of poly p in ideal I
  inv_basis(deg,G1,G2,...);      basis of space of homogeneous invariants of
                                 degree deg under the finite matrix group
                                 generated by G1,G2,...
  inv_basis_rey(RO,deg[,dim]);   basis of space of homogeneous invariants of
                                 degree deg and optionally dimension dim with
                                 help of Reynolds operator
  inv_ring_s(G1,G2,...[,intvec]); generators of the invariant ring (primary
                                 invariants according to Sturmfels)
  inv_ring_k(G1,G2,...[,intvec]); generators of the invariant ring (combination
                                 of algorithms by Kemper and Sturmfels for
                                 primary invariants)
  algebra_con(p,F);              check whether poly p is contained in invariant
                                 ring generated by entries in F
  module_con(f,P,S);             representing f in the Hironaka decomposition of
                                 the invariant ring into primary invariants P
                                 and secondary ones S
  orbit_var(F,s);                orbit variety of a finite matrix group whose
                                 invariant ring is generated by entries in F
  rel_orbit_var(I,F,s);          relative orbit variety with respect to
                                 invariant ideal I under finite matrix group,
                                 its invariant ring is generated by entries in F
  im_of_var(I,F);                image of variety defined by ideal I under
                                 finite matrix group whose invariant ring is
                                 generated by entries in F

////////////////////////////////////////////////////////////////////////////////
LIB "matrix.lib";
LIB "elim.lib";
LIB "general.lib";
////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////
// sign of integer a, returning 1 or -1 respectively
////////////////////////////////////////////////////////////////////////////////
proc sign(int i)
  USAGE:   sign(<int>);
  RETURN:  the sign of an integer (return type <int>)
  EXAMPLE: example sign; shows an example.
{ if (i>=0)
  { return(1);
  }
  else
  { return(-1);
  }
}
example
{ "  EXAMPLE:";
  echo=2;
           int i=-3;
           int j=3;
           sign(i);
           sign(j);
}

////////////////////////////////////////////////////////////////////////////////
// absolute value of integer a
////////////////////////////////////////////////////////////////////////////////
proc abs (int i)
  USAGE:   abs(<int>);
  RETURN:  the absolute value of an integer (return type <int>)
  EXAMPLE: example abs; shows an example.
{ return(i*sign(i));
}
example
{ "  EXAMPLE:";
  echo=2;
           int i=-3;
           int j=3;
           abs(i);
           abs(j);
}

////////////////////////////////////////////////////////////////////////////////
// Checks whether the last argument, being a matrix, is among the previous
// arguments, also being matrices
////////////////////////////////////////////////////////////////////////////////
proc unique (list #)
{ for (int i=1;i<size(#);i=i+1)
  { if (#[i]==#[size(#)])
    { return(0);
    }
  }
  return(1);
}

////////////////////////////////////////////////////////////////////////////////
// Computes the cyclotomic polynomial recursively, by dividing x^m-1 by the
// cyclotomic polynomial of proper divisors of m
////////////////////////////////////////////////////////////////////////////////
proc cycle (int m)
  USAGE:   cycle(<int>);
  RETURNS: the cyclotomic polynomial (type <poly>) as one in the first ring
           variable
  EXAMPLE: example cycle; shows an example
{ poly v1=var(1);
  if (m==1)
  { return(v1-1);                      // 1-st cyclotomic polynomial
  }
  poly min=v1^m-1;
  matrix s[1][2]=min,v1-1;             // dividing by the 1-st cyclotomic
  s=matrix(syz(ideal(s)));             // polynomial
  min=s[2,1];
  int i=2;
  int n;
  poly c;
  int flag=1;
  while(2*i<=m)                        // there are no proper divisors of m
  { if ((m%i)==0)                      // greater than m/2
    { if (flag==1)
      { n=i;                           // n stores the first proper divisor of
      }                                // m>1
      flag=0;
      c=cycle(i);                      // recursive computation
      s=min,c;
      s=matrix(syz(ideal(s)));         // dividing
      min=s[2,1];
    }
    if (n*i==m)                        // the earliest possible point to break
    { break;
    }
    i=i+1;
  }
  min=min/leadcoef(min);               // making sure that leading coefficient
  return(min);                         // is 1
}
example
{ echo=2;
          ring R=0,(x,y,z),dp;
          print(cycle(25));
}

////////////////////////////////////////////////////////////////////////////////
// Returns i such that root^i==n, i.e. it heavily relies on the right input.
////////////////////////////////////////////////////////////////////////////////
proc power(number n, number root)
{ int i=0;
   while((n/root^i)<>1)
   { i=i+1;
   }
   return(i);
}

////////////////////////////////////////////////////////////////////////////////
// Generates the Molien series when the characteristic of the base field is p>0
// and p does not divide the group order. Input is the entire group and a name
// for a new ring.
////////////////////////////////////////////////////////////////////////////////
proc p_molien(list #)
{ def br=basering;                     // keeping track of the base ring since
  int n=nvars(br);                     // we have to go into an extension of the
  int g=size(#)-2;                     // basefield -
  matrix G(1..g)=#[1..g];              // rewriting the group elements
  string newring=#[g+1];
  int flag=#[g+2];
  if (g<>1)
  { ring Q=0,x,dp;                     // we want to extend our ring as well as
                                       // the ring of rational numbers Q to
                                       // contain g-th primitive roots of unity
                                       // in order to factor characteristic
                                       // polynomials of group elements into
                                       // linear factors and lift eigenvalues to
                                       // characteristic 0 -
    poly minq=cycle(g);                // minq now contains the size-of-group-th
                                       // cyclotomic polynomial of Q, it is
                                       // irreducible there
    ring `newring`=(0,e),x,dp;
    map f=Q,ideal(e);
    minpoly=number(f(minq));           // e is now a g-th primitive root of
                                       // unity -
    kill Q, f;                         // no longer needed -
    poly p=1;                          // used to build the denominator of the
                                       // new term in the Molien series
    matrix s[1][2];                    // used for canceling -
    matrix M[1][2]=0,1;                // will contain Molien series -
    ring v1br=char(br),x,dp;           // we calculate the g-th cyclotomic
    poly minp=cycle(g);                // polynomial of the base field and pick
    minp=factorize(minp)[1][2];        // an irreducible factor of it -
    if (deg(minp)==1)                  // in this case the base field contains
    { ring bre=char(br),x,dp;          // g-th roots of unity already
      map f1=v1br,ideal(0);
      number e=-number((f1(minp)));    // e is a g-th primitive root of unity
    }
    else
    { ring bre=(char(br),e),x,dp;
      map f1=v1br,ideal(e);
      minpoly=number(f1(minp));        // e is a g-th primitive root of unity
    }
    map f2=br,ideal(0);                // we need f2 to map our group elements
                                       // to this new extension field bre
    matrix I=unitmat(n);
    poly p;                            // used for the characteristic polynomial
                                       // to factor -
    list L;                            // will contain the linear factors of the
    ideal F;                           // characteristic polynomial of the group
    intvec C;                          // elements and their powers
    int i, j, k;
    for (i=1;i<=g;i=i+1)
    { p=det(x*I-f2(G(i)));             // characteristic polynomial of G(i)
      L=factorize(p);
      F=L[1];
      C=L[2];
      for (j=2;j<=ncols(F);j=j+1)
      { F[j]=-1*(F[j]-x);              // F[j] is now an eigenvalue of G(i),
                                       // it is a power of a primitive g-th root
                                       // of unity -
        k=power(number(F[j]),e);       // F[j]==e^k
        setring `newring`;
        p=p*(1-x*(e^k))^C[j];          // building the denominator of the new
        setring bre;                   // term
      }
      setring `newring`;
      M[1,1]=M[1,1]*p+M[1,2];          // expanding M[1,1]/M[1,2] + 1/p
      M[1,2]=M[1,2]*p;
      p=1;
      s=matrix(syz(ideal(M)));         // canceling common terms of denominator
      M[1,1]=-s[2,1];                  // and enumerator
      M[1,2]=s[1,1];
      setring bre;
      if (flag)
      { "  Term "+string(i)+" has been computed.";
      }
    }
    if (flag)
    { "";
    }
    setring `newring`;
    map slead=`newring`,ideal(0);
    s=slead(M);                        // forcing the constant term of numerator
    M[1,1]=1/s[1,1]*M[1,1];            // and denominator to be 1
    M[1,2]=1/s[1,2]*M[1,2];
    kill slead;
    kill s;
    kill p;
  }
  else                                 // if the group only contains an identity
  { ring `newring`=0,x,dp;             // element, it is very easy to calculate
    matrix M[1][2]=1,(1-x)^n;          // the Molien series
  }
  // keepring `newring`;
  export `newring`;                    // TTO we keep the ring where we computed
                                       // the Molien series
  export M;                            // TTO so that we can keep the Molien
                                       // series
  setring br;
}

////////////////////////////////////////////////////////////////////////////////
// This procedure calculates all members of a finite matrix group in terms of
// the given generators. In one run trough the main loop, all left products of
// the generators with the new elements from the last run through the loop (or
// the generators themselves in the first run) will be formed. After that the
// newly generated elements will be added to the group and the loop starts over
// again unless no elements were added.
// Additionally, every time a new matrix is added to the group, its
// corresponding ring mapping in the Reynolds operator and if the
// characteristic is 0, its corresponding summand of the Molien series is
// calculated.
// When the characteristic of the basefield is p>0 such that it does not
// divide the group order, the Molien series is calculated at the end of the
// procedure.
// No matter when the Molien series is calculated, the procedure expands after
// every step to obtain a rational function.
// The first result of the procedure is the Reynolds operator, presented in
// form of a matrix; each row can be transformed into an ideal and from
// there can be used as a ring homomorphism via the command 'map'.
// If the characteristic is 0, the second result is a matrix, containing
// enumerator and denominator (with no common divisor) of the final
// rational function representing the Molien series.
// When the characteristic of the basefield is p>0 such that it does not
// divide the group order, the Molien series is returned in a ring of
// characteristic 0. It names was specified in the list of parameters.
////////////////////////////////////////////////////////////////////////////////
proc rey_mol (list #)
  USAGE:   rey_mol(<generators of a finite matrix group>[,<string>,<int>]);
           if the characteristic of the coefficient field is prime, <string>
           has to contain the name for a new polynomials ring with coefficient
           field of characteristic 0 that stores the Molien series - if <int> is
           not not equal to 0, some information will be printed during the run
  RETURNS: if the characteristic is 0: Reynolds operator (type <matrix>), Molien
           series (type <matrix> with two components, first being the numerator,
           second the denominator)
           if the characteristic is p>0 not dividing the group order: Reynolds
           operator (type <matrix>) - the Molien series will directly be stored
           under the name M (type <matrix>) in the ring `<string>`
           if the characteristic is p>0 dividing the group order: Reynolds
           operator (type <matrix>)
  EXAMPLE: example rey_mol; shows an example
{ def br=basering;                     // the Molien series depends on the
  int ch=char(br);                     // characteristic of the coefficient
  int flag;                            // field -
  if (ch<>0)                           // making sure the input is 'correct'...
  { if (typeof(#[size(#)])=="string")
    { flag=size(#)-1;
      string newring=#[size(#)];
      int v=0;                         // no information is default
    }
    else
    { if (typeof(#[size(#)-1])=="string")
      { flag=size(#)-2;
        string newring=#[size(#)-1];
        if (typeof(#[size(#)])<>"int")
        { "  ERROR:   if the second last parameter is <string>, the last must be";
          "           of type <int>";
          return();
        }
        int v=#[size(#)];
      }
      else
      { "  ERROR:   in characteristic p a <string> must be given for the name";
        "           of a new ring";
        return();
      }
    }
    if (newring=="")
    { "  ERROR:   <string> may not be empty";
      return();
    }
  }
  else
  { if (typeof(#[size(#)])=="int")
    { flag=size(#)-1;
      int v=#[size(#)];
    }
    else
    { flag=size(#);
      int v=0;                         // no information is default
    }
  }
  if (typeof(#[1])<>"matrix")
  { "  ERROR:   the parameters must be a list of matrices and optionally";
    "           a <string> and an <int>";
    return();
  }
  int n=nrows(#[1]);
  if (n<>nvars(br))
  { "  ERROR:   the number of variables of the basering needs to be the same";
    "           as the dimension of the matrices";
    return();
  }
  if (n<>ncols(#[1]))
  { "  ERROR:   matrices need to be square and of the same dimensions";
    return();
  }
  matrix vars=matrix(maxideal(1));     // creating an nx1-matrix containing the
  vars=transpose(vars);                // variables of the ring -
  matrix A(1)=#[1]*vars;               // calculating the first ring mapping -
                                       // A(1) will contain the Reynolds
                                       // operator -
  if (ch==0)                           // when ch==0 we can calculate the Molien
  { matrix I=diag(1,n);                // series in any case -
    poly v1=vars[1,1];                 // the Molien series will be in terms of
                                       // the first variable of the current
                                       // ring -
    matrix A(2)[1][2];                 // A(2) will contain the Molien series -
    A(2)[1,1]=1;                       // A(2)[1,1] will be the numerator
    A(2)[1,2]=det(I-v1*(#[1]));        // A(2)[1,2] will be the denominator -
    matrix s;                          // will help us canceling in the
                                       // fraction
  }
  matrix G(1)=#[1];                    // G(k) are elements of the group -
  poly p;                              // will contain the denominator of the
                                       // new term of the Molien series
  int i=1;
  for (int j=2;j<=flag;j=j+1)          // this loop adds the arguments to the
  {                                    // group, leaving out doubles and
                                       // checking whether the arguments are
                                       // compatible with the task of the
                                       // procedure
    if (not(typeof(#[j])=="matrix"))
    { "  ERROR:   the parameters must be a list of matrices and optionally";
      "           a <string> and an <int>";
      return();
    }
    if ((n!=nrows(#[j])) or (n!=ncols(#[j])))
    { "  ERROR:   matrices need to be square and of the same dimensions";
       return();
    }
    if (unique(G(1..i),#[j]))
    { i=i+1;
      matrix G(i)=#[j];
      A(1)=concat(A(1),#[j]*vars);     // adding ring homomorphisms to A(1)
      if (ch==0)
      { p=det(I-v1*#[j]);              // denominator of new term -
        A(2)[1,1]=A(2)[1,1]*p+A(2)[1,2]; // expanding A(2)[1,1]/A(2)[1,2] + 1/p
        A(2)[1,2]=A(2)[1,2]*p;
        s=matrix(syz(ideal(A(2))));    // canceling common factors
        A(2)[1,1]=-s[2,1];
        A(2)[1,2]=s[1,1];
      }
    }
  }
  int g=i;                             // G(1)..G(i) are generators without
                                       // doubles - g generally is the number
                                       // of elements in the group so far -
  j=i;                                 // j is the number of new elements that
                                       // we use as factors -
  int k, m, l;
  if (v)
  { if (ch==0)
    { "";
      "  Generating the entire matrix group, Reynolds operator and Molien series...";
      "";
    }
    else
    { "";
      "  Generating the entire matrix group and Reynolds operator...";
      "  If the characteristic of the basefield divides the order of the";
      "  group the result will be useless.";
      "";
    }
  }
  while (1)
  { l=0;   // l is the number of products we get in one going
    for (m=g-j+1;m<=g;m=m+1)
    { for (k=1;k<=i;k=k+1)
      { l=l+1;
        matrix P(l)=G(k)*G(m);         // possible new element
      }
    }
    j=0;
    for (k=1;k<=l;k=k+1)
    { if (unique(G(1..g),P(k)))
      { j=j+1;                         // a new factor for next run
        g=g+1;
        matrix G(g)=P(k);              // a new group element -
        A(1)=concat(A(1),P(k)*vars);   // adding new mapping to A(1)
        if (ch==0)
        { p=det(I-v1*P(k));            // denominator of new term -
          A(2)[1,1]=A(2)[1,1]*p+A(2)[1,2];
          A(2)[1,2]=A(2)[1,2]*p;       // expanding A(2)[1,1]/A(2)[1,2] + 1/p -
          s=matrix(syz(ideal(A(2))));  // canceling common factors
          A(2)[1,1]=-s[2,1];
          A(2)[1,2]=s[1,1];
        }
        if (v)
        { "  Group element "+string(g)+" has been found.";
        }
      }
      kill P(k);
    }
    if (j==0)                          // when we didn't add any new elements
    { break; }                         // in one run through the while loop
  }                                    // we are done -
  if (v)
  { if (g<=i)
    { "  There are only "+string(g)+" group elements.";
    }
    "";
  }
  A(1)=transpose(A(1));                // when we evaluate the Reynolds operator
                                       // later on, we actually want 1xn
                                       // matrices
  if (ch<>0 && minpoly==0)
  { if ((g%ch)<>0)
    { if (v)
      { "  Generating Molien series...";
        "";
      }
      p_molien(G(1..g),newring,v);     // the procedure that defines a ring of
                                       // characteristic 0 and calculates the
                                       // Molien series in it
      if (v)
      { "  Now we are done calculating Molien series and Reynolds operator.";
        "";
      }
      return(A(1));
    }
  }
  if (ch<>0 && minpoly<>0)
  { if ((g%ch)<>0)
    { if (voice==2)
      { "  WARNING: It is impossible for this program to calculate the Molien series";
        "           for finite groups over extension fields of prime characteristic.";
      }
      else
      { if (v)
        { "  Since it is impossible for this program to calculate the Molien series for";
          "  invariant rings over extension fields of prime characteristic, we have to";
          "  continue without it. The Reynolds operator is available, however.";
          "";
        }
      }
      return(A(1));
    }
  }
  if (ch<>0)
  { if ((g%ch)==0)
    { if (voice==2)
      { A(1)=0;
        "  WARNING: The characteristic of the coefficient field divides the group";
        "           order. Proceed without the Molien series or Reynolds operator!";
      }
      else
      { if (v)
        { "  The characteristic of the base field divides the group order.";
          "  We have to continue without Molien series and without Reynolds";
          "  operator..";
          "";
        }
      }
      return(A(1));
    }
  }
  if (ch==0)
  { map slead=br,ideal(0);
    s=slead(A(2));
    A(2)[1,1]=1/s[1,1]*A(2)[1,1];      // numerator and denominator have to have
    A(2)[1,2]=1/s[1,2]*A(2)[1,2];      // a constant term of 1
    if (v)
    { "  Now we are done calculating Molien series and Reynolds operator.";
      "";
    }
    return(A(1..2));
  }
}
example
{ "  EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
  "           note the case of prime characteristic";
  echo=2;
           ring R=0,(x,y,z),dp;
           matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
           matrix RM(1..2);
           RM(1..2)=rey_mol(A);
           print(RM(1..2));
           ring S=3,(x,y,z),dp;
           string newring="Qadjoint";
           matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
           matrix REY=rey_mol(A,newring);
           print(REY);
           setring Qadjoint;
           M;
           setring S;
           kill Qadjoint;
}

////////////////////////////////////////////////////////////////////////////////
// This procedure implements the following calculation:
// (1+a[1]x+a[2]x2+...+a[n]xn)/(1+b[1]x+b[2]x2+...+b[m]xm)=(1+(a[1]-b[1])x+...
// (1+b[1]x+b[2]x2+...+b[m]xm)
// ---------------------------
//    (a[1]-b[1])x+(a[2]-b[2])x2+...
//    (a[1]-b[1])x+b[1](a[1]-b[1])x2+...
////////////////////////////////////////////////////////////////////////////////
proc part_mol (matrix M, int n, list #)
  USAGE:   part_mol(M,n[,p]); M <matrix> (return value of 'rey_mol'), n <int>,
           indicating  number of terms in the expansion, p <poly> optionally, it
           ought to be the second return value of a previous run of 'part_mol'
           and avoids recalculating known terms
  RETURNS: n terms of partial expansion of the Molien series (type <poly>)
           (first n if there is no third argument given, otherwise the next n
           terms depending on a previous calculation) and an intermediate result
           (type <poly>) of the calculation to be used as third argument in a
           next run
  EXAMPLE: example part_mol; shows an example
{ poly A(2);                           // A(2) will contain the return value of
                                       // the intermediate result
  if (char(basering)<>0)
  { "  ERROR:   you have to change to a basering of characteristic 0, one in";
    "           which the Molien series is defined";
  }
  if (ncols(M)==2 && nrows(M)==1 && n>0 && size(#)<2)
  { def br=basering;                   // keeping track of the old ring
    map slead=br,ideal(0);
    matrix s=slead(M);
    if (s[1,1]<>1 || s[1,2]<>1)
    { "  ERROR:   the constant terms of enumerator and denominator are not 1";
      return();
    }

    if (size(#)==0)
    { A(2)=M[1,1];                     // if a third argument is not given, the
                                       // intermediate result from the last run
                                       // corresponds to the numerator - we need
    }                                  // its smallest term
    else
    { if (typeof(#[1])=="poly")
      { A(2)=#[1];                     // if a third term is given we 'start'
      }                                // with its smallest term
      else
      { "  ERROR:   <poly> as third argument expected";
        return();
      }
    }
    poly A(1)=M[1,2];                  // denominator of Molien series
                                       // (for now) -
    string mp=string(minpoly);
    execute "ring R=("+charstr(br)+"),("+varstr(br)+"),ds;";
    execute "minpoly=number("+mp+");";
    poly A(1)=0;                       // A(1) will contain the sum of n terms -
    poly min;                          // min will be our smallest term -
    poly A(2)=fetch(br,A(2));          // fetching A(2) and M[1,2] into R
    poly den=fetch(br,A(1));
    for (int i=1; i<=n; i=i+1)         // getting n terms and adding them up
    { min=lead(A(2));
      A(1)=A(1)+min;
      A(2)=A(2)-min*den;
    }
    setring br;                        // moving A(1) and A(2) back in the
    A(1)=fetch(R,A(1));                // actual ring for output
    A(2)=fetch(R,A(2));
    return(A(1..2));
  }
  else
  { "  ERROR:   the first argument has to be a 1x2-matrix, i.e. the matrix";
    "           returned by the procedure 'rey_mol', the second one";
    "           should be > 0 and there should be no more than 3 arguments;"
    return();
  }
}
example
{ "  EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
  echo=2;
           ring R=0,(x,y,z),dp;
           matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
           matrix B(1..2);
           B(1..2)=rey_mol(A);
           poly C(1..2);
           C(1..2)=part_mol(B(2),5);
           C(1);
           C(1..2)=part_mol(B(2),5,C(2));
           C(1);
}

////////////////////////////////////////////////////////////////////////////////
// RO will simply be cut into pieces and each row will act as a ring
// mapping of which the Reynolds operator is made up.
////////////////////////////////////////////////////////////////////////////////
proc eval_rey (matrix RO, poly f)
  USAGE:   eval_rey(RO,f); RO <matrix> (result of rey_mol),
           f <poly>
  RETURNS: image of f under the Reynolds operator (type <poly>)
  NOTE:    the characteristic of the coefficient field of the polynomial ring
           should not divide the order of the finite matrix group
  EXAMPLE: example eval_rey; shows an example
{ def br=basering;
  int n=nvars(br);
  if (ncols(RO)==n)
  { int m;                             // we need m to 'cut' the ring
                                       // homomorphisms 'out' of RO and to
    m=nrows(RO);                       // divide by the group order in the end
    poly p=0;
    map pRO;
    matrix RH[1][n];
    for (int i=1;i<=m;i=i+1)
    { RH=RO[i,1..n];
      pRO=br,ideal(RH);                // f is now the i-th ring homomorphism
      p=pRO(f)+p;
    }
    p=(1/poly(m))*p;
    return(p);
  }
  else
  { "  ERROR:   the number of columns in the matrix, being the first argument";
    "           should be the same as number of variables in the basering, in";
    "           fact it should be the matrix returned by 'rey_mol'";
    return();
  }
}
example
{ "  EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
  echo=2;
           ring R=0,(x,y,z),dp;
           matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
           matrix B(1..2);
           B(1..2)=rey_mol(A);
           poly p=x2;
           eval_rey(B(1),p);
}

////////////////////////////////////////////////////////////////////////////////
// Input: <ideal>=<f1,f2,...,fm> and <polynomial> g
// Question: Does g lie in the radical of <ideal>?
// Solution: Compute a standard basis G for <f1,f2,...,fm,gz-1> where z is a new
//           variable. Then g is contained in the radical of <ideal> <=> 1 is
//           generator in G.
////////////////////////////////////////////////////////////////////////////////
proc rad_con (poly g,ideal I)
  USAGE:   rad_con(<poly>,<ideal>);
  RETURNS: 1 (TRUE) (type <int>) if <poly> is contained in the radical of
           <ideal>, 0 (FALSE) (type <int>) otherwise
  EXAMPLE: example rad_con; shows an example
{ def br=basering;
  int n=nvars(br);
  int dB=degBound;
  degBound=0;
  string mp=string(minpoly);
  execute "ring R=("+charstr(br)+"),(x(1..n),z),dp;";
  execute "minpoly=number("+mp+");";
  ideal irrel=x(1..n);
  map f=br,irrel;
  poly p=f(g);
  ideal J=f(I)+ideal(p*z-1);
  J=std(J);
  degBound=dB;
  if (J[1]==1)
  { return(1);
  }
  else
  { return(0);
  }
}
example
{ "  EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7.";
  echo=2;
           ring R=0,(x,y,z),dp;
           ideal I=x2+y2,z2;
           poly f=x4+y4;
           rad_con(f,I);
           ideal J=x2+y2,z2,x4+y4;
           poly g=z;
           rad_con(g,I);
}

////////////////////////////////////////////////////////////////////////////////
// This procedure generates a basis of invariant polynomials in degree g. The
// way this works, is that we look how the generators act on a general
// polynomial of degree g - it turns out that one simply has to solve a system
// of linear equations.
////////////////////////////////////////////////////////////////////////////////
proc inv_basis (int g, list #)
  USAGE:   inv_basis(<int>,<generators of a finite matrix group>); <int>
           indicates in which degree (>0) we are looking for invariants
  RETURNS: the basis (type <ideal>) of the space of invariants of degree <int_1>
  EXAMPLE: example inv_basis; shows an example
{ if (g<=0)
  { "  ERROR:   the first argument should be > 0";
    return();
  }
  def br=basering;
  ideal mon=sort(maxideal(g))[1];      // needed for constructing a general
  int m=ncols(mon);                    // homogeneous polynomial of degree d
  int a=size(#);
  int i;
  int n=nvars(br);
  for (i=1;i<=a;i=i+1)                 // checking that input is ok
  { if (typeof(#[i])=="matrix")
    { if (nrows(#[i])==n && ncols(#[i])==n)
      { matrix G(i)=#[i];
      }
      else
      { "  ERROR:   the number of variables of the base ring needs to be the same";
        "           as the dimension of the square matrices";
        return();
      }
    }
    else
    { "  ERROR:   the last arguments should be a list of matrices";
      return();
    }
  }
  ideal vars_old=maxideal(1);
  execute "ring T=("+charstr(br)+"),("+varstr(br)+",p(1..m)),lp;";
  ideal vars=imap(br,vars_old);
  // p(1..m) are general coefficients of
  // the general polynomial
  map f;
  ideal mon=imap(br,mon);
  poly P=0;
  for (i=m;i>=1;i=i-1)
  { P=P+p(i)*mon[i];                   // P is the general polynomial
  }
  ideal I;                             // will help substituting variables in P
                                       // by linear combinations of variables -
  poly Pnew, temp;                     // Pnew is P with substitutions -
  matrix S[m*a][m];                    // will contain system of linear
                                       // equations
  int j, k;
  for (i=1;i<=a;i=i+1)                 // building system of linear equations
  { I=ideal(matrix(vars)*transpose(imap(br,G(i))));
    I=I,p(1..m);
    f=T,I;
    Pnew=f(P);
    for (j=1;j<=m;j=j+1)
    { temp=P/mon[j]-Pnew/mon[j];
      for (k=1;k<=m;k=k+1)
      { S[m*(i-1)+j,k]=temp/p(k);
      }
    }
  }
  setring br;
  map f=T,ideal(0);
  matrix S=f(S);
  matrix s=matrix(syz(S));             // s contains a basis of the space of
                                       // solutions -
  ideal I=ideal(matrix(mon)*s);        // I contains a basis of homogeneous
  if (I[1]<>0)                         // invariants of degree d
  { for (i=1;i<=ncols(I);i=i+1)
    { I[i]=I[i]/leadcoef(I[i]);        // setting leading coefficients to 1
    }
  }
  return(I);
}
example
{ "  EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
  echo=2;
             ring R=0,(x,y,z),dp;
             matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
             print(inv_basis(2,A));
}

////////////////////////////////////////////////////////////////////////////////
// This procedure generates invariant polynomials of degree g via the Reynolds
// operator and checks by calculating syzygies whether they are linearly
// independent. If they are the first column of syzygies does not contain any
// constant polynomials. If a third argument of type <int> is given, the
// program stopes once that many linearly independent polynomials have been
// found.
////////////////////////////////////////////////////////////////////////////////
proc inv_basis_rey (matrix RO, int g, list #)
  USAGE:   inv_basis_rey(<matrix>,<int_1>[,<int_2>]); <matrix> should be the
           Reynolds operator which is the first return value of rey_mol, <int_1>           indicates the degree of the invariants and <int_2> optionally the
           dimension of the space which is known from 'part_mol'
  RETURNS: the basis <ideal> of the space of invariants of degree <int_1>
  EXAMPLE: example inv_basis_rey; shows an example
{ if (g<=0)
  { "  ERROR:   the second argument should be > 0";
     return();
  }
  if (size(#)>0)
  { if (typeof(#[1])<>"int")
    { "  ERROR: the third argument should be of type <int>";
      return();
    }
    if (#[1]<0)
    { "  ERROR: the third argument should be and <int> >= 0";
      return();
    }
  }
  int i, k;
  ideal mon=sort(maxideal(g))[1];
  int j=ncols(mon);
  matrix S[ncols(mon)][1];             // will contain linear systems of
  int counter=0;                       // equations -
  degBound=g;                          // syzygies of higher degree need not be
                                       // computed -
  poly imRO;                           // image of Reynolds operator -
  ideal B;                             // will contain the basis
  for (i=j;i>0;i=i-1)
  { imRO=eval_rey(RO,mon[i]);
    if (imRO<>0)                       // the first candidate<>0 will definitely
    { if (counter==0)                  // be in the basis
      { B[1]=imRO;
        B[1]=B[1]/leadcoef(B[1]);
        counter=counter+1;
      }
      else                             // other candidates have to be checked
      { B=B,imRO;                      // for linear independence
        S=syz(B);
        k=1;
        while(k<>counter+2)
        { if (S[k,1]==0)               // checking whether there are constant
          { k=k+1;                     // entries <>0 in S
          }
          else
          { break;
          }
        }
        if (k==counter+2)              // this means that the loop was not
        { counter=counter+1;           // broken, we can keep B[counter]
          B[counter]=B[counter]/leadcoef(B[counter]);
        }
        else                           // we have to get rid of B[counter]
        { B[counter+1]=0;
          B=compress(B);
        }
      }
    }
    if (size(#)>0)
    { if (counter==#[1])               // we have found enough elements (if the
      { break;                         // user entered the right dim...
      }
    }
  }
  degBound=0;
  return(B);
}
example
{ "  EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
  echo=2;
             ring R=0,(x,y,z),dp;
             matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
             matrix B(1..2);
             B(1..2)=rey_mol(A);
             print(inv_basis_rey(B(1),6,8));
}

////////////////////////////////////////////////////////////////////////////////
// Procedure returning the succeeding vector after vec. It is used to list
// all the vectors of Z^n with first nonzero entry 1. They are listed by
// increasing sum of the absolute value of their entries.
////////////////////////////////////////////////////////////////////////////////
proc nextvec(intmat vec)
{ int n=ncols(vec);                    // p: >0, n: <0, p0: >=0, n0: <=0
  for (int i=1;i<=n;i=i+1)             // finding out which is the first
  { if (vec[1,i]<>0)                   // component <>0
    { break;
    }
  }
  intmat new[1][n];
  if (i>n)                             // 0,...,0 --> 1,0....,0
  { new[1,1]=1;
    return(new);
  }
  if (i==n)                            // 0,...,1 --> 1,1,0,...,0
  { new[1,1..2]=1,1;
    return(new);
  }
  if (i==n-1)
  { if (vec[1,n]==0)                   // 0,...,0,1,0 --> 0,...,0,1
    { new[1,n]=1;
      return(new);
    }
    if (vec[1,n]>0)                    // 0,..,0,1,p --> 0,...,0,1,-p
    { new[1,1..n]=vec[1,1..n-1],-vec[1,n];
      return(new);
    }
    new[1,1..2]=1,1-vec[1,n];          // 0,..,0,1,n --> 1,1-n,0,..,0
    return(new);
  }
  if (i>1)
  { intmat temp[1][n-i+1]=vec[1,i..n]; // 0,...,0,1,*,...,* --> 1,*,...,*
    temp=nextvec(temp);
    new[1,i..n]=temp[1,1..n-i+1];
    return(new);
  }                                    // case left: 1,*,...,*
  for (i=2;i<=n;i=i+1)
  { if (vec[1,i]>0)                    // make first positive negative and
    { vec[1,i]=-vec[1,i];              // return
      return(vec);
    }
    else
    { vec[1,i]=-vec[1,i];              // make all negatives before positives
    }                                  // positive
  }
  for (i=2;i<=n-1;i=i+1)               // case: 1,p,...,p after 1,n,...,n
  { if (vec[1,i]>0)
    { vec[1,2]=vec[1,i]-1;             // shuffleing things around...
      if (i>2)                         // same sum of absolute values of entries
      { vec[1,i]=0;
      }
      vec[1,i+1]=vec[1,i+1]+1;
      return(vec);
    }
  }                                    // case left: 1,0,...,0 --> 1,1,0,...,0
  new[1,2..3]=1,vec[1,n];              // and: 1,0,...,0,1 --> 0,1,1,0,...,0
  return(new);
}

////////////////////////////////////////////////////////////////////////////////
// Input is a list of nxm-matrices with n<m and rank n. Procedure checks whether
// the space generated by the rows of the last matrix lies in any of the spaces
// generated by other matrices' rows. Returns a boolean answer.
////////////////////////////////////////////////////////////////////////////////
proc space_con (list #)
{ matrix H;
  int n=nrows(#[1]);
  for (int i=1;i<size(#);i=i+1)
  { H=transpose(#[i]);
    H=concat(H,transpose(#[size(#)])); // concatenating works column-wise -
    H=bareiss(transpose(H));           // bareiss works row-wise -
    if (ncols(compress(transpose(H)))==n)  // means that the last rows of the
    { return(1);                       // matrix were in the span of the rows of
    }                                  // #[i]
  }
  return(0);
}

////////////////////////////////////////////////////////////////////////////////
// Maps integers to elements of the base field. It is only called if the base
// field is of prime characteristic. If the base field has q elements (depending
// on minpoly) 1..q is mapped to those q elements.
////////////////////////////////////////////////////////////////////////////////
proc intnumap (int i)
{ int p=char(basering);
  if (minpoly==0)                      // if no minpoly is given, we have p
  { i=i%p;                             // elements in the field
    return(number(i));
  }
  int d=pardeg(minpoly);
  if (i<0)
  { int bool=1;
    i=(-1)*i;
  }
  i=i%p^d;                             // base field has p^d elements
  number a=par(1);                     // a is the root of the minpoly, we have
  number out=0;                        // to construct a linear combination of
  int j=1;                             // a^k
  int k;
  while (1)
  { if (i<p^j)                         // finding an upper bound on i
    { for (k=0;k<j-1;k=k+1)
      { out=out+((i/p^k)%p)*a^k;       // finding how often p^k is contained in
      }                                // i
      out=out+(i/p^(j-1))*a^(j-1);
      if (defined(bool)=voice)
      { return((-1)*out);
      }
      return(out);
    }
    j=j+1;
  }
}

////////////////////////////////////////////////////////////////////////////////
// Attempting to construct n=[number of variables in the base ring] linear
// combinations of the m>n entries in Q such that the ideal generated by these
// combinations is of dimension 0. It is then a Noetherian normalization of the
// invariant ring. In characteristic 0 the existence of such a linear
// combination is ensured.
////////////////////////////////////////////////////////////////////////////////
proc noethernorm(ideal Q)
{ def br=basering;
  int lcm=deg(Q[1]);                   // will contain lowest common multiple of
  int ch=char(br);                     // degrees of polynomials in Q
  int n=nvars(br);
  int i, j;
  intvec degvec;
  int m=ncols(Q);
  degvec[1]=lcm;
  for (i=2;i<=m;i=i+1)
  { degvec[i]=deg(Q[i]);
    lcm=lcm*degvec[i]/gcd(lcm,degvec[i]); // lcm is now the least common
  }                                    // multiple of the first i elements of Q
  ideal A(1)=Q;
  for (i=1;i<=m;i=i+1)
  { A(1)[i]=(A(1)[i])^(lcm/degvec[i]); // now all elements in A(1) are of the
  }                                    // same degree, they are the elements of
                                       // Q raised to a power -
  matrix T[n][1];                      // will contain the n linear combinations
  matrix I[n][n]=unitmat(n);
  matrix H(1)[n][m];
  H(1)[1..n,1..n]=I[1..n,1..n];        // H(1) will be the first matrix, we try
  kill I;
  if ((n%2)==0)                        // H(1) ought to be of the form:
  { j=int(n)/int(2);                   // 1,0,...,0,0,1,0,...,0
  }                                    // 0,0,...,0,1,0,0,...,0
  else                                 //      .           .
  { j=int(n-1)/int(2);                 //      .           .
  }                                    //      .           .
  for (i=1;i<=j;i=i+1)                 // 1,0,...,0,0,0,0,...,0
  { H(1)=permcol(H(1),i,n-i+1);
  }
  H(1)[1,1]=1;
  int c=1;
  intmat vec[1][n*m];
  vec[1,1..n*m]=int(H(1)[1..n,1..m]);  // we rewrite H(1) as a vector
  while (1)
  { T=H(c)*transpose(matrix(A(1)));
    Q=ideal(T);
    attrib(Q,"isSB",1);
    if (dim(Q)>0)
    { if (dim(std(Q))==0)              // we found n linear combinations
      { A(1)=T;
        break;
      }
    }
    else                               // we found n linear combinations
    { A(1)=T;
      break;
    }
    matrix H(c+1)[n][m];               // we have to find a new matrix
    while(1)                           // generating n linear combinations
    { vec=nextvec(vec);
      if (ch==0)
      { H(c+1)[1..n,1..m]=vec[1,1..n*m];
      }
      else
      { for (i=1;i<=n;i=i+1)
        { for (j=1;j<=m;j=j+1)
          { H(c+1)[i,j]=intnumap(vec[1,(i-1)*m+j]); // mapping integers to the
          }                            // field
        }
      }
      if (minor(H(c+1),n)[1]<>0 && not(space_con(H(1..c+1)))) // if the ideal
      { c=c+1;                         // generated by the minors is not the 0
        break;                         // ideal and if the span of rows of
      }                                // H(c+1) is not in the span of rows
                                       // previously tried, then we found a new
                                       // interesting matrix
    }
  }
  if (ch==0)
  { poly p(1)=(1-var(1)^lcm)^n;        // since all elements are of degree
                                       // lcm, the denominator of the Hilbert
                                       // series of the ring generated by the
                                       // primary invariants equals p(1)
    return(A(1),p(1));
  }
  else
  { if (defined(Qa))                   // here is where we store Molien series
    { setring Qa;
      poly p(1)=(1-x^lcm)^n;           // since all elements are of degree
                                       // lcm, the denominator of the Hilbert
                                       // series of the ring generated by the
                                       // primary invariants equals p(1)
      setring br;
      return(A(1),p(1));
    }
    else
    { return(A(1));
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
// Computing the entire matrix group from generators and returning its
// cardinality.
////////////////////////////////////////////////////////////////////////////////
proc group (list #)
{ matrix G(1)=#[1];                    // first group element
  int i=1;
  for (int j=2;j<=size(#);j=j+1)       // throwing out doubles among the
  { if (unique(G(1..i),#[j]))          // generators
    { i=i+1;
      matrix G(i)=#[j];
    }
  }
  int g=i;                             // g: elements in the group so far, i:
  j=i;                                 // generators, j: new ones used as
  int m, k, l;                         // as factors, l: counting possible new
                                       // new elements
  while (1)
  { l=0;
    for (m=g-j+1;m<=g;m=m+1)
    { for (k=1;k<=i;k=k+1)
      { l=l+1;
        matrix P(l)=G(k)*G(m);         // possible new element
      }
    }
    j=0;
    for (k=1;k<=l;k=k+1)               // checking whether the P(k) are new
    { if (unique(G(1..g),P(k)))
      { j=j+1;
        g=g+1;
        matrix G(g)=P(k);              // adding new elements -
      }
      kill P(k);
    }
    if (j==0)                          // when we didn't add any new elements
    { break;                           // in one run through the while loop, we
    }                                  // are done
  }
  return(g);
}

////////////////////////////////////////////////////////////////////////////////
// If the characteristic of the base field is zero or prime not dividing the
// order of the group G, one can compute secondary invariants (free module
// generators) even without the Molien series. In other words, when the user
// enters a flag that tells the procedures inv_ring_s or inv_ring_k not to compute
// the Molien series, it the number of group elements will be computed (with
// group). If the characteristic is 0 or prime not dividing the order of the
// group, there are deg(P[1])*...*deg(P[n])/|G| free module generators where P
// contains the primary invariants. sec_minus_mol computes these secondary
// invariants by going through the various spaces of homogeneous invariants
// successively, starting with degree 1.
// list # is made of of various things. The last component is an integer, saying
// how many of the immediately preceding elements are bases of various vector
// spaces of homogeneous invariants. Before these bases, is a boolean variable.
// if it is 0, the preceding elements are group generators (we will use
// inv_basis), if it is 0, the Reynolds operator is passed on (and we can use
// inv_basis_rey).
// sP is the standard basis of the ideal generated by primary invariants in P. g
// is the cardinality of the group. v is the verbose-level.
////////////////////////////////////////////////////////////////////////////////
proc sec_minus_mol (ideal P, ideal sP, int g, int v, list#)
{ def br=basering;
  int n=nvars(br);
  int d=1;
  int r=size(#)-#[size(#)]-1;
  for (int i=r+1;i<size(#);i=i+1)
  { ideal B(i-r)=#[i];                 // rewriting the bases
  }
  for (i=1;i<=n;i=i+1)
  { d=d*deg(P[i]);                     // building the product of the degrees of
  }                                    // primary invariants -
  int bound=d/g;                       // number of secondary invariants
  if (v)
  { "  The invariant ring is Cohen-Macaulay.";
    "  We need to find "+string(d)+"/"+string(g)+"="+string(bound)+" secondary invariants.";
    "";
  }
  if (bound==1)                        // in this case, it is quick
  { if (v)
    { "  In degree 0 we have: 1";
      "";
      "  We're done!";
      "";
    }
    return(matrix(1));
  }
  qring Qring=sP;                      // secondary invariants are linearly
                                       // independent modulo the ideal generated
                                       // by primary invariants -
  ideal Smod;                          // stores secondary invariants modulo sP
                                       // that are homogeneous of the same
                                       // degree -
  ideal Bmod;                          // basis of homogeneous invariants modulo
                                       // sP -
  ideal sSmod;                         // standard basis of Smod modulo sP
  setring br;
  matrix S[1][bound]=1;                // stores all secondary invariants
  if (v)
  { "  In degree 0 we have: 1";
    "";
  }
  int counter=1;                       // counts secondary invariants -
  d=1;                                 // the degree of homogeneous invariants
  int degcounter=0;                    // counts secondary invariants of degree
                                       // d -
  int bool=1;                          // decides when std needs to be computed
  while (counter<>bound)
  { if (v)
    { "  Searching in degree "+string(d)+"...";
    }
    if (d>#[size(#)])                  // we need to compute basis of degree d
    {                                  // in this case -
      if (#[r])                        // in this case, we have the Reynolds
      { ideal B(d)=inv_basis_rey(#[r-1],d); // operator
      }
      else
      { ideal B(d)=inv_basis(d,#[1..r-1]);
      }
    }
    if (B(d)[1]<>0)                    // we only need to look for secondary
    { setring Qring;                   // invariants in this degre if B is not
      Smod=0;                          // the zero ideal
      Bmod=fetch(br,B(d));
      for (i=1;i<=ncols(Bmod);i=i+1)
      { if (degcounter<>0)
        { if (reduce(Bmod[i],std(ideal(0)))<>0) // in this case B[i] might be
          {                            // qualify as secondary invariant -
            if (bool)                  // compute a standard basis only if a new
            { sSmod=std(Smod);         // secondary invariant has been found in
            }                          // the last run -
            if (reduce(Bmod[i],sSmod)<>0) // if Bmod[i] is not contained in Smod
            { counter=counter+1;       // B[i] qualifies as secondary invariant
              degcounter=degcounter+1;
              Smod[degcounter]=Bmod[i];
              setring br;
              S[1,counter]=B(d)[i];
              if (v)
              { "           "+string(B(d)[i]);
              }
              bool=1;                  // we have to compute std next time
              setring Qring;
              if (counter==bound)      // in this case, we're done
              { break;
              }
            }
            else                       // next time, we don't need to compute
            { bool=0;                  // standard basis
            }
          }
        }
        else
        { if (reduce(Bmod[i],std(ideal(0)))<>0)
          { Smod[1]=Bmod[i];           // here we just add Bmod[i] without
            setring br;                // having to check linear independence
            counter=counter+1;
            degcounter=degcounter+1;
            S[1,counter]=B(d)[i];
            if (v)
            { "  We find: "+string(B(d)[i]);
            }
            setring Qring;
            bool=1;                    // next time, we have to compute std
            if (counter==bound)
            { break;
            }
          }
        }
      }
    }
    if (v and degcounter<>0)
    { "";
    }
    degcounter=0;
    setring br;
    d=d+1;                             // go to next degree
  }
  if (v)
  { "  We're done!";
  }
  return(S);
}

////////////////////////////////////////////////////////////////////////////////
// inv_ring_s calculates the primary and secondary invariants of the invariant
// ring with respect to a finite matrix group G. The primary invariants generate
// an invariant subring, lets say R, and the secondary invariants generate the
// invariant ring as an R-module. If the characteristic of the base field is
// zero or prime not dividing the group order, the secondary invariants are free
// generators and we have the Hironaka decomposition of the invariant ring.
// Otherwise the secondary invariants are possible not free generators.
// The procedure is based on the algorithms given by Sturmfels in "Algorithms
// in Invariant Theory" except for the one computing secondary invariants when
// the characteristic divides the group order which is based on Kemper's
// "Calculating Invariants Rings of Finite Groups over Arbitrary Fields".
////////////////////////////////////////////////////////////////////////////////
proc inv_ring_s (list #)
  USAGE:   inv_ring_s(<generators of a finite matrix group>[,<intvec>]);
           <intvec> has to contain 2 flags; if the first one equals 0, the
           program attempts to compute the Molien series and Reynolds operator,
           if it equals 1, the program is told that the characteristic of the
           base field divides the group order, if it is anything else the Molien
           series and Reynolds operator will not be computed; if the second flag
           does not equal 0, information about the various stages of the program
           will be printed while running
  RETURNS: generators of the invariant ring with respect to the matrix group
           generated by the matrices in the input; there are two return values
           of type <matrix>, the first containing primary invariants and the
           second secondary invariants, i.e. module generators over a Noetherian
           normalization
  EXAMPLE: example inv_ring_s; shows an example
{ def br=basering;
  int ch=char(br);                     // the algorithms depend very much on the
                                       // characteristic of the ground field
  int dB=degBound;
  degBound=0;
  int n=nvars(br);                     // n is the number of variables, as well
                                       // as the size of the matrices, as well
                                       // as the number of primary invariants,
                                       // we have to find
  if (typeof(#[size(#)])=="intvec")
  { if (size(#[size(#)])<>2)
    { "  ERROR:   <intvec> must have exactly two entires";
      return();
    }
    intvec flagvec=#[size(#)];
    if (flagvec[1]==0)
    { if (ch==0)
      { matrix R(1..2);                // one will contain Reynolds operator and
                                       // the other enumerator and denominator
                                       // of Molien series
        R(1..2)=rey_mol(#[1..size(#)-1],flagvec[2]);
      }
      else
      { string newring="Qa";
        matrix R(1)=rey_mol(#[1..size(#)-1],newring,flagvec[2]); // will contain
                                       // Reynolds operator, if Molien series
      }                                // can be computed, it will be stored in
                                       // the new ring Qa
    }
    else
    { for (int i=1;i<=size(#)-1;i=i+1) // checking whether the input is ok
      { if (not(typeof(#[i])=="matrix"))
        { "  ERROR:   the parameters must be a list of matrices and optionally";
          "           an <intvec>";
          return();
        }
        if (n<>ncols(#[i]) || n<>nrows(#[i]))
        { "  ERROR:   matrices need to be square and of the same dimensions as";
          "           the number of variables of the basering";
          return();
        }
      }
      kill i;
    }
  }
  else
  { if (typeof(#[size(#)])<>"matrix")
    { "  ERROR:   the parameters must be a list of matrices and optionally";
      "           an <intvec>";
      return();
    }
    if (ch==0)
    { matrix R(1..2);                  // will contain Reynolds operator and
                                       // enumerator and denominator of Molien
                                       // series
      R(1..2)=rey_mol(#[1..size(#)]);
    }
    else
    { string newring="Qa";             // we might need as a new ring of
                                       // characteristic 0 where we store the
                                       // Molien series -
      matrix R(1)=rey_mol(#[1..size(#)],newring); // will contain
                                       // Reynolds operator
    }
    intvec flagvec=0,0;                // default flags, no info
  }
  ideal Q=0;                           // will contain the candidates for
                                       // primary invariants -
  if (flagvec[1]==0 && flagvec[2])
  { "  We can start looking for primary invariants...";
    "";
  }
  else
  { if (flagvec[1] && flagvec[2])
    { "";
      "  We start by looking for primary invariants...";
      "";
    }
  }
  if ((ch==0 || defined(Qa)) && flagvec[1]==0) // i.e. we can use Molien series
  { if (ch==0)
    { poly p(1..2);                    // p(1) will be used for single terms of
                                       // the partial expansion, p(2) to store
      p(1..2)=part_mol(R(2),1);        // the intermediate result -
      poly v1=var(1);                  // we need v1 to split off coefficients
                                       // in the partial expansion of M (which
                                       // is in terms of the first variable) -
      poly d;                          // for splitting off the coefficient in
                                       // in one term of the partial expansion,
                                       // i.e. it stores the dimension of the
                                       // current homogeneous subspace
    }
    else
    { setring Qa;                      // Qa is where the Molien series is
                                       // stored -
      poly p(1..2);                    // p(1) will be used for single terms of
                                       // the partial expansion, p(2) to store
      p(1..2)=part_mol(M,1);           // the intermediate result -
      poly d;                          // stores the dimension of the current
                                       // homogeneous subspace
      setring br;
    }
    int g, di, counter, i, j, bool;    // g: current degree, di: d as integer,
                                       // counter: counts candidates in degree
                                       // g, i,j: going through monomials of
                                       // degree g, bool: indicating when the
                                       // ideal generated by the candidates
                                       // has dimension 0 -
    ideal mon;                         // will contain monomials of degree g -
    poly imRO;                         // the image of the Reynolds operator -
    while(1)                           // repeat until we reach dimension 0
    { if (ch==0)
      { p(1..2)=part_mol(R(2),1,p(2)); // 1 term of the partial expansion -
        g=deg(p(1));                   // current degree -
        d=coef(p(1),v1)[2,1];          // dimension of invariant space of degree
                                       // g -
        di=int(d);                     // just a type cast
      }
      else
      { setring Qa;
        p(1..2)=part_mol(M,1,p(2));    // 1 term of the partial expansion -
        g=deg(p(1));                   // current degree -
        d=coef(p(1),x)[2,1];           // dimension of invariant space of degree
                                       // g -
        di=int(d);                     // just a type cast
        setring br;
      }
      if (flagvec[2])
      { "  Searching for candidates in degree "+string(g)+":";
        "  There is/are "+string(di)+" linearly independent invariant(s) to choose from...";
      }
      mon=sort(maxideal(g))[1];        // all monomials of degree g -
      j=ncols(mon);
      counter=0;                       // we have 0 candidates of degree g so
                                       // far
      for (i=j;i>=1;i=i-1)
      { imRO=eval_rey(R(1),mon[i]);
        if (imRO<>0)
        { if (Q[1]==0)                 // if imRO is the first non-zero
          { counter=1;                 // invariant we find, the rad_con
            Q[1]=imRO/leadcoef(imRO);  // question is trivial and we just
            if (flagvec[2])            // include imRO
            { "  Found: "+string(Q[1]);
            }
            if (counter==di)           // if counter is up to di==d, we can
            { break;                   // leave the for-loop
            }
          }
          else
          { if (not(rad_con(imRO,Q)))  // if imRO is not contained in the
            { counter=counter+1;       // radical of Q, we add it to the
              Q=Q,imRO/leadcoef(imRO); // generators of Q
              if (flagvec[2])
              { "  Found: "+string(Q[ncols(Q)]);
              }
            }
              if (ncols(Q)>=n)         // when we have n or more candidates, we
              { attrib(Q,"isSB",1);    // test if dim(Q)==0, Singular might
              if (dim(Q)==0)           // recognize this property even if Q is
              { bool=1;                // no standard basis, but that is not
                break;                 // guaranteed -
              }                        // if dim(Q) is 0, we can construct a
              else                     // set of primary invariants from the
              { if (dim(std(Q))==0)    // generators of Q and we can leave both
                { bool=1;              // the for- and the while-loop
                  break;
                }
              }
              }
            if (counter==di)           // if counter is up to di, we can leave
            { break;                   // the for-loop
            }
          }
        }
      }
      if (n==1 or bool)                // if n=1, we're done when we've found
      { break;                         // the first
      }
    }
    if (flagvec[2])
    { "";
    }
    int m=ncols(Q);                    // m tells us if we found too many
                                       // candidates -
    ideal P=Q;                         // will eventually contain the primary
                                       // invariants -
    if (n<m)                           // the number of primary invariants
    { counter=m;                       // should be the same as the number of
      for (i=m-1;i>=1;i=i-1)           // variables in the basering; we are
      {                                // checking whether we can leave out some
        Q[i]=0;                        // candidates and still have full
                                       // radical -
        attrib(Q,"isSB",1);
        if (dim(Q)==0)                 // we're going backwards through the
        { P[i]=0;                      // candidates to throw out large degrees
          counter=counter-1;
        }
        else
        { if (dim(std(Q))==0)
          { P[i]=0;
            counter=counter-1;
          }
        }
        if (counter==n)
        { break;
        }
        Q=P;
      }
      P=compress(P);
      m=counter;
      if (m==n)
      { Q=std(P);                      // standard basis for computing secondary
                                       // invariants
      }
    }
    else                               // we need the standard basis of P to be
    { Q=std(P);                        // able to do calculations modulo primary
    }                                  // invariants
    intvec degvec;
    if (n<m)
    { if (flagvec[2] and ch==0)
      { "  We have too many candidates for primary invariants and have to find a";
        "  Noetherian normalization.";
        "";
      }
      if (ch<>0)
      { "  We have too many candidates for primary invariants and have to attempt";
        "  to construct a Noetherian normalization as linear combinations of powers";
        "  of the candidates. Careful! Termination is not guaranteed!";
        "";
      }
      P,p(1)=noethernorm(P);           // p(1) is the denominator of the Hilbert
                                       // series with respect to primary
                                       // invariants from P -
      Q=std(P);                        // we need to do calculations modulo
                                       // primary invariants -
      for (j=1;j<=n;j=j+1)             // we set the leading coefficients of the
      { P[j]=P[j]/leadcoef(P[j]);      // primary invariants to 1
      }
    }
    else                               // this is when m==n without Noetherian
    {                                  // normalization
      if (ch==0)
      { p(1)=1;
        for (j=1;j<=n;j=j+1)           // calculating the denominator of the
        { p(1)=p(1)*(1-v1^deg(P[j]));  // Hilbert series of the ring generated
        }                              // by the primary invariants
      }
      else
      { for (j=1;j<=n;j=j+1)           // degrees have to be taken in a ring
        { degvec[j]=deg(P[j]);         // of characteristic 0
          }
        setring Qa;
        p(1)=1;
          for (j=1;j<=n;j=j+1)         // calculating the denominator of the
          { p(1)=p(1)*(1-x^degvec[j]); // Hilbert series of the ring
          }                            // generated by the primary invariants
          setring br;
      }
    }
    if (flagvec[2])
    { "  These are the primary invariants: ";
      for (i=1;i<=n;i=i+1)
      { "   "+string(P[i]);
      }
      "";
    }
    if (ch==0)
    { matrix s[1][2]=R(2)[1,1]*p(1),R(2)[1,2]; // used for canceling
      s=matrix(syz(ideal(s)));
      p(1)=s[2,1];                     // the polynomial telling us where to
                                       // search for secondary invariants
      map slead=br,ideal(0);
      p(1)=1/slead(p(1))*p(1);         // smallest term of p(1) needs to be 1
      if (flagvec[2])
      { "  Polynomial telling us where to look for secondary invariants:";
        "   "+string(p(1));
        "";
      }
      matrix dimmat=coeffs(p(1),v1);   // dimmat will contain the number of
                                       // secondary invariants, we need to find
                                       // of a certain degree -
      m=nrows(dimmat);                 // m-1 is the highest degree
      degvec=0;
      for (j=1;j<=m;j=j+1)
      { if (dimmat[j,1]<>0)
        { degvec[j]=int(dimmat[j,1]);  // degvec contains the degrees of
        }                              // secondary invariants
      }
    }
    else
    { setring Qa;
      matrix s[1][2]=M[1,1]*p(1),M[1,2]; // used for canceling
      s=matrix(syz(ideal(s)));
      p(1)=s[2,1];                     // the polynomial telling us where to
                                       // search for secondary invariants
      map slead=Qa,ideal(0);
      p(1)=1/slead(p(1))*p(1);         // smallest term of p(1) needs to be 1
      if (flagvec[2])
      { "  Polynomial telling us where to look for secondary invariants:";
        "   "+string(p(1));
        "";
      }
      matrix dimmat=coeffs(p(1),x);    // dimmat will contain the number
                                       // of secondary invariants, we need
                                       // to find of a certain degree -
      m=nrows(dimmat);                 // m-1 is the highest
      degvec=0;
      for (j=1;j<=m;j=j+1)
      { if (dimmat[j,1]<>0)
        { degvec[j]=int(dimmat[j,1]);  // degvec[j] contains the number of
        }                              // secondary invariants of degree j-1
      }
      setring br;
      kill Qa;                         // all the information needed from Qa is
    }                                  // stored in dimmat -
    qring Qring=Q;                     // we need to do calculations modulo the
                                       // ideal generated by the primary
                                       // invariants, its standard basis is
                                       // stored in Q -
    poly imROmod;                      // imRO reduced -
    ideal Smod, sSmod;                 // secondary invariants of one degree
                                       // reduced and their standard basis
    setring br;
    kill Q;                            // Q might be big and isn't needed
                                       // anymore -
    ideal S=1;                         // secondary invariants, 1 definitely is
                                       // one
    if (flagvec[2])
    { "  Proceeding to look for secondary invariants...";
      "";
      "  In degree 0 we have: 1";
      "";
    }
    bool=0;                            // indicates when std-calculation is
                                       // necessary -
    for (i=2;i<=m;i=i+1)               // walking through degvec -
    { if (degvec[i]<>0)                // when it is == 0 we need to find 0
      {                                // elements of the degree i-1
        if (flagvec[2])
        { "  Searching in degree "+string(i-1)+", we need to find "+string(degvec[i])+" invariant(s)...";
        }
        mon=sort(maxideal(i-1))[1];    // all monomials of degree i-1 -
        counter=0;                     // we'll count up to degvec[i] -
        j=ncols(mon);                  // we'll go through mon from the end
        setring Qring;
        Smod=0;
        setring br;
        while (degvec[i]<>counter)     // we need to find degvec[i] linearly
        {                              // independent (in Qring) invariants -
          imRO=eval_rey(R(1),mon[j]);  // generating invariants
          setring Qring;
          imROmod=fetch(br,imRO);      // reducing the invariants
          if (reduce(imROmod,std(ideal(0)))<>poly(0) and counter<>0)
          {                            // if the first one is true and the
                                       // second false, imRO is the first
                                       // secondary invariant of that degree
                                       // that we want to add and we need not
                                       // check linear independence
            if (bool)
            { sSmod=std(Smod);
            }
            if (reduce(imROmod,sSmod)<>0)
            { Smod=Smod,imROmod;
              setring br;              // we make its leading coefficient to be
              imRO=imRO/leadcoef(imRO); // 1
              S=S,imRO;
              counter=counter+1;
              if (flagvec[2])
              { "           "+string(imRO);
              }
              bool=1;                  // next time we need to recalculate std
            }
            else
            { bool=0;                  // std-calculation is unnecessary
              setring br;
            }
          }
          else
          { if (reduce(imROmod,std(ideal(0)))<>poly(0) and counter==0)
            { Smod[1]=imROmod;         // here we just add imRO(mod) without
              setring br;              // having to check linear independence
              imRO=imRO/leadcoef(imRO);
              S=S,imRO;
              counter=counter+1;
              bool=1;                  // next time we need to calculate std
              if (flagvec[2])
              { "  We find: "+string(imRO);
              }
            }
            else
            { setring br;
            }
          }
          j=j-1;                       // going to next monomial
        }
        if (flagvec[2])
        { "";
        }
      }
    }
    degBound=dB;
    if (flagvec[2])
    { "  We're done!";
      "";
    }
    matrix FI(1)=matrix(P);
    matrix FI(2)=matrix(S);
    return(FI(1..2));
  }
                                       // this case is entered when either the
                                       // characteristic<>0 divides the group
                                       // order or when the Molien series could
                                       // not or has not been computed -
  if (flagvec[1]==0)                   // indicates that it has been attempted
  {                                    // to compute the Reynolds operator
                                       // etc. -
    int g=nrows(R(1));                 // order of the group -
    int flag=((g%ch)==0);              // flag is 1 if the characteristic
                                       // divides the order, it is 0 if it does
                                       // not -
    if (typeof(#[size(#)])=="intvec")  // getting a hold of the generators of
    { int gennum=size(#)-1;            // the group
    }
    else
    { int gennum=size(#);
    }
  }
  else
  { int flag=2;                        // flag is 2 if we don't know yet whether
    int gennum=size(#)-1;              // the group order is divisible by the
  }                                    // characteristic -
  int d=1;                             // d is set to the current degree, since
                                       // we know nothing about the finite
                                       // matrix group (via Molien series) we
                                       // have to start with degree 1 -
  int counter;                         // counts candidates for primary
                                       // invariants -
  int i, di, bool;
  while (1)
  { if (flagvec[2])
    { "  Searching for candidates in degree "+string(d)+":";
    }
    if (flag)                          // in this case we can not make use of
    {                                  // the Reynolds operator -
      ideal B(d)=inv_basis(d,#[1..gennum]); // we create a basis of the vector
                                       // space of all invariant polynomials of
    }                                  // degree d
    else
    {                                  // here the characteristic<>0 does not
      ideal B(d)=inv_basis_rey(R(1),d); // divide the group order, i.e. the
    }                                  // Reynolds operator can be used to
                                       // calculate a basis of the vector space
                                       // of all invariant polynomials of degree
                                       // d -
    di=ncols(B(d));                    // dimension of the homogeneous space -
    if (B(d)[1]<>0)                    // otherwise the space is empty
    { if (flagvec[2])
      { "  There is/are "+string(di)+" linearly independent invariant(s) to choose from...";
      }
      if (counter==0)                  // we have no candidates for primary
      {                                // invariants yet, i.e. don't have to
        Q[1]=B(d)[1];                  // check for radical containment
        if (flagvec[2])
        { "  Found: "+string(Q[1]);
        }
        i=2;                           // proceed with the second element of
        counter=1;                     // B(d)
        if (n==1)
        { break;
        }
      }
      else
      { i=1;                           // proceed with the first element of B(d)
      }
      while (i<=di)                    // goes through all polynomials in B(d) -
      { if (not(rad_con(B(d)[i],Q)))   // B(d)[i] is not in the radical of Q
        { counter=counter+1;
          Q=Q,B(d)[i];                 // including candidate
          if (flagvec[2])
          { "  Found: "+string(Q[counter]);
          }
          if (counter>=n)
          { attrib(Q,"isSB",1);
            if (dim(Q)==0)
            { bool=1;                  // when the dimension is 0, we're done
              break;                   // but this can only be when counter>=n
            }
            else
            { if (dim(std(Q))==0)
              { bool=1;                // bool indicates whether we are done
                break;
              }
            }
          }
        }
        i=i+1;                         // going to next element in basis
      }
      if (bool)
      { break;
      }
    }
    else
    { if (flagvec[2])
      { "  The space is 0-dimensional.";
      }
    }
    d=d+1;                             // up to the next degree
  }
  if (flagvec[2])
  { "";
  }
  int j;
  ideal P=Q;                           // P will contain primary invariants -
  if (n<counter)                       // we have too many candidates -
  { for (i=counter-1;i>=1;i=i-1)       // we take a look whether we can leave
    { Q[i]=0;                          // out some candidates, but have full
                                       // radical
      attrib(Q,"isSB",1);
      if (dim(Q)==0)                   // we're going backwards through the
      { P[i]=0;                        // candidates to throw out large degrees
        counter=counter-1;
      }
      else
      { if (dim(std(Q))==0)
        { P[i]=0;
          counter=counter-1;
        }
      }
      if (counter==n)
      { break;
      }
      Q=P;
    }
    P=compress(P);
    if (counter==n)
    { Q=std(P);
    }
  }
  else
  { Q=std(P);                          // we need to do calculations modulo
  }                                    // primary invariants
  if (n<counter)
  { if (flagvec[2] and ch==0)
    { "  We have too many candidates for primary invariants and have to find a";
      "  Noetherian normalization.";
      "";
    }
    if (ch<>0)
    { "  We have too many candidates for primary invariants and have to attempt";
      "  to construct a Noetherian normalization as linear combinations of powers";
      "  of the candidates. Careful! Termination is not guaranteed!";
      "";
    }
    P=noethernorm(P);
    for (j=1;j<=n;j=j+1)               // we set the lead coefficients of the
    { P[j]=P[j]/leadcoef(P[j]);        // primary invariants to be 1
    }
    Q=std(P);
  }
  if (flagvec[2])
  { "  These are the primary invariants: ";
    for (i=1;i<=n;i=i+1)
    { "   "+string(P[i]);
    }
    "";
    "  Proceeding to look for secondary invariants...";
  }
  // we can now proceed to calculate secondary invariants, we face the fact
  // that we can make no use of a Molien series - however, if the
  // characteristic does not divide the group order, we can make use of the
  // fact that the secondary invariants are free module generators and that we
  // need deg(P[1])*...*deg(P[n])/(cardinality of the group) of them
  if (flagvec[1]<>0 and flagvec[1]<>1)
  { int g=group(#[1..size(#)-1]);      // computing group order
    if (ch==0)
    { matrix FI(2)=sec_minus_mol(P,Q,g,flagvec[2],#[1..size(#)-1],0,B(1..d),d);
      matrix FI(1)=matrix(P);
      return(FI(1..2));
    }
    if (g%ch<>0)
    { matrix FI(2)=sec_minus_mol(P,Q,g,flagvec[2],#[1..size(#)-1],0,B(1..d),d);
      matrix FI(1)=matrix(P);
      return(FI(1..2));
    }
  }
  else
  { if (flag==0)                       // this is the case where we have a
    {                                  // nonzero minpoly, but the
                                       // characteristic does not divide the
                                       // group order
      matrix FI(2)=sec_minus_mol(P,Q,g,flagvec[2],R(1),1,B(1..d),d);
      matrix FI(1)=matrix(P);
      return(FI(1..2));
    }
  }
  if (flagvec[2])
  { "  Since the characteristic of the base field divides the group order, we do not";
    "  know whether the invariant ring is Cohen-Macaulay. We have to use Kemper's";
    "  algorithm and compute secondary invariants with respect to the trivial";
    "  subgroup of the given group.";
    "";

  }
  // we are using Kemper's algorithm with the trivial subgroup
  ring QQ=0,x,ds;                      // we lock at our primary invariants as
  ideal M=(1-x)^n;                     // such of the subgroup that only
                                       // contains the identity, this means that
                                       // ch does not divide the order anymore,
                                       // this means that we can make use of the
                                       // Molien series again - 1/M[1] is the
                                       // Molien series of that group, we now
                                       // calculate the secondary invariants of
                                       // this subgroup in the usual fashion
                                       // where the primary invariants are the
                                       // ones from the bigger group
  setring br;
  intvec degvec;                       // for the degrees of the primary
                                       // invariants -
  for (i=1;i<=n;i=i+1)                 // finding the degrees of these
  { degvec[i]=deg(P[i]);
  }
  setring QQ;                          // calculating the polynomial indicating
  M[2]=1;                              // where to search for secondary
  for (i=1;i<=n;i=i+1)                 // invariants (of the trivial subgroup)
  { M[2]=M[2]*(1-x^degvec[i]);
  }
  M=matrix(syz(M))[1,1];
  M[1]=M[1]/leadcoef(M[1]);
  if (flagvec[2])
  { "  Polynomial telling us where to look for these secondary invariants:";
    "   "+string(M[1]);
    "";
  }
  matrix dimmat=coeffs(M[1],x);        // storing the number of secondary
                                       // invariants we need in a certain
                                       // degree -
  int m=nrows(dimmat);                 // m-1 is the highest degree where we
                                       // need to search
  degvec=0;
  for (i=1;i<=m;i=i+1)                 // degvec will contain all the
  { if (dimmat[i,1]<>0)                // information about where to find
    { degvec[i]=int(dimmat[i,1]);      // secondary invariants, it is filled
    }                                  // with integers and therefore visible in
  }                                    // all rings
  kill QQ;
  setring br;
  ideal S=1;                           // 1 is a secondary invariant always -
  if (flagvec[2])
  { "  In degree 0 we have: 1";
    "";
  }
  ideal B;                             // basis of homogeneous invariants of a
                                       // certain degree with respect to the
                                       // trivial subgroup - i.e. all monomials
                                       // of that degree -
  qring Qring=Q;                       // need to do computations modulo primary
                                       // invariants -
  ideal Smod, sSmod, Bmod;             // Smod: secondary invariants of one
                                       // degree modulo Q, sSmod: standard basis
                                       // of the latter, Bmod: B modulo Q
  setring br;
  kill Q;                              // might be large
  int k;
  bool=0;                              // indicates when we need to do standard
                                       // basis computation -
  for (i=2;i<=m;i=i+1)                 // going through all entries of degvec
  { if (degvec[i]<>0)
    { B=sort(maxideal(i-1))[1];        // basis of the space of invariants (with
                                       // respect to the matrix subgroup
                                       // containing only the identity) of
                                       // degree i-1 -
      if (flagvec[2])
      { "  Searching in degree "+string(i-1)+", we need to find "+string(degvec[i])+" invariant(s)...";
      }
      counter=0;                       // we have 0 secondary invariants of
                                       // degree i-1
      setring Qring;
      Bmod=fetch(br,B);                // basis modulo primary invariants
      Smod=0;
      j=ncols(Bmod);                   // going backwards through Bmod
      while (degvec[i]<>counter)
      { if (reduce(Bmod[j],std(ideal(0)))<>0 && counter<>0)
        { if (bool)
          { sSmod=std(Smod);
          }
          if (reduce(Bmod[j],sSmod)<>0) // Bmod[j] qualifies as secondary
          { Smod=Smod,Bmod[j];         // invariant
            setring br;
            S=S,B[j];
            counter=counter+1;
            if (flagvec[2])
            { "           "+string(B[j]);
            }
            setring Qring;
            bool=1;                    // need to calculate std of Smod next
          }                            // time
          else
          { bool=0;                    // no std calculation necessary
          }
        }
        else
        { if (reduce(Bmod[j],std(ideal(0)))<>0 && counter==0)
          { Smod[1]=Bmod[j];           // in this case, we may just add B[j]
            setring br;
            S=S,B[j];
            if (flagvec[2])
            { "  We find: "+string(B[j]);
            }
            counter=counter+1;
            bool=1;                    // need to calculate std of Smod next
            setring Qring;             // time
          }
        }
        j=j-1;                         // next basis element
      }
      setring br;
    }
  }
  // now we have those secondary invariants
  k=ncols(S);                          // k: number of the secondary invariants,
                                       // we just calculated
  if (flagvec[2])
  { "";
    "  We calculate secondary invariants from the ones found for the trivial";
    "  subgroup.";
    "";
  }
  map f;                               // used to let generators act on
                                       // secondary invariants with respect to
                                       // the trivial group -
  matrix M(1)[gennum][k];              // M(1) will contain a module
  for (i=1;i<=gennum;i=i+1)
  { B=ideal(matrix(maxideal(1))*transpose(#[i])); // image of the various
                                       // variables under the i-th generator -
    f=br,B;                            // the corresponding mapping -
    B=f(S)-S;                          // these relations should be 0 -
    M(1)[i,1..k]=B[1..k];              // we will look for the syzygies of M(1)
  }
  module M(2)=res(M(1),2)[2];
  m=ncols(M(2));                       // number of generators of the module
                                       // M(2) -
  // the following steps calculates the intersection of the module M(2) with the
  // algebra A^k where A denote the subalgebra of the usual polynomial ring,
  // generated by the primary invariants
  string mp=string(minpoly);           // generating a ring where we can do
                                       // elimination
  execute "ring R=("+charstr(br)+"),(x(1..n),y(1..n),h),dp";
  execute "minpoly=number("+mp+");";
  map f=br,maxideal(1);                // canonical mapping
  matrix M[k][m+k*n];
  M[1..k,1..m]=matrix(f(M(2)));        // will contain a module -
  ideal P=f(P);                        // primary invariants in the new ring -
  for (i=1;i<=n;i=i+1)                 // constructing a module
  { for (j=1;j<=k;j=j+1)
    { M[j,m+(i-1)*k+j]=y(i)-P[i];
    }
  }
  M=elim(module(M),1,n);               // eliminating x(1..n), std-calculation
                                       // is done internally
  M=homog(module(M),h);                // homogenize for 'minbase'
  M=minbase(module(M));
  setring br;
  //execute "ideal v="+varstr(br)+",P,1"; // dehomogenizing -
  ideal v=maxideal(1),P,1;
  f=R,v;                               // replacing y(1..n) by primary
                                       // invariants -
  M(2)=f(M);                           // M(2) is the new module -
  matrix FI(1)=matrix(P);              // getting primary invariants ready for
                                       // output
  m=ncols(M(2));
  matrix FI(2)[1][m];
  FI(2)=matrix(S)*matrix(M(2));        // FI(2) contains the real secondary
                                       // invariants
  for (i=1; i<=m;i=i+1)
  { FI(2)[1,i]=FI(2)[1,i]/leadcoef(FI(2)[1,i]); // making elements nice
  }
  FI(2)=sort(ideal(FI(2)))[1];
  if (flagvec[2])
  { "  These are the secondary invariants: ";
    for (i=1;i<=m;i=i+1)
    { "   "+string(FI(2)[1,i]);
    }
    "";
    "  We're done!";
    "";
  }
  if ((flagvec[2] or (voice==2)) && flagvec[1]==1 && (m>1))
  { "  WARNING: The invariant ring might not have a Hironaka decomposition";
    "           if the characteristic of the coefficient field divides the";
    "           group order.";
  }
  else
  { if ((flagvec[2] or (voice==2)) and (m>1))
    { "  WARNING: The invariant ring might not have a Hironaka decomposition!";
      "           This is because the characteristic of the coefficient field";
      "           divides the group order.";
    }
  }
  degBound=dB;
  return(FI(1..2));
}
example
{ "  EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
  echo=2;
           ring R=0,(x,y,z),dp;
           matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
           matrix B(1..2);
           B(1..2)=inv_ring_s(A);
           print(B(1..2));
}

////////////////////////////////////////////////////////////////////////////////
// This procedure finds a linear combination of the generators in B such that it
// lowers the dimension of the ideal generated by the primary invariants found
// so far when added to the ideal. The coefficients lie in a field of
// characteristic 0.
////////////////////////////////////////////////////////////////////////////////
proc combi (ideal B,int b,ideal P,int d)
{
  intmat vec[1][b];                    // the zero vector -
  matrix t[1][1];                      // the linear combination
  while(1)
  { vec=nextvec(vec);                  // next vector
    t=vec*transpose(matrix(B));
    if (d-1==dim(std(P+ideal(t[1,1])))) // indicates that it was not necessary
    { return(t[1,1]);                  // to break out of the for-loop
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
// This procedure trys to find a linear combination of the generators in B such
// that it lowers the dimension of the ideal generated by the primary invariants
// found so far when added to the ideal. The coefficients lie in a finite field.
// It is not clear whether such a combination exists. In the worst case, all
// possibilities are tried.
////////////////////////////////////////////////////////////////////////////////
proc p_combi (ideal B, int b, ideal P, int di)
{ def br=basering;
  matrix vec(1)[1][b];                 // starting with 0-vector -
  intmat new[1][b];                    // new vector in characteristic 0 -
  matrix pnew[1][b];                   // new needs to be mapped into br -
  int counter=1;                       // we count how many vectors we try
  int i;
  int p=char(br);
  if (minpoly<>0)
  { int d=pardeg(minpoly);             // field has p^d elements
  }
  else
  { int d=1;                           // field has p^d elements
  }
  matrix t[1][1];                      // the linear combination
  ring R=0,x,dp;
  int bound=int((number(p)^(d*b)-1)/(number(p)^d-1)+1); // this is how many
                                       // linearly independent vectors of size
                                       // b exist having entries in the base
                                       // field of br
  setring br;
  while (counter<>bound)               // otherwise, we are done
  { new=nextvec(new);
    for (i=1;i<=b;i=i+1)
    { pnew[1,i]=intnumap(new[1,i]);    // mapping an integer into br
    }
    if (unique(vec(1..counter),pnew))  // checking whether we tried pnew before
    { counter=counter+1;
      matrix vec(counter)=pnew;        // keeping track of the ones we tried -
      t=vec(counter)*transpose(matrix(B)); // linear combination -
      if (di-1==dim(std(P+ideal(t[1,1])))) // indicates that it was not
      { return(t[1,1]);                // necessary to break out of the for-loop
      }
    }
  }
  return(0);
}

////////////////////////////////////////////////////////////////////////////////
// Finds out whether any basis element of the space of homogenous invariants of
// degree g (of dimension di) is not contained in the radical of P (of the ideal
// generated by the primary invariants found so far). It uses the Reynolds
// operator. It is used to indicate when we need to check whether nontrivial
// linear combinations of basis elements exists that lower the dimension of P
// when added.
////////////////////////////////////////////////////////////////////////////////
proc search (matrix RO, ideal P, int g, int di)
{ ideal B=inv_basis_rey(RO,g,di);      // basis of homogeneous invariants of
                                       // degree g
  int mdi=ncols(B);
  int bool=0;
  for (int i=1;i<=mdi;i=i+1)
  { if (not(rad_con(B[i],P)))          // indicating that we need to try and
    { bool=1;                          // find a linear combination of basis
    }                                  // elements in B
    else
    { B[i]=0;                          // getting rid of the ones that are fully
    }                                  // contained in the radical anyway
  }
  return(bool,compress(B));            // recycle B
}

////////////////////////////////////////////////////////////////////////////////
// Finds out whether any generator in B of some space of homogenous invariants
// is not contained in the radical of P (of the ideal generated by the primary
// invariants found so far). It is used to indicate when we need to check
// whether nontrivial linear combinations of basis elements exists that lower
// the dimension of P when added.
////////////////////////////////////////////////////////////////////////////////
proc searchalt (ideal B, ideal P)
{ int mdi=ncols(B);
  int bool=0;
  for (int i=1;i<=mdi;i=i+1)
  { if (not(rad_con(B[i],P)))          // indicating that we need to try and
    { bool=1;                          // find a linear combination of basis
    }                                  // elements in B
    else
    { B[i]=0;                          // getting rid of the ones that are fully
    }                                  // contained in the radical anyway
  }
  return(bool,compress(B));
}

////////////////////////////////////////////////////////////////////////////////
// 'inv_ring_k' and 'inv_ring_s' only differ in the way they are calculating the
// primary invariants. 'inv_ring_k' tries to find a set of primary invariants of
// possibly low degree. It does this by checking whether there is a linear
// combination of basis elements of a space of homogeneous invariants in a
// certain degree, such that the dimension of the variety generated by the
// primary invariant falls each time a new primary invariant is added. And this
// way we are done looking for primary invariants precisely when n (the number
// of variables of the basering) invariants are generated.
////////////////////////////////////////////////////////////////////////////////
proc inv_ring_k (list #)
  USAGE:   inv_ring_k(<generators of a finite matrix group>[,<intvec>]);
           <intvec> has to contain 2 flags; if the first one equals 0, the
           program attempts to compute the Molien series and Reynolds operator,
           if it equals 1, the program is told that the characteristic of the
           base field divides the group order, if it is anything else the Molien
           series and Reynolds operator will not be computed; if the second flag
           does not equal 0, information about the various stages of the program
           will be printed while running
  RETURNS: generators of the invariant ring with respect to the matrix group
           generated by the matrices in the input; there are two return values
           of type <matrix>, the first containing primary invariants and the
           second secondary invariants, i.e. module generators over a Noetherian
           normalization
  EXAMPLE: example inv_ring_k; shows an example
{ def br=basering;
  int ch=char(br);                     // the algorithms depend very much on the
                                       // characteristic of the ground field
  int dB=degBound;
  degBound=0;
  int n=nvars(br);                     // n is the number of variables, as well
                                       // as the size of the matrices, as well
                                       // as the number of primary invariants,
                                       // we should get
  if (typeof(#[size(#)])=="intvec")
  { if (size(#[size(#)])<>2)
    { "  ERROR:   <intvec> must have exactly two entires";
      return();
    }
    intvec flagvec=#[size(#)];
    if (flagvec[1]==0)
    { if (ch==0)
      { matrix R(1..2);                // one will contain Reynolds operator and
                                       // the other enumerator and denominator
                                       // of Molien series
        R(1..2)=rey_mol(#[1..size(#)-1],flagvec[2]);
      }
      else
      { string newring="Qa";
        matrix R(1)=rey_mol(#[1..size(#)-1],newring,flagvec[2]); // will contain
      }                                // Reynolds operator, if Molien series
    }                                  // can be computed, it will be stored in
                                       // the new ring Qa
    else
    { for (int i=1;i<=size(#)-1;i=i+1)
      { if (not(typeof(#[i])=="matrix"))
        { "  ERROR:   the parameters must be a list of matrices and optionally";
          "           an <intvec>";
          return();
        }
        if (n<>ncols(#[i]) || n<>nrows(#[i]))
        { "  ERROR:   matrices need to be square and of the same dimensions as";
          "           the number of variables of the basering";
          return();
        }
      }
      kill i;
    }
  }
  else
  { if (typeof(#[size(#)])<>"matrix")
    { "  ERROR:   the parameters must be a list of matrices and optionally";
      "           an <intvec>";
      return();
    }
    if (ch==0)
    { matrix R(1..2);                  // will contain Reynolds operator and
                                       // enumerator and denominator of Molien
                                       // series
      R(1..2)=rey_mol(#[1..size(#)]);
    }
    else
    { string newring="Qa";             // we might need as a new ring of
                                       // characteristic 0 where we store the
                                       // Molien series -
      matrix R(1)=rey_mol(#[1..size(#)],newring); // will contain
                                       // Reynolds operator
    }
    intvec flagvec=0,0;
  }
  ideal P=0;                           // will contain primary invariants
  if (flagvec[1]==0 && flagvec[2])
  { "  We can start looking for primary invariants...";
    "";
  }
  else
  { if (flagvec[1] && flagvec[2])
    { "";
      "  We start by looking for primary invariants...";
      "";
    }
  }
  if ((ch==0 || defined(Qa)) && flagvec[1]==0) // i.e. we can use Molien series
  { if (ch==0)
    { poly p(1..2);                    // p(1) will be used for single terms of
                                       // the partial expansion, p(2) to store
      p(1..2)=part_mol(R(2),1);        // the intermediate result -
      poly v1=var(1);                  // we need v1 to split off coefficients
                                       // in the partial expansion of M (which
                                       // is in terms of the first variable) -
      poly d;                          // for splitting off the coefficient in
                                       // in one term of the partial expansion,
                                       // i.e. it stores the dimension of the
                                       // current homogeneous subspace
    }
    else
    { setring Qa;                      // Qa is where the Molien series is
                                       // stored -
      poly p(1..2);                    // p(1) will be used for single terms of
                                       // the partial expansion, p(2) to store
      p(1..2)=part_mol(M,1);           // the intermediate result -
      poly d;                          // stores the dimension of the current
                                       // homogeneous subspace
      setring br;
    }
    int g, di, counter, i, j, m, bool; // g: current degree, di: d as integer,
                                       // counter: counts primary invariants in
                                       // degree g, i,j: going through monomials
                                       // of degree g, m: counting primary
                                       // invariants, bool: indicates whether
                                       // the case occurred that a new
                                       // polynomial did not lower the
                                       // dimension of the ideal generated by
                                       // previously found invariants -
    poly imRO;                         // the image of the Reynolds operator -
    ideal mon;                         // will contain monomials of degree g -
    while(1)                           // repeat until n polynomials are found
    { if (ch==0)
      { p(1..2)=part_mol(R(2),1,p(2)); // 1 term of the partial expansion -
        g=deg(p(1));                   // current degree -
        d=coef(p(1),v1)[2,1];          // dimension of invariant space of degree
                                       // g -
        di=int(d);                     // just a type cast
      }
      else
      {  setring Qa;
         p(1..2)=part_mol(M,1,p(2));   // 1 term of the partial expansion -
         g=deg(p(1));                  // current degree -
         d=coef(p(1),x)[2,1];          // dimension of invariant space of degree
                                       // g -
         di=int(d);                    // just a type cast
         setring br;
      }
      if (flagvec[2])
      { "  Searching for primary invariants in degree "+string(g)+":";
        "  There is/are "+string(di)+" linearly independent invariant(s) to choose from...";
      }
      mon=sort(maxideal(g))[1];        // all monomials of degree g -
      j=ncols(mon);
      counter=0;                       // we have 0 candidates of degree g so
                                       // far
      for (i=j;i>=1;i=i-1)
      { imRO=eval_rey(R(1),mon[i]);
        if (reduce(imRO,std(P))<>0)
        { if (P[1]==0)                 // if imRO is the first non-zero
          { counter=1;                 // invariant we find, the dim question is
            m=1;                       // trivial and we just include imRO
            P[1]=imRO/leadcoef(imRO);
            if (flagvec[2])
            { "  We find: "+string(P[1]);
            }
            if (counter==di)           // if counter is up to di==d, we can
            { break;                   // leave the for-loop
            }
          }
          else
          { P=P,imRO;                  // we add imRO to the generators of P
            attrib(P,"isSB",1);
            if (n-m-1<dim(P))          // here we are checking whether the
            { if (n-m-1<dim(std(P)))   // dimension is really going down with
              { P[m+1]=0;              // the new polynomial -
                P=compress(P);         // if the dimension does not go down
                                       // we get rid of imRO again -
                bool=1;                // we will have to go into the procedure
                                       // search later
              }
              else                     // we can keep imRO -
              { counter=counter+1;
                m=m+1;
                P[m]=P[m]/leadcoef(P[m]); // making m-th primary invariant
                if (flagvec[2])        // nice
                { if (counter<>1)
                  { "           "+string(P[m]);
                  }
                  else
                  { "  We find: "+string(P[m]);
                  }
                }
              }
            }
            else                       // we can keep imRO -
            { counter=counter+1;
              m=m+1;
              P[m]=P[m]/leadcoef(P[m]); // making m-th primary invariant
              if (flagvec[2])
              { if (counter<>1)
                { "           "+string(P[m]);
                }
                else
                { "  We find: "+string(P[m]);
                }
              }
            }
            if (n==m or (counter==di)) // if counter==di, we can leave the for
            { break;                   // loop; if n==m, we can leave both loops
            }
          }
        }
      }
      if (n==1 or n==m)
      { break;
      }
      if (bool)
      { if (not(defined(B)==voice))
        { ideal B;                     // will contain a subset of a basis of
          int T;                       // homogeneous invariants of degree g
          ideal Palt;                  // such that none is contained in the
          poly lin;                    // radical of P -
        }
        bool,B=search(R(1),P,g,di);    // checking whether we need to consider
                                       // nontrivial linear combinations of
                                       // basis elements of degree g
        di=ncols(B);
        counter=0;
      }
      if (bool && (di>1))              // indicates that some invariants are not
      {                                // in the radical, but don't lower the
                                       // dimension, if there is one element in
                                       // B, then there exists no linear
                                       // combination that lowers the dimension
        Palt=P,B;
        T=n-m-dim(std(Palt));
        while ((counter<>T) && (m<>n)) // runs until we are sure that there are
        {                              // no more primary invariant of this
                                       // degree -
                                       // otherwise we have to try and build a
                                       // sum of the basis elements of this
                                       // degree -
          if (ch==0)                   // we have to distinguish prime and non
          {                            // prime characteristic, in infinite
                                       // fields a (non-)solution is guaranteed
                                       // and here a systematic way of finding
                                       // such a solution is implemented -
            lin=combi(B,di,P,n-m);     // combi finds a combination
          }
          else
          { lin=p_combi(B,di,P,n-m);   // the subroutine p_combi finds out
                                       // whether there is a combination of the
                                       // basis elements at all such that it
                                       // lowers the dimension of P when added -
            if (lin==0)                // if the 0-polynomial is returned, it
            { break;                   // means that there was no combination -
            }
          }
          m=m+1;
          P[m]=lin;                    // we did find the combination lin
          if (flagvec[2])
          { "  We find: "+string(P[m]);
          }
          counter=counter+1;
        }
      }
      bool=0;
      if (m==n)                        // found all primary invariants
      { break;
      }
      if (flagvec[2])
      { "";
      }
    }
  }
  else
  {                                    // this case is entered when either the
                                       // characteristic<>0 divides the group
                                       // order or when the Molien series could
                                       // not or has not been computed -
    if (flagvec[1]==0)                 // indicates that the group order is
    { int g=nrows(R(1));               // known, here it is set to g -
      int flag=((g%ch)==0);            // flag is 1 if the characteristic
                                       // divides the order, it is 0 if it does
                                       // not -
      if (typeof(#[size(#)])=="intvec") // getting ahold of the generators of
      { int gennum=size(#)-1;          // the generators of the group
      }
      else
      { int gennum=size(#);
      }
    }
    else
    { int flag=2;                      // flag is 2 if we don't know yet whether
      int gennum=size(#)-1;            // the group order is divisible by the
    }                                  // characteristic -
    int d=1;                           // d is set to the current degree, since
                                       // we know nothing about the finite
                                       // matrix group (via Molien series) we
                                       // have to start with degree 1
    int j, counter, i, di, bool;       // counter: counts primary invariants,
                                       // i: goes through basis elements, di:
                                       // dimension of current space, bool:
                                       // indicates that the case occurred that
                                       // a basis element did not lower the
                                       // dimension, but was not in the radical
    while (1)
    { if (flagvec[2])
      { "  Searching for primary invariants in degree "+string(d)+":";
      }
      if (flag)                        // in this case we can not make use of
      {                                // the Reynolds operator -
        ideal B(d)=inv_basis(d,#[1..gennum]); // we create a basis of the vector
      }                                // space of all invariant polynomials of
                                       // degree d
      else
      {                                // here the characteristic<>0 does not
        ideal B(d)=inv_basis_rey(R(1),d); // divide the group order, i.e. the
      }                                // Reynolds operator can be used to
                                       // calculate a basis of the vector space
                                       // of all invariant polynomials of degree
                                       // d -
      di=ncols(B(d));
      if (B(d)[1]<>0)                  // otherwise the space is empty -
      { if (flagvec[2])
        { "  There is/are "+string(di)+" linearly independent invariant(s) to choose from...";
        }
        if (counter==0)                // we have no candidates for primary
        {                              // invariants yet, i.e. don't have to
          P[1]=B(d)[1];                // check for radical containment
          if (flagvec[2])
          { "  We find: "+string(P[1]);
          }
          i=2;                         // go to second basis element
          counter=1;
        }
        else
        { i=1;                         // go to first basis element
        }
        while (i<=di)                  // goes through all polynomials in B(d) -
        { P=P,B(d)[i];                 // adding candidate -
          attrib(P,"isSB",1);          // checking dimension -
          if (n-counter-1<dim(P))
          { if (n-counter-1<dim(std(P))) // in this case B(d)[i] would not lower
            { P[counter+1]=0;          // the dimension and we get rid of it
              P=compress(P);
              bool=1;
            }
            else                       // indicates that B(d)[i] qualifies
            { counter=counter+1;
              if (flagvec[2])
              { "  We find: "+string(P[counter]);
              }
              if (counter==n)          // in that case, we're done
              { break;
              }
            }
          }
          else                         // indicates that B(d)[i] qualifies
          { counter=counter+1;
            if (flagvec[2])
            { "  We find: "+string(P[counter]);
            }
            if (counter==n)            // in that case, we're done
            { break;
            }
          }
          i=i+1;                       // go to next basis element
        }
        if (counter==n)                // we're done
        { break;
        }
        if (bool)
        { if (not(defined(Ba)==voice))
          { ideal Ba;
            int T;
            ideal Palt;
            poly lin;
          }
          bool,Ba=searchalt(B(d),P);   // Ba will now contain a subset of
                                       // a basis of homogeneous invariants of
                                       // degree d such that none is contained
                                       // in the radical of P
          di=ncols(Ba);
        }
        if (bool && (di>1))            // this meant that we have to use
        {                              // Kemper's method, if there is one
                                       // element in Ba then there exists no
                                       // linear combination that lowers the
                                       // dimension
          Palt=P,Ba;
          T=n-counter-dim(std(Palt));
          while (counter<>n)           // runs until we are sure that there are
          {                            // no more primary invariant of this
                                       // degree -
                                       // otherwise we have to try and build a
                                       // sum of the basis elements of this
                                       // degree -
            if (ch==0)                 // we have to distinguish prime and non
            {                          // prime characteristic, in infinite
                                       // fields a (non-)solution is guaranteed
                                       // and here a systematic way of finding
                                       // such a solution is implemented -
              lin=combi(Ba,di,P,counter); // combi finds a combination
            }
            else
            { lin=p_combi(Ba,di,P,counter); // the subroutine p_combi finds out
                                       // whether there is a combination of the
                                       // basis elements at all such that it
                                       // lowers the dimension of P when added -
              if (lin==0)              // if the 0-polynomial is returned, it
              { break;
              }
            }
            counter=counter+1;         // otherwise, we did find a combination
            P[counter]=lin;
            if (flagvec[2])
            { "  We find: "+string(P[counter]);
            }
          }
        }
        bool=0;
        if (counter==n)                // found all primary invariants
        { break;
        }
        if (flagvec[2])
        { "";
        }
      }
      else
      { if (flagvec[2])
        { "  The space is 0-dimensional.";
        }
      }
      d=d+1;                           // up to the next degree
    }
  }
  if ((ch==0 || defined(Qa)) && flagvec[1]==0)
  { if (flagvec[2])
    { "";
    }
    ideal Q=std(P);                    // P contains the primary invariants -
    intvec degvec;                     // will contain the degrees of secondary
                                       // invariants -
    if (ch==0)                         // Molien series is stored in basering
    { p(1)=1;
      for (j=1;j<=n;j=j+1)             // calculating the denominator of the
      { p(1)=p(1)*(1-v1^deg(P[j]));    // Hilbert series of the ring generated
      }                                // generated by the primary invariants -
      matrix s[1][2]=R(2)[1,1]*p(1),R(2)[1,2]; // used for canceling
      s=matrix(syz(ideal(s)));
      p(1)=s[2,1];                     // the polynomial telling us where to
                                       // search for secondary invariants
      map slead=br,ideal(0);
      p(1)=1/slead(p(1))*p(1);         // smallest term of p(1) needs to be 1 -
      if (flagvec[2])
      { "  Polynomial telling us where to look for secondary invariants:";
        "   "+string(p(1));
        "";
      }
      matrix dimmat=coeffs(p(1),v1);   // dimmat will contain the number of
                                       // secondary invariants, we need to find
                                       // of a certain degree -
      m=nrows(dimmat);                 // m-1 is the highest degree
      degvec=0;
      for (j=1;j<=m;j=j+1)
      { if (dimmat[j,1]<>0)
        { degvec[j]=int(dimmat[j,1]);  // degvec[j] now contains the number of
        }                              // secondary invariants of degree j-1
      }
    }
    else
    { for (j=1;j<=n;j=j+1)             // degrees have to be taken in a ring of
      { degvec[j]=deg(P[j]);           // characteristic 0
      }
      setring Qa;
      p(1)=1;
      for (j=1;j<=n;j=j+1)             // calculating the denominator of the
      { p(1)=p(1)*(1-x^degvec[j]);     // Hilbert series of the ring generated
      }                                // by the primary invariants -
      matrix s[1][2]=M[1,1]*p(1),M[1,2]; // used for canceling
      s=matrix(syz(ideal(s)));
      p(1)=s[2,1];                     // the polynomial telling us where to
                                       // search for secondary invariants
      map slead=Qa,ideal(0);
      p(1)=1/slead(p(1))*p(1);         // smallest term of p(1) needs to be 1
      if (flagvec[2])
      { "  Polynomial telling us where to look for secondary invariants:";
        "   "+string(p(1));
        "";
      }
      matrix dimmat=coeffs(p(1),x);    // dimmat will contain the number of
                                       // secondary invariants, we need to find
                                       // find of a certain degree -
      m=nrows(dimmat);                 // m-1 actually is the highest degree
      degvec=0;
      for (j=1;j<=m;j=j+1)
      { if (dimmat[j,1]<>0)
        { degvec[j]=int(dimmat[j,1]);  // degvec[j-1] contains the number of
        }                              // secondary invariants of degree j-1
      }
      setring br;
      kill Qa;                         // all the information needed for Qa is
    }                                  // stored in degvec and dimmat -
    qring Qring=Q;                     // we need to do calculations modulo the
                                       // ideal generated by the elements of P,
                                       // its standard basis is stored in Q -
    poly imROmod;                      // imRO reduced -
    ideal Smod, sSmod;                 // secondary invariants of one degree
                                       // reduced and their standard basis
    setring br;
    kill Q;                            // Q might be big and isn't needed
                                       // anymore
    if (flagvec[2])
    { "  Proceeding to look for secondary invariants...";
      "";
      "  In degree 0 we have: 1";
      "";
    }
    bool=0;                            // indicates when standard basis
                                       // calculation is necessary -
    ideal S=1;                         // 1 definitely is a secondary invariant
    for (i=2;i<=m;i=i+1)               // walking through degvec -
    { if (degvec[i]<>0)                // when it is == 0 we need to find 0
      {                                // elements of the current degree being
                                       // i-1 -
        if (flagvec[2])
        { "  Searching in degree "+string(i-1)+", we need to find "+string(degvec[i,1])+" invariant(s)...";
        }
        mon=sort(maxideal(i-1))[1];    // all monomials of degree i-1 -
        counter=0;                     // we'll count up to degvec[i] -
        j=ncols(mon);                  // we'll go through mon from the end
        setring Qring;
        Smod=0;
        setring br;
        while (degvec[i]<>counter)     // need to find degvec[i] linearly
        {                              // independent (in Qring) invariants -
          imRO=eval_rey(R(1),mon[j]);  // generating invariants
          setring Qring;
          imROmod=fetch(br,imRO);      // reducing the invariants
          if (reduce(imROmod,std(ideal(0)))<>poly(0) and counter<>0)
          {                            // if the first condition is true and the
                                       // second false, imROmod is the first
                                       // secondary invariant of that degree
                                       // that we want to add and we need not
                                       // check linear independence
            if (bool)
            { sSmod=std(Smod);
            }
            if (reduce(imROmod,sSmod)<>0)
            { Smod=Smod,imROmod;
              setring br;              // we make its leading coefficient to be
              imRO=imRO/leadcoef(imRO); // 1
              S=S,imRO;
              if (flagvec[2])
              { "           "+string(imRO);
              }
              counter=counter+1;
              bool=1;                  // next time we need to recalculate std
            }
            else
            { bool=0;                  // std-calculation is unnecessary
              setring br;
            }
          }
          else
          { if (reduce(imROmod,std(ideal(0)))<>poly(0) and counter==0)
            { Smod[1]=imROmod;         // here we just add imRO(mod) without
              setring br;              // having to check linear independence
              imRO=imRO/leadcoef(imRO);
              S=S,imRO;
              counter=counter+1;
              bool=1;                  // next time we need to calculate std
              if (flagvec[2])
              { "  We find: "+string(imRO);
              }
            }
            else
            { setring br;
            }
          }
          j=j-1;                       // going to next monomial
        }
        if (flagvec[2])
        { "";
        }
      }
    }
    degBound=dB;
    if (flagvec[2])
    { "  We're done!";
      "";
    }
    matrix FI(1)=matrix(P);
    matrix FI(2)=matrix(S);
    return(FI(1..2));
  }
  else
  { if (flagvec[2])
    { "";
      "  Proceeding to look for secondary invariants...";
    }
    // we can now proceed to calculate secondary invariants, the problem
    // we face again is that we can make no use of a Molien series - however,
    // if the characteristic does not divide the group order, we can still make
    // use of the fact that the secondary invariants are free module generators
    // and that we need deg(P[1])*...*deg(P[n])/(cardinality of the group) of
    // them
    matrix FI(1)=matrix(P);            // primary invariants, ready for output -
    P=std(P);                          // for calculations module primary
                                       // invariants
    if (flagvec[1]<>0 and flagvec[1]<>1)
    { int g=group(#[1..size(#)-1]);    // computing group order
      if (ch==0)
      { matrix FI(2)=sec_minus_mol(ideal(FI(1)),P,g,flagvec[2],#[1..size(#)-1],0,B(1..d),d);
        return(FI(1..2));
      }
      if (g%ch<>0)
      { matrix FI(2)=sec_minus_mol(ideal(FI(1)),P,g,flagvec[2],#[1..size(#)-1],0,B(1..d),d);
          return(FI(1..2));
      }
    }
    else
    { if (flag==0)                     // this is the case where we have a
      {                                // nonzero minpoly, but the
                                       // characteristic does not divide the
                                       // group order
        matrix FI(2)=sec_minus_mol(ideal(FI(1)),P,g,flagvec[2],R(1),1,B(1..d),d);
        return(FI(1..2));
      }
    }
    if (flagvec[2])
    { "  Since the characteristic of the base field divides the group order, we do not";
      "  know whether the invariant ring is Cohen-Macaulay. We have to use Kemper's";
      "  algorithm and compute secondary invariants with respect to the trivial";
      "  subgroup of the given group.";
      "";

    }
    // are using Kemper's algorithm with the trivial subgroup
    ring QQ=0,x,dp;
    ideal M=(1-x)^n;                   // we look at our primary invariants as
                                       // such of the subgroup that only
                                       // contains the identity, this means that
                                       // ch does not divide the order anymore,
                                       // this means that we can make use of the
                                       // Molien series again - 1/M[1] is the
                                       // Molien series of that group, we now
                                       // calculate the secondary invariants of
                                       // this subgroup in the usual fashion
                                       // where the primary invariants are the
                                       // ones from the bigger group
    setring br;
    intvec degvec;                     // for the degrees of the primary
                                       // invariants -
    for (i=1;i<=n;i=i+1)               // finding the degrees of these
    { degvec[i]=deg(FI(1)[1,i]);
    }
    setring QQ;                        // calculating the polynomial indicating
    M[2]=1;                            // where to search for secondary
    for (i=1;i<=n;i=i+1)               // invariants (of the trivial subgroup)
    { M[2]=M[2]*(1-x^degvec[i]);
    }
    M=matrix(syz(M))[1,1];
    M[1]=M[1]/leadcoef(M[1]);
    if (flagvec[2])
    { "  Polynomial telling us where to look for these secondary invariants:";
      "   "+string(M[1]);
      "";
    }
    matrix dimmat=coeffs(M[1],x);      // storing the number of secondary
                                       // invariants we need in a certain
    int m=nrows(dimmat);               // m-1 is the highest degree where we
                                       // need to search
    degvec=0;
    for (i=1;i<=m;i=i+1)               // degvec will contain all the
    { if (dimmat[i,1]<>0)              // information about where to find
      { degvec[i]=int(dimmat[i,1]);    // secondary invariants, it is filled
      }                                // with integers and therefore visible in
    }                                  // all rings
    kill QQ;
    setring br;
    ideal B;
    ideal S=1;                         // 1 is a secondary invariant always
    if (flagvec[2])
    { "  In degree 0 we have: 1";
      "";
    }
    qring Qring=P;
    ideal Smod, Bmod, sSmod;           // Smod: secondary invariants of one
                                       // degree modulo P, sSmod: standard basis
                                       // of the latter, Bmod: B modulo P
    setring br;
    kill P;                            // might be large
    if (flagvec[1]==1)
    { int g;
    }
    for (i=2;i<=m;i=i+1)               // going through all entries of degvec
    { if (degvec[i]<>0)
      { B=sort(maxideal(i-1))[1];      // basis of the space of invariants (with
                                       // respect to the matrix subgroup
                                       // containing only the identity) of
                                       // degree i-1 -
        if (flagvec[2])
        { "  Searching in degree "+string(i-1)+", we need to find "+string(degvec[i])+" invariant(s)...";
        }
        counter=0;                     // we have 0 secondary invariants of
                                       // degree i-1 so far
        setring Qring;
        Bmod=fetch(br,B);              // basis modulo primary invariants
        Smod=0;
        j=ncols(Bmod);                 // going backwards through Bmod
        while (degvec[i]<>counter)
        { if (reduce(Bmod[j],std(ideal(0)))<>0 && counter<>0)
          { if (bool)
            { sSmod=std(Smod);
            }
            if (reduce(Bmod[j],sSmod)<>0) // Bmod[j] qualifies as secondary
            { Smod=Smod,Bmod[j];       // invariant
              setring br;
              S=S,B[j];
              counter=counter+1;
              if (flagvec[2])
              { "           "+string(B[j]);
              }
              setring Qring;
              bool=1;                  // need to calculate std of Smod next
            }                          // time
            else
            { bool=0;
            }
          }
          else
          { if (reduce(Bmod[j],std(ideal(0)))<>0 && counter==0)
            { Smod[1]=Bmod[j];         // in this case, we may just add B[j]
              setring br;
              S=S,B[j];
              if (flagvec[2])
              { "  We find: "+string(B[j]);
              }
              counter=counter+1;
              bool=1;                  // need to calculate std of Smod next
              setring Qring;           // time
            }
          }
          j=j-1;                       // next basis element
        }
        setring br;
      }
    }
    // now we have those secondary invariants
    int k=ncols(S);                    // k is the number of the secondary
                                       // invariants, we just calculated
    if (flagvec[2])
    { "";
      "  We calculate secondary invariants from the ones found for the trivial";
      "  subgroup.";
      "";
    }
    map f;                             // used to let generators act on
                                       // secondary invariants with respect to
                                       // the trivial group -
    matrix M(1)[gennum][k];            // M(1) will contain a module
    for (i=1;i<=gennum;i=i+1)
    { B=ideal(matrix(maxideal(1))*transpose(#[i])); // image of the various
                                       // variables under the i-th generator -
      f=br,B;                          // the corresponding mapping -
      B=f(S)-S;                        // these relations should be 0 -
      M(1)[i,1..k]=B[1..k];            // we will look for the syzygies of M(1)
    }
    module M(2)=res(M(1),2)[2];
    m=ncols(M(2));                     // number of generators of the module
                                       // M(2) -
    // the following steps calculates the intersection of the module M(2) with
    // the algebra A^k where A denote the subalgebra of the usual polynomial
    // ring, generated by the primary invariants
    string mp=string(minpoly);         // generating a ring where we can do
                                       // elimination
    execute "ring R=("+charstr(br)+"),(x(1..n),y(1..n),h),dp;";
    execute "minpoly=number("+mp+");";
    map f=br,maxideal(1);              // canonical mapping
    matrix M[k][m+k*n];
    M[1..k,1..m]=matrix(f(M(2)));      // will contain a module -
    matrix FI(1)=f(FI(1));             // primary invariants in the new ring
    for (i=1;i<=n;i=i+1)
    { for (j=1;j<=k;j=j+1)
      { M[j,m+(i-1)*k+j]=y(i)-FI(1)[1,i];
      }
    }
    M=elim(module(M),1,n);             // eliminating x(1..n), std-calculation
                                       // is done internally -
    M=homog(module(M),h);              // homogenize for 'minbase'
    M=minbase(module(M));
    setring br;
    //execute "ideal v="+varstr(br)+",ideal(FI(1)),1";
    ideal v=maxideal(1),ideal(FI(1)),1;
    f=R,v;                             // replacing y(1..n) by primary
                                       // invariants -
    M(2)=f(M);                         // M(2) is the new module -
    m=ncols(M(2));
    matrix FI(2)[1][m];
    FI(2)=matrix(S)*matrix(M(2));      // FI(2) now contains the secondary
                                       // invariants
    for (i=1; i<=m;i=i+1)
    { FI(2)[1,i]=FI(2)[1,i]/leadcoef(FI(2)[1,i]); // making elements nice
    }
    FI(2)=sort(ideal(FI(2)))[1];
    if (flagvec[2])
    { "  These are the secondary invariants: ";
      for (i=1;i<=m;i=i+1)
      { "   "+string(FI(2)[1,i]);
      }
      "";
      "  We're done!";
      "";
    }
    if ((flagvec[2] or (voice==2)) && flagvec[1]==1 && (m>1))
    { "  WARNING: The invariant ring might not have a Hironaka decomposition";
      "           if the characteristic of the coefficient field divides the";
      "           group order.";
    }
    else
    { if ((flagvec[2] or (voice==2)) and (m>1))
      { "  WARNING: The invariant ring might not have a Hironaka decomposition!"
;
        "           This is because the characteristic of the coefficient field"
;
        "           divides the group order.";
      }
    }
    degBound=dB;
    return(FI(1..2));
  }
}
example
{ "  EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
  echo=2;
           ring R=0,(x,y,z),dp;
           matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
           matrix B(1..2);
           B(1..2)=inv_ring_k(A);
           print(B(1..2));
}

////////////////////////////////////////////////////////////////////////////////
// The procedure introduces m new variables y(i), m being the number of
// generators {f_1,...,f_m} of the subring in the variables x(1),...,x(n).
// A Groebner basis of {f_1-y(1),...,f_m-y(m)} is computed with respect to
// the product ordering of x^a*y^b > y^d*y^e if x^a > x^d or else if y^b > y^e.
// f reduces to a polynomial only in the y(i) <=> p is contained in the subring
// generated by the polynomials in F.
////////////////////////////////////////////////////////////////////////////////
proc algebra_con (poly p, matrix F)
  USAGE:   algebra_con(<poly>,<matrix>); <poly> is arbitrary in the basering,
           <matrix> defines a subring of the basering
  RETURNS: if <poly> is contained in the ring, 1 (TRUE) (type <int>) is
           returned as well as a comment showing a representation of <poly>
           where y(i) represents the i-th element in <matrix>. 0 (type <int>)
           is returned if <poly> is not contained
  EXAMPLE: example algebra_con; shows an example
{ if (nrows(F)==1)
  { def br=basering;
    int n=nvars(br);
    int m=ncols(F);
    ring R=0,(x(1..n),y(1..m)),(dp(n),dp(m));
    ideal vars=x(1..n);
    map emb=br,vars;
    ideal F=ideal(emb(F));
    ideal check=emb(p);
    for (int i=1;i<=m;i=i+1)
    { F[i]=F[i]-y(i);
    }
    F=std(F);
    check[1]=reduce(check[1],F);
    F=elim(check,1,n);
    if (F[1]<>0)
    { "\/\/ "+string(check);
      return(1);
    }
    else
    { return(0);
    }
  }
  else
  { "  ERROR:   <matrix> may only have one row";
    return();
  }
}
example
{ "  EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
  echo=2;
           ring R=0,(x,y,z),dp;
           matrix F[1][7]=x2+y2,z2,x4+y4,x4+y4,1,x2z-1y2z,xyz,x3y-1xy3;
           poly p1=(x2z-1y2z)*z2;
           algebra_con(p1,F);
           poly p2=z;
           algebra_con(p2,F);
}

////////////////////////////////////////////////////////////////////////////////
// The procedure introduces n+m new variables y(i) and z(j), n being the number
// of primary generators {p_1,...,p_n} and m the number of secondary ones
// {s_1,...,s_m} in the variables x(1),...,x(n). A Groebner basis of
// {p_1-y(1),...,p_n-y(n),s_1-z(1),...,s_m-z(m)} is computed with respect to the
// product ordering of x^a*y^b*z^c > x^d*y^e*z^f if x^a > x^d with respect
// to the purely lexicographical ordering or else if z^c > z^f with respect
// to the degree lexicographical ordering or else if y^b > y^e with respect
// to the purely lexicographical ordering again. f reduces to a polynomial
// only in y(i) and z(j) (more specifically, linear in the z(j)) <=> f is
// contained in the Cohen-Macaulay ring.
////////////////////////////////////////////////////////////////////////////////
proc module_con(poly f, matrix P, matrix S)
  USAGE:   module_con(<poly>,<matrix_1>,<matrix_2>); <poly> is arbitrary in
           the basering, <matrix_1> should represent the primary generators of
           a Cohen-Macaulay ring, <matrix_2> the secondary ones
  RETURNS: if <poly> is contained in the ring, 1 (TRUE) (type <int>) is
           returned as well as a comment showing the unique representation
           of <poly> with respect to a Hironaka decomposition; y(i) represents
           the i-th element in <matrix_2> and z(j) represents the j-th element
           in <matrix_1>. 0 (type <int>) is returned if <poly> is not contained.
  EXAMPLE: example module_con; shows an example
{ def br=basering;
  int n=nvars(br);
  if (ncols(P)==n and nrows(P)==1 and nrows(S)==1)
  { int m=ncols(S);
    ring R=0,(x(1..n),y(1..m),z(1..n)),(lp(n),dp(m),lp(n));
    ideal vars=x(1..n);
    map emb=br,vars;
    matrix P=emb(P);
    matrix S=emb(S);
    ideal check=emb(f);
    ideal I;
    for (int i=1;i<=m;i=i+1)
    { I[i]=S[1,i]-y(i);
    }
    for (i=1;i<=n;i=i+1)
    { I[n+i]=P[1,i]-z(i);
    }
    I=std(I);
    check[1]=reduce(check[1],I);
    I=elim(check,1,n);                 // checking whether all variables from
    if (I[1]<>0)                       // the former ring have disappeared
    { "\/\/ "+string(check);
      return(1);
    }
    else
    { return(0);
    }
  }
  else
  { "  ERROR:   <matrix_1> must have the same number of columns as the basering";
    "           and both <matrix_1> and <matrix_2> may only have one row";
    return();
  }
}
example
{ "  EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
  echo=2;
           ring R=0,(x,y,z),dp;
           matrix P[1][3]=x2+y2,z2,x4+y4;
           matrix S[1][4]=1,x2z-1y2z,xyz,x3y-1xy3;
           poly p1=(x2z-1y2z)*xyz;
           module_con(p1,P,S);
           poly p2=z;
           module_con(p2,P,S);
}

////////////////////////////////////////////////////////////////////////////////
// 'orbit_var' calculates the syzygy ideal of the generators of the
// invariant ring, then eliminates the variables of the original ring and
// the polynomials that are left over define the orbit variety
////////////////////////////////////////////////////////////////////////////////
proc orbit_var (matrix F,string newring)
  USAGE:   orbit_var(<matrix>,<string>); <matrix> defines an invariant ring,
           <string> is the name for a new ring
  RETURN:  a Groebner basis (type <ideal>, named G) for the ideal defining the
           orbit variety (i.e. the syzygy ideal) in the new ring (named
           <string>)
  EXAMPLE: example orbit_var; shows an example
{ if (newring=="")
  { "  ERROR:   the second argument may not be an empty <string>";
    return();
  }
  if (nrows(F)==1)
  { def br=basering;
    int n=nvars(br);
    int m=ncols(F);
    string mp=string(minpoly);
    execute "ring R=("+charstr(br)+"),("+varstr(br)+",y(1..m)),dp;";
    execute "minpoly=number("+mp+");";
    ideal I=ideal(imap(br,F));
    for (int i=1;i<=m;i=i+1)
    { I[i]=I[i]-y(i);
    }
    I=elim(I,1,n);
    execute "ring "+newring+"=("+charstr(br)+"),(y(1..m)),dp(m);";
    execute "minpoly=number("+mp+");";
    ideal vars;
    for (i=2;i<=n;i=i+1)
    { vars[i]=0;
    }
    vars=vars,y(1..m);
    map emb=R,vars;
    ideal G=emb(I);
    kill emb, vars, R;
    keepring `newring`;
    // execute "keepring "+newring+";";
    return();
  }
  else
  { "  ERROR:   the <matrix> may only have one row";
    return();
  }
}
example
{ "  EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
  echo=2;
           ring R=0,(x,y,z),dp;
           matrix F[1][7]=x2+y2,z2,x4+y4,1,x2z-1y2z,xyz,x3y-1xy3;
           string newring="E";
           orbit_var(F,newring);
           print(G);
           basering;
}

////////////////////////////////////////////////////////////////////////////////
// Let f1,...,fm be generators of the invariant ring, y1,...,ym new variables
// and h1,...,hk generators of I. 'rel_orbit_var'  computes a standard basis of
// the ideal generated by f1-y1,...,fm-ym with respect to a pure lexicographic
// order. Further, a standard basis of the the ideal generated by the elements
// of the previously found standard basis and h1,...,hk is found. Eliminating
// the original variables yields generators of the relative orbit variety with
// respect to I.
////////////////////////////////////////////////////////////////////////////////
proc rel_orbit_var(ideal I,matrix F, string newring)
  USAGE:   rel_orbit_var(<ideal>,<matrix>,<string>); <ideal> defines an
           ideal invariant under the action of a group, <matrix> defines the
           invariant ring of this group, <string> is a name for a new ring
  RETURN:  a Groebner basis (type <ideal>, named G) for the ideal defining the
           relative orbit variety with respect to <ideal> in the new ring (named
           <string>)
  EXAMPLE: example rel_orbit_var; shows an example
{ if (newring=="")
  { "  ERROR:   the third argument may not be empty a <string>";
    return();
  }
  if (nrows(F)==1)
  { def br=basering;
    int n=nvars(br);
    int m=ncols(F);
    string mp=string(minpoly);
    execute "ring R=("+charstr(br)+"),("+varstr(br)+",y(1..m)),lp;";
    execute "minpoly=number("+mp+");";
    ideal J=ideal(imap(br,F));
    ideal I=imap(br,I);
    for (int i=1;i<=m;i=i+1)
    { J[i]=J[i]-y(i);
    }
    J=std(J);
    J=J,I;
    option(redSB);
    J=std(J);
    ideal vars;
    //for (i=1;i<=n;i=i+1)
    //{ vars[i]=0;
    //}
    vars[n]=0;
    vars=vars,y(1..m);
    map emb=R,vars;
    ideal G=emb(J);
    J=J-G;
    for (i=1;i<=ncols(G);i=i+1)
    { if (J[i]<>0)
      { G[i]=0;
      }
    }
    G=compress(G);
    execute "ring "+newring+"=("+charstr(br)+"),(y(1..m)),lp;";
    execute "minpoly=number("+mp+");";
    ideal vars;
    for (i=2;i<=n;i=i+1)
    { vars[i]=0;
    }
    vars=vars,y(1..m);
    map emb=R,vars;
    ideal G=emb(G);
    kill vars, emb;
    keepring `newring`;
    // execute "keepring "+newring+";";
    return();
  }
  else
  { "  ERROR:   the <matrix> may only have one row";
    return();
  }
}
example
{ "  EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.6.3:";
  echo=2;
           ring R=0,(x,y,z),dp;
           matrix F[1][3]=x+y+z,xy+xz+yz,xyz;
           ideal I=x2+y2+z2-1,x2y+y2z+z2x-2x-2y-2z,xy2+yz2+zx2-2x-2y-2z;
           string newring="E";
           rel_orbit_var(I,F,newring);
           print(G);
           basering;
}

////////////////////////////////////////////////////////////////////////////////
// Let f1,...,fm be generators of the invariant ring, y1,...,ym new variables
// and h1,...,hk generators of I. 'im_of_var' calls 'rel_orbit_var' with input
// I, F and the string newring. In the output the variables y1,...,ym are
// replaced by f1,...,fm. The result is the output of 'im_of_var' and defines
// the variety under the matrix group.
////////////////////////////////////////////////////////////////////////////////
proc im_of_var(ideal I,matrix F)
  USAGE:   im_of_var(<ideal>,<matrix>); <ideal> is arbitrary, <matrix>
           defines an invariant ring of a certain matrix group
  RETURN:  the <ideal> defining the image of the variety defined by the input
           ideal with respect to that group
  EXAMPLE: example im_of_var; shows an example
{ if (nrows(F)==1)
  { def br=basering;
    int n=nvars(br);
    string newring="E";
    rel_orbit_var(I,F,newring);
    execute "ring R=("+charstr(br)+"),("+varstr(br)+","+varstr(E)+"),lp;";
    ideal F=imap(br,F);
    for (int i=1;i<=n;i=i+1)
    { F=0,F;
    }
    setring br;
    map emb2=E,F;
    return(compress(emb2(G)));
  }
  else
  { "  ERROR:   the <matrix> may only have one row";
    return();
  }
}
example
{ "  EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.6.8:";
  echo=2;
           ring R=0,(x,y,z),dp;
           matrix F[1][3]=x+y+z,xy+xz+yz,xyz;
           ideal I=xy;
           print(im_of_var(I,F));
}