// $Id: finvar.lib,v 1.12 1998-05-29 15:01:57 Singular Exp $
// author: Agnes Eileen Heydtmann, email:agnes@math.uni-sb.de
// last change: 3.5.98
//////////////////////////////////////////////////////////////////////////////
version="$Id: finvar.lib,v 1.12 1998-05-29 15:01:57 Singular Exp $"
info="
LIBRARY:  finvar.lib       LIBRARY TO CALCULATE INVARIANT RINGS & MORE
                           (c) Agnes Eileen Heydtmann,
                           send bugs and comments to agnes@math.uni-sb.de

  cyclotomic(...)                        cyclotomic polynomial
  group_reynolds(...)                    finite group and Reynolds operator
  molien(...)                            Molien series
  reynolds_molien(...)                   Reynolds operator and Molien series
  partial_molien(...)                    partial expansion of Molien series
  evaluate_reynolds(...)                 image under the Reynolds operator
  invariant_basis(...)                   basis of homogeneous invariants
  invariant_basis_reynolds(...)          basis of homogeneous invariants
  primary_char0(...)                     primary invariants
  primary_charp(...)                     primary invariants
  primary_char0_no_molien(...)           primary invariants
  primary_charp_no_molien(...)           primary invariants
  primary_charp_without(...)             primary invariants
  primary_invariants(...)                primary invariants
  primary_char0_random(...)              primary invariants
  primary_charp_random(...)              primary invariants
  primary_char0_no_molien_random(...)    primary invariants
  primary_charp_no_molien_random(...)    primary invariants
  primary_charp_without_random(...)      primary invariants
  primary_invariants_random(...)         primary invariants
  power_products(...)                    exponents for power products
  secondary_char0(...)                   secondary invariants
  secondary_charp(...)                   secondary invariants
  secondary_no_molien(...)               secondary invariants
  secondary_with_irreducible_ones_no_molien(...)
                                         secondary invariants
  secondary_not_cohen_macaulay(...)      secondary invariants
  invariant_ring(...)                    primary and secondary invariants
  invariant_ring_random(...)             primary and secondary invariants
  algebra_containment(...)               answers query of algebra containment
  module_containment(...)                answers query of module containment
  orbit_variety(...)                     ideal of the orbit variety
  relative_orbit_variety(...)            ideal of a relative orbit variety
  image_of_variety(...)                  ideal of the image of a variety
";

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

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

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

proc group_reynolds (list #)
"USAGE:   group_reynolds(G1,G2,...[,v]);
         G1,G2,...: nxn <matrices> generating a finite matrix group, v: an
         optional <int>
ASSUME:  n is the number of variables of the basering, g the number of group
         elements
RETURN:  a <list>, the first list element will be a gxn <matrix> representing
         the Reynolds operator if we are in the non-modular case; if the
         characteristic is >0, minpoly==0 and the finite group non-cyclic the
         second list element is an <int> giving the lowest common multiple of
         the matrix group elements (used in molien); in general all other list
         elements are nxn <matrices> listing all elements of the finite group
DISPLAY: information if v does not equal 0
EXAMPLE: example group_reynolds; shows an example
THEORY:  The entire matrix group is generated by getting all left products of
         the generators with the new elements from the last run through the loop
         (or the generators themselves during the first run). All the ones that
         have been generated before are thrown out and the program terminates
         when there are no new elements found in one run. Additionally each time
         a new group element is found the corresponding ring mapping of which
         the Reynolds operator is made up is generated. They are stored in the
         rows of the first return value.
"
{ int ch=char(basering);               // the existance of the Reynolds operator
                                       // is dependent on the characteristic of
                                       // the base field
  int gen_num;                         // number of generators
 //------------------------ making sure the input is okay ---------------------
  if (typeof(#[size(#)])=="int")
  { if (size(#)==1)
    { "ERROR:   there are no matrices given among the parameters";
      return();
    }
    int v=#[size(#)];
    gen_num=size(#)-1;
  }
  else                                 // last parameter is not <int>
  { int v=0;                           // no information is default
    gen_num=size(#);
  }
  if (typeof(#[1])<>"matrix")
  { "ERROR:   the parameters must be a list of matrices and maybe an <int>";
    return();
  }
  int n=nrows(#[1]);
  if (n<>nvars(basering))
  { "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 REY=#[1]*vars;                // calculating the first ring mapping -
                                       // REY will contain the Reynolds
                                       // operator -
  matrix G(1)=#[1];                    // G(k) are elements of the group -
  if (ch<>0 && minpoly==0 && gen_num<>1) // finding out of which order the group
  { matrix I=diag(1,n);                // element is
    matrix TEST=G(1);
    int o1=1;
    int o2;
    while (TEST<>I)
    { TEST=TEST*G(1);
      o1=o1+1;
    }
  }
  int i=1;
 // -------------- doubles among the generators should be avoided -------------
  for (int j=2;j<=gen_num;j=j+1)       // this loop adds the parameters to the
  {                                    // group, leaving out doubles and
                                       // checking whether the parameters are
                                       // compatible with the task of the
                                       // procedure
    if (not(typeof(#[j])=="matrix"))
    { "ERROR:   the parameters must be a list of matrices and maybe 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];
      if (ch<>0 && minpoly==0)         // finding out of which order the group
      { TEST=G(i);                     // element is
        o2=1;
        while (TEST<>I)
        { TEST=TEST*G(i);
          o2=o2+1;
        }
        o1=o1*o2/gcd(o1,o2);           // lowest common multiple of the element
      }                                // orders -
      REY=concat(REY,#[j]*vars);       // adding ring homomorphisms to REY
    }
  }
  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)
  { "";
    "  Generating the entire matrix group and the Reynolds operator...";
    "";
  }
 // -------------- main loop that finds all the group elements ----------------
  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 -
        if (ch<>0 && minpoly==0 && i<>1) // finding out of which order the group
        { TEST=G(g);                   // element is
          o2=1;
          while (TEST<>I)
          { TEST=TEST*G(g);
            o2=o2+1;
          }
          o1=o1*o2/gcd(o1,o2);         // lowest common multiple of the element
        }                              // orders -
        REY=concat(REY,P(k)*vars);     // adding new mapping to REY
        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.";
    }
    "";
  }
  REY=transpose(REY);                  // when we evaluate the Reynolds operator
                                       // later on, we actually want 1xn
                                       // matrices
  if (ch<>0)
  { if ((g%ch)==0)
    { if (voice==2)
      { "WARNING: The characteristic of the coefficient field divides the group order.";
        "         Proceed without the Reynolds operator!";
      }
      else
      { if (v)
        { "  The characteristic of the base field divides the group order.";
          "  We have to continue without Reynolds operator...";
          "";
        }
      }
      kill REY;
      matrix REY[1][1]=0;
      return(REY,G(1..g));
    }
    if (minpoly==0)
    { if (i>1)
      { return(REY,o1,G(1..g));
      }
      return(REY,G(1..g));
    }
  }
  if (v)
  { "  Done generating the group and the Reynolds operator.";
    "";
  }
  return(REY,G(1..g));
}
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;
         list L=group_reynolds(A);
         print(L[1]);
         print(L[2..size(L)]);
}

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

proc molien (list #)
"USAGE:   molien(G1,G2,...[,ringname,lcm,flags]);
         G1,G2,...: nxn <matrices> generating a finite matrix group, ringname:
         a <string> giving a name for a new ring of characteristic 0 for the
         Molien series in case of prime characteristic, lcm: an <int> giving the
         lowest common multiple of the elements' orders in case of prime
         characteristic, minpoly==0 and a non-cyclic group, flags: an optional
         <intvec> with three components: if the first element is not equal to 0
         characteristic 0 is simulated, i.e. the Molien series is computed as
         if the base field were characteristic 0 (the user must choose a field
         of large prime characteristic, e.g. 32003), the second component should
         give the size of intervals between canceling common factors in the
         expansion of the Molien series, 0 (the default) means only once after
         generating all terms, in prime characteristic also a negative number
         can be given to indicate that common factors should always be canceled
         when the expansion is simple (the root of the extension field does not
         occur among the coefficients)
ASSUME:  n is the number of variables of the basering, G1,G2... are the group
         elements generated by group_reynolds(), lcm is the second return value
         of group_reynolds()
RETURN:  in case of characteristic 0 a 1x2 <matrix> giving enumerator and
         denominator of Molien series; in case of prime characteristic a ring
         with the name `ringname` of characteristic 0 is created where the same
         Molien series (named M) is stored
DISPLAY: information if the third component of flags does not equal 0
EXAMPLE: example molien; shows an example
THEORY:  In characteristic 0 the terms 1/det(1-xE) for all group elements of the
         Molien series are computed in a straight forward way. In prime
         characteristic a Brauer lift is involved. The returned matrix gives
         enumerator and denominator of the expanded version where common factors
         have been canceled.
"
{ def br=basering;                     // the Molien series depends on the
  int ch=char(br);                     // characteristic of the coefficient
                                       // field -
  int g;                               // size of the group
 //---------------------- making sure the input is okay -----------------------
  if (typeof(#[size(#)])=="intvec")
  { if (size(#[size(#)])==3)
    { int mol_flag=#[size(#)][1];
      if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0)))
      { "ERROR:   the second component of <intvec> should be >=0"
        return();
      }
      int interval=#[size(#)][2];
      int v=#[size(#)][3];
    }
    else
    { "ERROR:   <intvec> should have three components";
      return();
    }
    if (ch<>0)
    { if (typeof(#[size(#)-1])=="int")
      { int r=#[size(#)-1];
        if (typeof(#[size(#)-2])<>"string")
        { "ERROR:   in characteristic p>0 a <string> must be given for the name of a new";
          "         ring where the Molien series can be stored";
          return();
        }
        else
        { if (#[size(#)-2]=="")
          { "ERROR:   <string> may not be empty";
            return();
          }
          string newring=#[size(#)-2];
          g=size(#)-3;
        }
      }
      else
      { if (typeof(#[size(#)-1])<>"string")
        { "ERROR:   in characteristic p>0 a <string> must be given for the name of a new";
          "         ring where the Molien series can be stored";
          return();
        }
        else
        { if (#[size(#)-1]=="")
          { "ERROR:   <string> may not be empty";
            return();
          }
          string newring=#[size(#)-1];
          g=size(#)-2;
          int r=g;
        }
      }
    }
    else                               // then <string> ist not needed
    { g=size(#)-1;
    }
  }
  else                                 // last parameter is not <intvec>
  { int v=0;                           // no information is default
    int mol_flag=0;                    // computing of Molien series is default
    int interval=0;
    if (ch<>0)
    { if (typeof(#[size(#)])=="int")
      { int r=#[size(#)];
        if (typeof(#[size(#)-1])<>"string")
        { "ERROR:   in characteristic p>0 a <string> must be given for the name of a new";
          "         ring where the Molien series can be stored";
            return();
        }
        else
        { if (#[size(#)-1]=="")
          { "ERROR:   <string> may not be empty";
            return();
          }
          string newring=#[size(#)-1];
          g=size(#)-2;
        }
      }
      else
      { if (typeof(#[size(#)])<>"string")
        { "ERROR:   in characteristic p>0 a <string> must be given for the name of a new";
          "         ring where the Molien series can be stored";
          return();
        }
        else
        { if (#[size(#)]=="")
          { "ERROR:   <string> may not be empty";
            return();
          }
          string newring=#[size(#)];
          g=size(#)-1;
          int r=g;
        }
      }
    }
    else
    { g=size(#);
    }
  }
  if (ch<>0)
  { if ((g/r)*r<>g)
   { "ERROR:   <int> should divide the group order."
      return();
    }
  }
  if (ch<>0)
  { if ((g%ch)==0)
    { if (voice==2)
      { "WARNING: The characteristic of the coefficient field divides the group";
        "         order. Proceed without the Molien series!";
      }
      else
      { if (v)
        { "  The characteristic of the base field divides the group order.";
          "  We have to continue without Molien series...";
          "";
        }
      }
    }
    if (minpoly<>0 && mol_flag==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.";
          "";
        }
      }
      return();
    }
  }
 //----------------------------------------------------------------------------
  if (not(typeof(#[1])=="matrix"))
  { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
    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 square matrices";
    return();
  }
  if (v && voice<>2)
  { "";
    "  Generating the Molien series...";
    "";
  }
  if (v && voice==2)
  { "";
  }
 //------------- calculating Molien series in characteristic 0 ----------------
  if (ch==0)                           // when ch==0 we can calculate the Molien
  { matrix I=diag(1,n);                // series in any case -
    poly v1=maxideal(1)[1];            // the Molien series will be in terms of
                                       // the first variable of the current
                                       // ring -
    matrix M[1][2];                    // M will contain the Molien series -
    M[1,1]=0;                          // M[1,1] will be the numerator -
    M[1,2]=1;                          // M[1,2] will be the denominator -
    matrix s;                          // will help us canceling in the
                                       // fraction
    poly p;                            // will contain the denominator of the
                                       // new term of the Molien series
 //------------ computing 1/det(1+xE) for all E in the group ------------------
    for (int j=1;j<=g;j=j+1)
    { if (not(typeof(#[j])=="matrix"))
      { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
        return();
      }
      if ((n<>nrows(#[j])) or (n<>ncols(#[j])))
      { "ERROR:   matrices need to be square and of the same dimensions";
         return();
      }
      p=det(I-v1*#[j]);                // denominator of new term -
      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;
      if (interval<>0)                 // canceling common terms of denominator
      { if ((j/interval)*interval==j or j==g) // and enumerator -
        { s=matrix(syz(ideal(M)));     // once gcd() is faster than syz() these
          M[1,1]=-s[2,1];              // three lines should be replaced by the
          M[1,2]=s[1,1];               // following three
          // p=gcd(M[1,1],M[1,2]);
          // M[1,1]=M[1,1]/p;
          // M[1,2]=M[1,2]/p;
        }
      }
      if (v)
      { "  Term "+string(j)+" of the Molien series has been computed.";
      }
    }
    if (interval==0)                   // canceling common terms of denominator
    {                                  // and enumerator -
      s=matrix(syz(ideal(M)));         // once gcd() is faster than syz() these
      M[1,1]=-s[2,1];                  // three lines should be replaced by the
      M[1,2]=s[1,1];                   // following three
      // p=gcd(M[1,1],M[1,2]);
      // M[1,1]=M[1,1]/p;
      // M[1,2]=M[1,2]/p;
    }
    map slead=br,ideal(0);
    s=slead(M);
    M[1,1]=1/s[1,1]*M[1,1];            // numerator and denominator have to have
    M[1,2]=1/s[1,2]*M[1,2];            // a constant term of 1
    if (v)
    { "";
      "  We are done calculating the Molien series.";
      "";
    }
    return(M);
  }
 //---- calculating Molien series in prime characteristic with Brauer lift ----
  if (ch<>0 && mol_flag==0)
  { if (g<>1)
    { matrix G(1..g)=#[1..g];
      if (interval<0)
      { string Mstring;
      }
 //------ preparing everything for Brauer lifts into characteristic 0 ---------
      ring Q=0,x,dp;                   // we want to extend our ring as well as
                                       // the ring of rational numbers Q to
                                       // contain r-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=cyclotomic(r);         // 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 r-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 r-th cyclotomic
      poly minp=cyclotomic(r);         // 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;        // r-th roots of unity already
        map f1=v1br,ideal(0);
        number e=-number((f1(minp)));  // e is a r-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 r-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 xI=diag(x,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;
 // -------------- finding all the terms of the Molien series -----------------
      for (i=1;i<=g;i=i+1)
      { setring br;
        if (not(typeof(#[i])=="matrix"))
        { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
          return();
        }
        if ((n<>nrows(#[i])) or (n<>ncols(#[i])))
        { "ERROR:   matrices need to be square and of the same dimensions";
           return();
        }
        setring bre;
        p=det(xI-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 r-th root
                                       // of unity -
          k=exponent(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
        }
//         -----------
//         k=0;
//         while(k<r)
//         { map f3=basering,ideal(e^k);
//           while (f3(p)==0)
//           { p=p/(x-e^k);
//             setring `newring`;
//             p=p*(1-x*(e^k));        // building the denominator of the new
//             setring bre;
//           }
//           kill f3;
//           if (p==1)
//           { break;
//           }
//           k=k+1;
//         }
        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;
        if (interval<0)
        { if (i<>g)
          { Mstring=string(M);
            for (j=1;j<=size(Mstring);j=j+1)
            { if (Mstring[j]=="e")
              { interval=0;
                break;
              }
            }
          }
          if (interval<>0)
          { s=matrix(syz(ideal(M)));   // once gcd() is faster than syz()
            M[1,1]=-s[2,1];            // these three lines should be
            M[1,2]=s[1,1];             // replaced by the following three
            // p=gcd(M[1,1],M[1,2]);
            // M[1,1]=M[1,1]/p;
            // M[1,2]=M[1,2]/p;
          }
          else
          { interval=-1;
          }
        }
        else
        { if (interval<>0)             // canceling common terms of denominator
          { if ((i/interval)*interval==i or i==g) // and enumerator
            { s=matrix(syz(ideal(M))); // once gcd() is faster than syz()
              M[1,1]=-s[2,1];          // these three lines should be
              M[1,2]=s[1,1];           // replaced by the following three
              // p=gcd(M[1,1],M[1,2]);
              // M[1,1]=M[1,1]/p;
              // M[1,2]=M[1,2]/p;
            }
          }
        }
        p=1;
        setring bre;
        if (v)
        { "  Term "+string(i)+" of the Molien series has been computed.";
        }
      }
      if (v)
      { "";
      }
      setring `newring`;
      if (interval==0)                 // canceling common terms of denominator
      {                                // and enumerator -
        s=matrix(syz(ideal(M)));       // once gcd() is faster than syz() these
        M[1,1]=-s[2,1];                // three lines should be replaced by the
        M[1,2]=s[1,1];                 // following three
        // p=gcd(M[1,1],M[1,2]);
        // M[1,1]=M[1,1]/p;
        // M[1,2]=M[1,2]/p;
      }
      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
    }
    export `newring`;                  // we keep the ring where we computed the
    export M;                          // Molien series in such that we can
    setring br;                        // keep it
    if (v)
    { "  We are done calculating the Molien series.";
      "";
    }
  }
  else                                 // i.e. char<>0 and mol_flag<>0, the user
  {                                    // has specified that we are dealing with
                                       // a ring of large characteristic which
                                       // can be treated like a ring of
                                       // characteristic 0; we'll avoid the
                                       // Brauer lifts
 //----------------------- simulating characteristic 0 ------------------------
    string chst=charstr(br);
    for (int i=1;i<=size(chst);i=i+1)
    { if (chst[i]==",")
      { break;
      }
    }
 //----------------- generating ring of characteristic 0 ----------------------
    if (minpoly==0)
    { if (i>size(chst))
      { execute "ring "+newring+"=0,("+varstr(br)+"),("+ordstr(br)+")";
      }
      else
      { chst=chst[i..size(chst)];
        execute "ring "+newring+"=(0"+chst+"),("+varstr(br)+"),("+ordstr(br)+")";
      }
    }
    else
    { string minp=string(minpoly);
      minp=minp[2..size(minp)-1];
      chst=chst[i..size(chst)];
      execute "ring "+newring+"=(0"+chst+"),("+varstr(br)+"),("+ordstr(br)+")";
      execute "minpoly="+minp;
    }
    matrix I=diag(1,n);
    poly v1=maxideal(1)[1];            // the Molien series will be in terms of
                                       // the first variable of the current
                                       // ring -
    matrix M[1][2];                    // M will contain the Molien series -
    M[1,1]=0;                          // M[1,1] will be the numerator -
    M[1,2]=1;                          // M[1,2] will be the denominator -
    matrix s;                          // will help us canceling in the
                                       // fraction
    poly p;                            // will contain the denominator of the
                                       // new term of the Molien series
    int j;
    string links, rechts;
 //----------------- finding all terms of the Molien series -------------------
    for (i=1;i<=g;i=i+1)
    { setring br;
      if (not(typeof(#[i])=="matrix"))
      { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
        return();
      }
      if ((n<>nrows(#[i])) or (n<>ncols(#[i])))
      { "ERROR:   matrices need to be square and of the same dimensions";
         return();
      }
      string stM(i)=string(#[i]);
      for (j=1;j<=size(stM(i));j=j+1)
      { if (stM(i)[j]=="
")
        { links=stM(i)[1..j-1];
          rechts=stM(i)[j+1..size(stM(i))];
          stM(i)=links+rechts;
        }
      }
      setring `newring`;
      execute "matrix G(i)["+string(n)+"]["+string(n)+"]="+stM(i);
      p=det(I-v1*G(i));                // denominator of new term -
      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;
      if (interval<>0)                 // canceling common terms of denominator
      { if ((i/interval)*interval==i or i==g) // and enumerator
        {
          s=matrix(syz(ideal(M)));     // once gcd() is faster than syz() these
          M[1,1]=-s[2,1];              // three lines should be replaced by the
          M[1,2]=s[1,1];               // following three
          // p=gcd(M[1,1],M[1,2]);
          // M[1,1]=M[1,1]/p;
          // M[1,2]=M[1,2]/p;
        }
      }
      if (v)
      { "  Term "+string(i)+" of the Molien series has been computed.";
      }
    }
    if (interval==0)                   // canceling common terms of denominator
    {                                  // and enumerator -
      s=matrix(syz(ideal(M)));         // once gcd() is faster than syz() these
      M[1,1]=-s[2,1];                  // three lines should be replaced by the
      M[1,2]=s[1,1];                   // following three
      // p=gcd(M[1,1],M[1,2]);
      // M[1,1]=M[1,1]/p;
      // M[1,2]=M[1,2]/p;
    }
    map slead=`newring`,ideal(0);
    s=slead(M);
    M[1,1]=1/s[1,1]*M[1,1];            // numerator and denominator have to have
    M[1,2]=1/s[1,2]*M[1,2];            // a constant term of 1
    if (v)
    { "";
      "  We are done calculating the Molien series.";
      "";
    }
    kill G(1..g), s, slead, p, v1, I;
    export `newring`;                  // we keep the ring where we computed the
    export M;                          // the Molien series such that we can
    setring br;                        // keep it
  }
}
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;
         list L=group_reynolds(A);
         matrix M=molien(L[2..size(L)]);
         print(M);
         ring S=3,(x,y,z),dp;
         string newring="alksdfjlaskdjf";
         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
         list L=group_reynolds(A);
         molien(L[2..size(L)],newring);
         setring alksdfjlaskdjf;
         print(M);
         setring S;
         kill alksdfjlaskdjf;
}

proc reynolds_molien (list #)
"USAGE:   reynolds_molien(G1,G2,...[,ringname,flags]);
         G1,G2,...: nxn <matrices> generating a finite matrix group, ringname:
         a <string> giving a name for a new ring of characteristic 0 for the
         Molien series in case of prime characteristic, flags: an optional
         <intvec> with three components: if the first element is not equal to 0
         characteristic 0 is simulated, i.e. the Molien series is computed as
         if the base field were characteristic 0 (the user must choose a field
         of large prime characteristic, e.g. 32003) the second component should
         give the size of intervals between canceling common factors in the
         expansion of the Molien series, 0 (the default) means only once after
         generating all terms, in prime characteristic also a negative number
         can be given to indicate that common factors should always be canceled
         when the expansion is simple (the root of the extension field does not
         occur among the coefficients)
ASSUME:  n is the number of variables of the basering, G1,G2... are the group
         elements generated by group_reynolds(), g is the size of the group
RETURN:  a gxn <matrix> representing the Reynolds operator is the first return
         value and in case of characteristic 0 a 1x2 <matrix> giving enumerator
         and denominator of Molien series is the second one; in case of prime
         characteristic a ring with the name `ringname` of characteristic 0 is
         created where the same Molien series (named M) is stored
DISPLAY: information if the third component of flags does not equal 0
EXAMPLE: example reynolds_molien; shows an example
THEORY:  The entire matrix group is generated by getting all left products of
         the generators with the new elements from the last run through the loop
         (or the generators themselves during the first run). All the ones that
         have been generated before are thrown out and the program terminates
         when there are no new elements found in one run. Additionally each time
         a new group element is found the corresponding ring mapping of which
         the Reynolds operator is made up is generated. They are stored in the
         rows of the first return value. In characteristic 0 the terms
         1/det(1-xE) is computed whenever a new element E is found. In prime
         characteristic a Brauer lift is involved and the terms are only
         computed after the entire matrix group is generated (to avoid the
         modular case). The returned matrix gives enumerator and denominator of
         the expanded version where common factors have been canceled.
"
{ def br=basering;                     // the Molien series depends on the
  int ch=char(br);                     // characteristic of the coefficient
                                       // field
  int gen_num;
 //------------------- making sure the input is okay --------------------------
  if (typeof(#[size(#)])=="intvec")
  { if (size(#[size(#)])==3)
    { int mol_flag=#[size(#)][1];
      if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0)))
      { "ERROR:   the second component of the <intvec> should be >=0";
        return();
      }
      int interval=#[size(#)][2];
      int v=#[size(#)][3];
    }
    else
    { "ERROR:   <intvec> should have three components";
      return();
    }
    if (ch<>0)
    { if (typeof(#[size(#)-1])<>"string")
      { "ERROR:   in characteristic p a <string> must be given for the name";
        "         of a new ring where the Molien series can be stored";
        return();
      }
      else
      { if (#[size(#)-1]=="")
        { "ERROR:   <string> may not be empty";
          return();
        }
        string newring=#[size(#)-1];
        gen_num=size(#)-2;
      }
    }
    else                               // then <string> ist not needed
    { gen_num=size(#)-1;
    }
  }
  else                                 // last parameter is not <intvec>
  { int v=0;                           // no information is default
    int interval;
    int mol_flag=0;                    // computing of Molien series is default
    if (ch<>0)
    { if (typeof(#[size(#)])<>"string")
      { "ERROR:   in characteristic p a <string> must be given for the name";
        "         of a new ring where the Molien series can be stored";
        return();
      }
      else
      { if (#[size(#)]=="")
        { "ERROR:   <string> may not be empty";
          return();
        }
        string newring=#[size(#)];
        gen_num=size(#)-1;
      }
    }
    else
    { gen_num=size(#);
    }
  }
 // ----------------- computing the terms with Brauer lift --------------------
  if (ch<>0 && mol_flag==0)
  { list L=group_reynolds(#[1..gen_num],v);
    if (L[1]==0)
    { if (voice==2)
      { "WARNING: The characteristic of the coefficient field divides the group order.";
        "         Proceed without the Reynolds operator or the Molien series!";
        return();
      }
      if (v)
      { "  The characteristic of the base field divides the group order.";
        "  We have to continue without Reynolds operator or the Molien series...";
        return();
      }
    }
    if (minpoly<>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.";
        return(L[1]);
      }
      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.";
          return(L[1]);
        }
      }
    }
    if (typeof(L[2])=="int")
    { molien(L[3..size(L)],newring,L[2],intvec(mol_flag,interval,v));
    }
    else
    { molien(L[2..size(L)],newring,intvec(mol_flag,interval,v));
    }
    return(L[1]);
  }
 //----------- computing Molien series in the straight forward way ------------
  if (ch==0)
  { if (typeof(#[1])<>"matrix")
    { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
      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 -
    poly v1=vars[1,1];                 // the Molien series will be in terms of
                                       // the first variable of the current
                                       // ring
    matrix I=diag(1,n);
    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
    matrix G(1)=#[1];                  // G(k) are elements of the group -
    A(2)[1,2]=det(I-v1*(G(1)));        // A(2)[1,2] will be the denominator -
    matrix s;                          // will help us canceling in the
                                       // fraction
    poly p;                            // will contain the denominator of the
                                       // new term of the Molien series
    int i=1;
    for (int j=2;j<=gen_num;j=j+1)     // this loop adds the parameters to the
    {                                  // group, leaving out doubles and
                                       // checking whether the parameters are
                                       // compatible with the task of the
                                       // procedure
      if (not(typeof(#[j])=="matrix"))
      { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
        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) -
        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;
        if (interval<>0)               // canceling common terms of denominator
        { if ((i/interval)*interval==i) // and enumerator
          {
            s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz() these
            A(2)[1,1]=-s[2,1];         // three lines should be replaced by the
            A(2)[1,2]=s[1,1];          // following three
            // p=gcd(A(2)[1,1],A(2)[1,2]);
            // A(2)[1,1]=A(2)[1,1]/p;
            // A(2)[1,2]=A(2)[1,2]/p;
          }
        }
      }
    }
    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)
    { "";
      "  Generating the entire matrix group. Whenever a new group element is found,";
      "  the coressponding ring homomorphism of the Reynolds operator and the";
      "  corresponding term of the Molien series is generated.";
      "";
    }
 //----------- computing 1/det(I-xE) whenever a new element E is found --------
    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)
          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 -
          if (interval<>0)             // canceling common terms of denominator
          { if ((g/interval)*interval==g) // and enumerator
            {
              s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz()
              A(2)[1,1]=-s[2,1];       // these three lines should be replaced
              A(2)[1,2]=s[1,1];        // by the following three
              // p=gcd(A(2)[1,1],A(2)[1,2]);
              // A(2)[1,1]=A(2)[1,1]/p;
              // A(2)[1,2]=A(2)[1,2]/p;
            }
          }
          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 (interval==0)                   // canceling common terms of denominator
    {                                  // and enumerator -
      s=matrix(syz(ideal(A(2))));      // once gcd() is faster than syz()
      A(2)[1,1]=-s[2,1];               // these three lines should be replaced
      A(2)[1,2]=s[1,1];                // by the following three
      // p=gcd(A(2)[1,1],A(2)[1,2]);
      // A(2)[1,1]=A(2)[1,1]/p;
      // A(2)[1,2]=A(2)[1,2]/p;
    }
    if (interval<>0)                   // canceling common terms of denominator
    { if ((g/interval)*interval<>g)    // and enumerator
      {
        s=matrix(syz(ideal(A(2))));    // once gcd() is faster than syz()
        A(2)[1,1]=-s[2,1];             // these three lines should be replaced
        A(2)[1,2]=s[1,1];              // by the following three
        // p=gcd(A(2)[1,1],A(2)[1,2]);
        // A(2)[1,1]=A(2)[1,1]/p;
        // A(2)[1,2]=A(2)[1,2]/p;
      }
    }
    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));
  }
 //------------------------ simulating characteristic 0 -----------------------
  else                                 // if ch<>0 and mol_flag<>0
  { if (typeof(#[1])<>"matrix")
    { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
      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
    string chst=charstr(br);
    for (int i=1;i<=size(chst);i=i+1)
    { if (chst[i]==",")
      { break;
      }
    }
    if (minpoly==0)
    { if (i>size(chst))
      { execute "ring "+newring+"=0,("+varstr(br)+"),("+ordstr(br)+")";
      }
      else
      { chst=chst[i..size(chst)];
        execute "ring "+newring+"=(0"+chst+"),("+varstr(br)+"),("+ordstr(br)+")";
      }
    }
    else
    { string minp=string(minpoly);
      minp=minp[2..size(minp)-1];
      chst=chst[i..size(chst)];
      execute "ring "+newring+"=(0"+chst+"),("+varstr(br)+"),("+ordstr(br)+")";
      execute "minpoly="+minp;
    }
    poly v1=var(1);                    // the Molien series will be in terms of
                                       // the first variable of the current
                                       // ring
    matrix I=diag(1,n);
    int o;
    setring br;
    matrix G(1)=#[1];
    string links, rechts;
    string stM(1)=string(#[1]);
    for (o=1;o<=size(stM(1));o=o+1)
    { if (stM(1)[o]=="
")
      { links=stM(1)[1..o-1];
        rechts=stM(1)[o+1..size(stM(1))];
        stM(1)=links+rechts;
      }
    }
    setring `newring`;
    execute "matrix G(1)["+string(n)+"]["+string(n)+"]="+stM(1);
    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*(G(1)));        // A(2)[1,2] will be the denominator -
    matrix s;                          // will help us canceling in the
                                       // fraction
    poly p;                            // will contain the denominator of the
                                       // new term of the Molien series
    i=1;
    for (int j=2;j<=gen_num;j=j+1)     // this loop adds the parameters to the
    {                                  // group, leaving out doubles and
                                       // checking whether the parameters are
                                       // compatible with the task of the
                                       // procedure
      setring br;
      if (not(typeof(#[j])=="matrix"))
      { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
        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),G(i)*vars);   // adding ring homomorphisms to A(1)
        string stM(i)=string(G(i));
        for (o=1;o<=size(stM(i));o=o+1)
        { if (stM(i)[o]=="
")
          { links=stM(i)[1..o-1];
            rechts=stM(i)[o+1..size(stM(i))];
            stM(i)=links+rechts;
          }
        }
        setring `newring`;
        execute "matrix G(i)["+string(n)+"]["+string(n)+"]="+stM(i);
        p=det(I-v1*G(i));              // 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;
        if (interval<>0)               // canceling common terms of denominator
        { if ((i/interval)*interval==i) // and enumerator
          {
            s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz() these
            A(2)[1,1]=-s[2,1];         // three lines should be replaced by the
            A(2)[1,2]=s[1,1];          // following three
            // p=gcd(A(2)[1,1],A(2)[1,2]);
            // A(2)[1,1]=A(2)[1,1]/p;
            // A(2)[1,2]=A(2)[1,2]/p;
          }
        }
        setring br;
      }
    }
    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)
    { "";
      "  Generating the entire matrix group. Whenever a new group element is found,";
      "  the coressponding ring homomorphism of the Reynolds operator and the";
      "  corresponding term of the Molien series is generated.";
      "";
    }
 // taking all elements in a ring of characteristic 0 and computing the terms
 // of the Molien series there
    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),G(g)*vars); // adding new mapping to A(1)
          string stM(g)=string(G(g));
          for (o=1;o<=size(stM(g));o=o+1)
          { if (stM(g)[o]=="
")
            { links=stM(g)[1..o-1];
              rechts=stM(g)[o+1..size(stM(g))];
              stM(g)=links+rechts;
            }
          }
          setring `newring`;
          execute "matrix G(g)["+string(n)+"]["+string(n)+"]="+stM(g);
          p=det(I-v1*G(g));            // 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 -
          if (interval<>0)             // canceling common terms of denominator
          { if ((g/interval)*interval==g) // and enumerator
            {
              s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz()
              A(2)[1,1]=-s[2,1];       // these three lines should be replaced
              A(2)[1,2]=s[1,1];        // by the following three
              // p=gcd(A(2)[1,1],A(2)[1,2]);
              // A(2)[1,1]=A(2)[1,1]/p;
              // A(2)[1,2]=A(2)[1,2]/p;
            }
          }
          if (v)
          { "  Group element "+string(g)+" has been found.";
          }
          setring br;
        }
        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
    setring `newring`;
    if (interval==0)                   // canceling common terms of denominator
    {                                  // and enumerator -
      s=matrix(syz(ideal(A(2))));      // once gcd() is faster than syz()
      A(2)[1,1]=-s[2,1];               // these three lines should be replaced
      A(2)[1,2]=s[1,1];                // by the following three
      // p=gcd(A(2)[1,1],A(2)[1,2]);
      // A(2)[1,1]=A(2)[1,1]/p;
      // A(2)[1,2]=A(2)[1,2]/p;
    }
    if (interval<>0)                   // canceling common terms of denominator
    { if ((g/interval)*interval<>g)    // and enumerator
      {
        s=matrix(syz(ideal(A(2))));    // once gcd() is faster than syz()
        A(2)[1,1]=-s[2,1];             // these three lines should be replaced
        A(2)[1,2]=s[1,1];              // by the following three
        // p=gcd(A(2)[1,1],A(2)[1,2]);
        // A(2)[1,1]=A(2)[1,1]/p;
        // A(2)[1,2]=A(2)[1,2]/p;
      }
    }
    map slead=`newring`,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.";
      "";
    }
    matrix M=A(2);
    kill G(1..g), s, slead, p, v1, I, A(2);
    export `newring`;                  // we keep the ring where we computed the
    export M;                          // the Molien series such that we can
    setring br;                        // keep it
    return(A(1));
  }
}
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 REY,M=reynolds_molien(A);
         print(REY);
         print(M);
         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=reynolds_molien(A,newring);
         print(REY);
         setring Qadjoint;
         print(M);
         setring S;
         kill Qadjoint;
}

proc partial_molien (matrix M, int n, list #)
"USAGE:   partial_molien(M,n[,p]);
         M: a 1x2 <matrix>, n: an <int> indicating  number of terms in the
         expansion, p: an optional <poly>
ASSUME:  M is the return value of molien or the second return value of
         reynolds_molien, p ought to be the second return value of a previous
         run of partial_molien and avoids recalculating known terms
RETURN:  n terms (type <poly>) of the partial expansion of the Molien series
         (first n if there is no third parameter 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 parameter in a next
         run of partial_molien
THEORY:  The following calculation is implemented:
         (1+a1x+a2x^2+...+anx^n)/(1+b1x+b2x^2+...+bmx^m)=(1+(a1-b1)x+...
         (1+b1x+b2x^2+...+bmx^m)
         -----------------------
            (a1-b1)x+(a2-b2)x^2+...
            (a1-b1)x+b1(a1-b1)x^2+...
EXAMPLE: example partial_molien; 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 parameter 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 parameter 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 parameter has to be a 1x2-matrix, i.e. the matrix";
    "         returned by the procedure 'reynolds_molien', the second one";
    "         should be > 0 and there should be no more than 3 parameters;"
    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 REY,M=reynolds_molien(A);
         poly p(1..2);
         p(1..2)=partial_molien(M,5);
         p(1);
         p(1..2)=partial_molien(M,5,p(2));
         p(1);
}

proc evaluate_reynolds (matrix REY, ideal I)
"USAGE:   evaluate_reynolds(REY,I);
         REY: a <matrix> representing the Reynolds operator, I: an arbitrary
         <ideal>
ASSUME:  REY is the first return value of group_reynolds() or reynolds_molien()
RETURNS: image of the polynomials defining I under the Reynolds operator
         (type <ideal>)
NOTE:    the characteristic of the coefficient field of the polynomial ring
         should not divide the order of the finite matrix group
EXAMPLE: example evaluate_reynolds; shows an example
THEORY:  REY has been constructed in such a way that each row serves as a ring
         mapping of which the Reynolds operator is made up.
"
{ def br=basering;
  int n=nvars(br);
  if (ncols(REY)==n)
  { int m=nrows(REY);                  // we need m to 'cut' the ring
                                       // homomorphisms 'out' of REY and to
                                       // divide by the group order in the end
    int num_poly=ncols(I);
    matrix MI=matrix(I);
    matrix MiI[1][num_poly];
    map pREY;
    matrix rowREY[1][n];
    for (int i=1;i<=m;i=i+1)
    { rowREY=REY[i,1..n];
      pREY=br,ideal(rowREY);           // f is now the i-th ring homomorphism
      MiI=pREY(MI)+MiI;
    }
    MiI=(1/number(m))*MiI;
    return(ideal(MiI));
  }
  else
  { "ERROR:   the number of columns in the <matrix> should be the same as the";
    "         number of variables in the basering; in fact it should be first";
    "         return value of group_reynolds() or reynolds_molien().";
    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;
         list L=group_reynolds(A);
         ideal I=x2,y2,z2;
         print(evaluate_reynolds(L[1],I));
}

proc invariant_basis (int g,list #)
"USAGE:   invariant_basis(g,G1,G2,...);
         g: an <int> indicating of which degree (>0) the homogeneous basis
         shoud be, G1,G2,...: <matrices> generating a finite matrix group
RETURNS: the basis (type <ideal>) of the space of invariants of degree g
EXAMPLE: example invariant_basis; shows an example
THEORY:  A general polynomial of degree g is generated and the generators of the
         matrix group applied. The difference ought to be 0 and this way a
         system of linear equations is created. It is solved by computing
         syzygies.
"
{ if (g<=0)
  { "ERROR:   the first parameter 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 g
  mon=sort(mon,intvec(m..1))[1];
  int a=size(#);
  int i;
  int n=nvars(br);
 //---------------------- checking that the input is ok -----------------------
  for (i=1;i<=a;i=i+1)
  { 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 parameters should be a list of matrices";
      return();
    }
  }
 //----------------------------------------------------------------------------
  execute "ring T=("+charstr(br)+"),("+varstr(br)+",p(1..m)),lp;";
  // p(1..m) are the general coefficients of the general polynomial of degree g
  execute "ideal vars="+varstr(br)+";";
  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;
 //------------------- building the system of linear equations ----------------
  for (i=1;i<=a;i=i+1)
  { 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(invariant_basis(2,A));
}

proc invariant_basis_reynolds (matrix REY,int d,list #)
"USAGE:   invariant_basis_reynolds(REY,d[,flags]);
         REY: a <matrix> representing the Reynolds operator, d: an <int>
         indicating of which degree (>0) the homogeneous basis shoud be, flags:
         an optional <intvec> with two entries: its first component gives the
         dimension of the space (default <0 meaning unknown) and its second
         component is used as the number of polynomials that should be mapped
         to invariants during one call of evaluate_reynolds if the dimension of
         the space is unknown or the number such that number x dimension
         polynomials are mapped to invariants during one call of
         evaluate_reynolds
ASSUME:  REY is the first return value of group_reynolds() or reynolds_molien()
         and flags[1] given by partial_molien
RETURN:  the basis (type <ideal>) of the space of invariants of degree d
EXAMPLE: example invariant_basis_reynolds; shows an example
THEORY:  Monomials of degree d are mapped to invariants with the Reynolds
         operator. A linearly independent set is generated with the help of
         minbase.
"
{
 //---------------------- checking that the input is ok -----------------------
  if (d<=0)
  { "  ERROR:   the second parameter should be > 0";
     return();
  }
  if (size(#)>1)
  { "  ERROR:   there should be at most three parameters";
    return();
  }
  if (size(#)==1)
  { if (typeof(#[1])<>"intvec")
    { "  ERROR: the third parameter should be of type <intvec>";
      return();
    }
    if (size(#[1])<>2)
    { "  ERROR: there should be two components in <intvec>";
      return();
    }
    else
    { int cd=#[1][1];
      int step_fac=#[1][2];
    }
    if (step_fac<=0)
    { "  ERROR: the second component of <intvec> should be > 0";
      return();
    }
    if (cd==0)
    { return(ideal(0));
    }
  }
  else
  { int step_fac=1;
    int cd=-1;
  }
  if (ncols(REY)<>nvars(basering))
  { "ERROR:   the number of columns in the <matrix> should be the same as the";
    "         number of variables in the basering; in fact it should be first";
    "         return value of group_reynolds() or reynolds_molien().";
    return();
  }
 //----------------------------------------------------------------------------
  ideal mon=sort(maxideal(d))[1];
  degBound=d;
  int j=ncols(mon);
  mon=sort(mon,intvec(j..1))[1];
  ideal B;                             // will contain the basis
  if (cd<0)
  { if (step_fac>j)                    // all of mon will be mapped to
    { B=evaluate_reynolds(REY,mon);    // invariants at once
      B=minbase(B);
      degBound=0;
      return(B);
    }
  }
  else
  { if (step_fac*cd>j)                 // all of mon will be mapped to
    { B=evaluate_reynolds(REY,mon);    // invariants at once
      B=minbase(B);
      degBound=0;
      return(B);
    }
  }
  int i,k;
  int upper_bound=0;
  int lower_bound=0;
  ideal part_mon;                      // a part of mon of size step_fac*cd
  while (1)
  { lower_bound=upper_bound+1;
    if (cd<0)
    { upper_bound=upper_bound+step_fac;
    }
    else
    { upper_bound=upper_bound+step_fac*cd;
    }
    if (upper_bound>j)
    { upper_bound=j;
    }
    part_mon=mon[lower_bound..upper_bound];
    B=minbase(B+evaluate_reynolds(REY,part_mon));
    if ((ncols(B)==cd and B[1]<>0) or upper_bound==j)
    { 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;
           intvec flags=0,1,0;
           matrix REY,M=reynolds_molien(A,flags);
           flags=8,6;
           print(invariant_basis_reynolds(REY,6,flags));
}

////////////////////////////////////////////////////////////////////////////////
// This procedure generates linearly independent invariant polynomials of degree
// d that do not reduce to 0 modulo the primary invariants. It does this by
// applying the Reynolds operator to the monomials returned by kbase(sP,d). The
// result is used when computing secondary invariants.
////////////////////////////////////////////////////////////////////////////////
proc sort_of_invariant_basis (ideal sP,matrix REY,int d,int step_fac)
{ ideal mon=kbase(sP,d);
  degBound=d;
  int j=ncols(mon);
  int i;
  mon=sort(mon,intvec(j..1))[1];
  ideal B;                             // will contain the "sort of basis"
  if (step_fac>j)
  { B=compress(evaluate_reynolds(REY,mon));
    for (i=1;i<=ncols(B);i=i+1)        // those are taken our that are o mod sP
    { if (reduce(B[i],sP)==0)
      { B[i]=0;
      }
    }
    B=minbase(B);                      // here are the linearly independent ones
    degBound=0;
    return(B);
  }
  int upper_bound=0;
  int lower_bound=0;
  ideal part_mon;                      // parts of mon
  while (1)
  { lower_bound=upper_bound+1;
    upper_bound=upper_bound+step_fac;
    if (upper_bound>j)
    { upper_bound=j;
    }
    part_mon=mon[lower_bound..upper_bound];
    part_mon=compress(evaluate_reynolds(REY,part_mon));
    for (i=1;i<=ncols(part_mon);i=i+1)
    { if (reduce(part_mon[i],sP)==0)
      { part_mon[i]=0;
      }
    }
    B=minbase(B+part_mon);             // here are the linearly independent ones
    if (upper_bound==j)
    { degBound=0;
      return(B);
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
// 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 next_vector(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=next_vector(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);
}

////////////////////////////////////////////////////////////////////////////////
// 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 int_number_map (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;
  }
}

////////////////////////////////////////////////////////////////////////////////
// This procedure finds dif primary invariants in degree d. It returns all
// primary invariants found so far. The coefficients lie in a field of
// characteristic 0.
////////////////////////////////////////////////////////////////////////////////
proc search (int n,int d,ideal B,int cd,ideal P,ideal sP,int i,int dif,int dB,ideal CI)
{ intmat vec[1][cd];                   // the coefficients for the next
                                       // combination -
  degBound=0;
  poly test_poly;                      // the linear combination to test
  int test_dim;
  intvec h;                            // Hilbert series
  int j=i+1;
  matrix tB=transpose(B);
  ideal TEST;
  while(j<=i+dif)
  { CI=CI+ideal(var(j)^d);             // homogeneous polynomial of the same
                                       // degree as the one we're looking for is
                                       // added
    // h=hilb(std(CI),1);
    dB=dB+d-1;                         // used as degBound
    while(1)
    { vec=next_vector(vec);            // next vector
      test_poly=(vec*tB)[1,1];
      // degBound=dB;
      TEST=sP+ideal(test_poly);
      attrib(TEST,"isSB",1);
      test_dim=dim(TEST);
      // degBound=0;
      if (n-test_dim==j)               // the dimension has been lowered by one
      { sP=TEST;
        break;
      }
      // degBound=dB;
      TEST=std(sP+ideal(test_poly));   // should soon be replaced by next line
      // TEST=std(sP,test_poly,h);        // Hilbert driven std-calculation
      test_dim=dim(TEST);
      // degBound=0;
      if (n-test_dim==j)               // the dimension has been lowered by one
      { sP=TEST;
        break;
      }
    }
    P[j]=test_poly;                    // test_poly ist added to primary
    j=j+1;                             // invariants
  }
  return(P,sP,CI,dB);
}

////////////////////////////////////////////////////////////////////////////////
// This procedure finds at most dif primary invariants in degree d. It returns
// all primary invariants found so far. The coefficients lie in the field of
// characteristic p>0.
////////////////////////////////////////////////////////////////////////////////
proc p_search (int n,int d,ideal B,int cd,ideal P,ideal sP,int i,int dif,int dB,ideal CI)
{ def br=basering;
  degBound=0;
  matrix vec(1)[1][cd];                // starting with 0-vector -
  intmat new[1][cd];                   // the coefficients for the next
                                       // combination -
  matrix pnew[1][cd];                  // new needs to be mapped into br -
  int counter=1;                       // counts the vectors
  int j;
  int p=char(br);
  if (minpoly<>0)
  { int ext_deg=pardeg(minpoly);       // field has p^d elements
  }
  else
  { int ext_deg=1;                     // field has p^d elements
  }
  poly test_poly;                      // the linear combination to test
  int test_dim;
  ring R=0,x,dp;                       // just to calculate next variable
                                       // bound -
  number bound=(number(p)^(ext_deg*cd)-1)/(number(p)^ext_deg-1)+1; // this is
                                       // how many linearly independent vectors
                                       // of size cd exist having entries in the
                                       // base field of br
  setring br;
  intvec h;                            // Hilbert series
  int k=i+1;
  matrix tB=transpose(B);
  ideal TEST;
  while (k<=i+dif)
  { CI=CI+ideal(var(k)^d);             // homogeneous polynomial of the same
                                       // degree as the one we're looking for is
                                       // added
    // h=hilb(std(CI),1);
    dB=dB+d-1;                         // used as degBound
    setring R;
    while (number(counter)<>bound)     // otherwise, we are done
    { setring br;
      new=next_vector(new);
      for (j=1;j<=cd;j=j+1)
      { pnew[1,j]=int_number_map(new[1,j]); // 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 -
        test_poly=(vec(counter)*tB)[1,1]; // linear combination -
        // degBound=dB;
        TEST=sP+ideal(test_poly);
        attrib(TEST,"isSB",1);
        test_dim=dim(TEST);
        // degBound=0;
        if (n-test_dim==k)             // the dimension has been lowered by one
        { sP=TEST;
          setring R;
          break;
        }
        // degBound=dB;
        TEST=std(sP+ideal(test_poly)); // should soon to be replaced by next
                                       // line
        // TEST=std(sP,test_poly,h);      // Hilbert driven std-calculation
        test_dim=dim(TEST);
        // degBound=0;
        if (n-test_dim==k)             // the dimension has been lowered by one
        { sP=TEST;
          setring R;
          break;
        }
      }
      setring R;
    }
    if (number(counter)<=bound)
    { setring br;
      P[k]=test_poly;                  // test_poly ist added to primary
    }                                  // invariants
    else
    { setring br;
      CI=CI[1..size(CI)-1];
      return(P,sP,CI,dB-d+1);
    }
    k=k+1;
  }
  return(P,sP,CI,dB);
}

proc primary_char0 (matrix REY,matrix M,list #)
"USAGE:   primary_char0(REY,M[,v]);
         REY: a <matrix> representing the Reynolds operator, M: a 1x2 <matrix>
         representing the Molien series, v: an optional <int>
ASSUME:  REY is the first return value of group_reynolds or reynolds_molien and
         M the one of molien or the second one of reynolds_molien
DISPLAY: information about the various stages of the programme if v does not
         equal 0
RETURN:  primary invariants (type <matrix>) of the invariant ring
EXAMPLE: example primary_char0; shows an example
THEORY:  Bases of homogeneous invariants are generated successively and those
         are chosen as primary invariants that lower the dimension of the ideal
         generated by the previously found invariants (see paper \"Generating a
         Noetherian Normalization of the Invariant Ring of a Finite Group\" by
         Decker, Heydtmann, Schreyer (1997) to appear in JSC).
"
{ degBound=0;
  if (char(basering)<>0)
  { "ERROR:   primary_char0 should only be used with rings of characteristic 0.";
    return();
  }
 //----------------- checking input and setting verbose mode ------------------
  if (size(#)>1)
  { "ERROR:   primary_char0 can only have three parameters.";
    return();
  }
  if (size(#)==1)
  { if (typeof(#[1])<>"int")
    { "ERROR:   The third parameter should be of type <int>.";
      return();
    }
    else
    { int v=#[1];
    }
  }
  else
  { int v=0;
  }
  int n=nvars(basering);               // 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 (ncols(REY)<>n)
  { "ERROR:   First parameter ought to be the Reynolds operator."
    return();
  }
  if (ncols(M)<>2 or nrows(M)<>1)
  { "ERROR:   Second parameter ought to be the Molien series."
    return();
  }
 //----------------------------------------------------------------------------
  if (v && voice<>2)
  { "  We can start looking for primary invariants...";
    "";
  }
  if (v && voice==2)
  { "";
  }
 //------------------------- initializing variables ---------------------------
  int dB;
  poly p(1..2);                        // p(1) will be used for single terms of
                                       // the partial expansion, p(2) to store
  p(1..2)=partial_molien(M,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) -
  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
                                       // space of invariants of degree d,
                                       // newdim: dimension the ideal generated
                                       // the primary invariants plus basis
                                       // elements, dif=n-i-newdim, i.e. the
                                       // number of new primary invairants that
                                       // should be added in this degree -
  ideal P,Pplus,sPplus,CI,B;           // P: will contain primary invariants,
                                       // Pplus: P+B, CI: a complete
                                       // intersection with the same Hilbert
                                       // function as P
  ideal sP=std(P);
  dB=1;                                // used as degree bound
  int i=0;
 //-------------- loop that searches for primary invariants  ------------------
  while(1)                             // repeat until n primary invariants are
  {                                    // found -
    p(1..2)=partial_molien(M,1,p(2));  // next term of the partial expansion -
    d=deg(p(1));                       // degree where we'll search -
    cd=int(coef(p(1),v1)[2,1]);        // dimension of the homogeneous space of
                                       // inviarants of degree d
    if (v)
    { "  Computing primary invariants in degree "+string(d)+":";
    }
    B=invariant_basis_reynolds(REY,d,intvec(cd,6)); // basis of invariants of
                                       // degree d
    if (B[1]<>0)
    { Pplus=P+B;
      sPplus=std(Pplus);
      newdim=dim(sPplus);
      dif=n-i-newdim;
    }
    else
    { dif=0;
    }
    if (dif<>0)                        // we have to find dif new primary
    {                                  // invariants
      if (cd<>dif)
      { P,sP,CI,dB=search(n,d,B,cd,P,sP,i,dif,dB,CI); // searching for dif invariants
      }                                // i.e. we can take all of B
      else
      { for(j=i+1;j>i+dif;j=j+1)
        { CI=CI+ideal(var(j)^d);
        }
        dB=dB+dif*(d-1);
        P=Pplus;
        sP=sPplus;
      }
      if (v)
      { for (j=1;j<=dif;j=j+1)
        { "  We find: "+string(P[i+j]);
        }
      }
      i=i+dif;
      if (i==n)                        // found all primary invariants
      { if (v)
        { "";
          "  We found all primary invariants.";
          "";
        }
        return(matrix(P));
      }
    }                                  // done with degree d
  }
}
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 REY,M=reynolds_molien(A);
         matrix P=primary_char0(REY,M);
         print(P);
}

proc primary_charp (matrix REY,string ring_name,list #)
"USAGE:   primary_charp(REY,ringname[,v]);
         REY: a <matrix> representing the Reynolds operator, ringname: a
         <string> giving the name of a ring where the Molien series is stored,
         v: an optional <int>
ASSUME:  REY is the first return value of group_reynolds or reynolds_molien and
         ringname gives the name of a ring of characteristic 0 that has been
         created by molien or reynolds_molien
DISPLAY: information about the various stages of the programme if v does not
         equal 0
RETURN:  primary invariants (type <matrix>) of the invariant ring
EXAMPLE: example primary_charp; shows an example
THEORY:  Bases of homogeneous invariants are generated successively and those
         are chosen as primary invariants that lower the dimension of the ideal
         generated by the previously found invariants (see paper \"Generating a
         Noetherian Normalization of the Invariant Ring of a Finite Group\" by
         Decker, Heydtmann, Schreyer (1997) to appear in JSC).
"
{ degBound=0;
// ---------------- checking input and setting verbose mode -------------------
  if (char(basering)==0)
  { "ERROR:   primary_charp should only be used with rings of characteristic p>0.";
    return();
  }
  if (size(#)>1)
  { "ERROR:   primary_charp can only have three parameters.";
    return();
  }
  if (size(#)==1)
  { if (typeof(#[1])<>"int")
    { "ERROR:   The third parameter should be of type <int>.";
      return();
    }
    else
    { int v=#[1];
    }
  }
  else
  { int v=0;
  }
  def br=basering;
  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 (ncols(REY)<>n)
  { "ERROR:   First parameter ought to be the Reynolds operator."
    return();
  }
  if (typeof(`ring_name`)<>"ring")
  { "ERROR:   Second parameter ought to the name of a ring where the Molien";
    "         is stored.";
    return();
  }
 //----------------------------------------------------------------------------
  if (v && voice<>2)
  { "  We can start looking for primary invariants...";
    "";
  }
  if (v && voice==2)
  { "";
  }
 //----------------------- initializing variables -----------------------------
  int dB;
  setring `ring_name`;                 // the Molien series is stores here -
  poly p(1..2);                        // p(1) will be used for single terms of
                                       // the partial expansion, p(2) to store
  p(1..2)=partial_molien(M,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)
  setring br;
  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
                                       // space of invariants of degree d,
                                       // newdim: dimension the ideal generated
                                       // the primary invariants plus basis
                                       // elements, dif=n-i-newdim, i.e. the
                                       // number of new primary invairants that
                                       // should be added in this degree -
  ideal P,Pplus,sPplus,CI,B;           // P: will contain primary invariants,
                                       // Pplus: P+B, CI: a complete
                                       // intersection with the same Hilbert
                                       // function as P
  ideal sP=std(P);
  dB=1;                                // used as degree bound
  int i=0;
 //---------------- loop that searches for primary invariants -----------------
  while(1)                             // repeat until n primary invariants are
  {                                    // found
    setring `ring_name`;
    p(1..2)=partial_molien(M,1,p(2));  // next term of the partial expansion -
    d=deg(p(1));                       // degree where we'll search -
    cd=int(coef(p(1),v1)[2,1]);        // dimension of the homogeneous space of
                                       // inviarants of degree d
    setring br;
    if (v)
    { "  Computing primary invariants in degree "+string(d)+":";
    }
    B=invariant_basis_reynolds(REY,d,intvec(cd,6)); // basis of invariants of
                                       // degree d
    if (B[1]<>0)
    { Pplus=P+B;
      sPplus=std(Pplus);
      newdim=dim(sPplus);
      dif=n-i-newdim;
    }
    else
    { dif=0;
    }
    if (dif<>0)                        // we have to find dif new primary
    {                                  // invariants
      if (cd<>dif)
      { P,sP,CI,dB=p_search(n,d,B,cd,P,sP,i,dif,dB,CI);
      }
      else                             // i.e. we can take all of B
      { for(j=i+1;j>i+dif;j=j+1)
        { CI=CI+ideal(var(j)^d);
        }
        dB=dB+dif*(d-1);
        P=Pplus;
        sP=sPplus;
      }
      if (v)
      { for (j=1;j<=size(P)-i;j=j+1)
        { "  We find: "+string(P[i+j]);
        }
      }
      i=size(P);
      if (i==n)                        // found all primary invariants
      { if (v)
        { "";
          "  We found all primary invariants.";
          "";
        }
        return(matrix(P));
      }
    }                                  // done with degree d
  }
}
example
{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed into";
  "         characteristic 3)";
  echo=2;
         ring R=3,(x,y,z),dp;
         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
         list L=group_reynolds(A);
         string newring="alskdfj";
         molien(L[2..size(L)],newring);
         matrix P=primary_charp(L[1],newring);
         kill `newring`;
         print(P);
}

proc primary_char0_no_molien (matrix REY, list #)
"USAGE:   primary_char0_no_molien(REY[,v]);
         REY: a <matrix> representing the Reynolds operator, v: an optional
         <int>
ASSUME:  REY is the first return value of group_reynolds or reynolds_molien
DISPLAY: information about the various stages of the programme if v does not
         equal 0
RETURN:  primary invariants (type <matrix>) of the invariant ring and an
         <intvec> listing some of the degrees where no non-trivial homogeneous
         invariants are to be found
EXAMPLE: example primary_char0_no_molien; shows an example
THEORY:  Bases of homogeneous invariants are generated successively and those
         are chosen as primary invariants that lower the dimension of the ideal
         generated by the previously found invariants (see paper \"Generating a
         Noetherian Normalization of the Invariant Ring of a Finite Group\" by
         Decker, Heydtmann, Schreyer (1997) to appear in JSC).
"
{ degBound=0;
 //-------------- checking input and setting verbose mode ---------------------
  if (char(basering)<>0)
  { "ERROR:   primary_char0_no_molien should only be used with rings of";
    "         characteristic 0.";
    return();
  }
  if (size(#)>1)
  { "ERROR:   primary_char0_no_molien can only have two parameters.";
    return();
  }
  if (size(#)==1)
  { if (typeof(#[1])<>"int")
    { "ERROR:   The second parameter should be of type <int>.";
      return();
    }
    else
    { int v=#[1];
    }
  }
  else
  { int v=0;
  }
  int n=nvars(basering);               // 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 (ncols(REY)<>n)
  { "ERROR:   First parameter ought to be the Reynolds operator."
    return();
  }
 //----------------------------------------------------------------------------
  if (v && voice<>2)
  { "  We can start looking for primary invariants...";
    "";
  }
  if (v && voice==2)
  { "";
  }
 //----------------------- initializing variables -----------------------------
  int dB;
  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
                                       // space of invariants of degree d,
                                       // newdim: dimension the ideal generated
                                       // the primary invariants plus basis
                                       // elements, dif=n-i-newdim, i.e. the
                                       // number of new primary invairants that
                                       // should be added in this degree -
  ideal P,Pplus,CI,B;                  // P: will contain primary invariants,
                                       // Pplus: P+B, CI: a complete
                                       // intersection with the same Hilbert
                                       // function as P
  ideal sP=std(P);
  dB=1;                                // used as degree bound -
  d=0;                                 // initializing
  int i=0;
  intvec deg_vector;
 //------------------ loop that searches for primary invariants ---------------
  while(1)                             // repeat until n primary invariants are
  {                                    // found -
    d=d+1;                             // degree where we'll search
    if (v)
    { "  Computing primary invariants in degree "+string(d)+":";
    }
    B=invariant_basis_reynolds(REY,d,intvec(-1,6)); // basis of invariants of
                                       // degree d
    if (B[1]<>0)
    { Pplus=P+B;
      newdim=dim(std(Pplus));
      dif=n-i-newdim;
    }
    else
    { dif=0;
      deg_vector=deg_vector,d;
    }
    if (dif<>0)                        // we have to find dif new primary
    {                                  // invariants
      cd=size(B);
      if (cd<>dif)
      { P,sP,CI,dB=search(n,d,B,cd,P,sP,i,dif,dB,CI);
      }
      else                             // i.e. we can take all of B
      { for(j=i+1;j<=i+dif;j=j+1)
        { CI=CI+ideal(var(j)^d);
        }
        dB=dB+dif*(d-1);
        P=Pplus;
        sP=std(P);
      }
      if (v)
      { for (j=1;j<=dif;j=j+1)
        { "  We find: "+string(P[i+j]);
        }
      }
      i=i+dif;
      if (i==n)                        // found all primary invariants
      { if (v)
        { "";
          "  We found all primary invariants.";
          "";
        }
        if (deg_vector==0)
        { return(matrix(P));
        }
        else
        { return(matrix(P),compress(deg_vector));
        }
      }
    }                                  // done with degree d
    else
    { if (v)
      { "  None here...";
      }
    }
  }
}
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;
         list L=group_reynolds(A);
         list l=primary_char0_no_molien(L[1]);
         print(l[1]);
}

proc primary_charp_no_molien (matrix REY, list #)
"USAGE:   primary_charp_no_molien(REY[,v]);
         REY: a <matrix> representing the Reynolds operator, v: an optional
         <int>
ASSUME:  REY is the first return value of group_reynolds or reynolds_molien
DISPLAY: information about the various stages of the programme if v does not
         equal 0
RETURN:  primary invariants (type <matrix>) of the invariant ring  and an
         <intvec> listing some of the degrees where no non-trivial homogeneous
         invariants are to be found
EXAMPLE: example primary_charp_no_molien; shows an example
THEORY:  Bases of homogeneous invariants are generated successively and those
         are chosen as primary invariants that lower the dimension of the ideal
         generated by the previously found invariants (see paper \"Generating a
         Noetherian Normalization of the Invariant Ring of a Finite Group\" by
         Decker, Heydtmann, Schreyer (1997) to appear in JSC).
"
{ degBound=0;
 //----------------- checking input and setting verbose mode ------------------
  if (char(basering)==0)
  { "ERROR:   primary_charp_no_molien should only be used with rings of";
    "         characteristic p>0.";
    return();
  }
  if (size(#)>1)
  { "ERROR:   primary_charp_no_molien can only have two parameters.";
    return();
  }
  if (size(#)==1)
  { if (typeof(#[1])<>"int")
    { "ERROR:   The second parameter should be of type <int>.";
      return();
    }
    else
    { int v=#[1];
    }
  }
  else
  { int v=0;
  }
  int n=nvars(basering);               // 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 (ncols(REY)<>n)
  { "ERROR:   First parameter ought to be the Reynolds operator."
    return();
  }
 //----------------------------------------------------------------------------
  if (v && voice<>2)
  { "  We can start looking for primary invariants...";
    "";
  }
  if (v && voice==2)
  { "";
  }
 //-------------------- initializing variables --------------------------------
  int dB;
  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
                                       // space of invariants of degree d,
                                       // newdim: dimension the ideal generated
                                       // the primary invariants plus basis
                                       // elements, dif=n-i-newdim, i.e. the
                                       // number of new primary invairants that
                                       // should be added in this degree -
  ideal P,Pplus,sPplus,CI,B;           // P: will contain primary invariants,
                                       // Pplus: P+B, CI: a complete
                                       // intersection with the same Hilbert
                                       // function as P
  ideal sP=std(P);
  dB=1;                                // used as degree bound -
  d=0;                                 // initializing
  int i=0;
  intvec deg_vector;
 //------------------ loop that searches for primary invariants ---------------
  while(1)                             // repeat until n primary invariants are
  {                                    // found -
    d=d+1;                             // degree where we'll search
    if (v)
    { "  Computing primary invariants in degree "+string(d)+":";
    }
    B=invariant_basis_reynolds(REY,d,intvec(-1,6)); // basis of invariants of
                                       // degree d
    if (B[1]<>0)
    { Pplus=P+B;
      sPplus=std(Pplus);
      newdim=dim(sPplus);
      dif=n-i-newdim;
    }
    else
    { dif=0;
      deg_vector=deg_vector,d;
    }
    if (dif<>0)                        // we have to find dif new primary
    {                                  // invariants
      cd=size(B);
      if (cd<>dif)
      { P,sP,CI,dB=p_search(n,d,B,cd,P,sP,i,dif,dB,CI);
      }
      else                             // i.e. we can take all of B
      { for(j=i+1;j<=i+dif;j=j+1)
        { CI=CI+ideal(var(j)^d);
        }
        dB=dB+dif*(d-1);
        P=Pplus;
        sP=sPplus;
      }
      if (v)
      { for (j=1;j<=size(P)-i;j=j+1)
        { "  We find: "+string(P[i+j]);
        }
      }
      i=size(P);
      if (i==n)                        // found all primary invariants
      { if (v)
        { "";
          "  We found all primary invariants.";
          "";
        }
        if (deg_vector==0)
        { return(matrix(P));
        }
        else
        { return(matrix(P),compress(deg_vector));
        }
      }
    }                                  // done with degree d
    else
    { if (v)
      { "  None here...";
      }
    }
  }
}
example
{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed into";
  "         characteristic 3)";
  echo=2;
         ring R=3,(x,y,z),dp;
         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
         list L=group_reynolds(A);
         list l=primary_charp_no_molien(L[1]);
         print(l[1]);
}

proc primary_charp_without (list #)
"USAGE:   primary_charp_without(G1,G2,...[,v]);
         G1,G2,...: <matrices> generating a finite matrix group, v: an optional
         <int>
DISPLAY: information about the various stages of the programme if v does not
         equal 0
RETURN:  primary invariants (type <matrix>) of the invariant ring
EXAMPLE: example primary_charp_without; shows an example
THEORY:  Bases of homogeneous invariants are generated successively and those
         are chosen as primary invariants that lower the dimension of the ideal
         generated by the previously found invariants (see paper \"Generating a
         Noetherian Normalization of the Invariant Ring of a Finite Group\" by
         Decker, Heydtmann, Schreyer (1997) to appear in JSC). No Reynolds
         operator or Molien series is used.
"
{ degBound=0;
 //--------------------- checking input and setting verbose mode --------------
  if (char(basering)==0)
  { "ERROR:   primary_charp_without should only be used with rings of";
    "         characteristic 0.";
    return();
  }
  if (size(#)==0)
  { "ERROR:   There are no parameters.";
    return();
  }
  if (typeof(#[size(#)])=="int")
  { int v=#[size(#)];
    int gen_num=size(#)-1;
    if (gen_num==0)
    { "ERROR:   There are no generators of a finite matrix group given.";
      return();
    }
  }
  else
  { int v=0;
    int gen_num=size(#);
  }
  int n=nvars(basering);               // 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
  for (int i=1;i<=gen_num;i=i+1)
  { if (typeof(#[i])=="matrix")
    { if (nrows(#[i])<>n or ncols(#[i])<>n)
      { "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 first parameters should be a list of matrices";
      return();
    }
  }
 //----------------------------------------------------------------------------
  if (v && voice==2)
  { "";
  }
 //---------------------------- initializing variables ------------------------
  int dB;
  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
                                       // space of invariants of degree d,
                                       // newdim: dimension the ideal generated
                                       // the primary invariants plus basis
                                       // elements, dif=n-i-newdim, i.e. the
                                       // number of new primary invairants that
                                       // should be added in this degree -
  ideal P,Pplus,sPplus,CI,B;           // P: will contain primary invariants,
                                       // Pplus: P+B, CI: a complete
                                       // intersection with the same Hilbert
                                       // function as P
  ideal sP=std(P);
  dB=1;                                // used as degree bound -
  d=0;                                 // initializing
  i=0;
  intvec deg_vector;
 //-------------------- loop that searches for primary invariants -------------
  while(1)                             // repeat until n primary invariants are
  {                                    // found -
    d=d+1;                             // degree where we'll search
    if (v)
    { "  Computing primary invariants in degree "+string(d)+":";
    }
    B=invariant_basis(d,#[1..gen_num]); // basis of invariants of degree d
    if (B[1]<>0)
    { Pplus=P+B;
      sPplus=std(Pplus);
      newdim=dim(sPplus);
      dif=n-i-newdim;
    }
    else
    { dif=0;
      deg_vector=deg_vector,d;
    }
    if (dif<>0)                        // we have to find dif new primary
    {                                  // invariants
      cd=size(B);
      if (cd<>dif)
      { P,sP,CI,dB=p_search(n,d,B,cd,P,sP,i,dif,dB,CI);
      }
      else                             // i.e. we can take all of B
      { for(j=i+1;j<=i+dif;j=j+1)
        { CI=CI+ideal(var(j)^d);
        }
        dB=dB+dif*(d-1);
        P=Pplus;
        sP=sPplus;
      }
      if (v)
      { for (j=1;j<=size(P)-i;j=j+1)
        { "  We find: "+string(P[i+j]);
        }
      }
      i=size(P);
      if (i==n)                        // found all primary invariants
      { if (v)
        { "";
          "  We found all primary invariants.";
          "";
        }
        return(matrix(P));
      }
    }                                  // done with degree d
    else
    { if (v)
      { "  None here...";
      }
    }
  }
}
example
{ echo=2;
         ring R=2,(x,y,z),dp;
         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
         matrix P=primary_charp_without(A);
         print(P);
}

proc primary_invariants (list #)
"USAGE:   primary_invariants(G1,G2,...[,flags]);
         G1,G2,...: <matrices> generating a finite matrix group, flags: an
         optional <intvec> with three entries, if the first one equals 0 (also
         the default), the programme attempts to compute the Molien series and
         Reynolds operator, if it equals 1, the programme is told that the
         Molien series should not be computed, if it equals -1 characteristic 0
         is simulated, i.e. the Molien series is computed as if the base field
         were characteristic 0 (the user must choose a field of large prime
         characteristic, e.g. 32003) and if the first one is anything else, it
         means that the characteristic of the base field divides the group
         order, the second component should give the size of intervals between
         canceling common factors in the expansion of the Molien series, 0 (the
         default) means only once after generating all terms, in prime
         characteristic also a negative number can be given to indicate that
         common factors should always be canceled when the expansion is simple
         (the root of the extension field does not occur among the coefficients)
DISPLAY: information about the various stages of the programme if the third
         flag does not equal 0
RETURN:  primary invariants (type <matrix>) of the invariant ring and if
         computable Reynolds operator (type <matrix>) and Molien series (type
         <matrix>), if the first flag is 1 and we are in the non-modular case
         then an <intvec> is returned giving some of the degrees where no
         non-trivial homogeneous invariants can be found
EXAMPLE: example primary_invariants; shows an example
THEORY:  Bases of homogeneous invariants are generated successively and those
         are chosen as primary invariants that lower the dimension of the ideal
         generated by the previously found invariants (see paper \"Generating a
         Noetherian Normalization of the Invariant Ring of a Finite Group\" by
         Decker, Heydtmann, Schreyer (1997) to appear in JSC).
"
{
 // ----------------- checking input and setting flags ------------------------
  if (size(#)==0)
  { "ERROR:   There are no parameters.";
    return();
  }
  int ch=char(basering);               // the algorithms depend very much on the
                                       // characteristic of the ground field
  int n=nvars(basering);               // 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
  int gen_num;
  int mol_flag,v;
  if (typeof(#[size(#)])=="intvec")
  { if (size(#[size(#)])<>3)
    { "ERROR:   <intvec> should have three entries.";
      return();
    }
    gen_num=size(#)-1;
    mol_flag=#[size(#)][1];
    if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag==-1)))
    { "ERROR:   the second component of <intvec> should be >=0";
      return();
    }
    int interval=#[size(#)][2];
    v=#[size(#)][3];
    if (gen_num==0)
    { "ERROR:   There are no generators of a finite matrix group given.";
      return();
    }
  }
  else
  { gen_num=size(#);
    mol_flag=0;
    int interval=0;
    v=0;
  }
  for (int i=1;i<=gen_num;i=i+1)
  { if (typeof(#[i])=="matrix")
    { if (nrows(#[i])<>n or ncols(#[i])<>n)
      { "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 first parameters should be a list of matrices";
      return();
    }
  }
 //----------------------------------------------------------------------------
  if (mol_flag==0)
  { if (ch==0)
    { matrix REY,M=reynolds_molien(#[1..gen_num],intvec(mol_flag,interval,v));
                                       // one will contain Reynolds operator and
                                       // the other enumerator and denominator
                                       // of Molien series
      matrix P=primary_char0(REY,M,v);
      return(P,REY,M);
    }
    else
    { list L=group_reynolds(#[1..gen_num],v);
      if (L[1]<>0)                     // testing whether we are in the modular
      { string newring="aksldfalkdsflkj"; // case
        if (minpoly==0)
        { if (v)
          { "  We are dealing with the non-modular case.";
          }
          if (typeof(L[2])=="int")
          { molien(L[3..size(L)],newring,L[2],intvec(mol_flag,interval,v));
          }
          else
          { molien(L[2..size(L)],newring,intvec(mol_flag,interval,v));
          }
          matrix P=primary_charp(L[1],newring,v);
          return(P,L[1],newring);
        }
        else
        { if (v)
          { "  Since it is impossible for this programme to calculate the Molien series for";
            "  invariant rings over extension fields of prime characteristic, we have to";
            "  continue without it.";
            "";

          }
          list l=primary_charp_no_molien(L[1],v);
          if (size(l)==2)
          { return(l[1],L[1],l[2]);
          }
          else
          { return(l[1],L[1]);
          }
        }
      }
      else                             // the modular case
      { if (v)
        { "  There is also no Molien series, we can make use of...";
          "";
          "  We can start looking for primary invariants...";
          "";
        }
        return(primary_charp_without(#[1..gen_num],v));
      }
    }
  }
  if (mol_flag==1)                     // the user wants no calculation of the
  { list L=group_reynolds(#[1..gen_num],v); // Molien series
    if (ch==0)
    { list l=primary_char0_no_molien(L[1],v);
      if (size(l)==2)
      { return(l[1],L[1],l[2]);
      }
      else
      { return(l[1],L[1]);
      }
    }
    else
    { if (L[1]<>0)                     // testing whether we are in the modular
      { list l=primary_charp_no_molien(L[1],v); // case
        if (size(l)==2)
        { return(l[1],L[1],l[2]);
        }
        else
        { return(l[1],L[1]);
        }
      }
      else                             // the modular case
      { if (v)
        { "  We can start looking for primary invariants...";
          "";
        }
        return(primary_charp_without(#[1..gen_num],v));
      }
    }
  }
  if (mol_flag==-1)
  { if (ch==0)
    { "ERROR:   Characteristic 0 can only be simulated in characteristic p>>0.";
      return();
    }
    list L=group_reynolds(#[1..gen_num],v);
    string newring="aksldfalkdsflkj";
    molien(L[2..size(L)],newring,intvec(1,interval,v));
    matrix P=primary_charp(L[1],newring,v);
    return(P,L[1],newring);
  }
  else                                 // the user specified that the
  { if (ch==0)                         // characteristic divides the group order
    { "ERROR:   The characteristic cannot divide the group order when it is 0.";
      return();
    }
    if (v)
    { "";
    }
    return(primary_charp_without(#[1..gen_num],v));
  }
}
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;
         list L=primary_invariants(A);
         print(L[1]);
}

////////////////////////////////////////////////////////////////////////////////
// This procedure finds dif primary invariants in degree d. It returns all
// primary invariants found so far. The coefficients lie in a field of
// characteristic 0.
////////////////////////////////////////////////////////////////////////////////
proc search_random (int n,int d,ideal B,int cd,ideal P,int i,int dif,int dB,ideal CI,int max)
{ string answer;
  degBound=0;
  int j,k,test_dim,flag;
  matrix test_matrix[1][dif];          // the linear combination to test
  intvec h;                            // Hilbert series
  for (j=i+1;j<=i+dif;j=j+1)
  { CI=CI+ideal(var(j)^d);             // homogeneous polynomial of the same
                                       // degree as the one we're looking for
                                       // is added
  }
  ideal TEST;
  // h=hilb(std(CI),1);
  dB=dB+dif*(d-1);                     // used as degBound
  while (1)
  { test_matrix=matrix(B)*random(max,cd,dif);
    // degBound=dB;
    TEST=P+ideal(test_matrix);
    attrib(TEST,"isSB",1);
    test_dim=dim(TEST);
    // degBound=0;
    if (n-test_dim==i+dif)
    { break;
    }
    // degBound=dB;
    test_dim=dim(std(TEST));
    // test_dim=dim(std(TEST,h));         // Hilbert driven std-calculation
    // degBound=0;
    if (n-test_dim==i+dif)
    { break;
    }
    else
    { "HELP:    The "+string(dif)+" random combination(s) of the "+string(cd)+" basis elements with";
      "         coefficients in the range from -"+string(max)+" to "+string(max)+" did not lower the";
      "         dimension by "+string(dif)+". You can abort, try again or give a new range:";
      answer="";
      while (answer<>"n
" && answer<>"y
")
      { "         Do you want to abort (y/n)?";
        answer=read("");
      }
      if (answer=="y
")
      { flag=1;
        break;
      }
      answer="";
      while (answer<>"n
" && answer<>"y
")
      { "         Do you want to try again (y/n)?";
        answer=read("");
      }
      if (answer=="n
")
      { flag=1;
        while (flag)
        { "         Give a new <int> > "+string(max)+" that bounds the range of coefficients:";
          answer=read("");
          for (j=1;j<=size(answer)-1;j=j+1)
          { for (k=0;k<=9;k=k+1)
            { if (answer[j]==string(k))
              { break;
              }
            }
            if (k>9)
            { flag=1;
              break;
            }
            flag=0;
          }
          if (not(flag))
          { execute "test_dim="+string(answer[1..size(answer)]);
            if (test_dim<=max)
            { flag=1;
            }
            else
            { max=test_dim;
            }
          }
        }
      }
    }
  }
  if (not(flag))
  { P[(i+1)..(i+dif)]=test_matrix[1,1..dif];
  }
  return(P,CI,dB);
}

////////////////////////////////////////////////////////////////////////////////
// This procedure finds at most dif primary invariants in degree d. It returns
// all primary invariants found so far. The coefficients lie in the field of
// characteristic p>0.
////////////////////////////////////////////////////////////////////////////////
proc p_search_random (int n,int d,ideal B,int cd,ideal P,int i,int dif,int dB,ideal CI,int max)
{ string answer;
  degBound=0;
  int j,k,test_dim,flag;
  matrix test_matrix[1][dif];          // the linear combination to test
  intvec h;                            // Hilbert series
  ideal TEST;
  while (dif>0)
  { for (j=i+1;j<=i+dif;j=j+1)
    { CI=CI+ideal(var(j)^d);           // homogeneous polynomial of the same
                                       // degree as the one we're looking for
                                       // is added
    }
    // h=hilb(std(CI),1);
    dB=dB+dif*(d-1);                   // used as degBound
    test_matrix=matrix(B)*random(max,cd,dif);
    // degBound=dB;
    TEST=P+ideal(test_matrix);
    attrib(TEST,"isSB",1);
    test_dim=dim(TEST);
    // degBound=0;
    if (n-test_dim==i+dif)
    { break;
    }
    // degBound=dB;
    test_dim=dim(std(TEST));
    // test_dim=dim(std(TEST,h));         // Hilbert driven std-calculation
    // degBound=0;
    if (n-test_dim==i+dif)
    { break;
    }
    else
    { "HELP:    The "+string(dif)+" random combination(s) of the "+string(cd)+" basis elements with";
      "         coefficients in the range from -"+string(max)+" to "+string(max)+" did not lower the";
      "         dimension by "+string(dif)+". You can abort, try again, lower the number of";
      "         combinations searched for by 1 or give a larger coefficient range:";
      answer="";
      while (answer<>"n
" && answer<>"y
")
      { "         Do you want to abort (y/n)?";
        answer=read("");
      }
      if (answer=="y
")
      { flag=1;
        break;
      }
      answer="";
      while (answer<>"n
" && answer<>"y
")
      { "         Do you want to try again (y/n)?";
        answer=read("");
      }
      if (answer=="n
")
      { answer="";
        while (answer<>"n
" && answer<>"y
")
        { "         Do you want to lower the number of combinations by 1 (y/n)?";
          answer=read("");
        }
        if (answer=="y
")
        { dif=dif-1;
        }
        else
        { flag=1;
          while (flag)
          { "         Give a new <int> > "+string(max)+" that bounds the range of coefficients:";
            answer=read("");
            for (j=1;j<=size(answer)-1;j=j+1)
            { for (k=0;k<=9;k=k+1)
              { if (answer[j]==string(k))
                { break;
                }
              }
              if (k>9)
              { flag=1;
                break;
              }
              flag=0;
            }
            if (not(flag))
            { execute "test_dim="+string(answer[1..size(answer)]);
              if (test_dim<=max)
              { flag=1;
              }
              else
              { max=test_dim;
              }
            }
          }
        }
      }
    }
    CI=CI[1..i];
    dB=dB-dif*(d-1);
  }
  if (dif && not(flag))
  { P[(i+1)..(i+dif)]=test_matrix[1,1..dif];
  }
  if (dif && flag)
  { P[n+1]=0;
  }
  return(P,CI,dB);
}

proc primary_char0_random (matrix REY,matrix M,int max,list #)
"USAGE:   primary_char0_random(REY,M,r[,v]);
         REY: a <matrix> representing the Reynolds operator, M: a 1x2 <matrix>
         representing the Molien series, r: an <int> where -|r| to |r| is the
         range of coefficients of the random combinations of bases elements,
         v: an optional <int>
ASSUME:  REY is the first return value of group_reynolds or reynolds_molien and
         M the one of molien or the second one of reynolds_molien
DISPLAY: information about the various stages of the programme if v does not
         equal 0
RETURN:  primary invariants (type <matrix>) of the invariant ring
EXAMPLE: example primary_char0_random; shows an example
THEORY:  Bases of homogeneous invariants are generated successively and random
         linear combinations are chosen as primary invariants that lower the
         dimension of the ideal generated by the previously found invariants
         (see paper \"Generating a Noetherian Normalization of the Invariant Ring
         of a Finite Group\" by Decker, Heydtmann, Schreyer (1997) to appear in
         JSC).
"
{ degBound=0;
  if (char(basering)<>0)
  { "ERROR:   primary_char0_random should only be used with rings of";
    "         characteristic 0.";
    return();
  }
 //----------------- checking input and setting verbose mode ------------------
  if (size(#)>1)
  { "ERROR:   primary_char0_random can only have four parameters.";
    return();
  }
  if (size(#)==1)
  { if (typeof(#[1])<>"int")
    { "ERROR:   The fourth parameter should be of type <int>.";
      return();
    }
    else
    { int v=#[1];
    }
  }
  else
  { int v=0;
  }
  int n=nvars(basering);               // 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 (ncols(REY)<>n)
  { "ERROR:   First parameter ought to be the Reynolds operator."
    return();
  }
  if (ncols(M)<>2 or nrows(M)<>1)
  { "ERROR:   Second parameter ought to be the Molien series."
    return();
  }
 //----------------------------------------------------------------------------
  if (v && voice<>2)
  { "  We can start looking for primary invariants...";
    "";
  }
  if (v && voice==2)
  { "";
  }
 //------------------------- initializing variables ---------------------------
  int dB;
  poly p(1..2);                        // p(1) will be used for single terms of
                                       // the partial expansion, p(2) to store
  p(1..2)=partial_molien(M,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) -
  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
                                       // space of invariants of degree d,
                                       // newdim: dimension the ideal generated
                                       // the primary invariants plus basis
                                       // elements, dif=n-i-newdim, i.e. the
                                       // number of new primary invairants that
                                       // should be added in this degree -
  ideal P,Pplus,CI,B;                  // P: will contain primary invariants,
                                       // Pplus: P+B,CI: a complete
                                       // intersection with the same Hilbert
                                       // function as P -
  dB=1;                                // used as degree bound
  int i=0;
 //-------------- loop that searches for primary invariants  ------------------
  while(1)                             // repeat until n primary invariants are
  {                                    // found -
    p(1..2)=partial_molien(M,1,p(2));  // next term of the partial expansion -
    d=deg(p(1));                       // degree where we'll search -
    cd=int(coef(p(1),v1)[2,1]);        // dimension of the homogeneous space of
                                       // inviarants of degree d
    if (v)
    { "  Computing primary invariants in degree "+string(d)+":";
    }
    B=invariant_basis_reynolds(REY,d,intvec(cd,6)); // basis of invariants of
                                       // degree d
    if (B[1]<>0)
    { Pplus=P+B;
      newdim=dim(std(Pplus));
      dif=n-i-newdim;
    }
    else
    { dif=0;
    }
    if (dif<>0)                        // we have to find dif new primary
    {                                  // invariants
      if (cd<>dif)
      { P,CI,dB=search_random(n,d,B,cd,P,i,dif,dB,CI,max); // searching for
      }                                // dif invariants -
      else                             // i.e. we can take all of B
      { for(j=i+1;j>i+dif;j=j+1)
        { CI=CI+ideal(var(j)^d);
        }
        dB=dB+dif*(d-1);
        P=Pplus;
      }
      if (ncols(P)==i)
      { "WARNING: The return value is not a set of primary invariants, but";
        "         polynomials qualifying as the first "+string(i)+" primary invariants.";
        return(matrix(P));
      }
      if (v)
      { for (j=1;j<=dif;j=j+1)
        { "  We find: "+string(P[i+j]);
        }
      }
      i=i+dif;
      if (i==n)                        // found all primary invariants
      { if (v)
        { "";
          "  We found all primary invariants.";
          "";
        }
        return(matrix(P));
      }
    }                                  // done with degree d
  }
}
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 REY,M=reynolds_molien(A);
         matrix P=primary_char0_random(REY,M,1);
         print(P);
}

proc primary_charp_random (matrix REY,string ring_name,int max,list #)
"USAGE:   primary_charp_random(REY,ringname,r[,v]);
         REY: a <matrix> representing the Reynolds operator, ringname: a
         <string> giving the name of a ring where the Molien series is stored,
         r: an <int> where -|r| to |r| is the range of coefficients of the
         random combinations of bases elements, v: an optional <int>
ASSUME:  REY is the first return value of group_reynolds or reynolds_molien and
         ringname gives the name of a ring of characteristic 0 that has been
         created by molien or reynolds_molien
DISPLAY: information about the various stages of the programme if v does not
         equal 0
RETURN:  primary invariants (type <matrix>) of the invariant ring
EXAMPLE: example primary_charp_random; shows an example
THEORY:  Bases of homogeneous invariants are generated successively and random
         linear combinations are chosen as primary invariants that lower the
         dimension of the ideal generated by the previously found invariants
         (see paper \"Generating a Noetherian Normalization of the Invariant Ring
         of a Finite Group\" by Decker, Heydtmann, Schreyer (1997) to appear in
         JSC).
"
{ degBound=0;
 // ---------------- checking input and setting verbose mode ------------------
  if (char(basering)==0)
  { "ERROR:   primary_charp_random should only be used with rings of";
    "         characteristic p>0.";
    return();
  }
  if (size(#)>1)
  { "ERROR:   primary_charp_random can only have four parameters.";
    return();
  }
  if (size(#)==1)
  { if (typeof(#[1])<>"int")
    { "ERROR:   The fourth parameter should be of type <int>.";
      return();
    }
    else
    { int v=#[1];
    }
  }
  else
  { int v=0;
  }
  def br=basering;
  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 (ncols(REY)<>n)
  { "ERROR:   First parameter ought to be the Reynolds operator."
    return();
  }
  if (typeof(`ring_name`)<>"ring")
  { "ERROR:   Second parameter ought to the name of a ring where the Molien";
    "         is stored.";
    return();
  }
 //----------------------------------------------------------------------------
  if (v && voice<>2)
  { "  We can start looking for primary invariants...";
    "";
  }
  if (v && voice==2)
  { "";
  }
 //----------------------- initializing variables -----------------------------
  int dB;
  setring `ring_name`;                 // the Molien series is stores here -
  poly p(1..2);                        // p(1) will be used for single terms of
                                       // the partial expansion, p(2) to store
  p(1..2)=partial_molien(M,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)
  setring br;
  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
                                       // space of invariants of degree d,
                                       // newdim: dimension the ideal generated
                                       // the primary invariants plus basis
                                       // elements, dif=n-i-newdim, i.e. the
                                       // number of new primary invairants that
                                       // should be added in this degree -
  ideal P,Pplus,CI,B;                  // P: will contain primary invariants,
                                       // Pplus: P+B, CI: a complete
                                       // intersection with the same Hilbert
                                       // function as P -
  dB=1;                                // used as degree bound
  int i=0;
 //---------------- loop that searches for primary invariants -----------------
  while(1)                             // repeat until n primary invariants are
  {                                    // found
    setring `ring_name`;
    p(1..2)=partial_molien(M,1,p(2));  // next term of the partial expansion -
    d=deg(p(1));                       // degree where we'll search -
    cd=int(coef(p(1),v1)[2,1]);        // dimension of the homogeneous space of
                                       // inviarants of degree d
    setring br;
    if (v)
    { "  Computing primary invariants in degree "+string(d)+":";
    }
    B=invariant_basis_reynolds(REY,d,intvec(cd,6)); // basis of invariants of
                                       // degree d
    if (B[1]<>0)
    { Pplus=P+B;
      newdim=dim(std(Pplus));
      dif=n-i-newdim;
    }
    else
    { dif=0;
    }
    if (dif<>0)                        // we have to find dif new primary
    {                                  // invariants
      if (cd<>dif)
      { P,CI,dB=p_search_random(n,d,B,cd,P,i,dif,dB,CI,max);
      }
      else                             // i.e. we can take all of B
      { for(j=i+1;j>i+dif;j=j+1)
        { CI=CI+ideal(var(j)^d);
        }
        dB=dB+dif*(d-1);
        P=Pplus;
      }
      if (ncols(P)==n+1)
      { "WARNING: The first return value is not a set of primary invariants,";
        "         but polynomials qualifying as the first "+string(i)+" primary invariants.";
        return(matrix(P));
      }
      if (v)
      { for (j=1;j<=size(P)-i;j=j+1)
        { "  We find: "+string(P[i+j]);
        }
      }
      i=size(P);
      if (i==n)                        // found all primary invariants
      { if (v)
        { "";
          "  We found all primary invariants.";
          "";
        }
        return(matrix(P));
      }
    }                                  // done with degree d
  }
}
example
{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed into";
  "         characteristic 3)";
  echo=2;
         ring R=3,(x,y,z),dp;
         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
         list L=group_reynolds(A);
         string newring="alskdfj";
         molien(L[2..size(L)],newring);
         matrix P=primary_charp_random(L[1],newring,1);
         kill `newring`;
         print(P);
}

proc primary_char0_no_molien_random (matrix REY, int max, list #)
"USAGE:   primary_char0_no_molien_random(REY,r[,v]);
         REY: a <matrix> representing the Reynolds operator, r: an <int> where
         -|r| to |r| is the range of coefficients of the random combinations of
         bases elements, v: an optional <int>
ASSUME:  REY is the first return value of group_reynolds or reynolds_molien
DISPLAY: information about the various stages of the programme if v does not
         equal 0
RETURN:  primary invariants (type <matrix>) of the invariant ring  and an
         <intvec> listing some of the degrees where no non-trivial homogeneous
         invariants are to be found
EXAMPLE: example primary_char0_no_molien_random; shows an example
THEORY:  Bases of homogeneous invariants are generated successively and random
         linear combinations are chosen as primary invariants that lower the
         dimension of the ideal generated by the previously found invariants
         (see paper \"Generating a Noetherian Normalization of the Invariant Ring
         of a Finite Group\" by Decker, Heydtmann, Schreyer (1997) to appear in
         JSC).
"
{ degBound=0;
 //-------------- checking input and setting verbose mode ---------------------
  if (char(basering)<>0)
  { "ERROR:   primary_char0_no_molien_random should only be used with rings of";
    "         characteristic 0.";
    return();
  }
  if (size(#)>1)
  { "ERROR:   primary_char0_no_molien_random can only have three parameters.";
    return();
  }
  if (size(#)==1)
  { if (typeof(#[1])<>"int")
    { "ERROR:   The third parameter should be of type <int>.";
      return();
    }
    else
    { int v=#[1];
    }
  }
  else
  { int v=0;
  }
  int n=nvars(basering);               // 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 (ncols(REY)<>n)
  { "ERROR:   First parameter ought to be the Reynolds operator."
    return();
  }
 //----------------------------------------------------------------------------
  if (v && voice<>2)
  { "  We can start looking for primary invariants...";
    "";
  }
  if (v && voice==2)
  { "";
  }
 //----------------------- initializing variables -----------------------------
  int dB;
  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
                                       // space of invariants of degree d,
                                       // newdim: dimension the ideal generated
                                       // the primary invariants plus basis
                                       // elements, dif=n-i-newdim, i.e. the
                                       // number of new primary invairants that
                                       // should be added in this degree -
  ideal P,Pplus,CI,B;                  // P: will contain primary invariants,
                                       // Pplus: P+B, CI: a complete
                                       // intersection with the same Hilbert
                                       // function as P -
  dB=1;                                // used as degree bound -
  d=0;                                 // initializing
  int i=0;
  intvec deg_vector;
 //------------------ loop that searches for primary invariants ---------------
  while(1)                             // repeat until n primary invariants are
  {                                    // found -
    d=d+1;                             // degree where we'll search
    if (v)
    { "  Computing primary invariants in degree "+string(d)+":";
    }
    B=invariant_basis_reynolds(REY,d,intvec(-1,6)); // basis of invariants of
                                       // degree d
    if (B[1]<>0)
    { Pplus=P+B;
      newdim=dim(std(Pplus));
      dif=n-i-newdim;
    }
    else
    { dif=0;
      deg_vector=deg_vector,d;
    }
    if (dif<>0)                        // we have to find dif new primary
    {                                  // invariants
      cd=size(B);
      if (cd<>dif)
      { P,CI,dB=search_random(n,d,B,cd,P,i,dif,dB,CI,max);
      }
      else                             // i.e. we can take all of B
      { for(j=i+1;j<=i+dif;j=j+1)
        { CI=CI+ideal(var(j)^d);
        }
        dB=dB+dif*(d-1);
        P=Pplus;
      }
      if (ncols(P)==i)
      { "WARNING: The first return value is not a set of primary invariants,";
        "         but polynomials qualifying as the first "+string(i)+" primary invariants.";
        return(matrix(P));
      }
      if (v)
      { for (j=1;j<=dif;j=j+1)
        { "  We find: "+string(P[i+j]);
        }
      }
      i=i+dif;
      if (i==n)                        // found all primary invariants
      { if (v)
        { "";
          "  We found all primary invariants.";
          "";
        }
        if (deg_vector==0)
        { return(matrix(P));
        }
        else
        { return(matrix(P),compress(deg_vector));
        }
      }
    }                                  // done with degree d
    else
    { if (v)
      { "  None here...";
      }
    }
  }
}
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;
         list L=group_reynolds(A);
         list l=primary_char0_no_molien_random(L[1],1);
         print(l[1]);
}

proc primary_charp_no_molien_random (matrix REY, int max, list #)
"USAGE:   primary_charp_no_molien_random(REY,r[,v]);
         REY: a <matrix> representing the Reynolds operator, r: an <int> where
         -|r| to |r| is the range of coefficients of the random combinations of
         bases elements, v: an optional <int>
ASSUME:  REY is the first return value of group_reynolds or reynolds_molien
DISPLAY: information about the various stages of the programme if v does not
         equal 0
RETURN:  primary invariants (type <matrix>) of the invariant ring  and an
         <intvec> listing some of the degrees where no non-trivial homogeneous
         invariants are to be found
EXAMPLE: example primary_charp_no_molien_random; shows an example
THEORY:  Bases of homogeneous invariants are generated successively and random
         linear combinations are chosen as primary invariants that lower the
         dimension of the ideal generated by the previously found invariants
         (see paper \"Generating a Noetherian Normalization of the Invariant Ring
         of a Finite Group\" by Decker, Heydtmann, Schreyer (1997) to appear in
         JSC).
"
{ degBound=0;
 //----------------- checking input and setting verbose mode ------------------
  if (char(basering)==0)
  { "ERROR:   primary_charp_no_molien_random should only be used with rings of";
    "         characteristic p>0.";
    return();
  }
  if (size(#)>1)
  { "ERROR:   primary_charp_no_molien_random can only have three parameters.";
    return();
  }
  if (size(#)==1)
  { if (typeof(#[1])<>"int")
    { "ERROR:   The third parameter should be of type <int>.";
      return();
    }
    else
    { int v=#[1];
    }
  }
  else
  { int v=0;
  }
  int n=nvars(basering);               // 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 (ncols(REY)<>n)
  { "ERROR:   First parameter ought to be the Reynolds operator."
    return();
  }
 //----------------------------------------------------------------------------
  if (v && voice<>2)
  { "  We can start looking for primary invariants...";
    "";
  }
  if (v && voice==2)
  { "";
  }
 //-------------------- initializing variables --------------------------------
  int dB;
  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
                                       // space of invariants of degree d,
                                       // newdim: dimension the ideal generated
                                       // the primary invariants plus basis
                                       // elements, dif=n-i-newdim, i.e. the
                                       // number of new primary invairants that
                                       // should be added in this degree -
  ideal P,Pplus,CI,B;                  // P: will contain primary invariants,
                                       // Pplus: P+B, CI: a complete
                                       // intersection with the same Hilbert
                                       // function as P -
  dB=1;                                // used as degree bound -
  d=0;                                 // initializing
  int i=0;
  intvec deg_vector;
 //------------------ loop that searches for primary invariants ---------------
  while(1)                             // repeat until n primary invariants are
  {                                    // found -
    d=d+1;                             // degree where we'll search
    if (v)
    { "  Computing primary invariants in degree "+string(d)+":";
    }
    B=invariant_basis_reynolds(REY,d,intvec(-1,6)); // basis of invariants of
                                       // degree d
    if (B[1]<>0)
    { Pplus=P+B;
      newdim=dim(std(Pplus));
      dif=n-i-newdim;
    }
    else
    { dif=0;
      deg_vector=deg_vector,d;
    }
    if (dif<>0)                        // we have to find dif new primary
    {                                  // invariants
      cd=size(B);
      if (cd<>dif)
      { P,CI,dB=p_search_random(n,d,B,cd,P,i,dif,dB,CI,max);
      }
      else                             // i.e. we can take all of B
      { for(j=i+1;j<=i+dif;j=j+1)
        { CI=CI+ideal(var(j)^d);
        }
        dB=dB+dif*(d-1);
        P=Pplus;
      }
      if (ncols(P)==n+1)
      { "WARNING: The first return value is not a set of primary invariants,";
        "         but polynomials qualifying as the first "+string(i)+" primary invariants.";
        return(matrix(P));
      }
      if (v)
      { for (j=1;j<=size(P)-i;j=j+1)
        { "  We find: "+string(P[i+j]);
        }
      }
      i=size(P);
      if (i==n)                        // found all primary invariants
      { if (v)
        { "";
          "  We found all primary invariants.";
          "";
        }
        if (deg_vector==0)
        { return(matrix(P));
        }
        else
        { return(matrix(P),compress(deg_vector));
        }
      }
    }                                  // done with degree d
    else
    { if (v)
      { "  None here...";
      }
    }
  }
}
example
{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed into";
  "         characteristic 3)";
  echo=2;
         ring R=3,(x,y,z),dp;
         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
         list L=group_reynolds(A);
         list l=primary_charp_no_molien_random(L[1],1);
         print(l[1]);
}

proc primary_charp_without_random (list #)
"USAGE:   primary_charp_without_random(G1,G2,...,r[,v]);
         G1,G2,...: <matrices> generating a finite matrix group, r: an <int>
         where -|r| to |r| is the range of coefficients of the random
         combinations of bases elements, v: an optional <int>
DISPLAY: information about the various stages of the programme if v does not
         equal 0
RETURN:  primary invariants (type <matrix>) of the invariant ring
EXAMPLE: example primary_charp_without_random; shows an example
THEORY:  Bases of homogeneous invariants are generated successively and random
         linear combinations are chosen as primary invariants that lower the
         dimension of the ideal generated by the previously found invariants
         (see paper \"Generating a Noetherian Normalization of the Invariant Ring
         of a Finite Group\" by Decker, Heydtmann, Schreyer (1997) to appear in
         JSC). No Reynolds operator or Molien series is used.
"
{ degBound=0;
 //--------------------- checking input and setting verbose mode --------------
  if (char(basering)==0)
  { "ERROR:   primary_charp_without_random should only be used with rings of";
    "         characteristic 0.";
    return();
  }
  if (size(#)<2)
  { "ERROR:   There are too few parameters.";
    return();
  }
  if (typeof(#[size(#)])=="int" && typeof(#[size(#)-1])=="int")
  { int v=#[size(#)];
    int max=#[size(#)-1];
    int gen_num=size(#)-2;
    if (gen_num==0)
    { "ERROR:   There are no generators of a finite matrix group given.";
      return();
    }
  }
  else
  { if (typeof(#[size(#)])=="int")
    { int max=#[size(#)];
      int v=0;
      int gen_num=size(#)-1;
    }
    else
    { "ERROR:   The last parameter should be an <int>.";
      return();
    }
  }
  int n=nvars(basering);               // 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
  for (int i=1;i<=gen_num;i=i+1)
  { if (typeof(#[i])=="matrix")
    { if (nrows(#[i])<>n or ncols(#[i])<>n)
      { "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 first parameters should be a list of matrices";
      return();
    }
  }
 //----------------------------------------------------------------------------
  if (v && voice==2)
  { "";
  }
 //---------------------------- initializing variables ------------------------
  int dB;
  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
                                       // space of invariants of degree d,
                                       // newdim: dimension the ideal generated
                                       // the primary invariants plus basis
                                       // elements, dif=n-i-newdim, i.e. the
                                       // number of new primary invairants that
                                       // should be added in this degree -
  ideal P,Pplus,CI,B;                  // P: will contain primary invariants,
                                       // Pplus: P+B, CI: a complete
                                       // intersection with the same Hilbert
                                       // function as P -
  dB=1;                                // used as degree bound -
  d=0;                                 // initializing
  i=0;
  intvec deg_vector;
 //-------------------- loop that searches for primary invariants -------------
  while(1)                             // repeat until n primary invariants are
  {                                    // found -
    d=d+1;                             // degree where we'll search
    if (v)
    { "  Computing primary invariants in degree "+string(d)+":";
    }
    B=invariant_basis(d,#[1..gen_num]); // basis of invariants of degree d
    if (B[1]<>0)
    { Pplus=P+B;
      newdim=dim(std(Pplus));
      dif=n-i-newdim;
    }
    else
    { dif=0;
      deg_vector=deg_vector,d;
    }
    if (dif<>0)                        // we have to find dif new primary
    {                                  // invariants
      cd=size(B);
      if (cd<>dif)
      { P,CI,dB=p_search_random(n,d,B,cd,P,i,dif,dB,CI,max);
      }
      else                             // i.e. we can take all of B
      { for(j=i+1;j<=i+dif;j=j+1)
        { CI=CI+ideal(var(j)^d);
        }
        dB=dB+dif*(d-1);
        P=Pplus;
      }
      if (ncols(P)==n+1)
      { "WARNING: The first return value is not a set of primary invariants,";
        "         but polynomials qualifying as the first "+string(i)+" primary invariants.";
        return(matrix(P));
      }
      if (v)
      { for (j=1;j<=size(P)-i;j=j+1)
        { "  We find: "+string(P[i+j]);
        }
      }
      i=size(P);
      if (i==n)                        // found all primary invariants
      { if (v)
        { "";
          "  We found all primary invariants.";
          "";
        }
        return(matrix(P));
      }
    }                                  // done with degree d
    else
    { if (v)
      { "  None here...";
      }
    }
  }
}
example
{ echo=2;
         ring R=2,(x,y,z),dp;
         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
         matrix P=primary_charp_without_random(A,1);
         print(P);
}

proc primary_invariants_random (list #)
"USAGE:   primary_invariants_random(G1,G2,...,r[,flags]);
         G1,G2,...: <matrices> generating a finite matrix group, r: an <int>
         where -|r| to |r| is the range of coefficients of the random
         combinations of bases elements, flags: an optional <intvec> with three
         entries, if the first one equals 0 (also the default), the programme
         attempts to compute the Molien series and Reynolds operator, if it
         equals 1, the programme is told that the Molien series should not be
         computed, if it equals -1 characteristic 0 is simulated, i.e. the
         Molien series is computed as if the base field were characteristic 0
         (the user must choose a field of large prime characteristic, e.g.
         32003) and if the first one is anything else, it means that the
         characteristic of the base field divides the group order, the second
         component should give the size of intervals between canceling common
         factors in the expansion of the Molien series, 0 (the default) means
         only once after generating all terms, in prime characteristic also a
         negative number can be given to indicate that common factors should
         always be canceled when the expansion is simple (the root of the
         extension field does not occur among the coefficients)
DISPLAY: information about the various stages of the programme if the third
         flag does not equal 0
RETURN:  primary invariants (type <matrix>) of the invariant ring and if
         computable Reynolds operator (type <matrix>) and Molien series (type
         <matrix>), if the first flag is 1 and we are in the non-modular case
         then an <intvec> is returned giving some of the degrees where no
         non-trivial homogeneous invariants can be found
EXAMPLE: example primary_invariants_random; shows an example
THEORY:  Bases of homogeneous invariants are generated successively and random
         linear combinations are chosen as primary invariants that lower the
         dimension of the ideal generated by the previously found invariants
         (see paper \"Generating a Noetherian Normalization of the Invariant Ring
         of a Finite Group\" by Decker, Heydtmann, Schreyer (1997) to appear in
         JSC).
"
{
 // ----------------- checking input and setting flags ------------------------
  if (size(#)<2)
  { "ERROR:   There are too few parameters.";
    return();
  }
  int ch=char(basering);               // the algorithms depend very much on the
                                       // characteristic of the ground field
  int n=nvars(basering);               // 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
  int gen_num;
  int mol_flag,v;
  if (typeof(#[size(#)])=="intvec" && typeof(#[size(#)-1])=="int")
  { if (size(#[size(#)])<>3)
    { "ERROR:   <intvec> should have three entries.";
      return();
    }
    gen_num=size(#)-2;
    mol_flag=#[size(#)][1];
    if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0)))
    { "ERROR:   the second component of <intvec> should be >=0";
      return();
    }
    int interval=#[size(#)][2];
    v=#[size(#)][3];
    int max=#[size(#)-1];
    if (gen_num==0)
    { "ERROR:   There are no generators of a finite matrix group given.";
      return();
    }
  }
  else
  { if (typeof(#[size(#)])=="int")
    { gen_num=size(#)-1;
      mol_flag=0;
      int interval=0;
      v=0;
      int max=#[size(#)];
    }
    else
    { "ERROR:   If the two last parameters are not <int> and <intvec>, the last";
      "         parameter should be an <int>.";
      return();
    }
  }
  for (int i=1;i<=gen_num;i=i+1)
  { if (typeof(#[i])=="matrix")
    { if (nrows(#[i])<>n or ncols(#[i])<>n)
      { "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 first parameters should be a list of matrices";
      return();
    }
  }
 //----------------------------------------------------------------------------
  if (mol_flag==0)
  { if (ch==0)
    { matrix REY,M=reynolds_molien(#[1..gen_num],intvec(0,interval,v));
                                       // one will contain Reynolds operator and
                                       // the other enumerator and denominator
                                       // of Molien series
      matrix P=primary_char0_random(REY,M,max,v);
      return(P,REY,M);
    }
    else
    { list L=group_reynolds(#[1..gen_num],v);
      if (L[1]<>0)                     // testing whether we are in the modular
      { string newring="aksldfalkdsflkj"; // case
        if (minpoly==0)
        { if (v)
          { "  We are dealing with the non-modular case.";
          }
          molien(L[2..size(L)],newring,intvec(0,interval,v));
          matrix P=primary_charp_random(L[1],newring,max,v);
          return(P,L[1],newring);
        }
        else
        { if (v)
          { "  Since it is impossible for this programme to calculate the Molien series for";
            "  invariant rings over extension fields of prime characteristic, we have to";
            "  continue without it.";
            "";

          }
          list l=primary_charp_no_molien_random(L[1],max,v);
          if (size(l)==2)
          { return(l[1],L[1],l[2]);
          }
          else
          { return(l[1],L[1]);
          }
        }
      }
      else                             // the modular case
      { if (v)
        { "  There is also no Molien series, we can make use of...";
          "";
          "  We can start looking for primary invariants...";
          "";
        }
        return(primary_charp_without_random(#[1..gen_num],max,v));
      }
    }
  }
  if (mol_flag==1)                     // the user wants no calculation of the
  { list L=group_reynolds(#[1..gen_num],v); // Molien series
    if (ch==0)
    { list l=primary_char0_no_molien_random(L[1],max,v);
      if (size(l)==2)
      { return(l[1],L[1],l[2]);
      }
      else
      { return(l[1],L[1]);
      }
    }
    else
    { if (L[1]<>0)                     // testing whether we are in the modular
      { list l=primary_charp_no_molien_random(L[1],max,v); // case
        if (size(l)==2)
        { return(l[1],L[1],l[2]);
        }
        else
        { return(l[1],L[1]);
        }
      }
      else                             // the modular case
      { if (v)
        { "  We can start looking for primary invariants...";
          "";
        }
        return(primary_charp_without_random(#[1..gen_num],max,v));
      }
    }
  }
  if (mol_flag==-1)
  { if (ch==0)
    { "ERROR:   Characteristic 0 can only be simulated in characteristic p>>0.";
      return();
    }
    list L=group_reynolds(#[1..gen_num],v);
    string newring="aksldfalkdsflkj";
    molien(L[2..size(L)],newring,intvec(1,interval,v));
    matrix P=primary_charp_random(L[1],newring,max,v);
    return(P,L[1],newring);
  }
  else                                 // the user specified that the
  { if (ch==0)                         // characteristic divides the group order
    { "ERROR:   The characteristic cannot divide the group order when it is 0.";
      return();
    }
    if (v)
    { "";
    }
    return(primary_charp_without_random(#[1..gen_num],max,v));
  }
}
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;
         list L=primary_invariants_random(A,1);
         print(L[1]);
}

proc concat_intmat(intmat A,intmat B)
{ int n=nrows(A);
  int m1=ncols(A);
  int m2=ncols(B);
  intmat C[n][m1+m2];
  C[1..n,1..m1]=A[1..n,1..m1];
  C[1..n,m1+1..m1+m2]=B[1..n,1..m2];
  return(C);
}

proc power_products(intvec deg_vec,int d)
"USAGE:   power_products(dv,d);
         dv: an <intvec> giving the degrees of homogeneous polynomials, d: the
         degree of the desired power products
RETURN:  a size(dv)*m <intmat> where each column ought to be interpreted as
         containing the exponents of the corresponding polynomials. The product
         of the powers is then homogeneous of degree d.
EXAMPLE: example power_products; gives an example
"
{ ring R=0,x,dp;
  if (d<=0)
  { "ERROR:   The <int> may not be <= 0";
    return();
  }
  int d_neu,j,nc;
  int s=size(deg_vec);
  intmat PP[s][1];
  intmat TEST[s][1];
  for (int i=1;i<=s;i=i+1)
  { if (i<0)
    { "ERROR:   The entries of <intvec> may not be <= 0";
      return();
    }
    d_neu=d-deg_vec[i];
    if (d_neu>0)
    { intmat PPd_neu=power_products(intvec(deg_vec[i..s]),d_neu);
      if (size(ideal(PPd_neu))<>0)
      { nc=ncols(PPd_neu);
        intmat PPd_neu_gross[s][nc];
        PPd_neu_gross[i..s,1..nc]=PPd_neu[1..s-i+1,1..nc];
        for (j=1;j<=nc;j=j+1)
        { PPd_neu_gross[i,j]=PPd_neu_gross[i,j]+1;
        }
        PP=concat_intmat(PP,PPd_neu_gross);
        kill PPd_neu_gross;
      }
      kill PPd_neu;
    }
    if (d_neu==0)
    { intmat PPd_neu[s][1];
      PPd_neu[i,1]=1;
      PP=concat_intmat(PP,PPd_neu);
      kill PPd_neu;
    }
  }
  if (matrix(PP)<>matrix(TEST))
  { PP=compress(PP);
  }
  return(PP);
}
example
{ echo=2;
         intvec dv=5,5,5,10,10;
         print(power_products(dv,10));
         print(power_products(dv,7));
}

proc secondary_char0 (matrix P, matrix REY, matrix M, list #)
"USAGE:   secondary_char0(P,REY,M[,v]);
         P: a 1xn <matrix> with primary invariants, REY: a gxn <matrix>
         representing the Reynolds operator, M: a 1x2 <matrix> giving enumerator
         and denominator of the Molien series, v: an optional <int>
ASSUME:  n is the number of variables of the basering, g the size of the group,
         REY is the 1st return value of group_reynolds(), reynolds_molien() or
         the second one of primary_invariants(), M the return value of molien()
         or the second one of reynolds_molien() or the third one of
         primary_invariants()
RETURN:  secondary invariants of the invariant ring (type <matrix>) and
         irreducible secondary invariants (type <matrix>)
DISPLAY: information if v does not equal 0
EXAMPLE: example secondary_char0; shows an example
THEORY:  The secondary invariants are calculated by finding a basis (in terms of
         monomials) of the basering modulo the primary invariants, mapping those
         to invariants with the Reynolds operator and using these images or
         their power products such that they are linearly independent modulo the
         primary invariants (see paper \"Some Algorithms in Invariant Theory of
         Finite Groups\" by Kemper and Steel (1997)).
"
{ def br=basering;
  degBound=0;
 //----------------- checking input and setting verbose mode ------------------
  if (char(br)<>0)
  { "ERROR:   secondary_char0 should only be used with rings of characteristic 0.";
    return();
  }
  int i;
  if (size(#)>0)
  { if (typeof(#[size(#)])=="int")
    { int v=#[size(#)];
    }
    else
    { int v=0;
    }
  }
  else
  { int v=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 (ncols(P)<>n)
  { "ERROR:   The first parameter ought to be the matrix of the primary";
    "         invariants."
    return();
  }
  if (ncols(REY)<>n)
  { "ERROR:   The second parameter ought to be the Reynolds operator."
    return();
  }
  if (ncols(M)<>2 or nrows(M)<>1)
  { "ERROR:   The third parameter ought to be the Molien series."
    return();
  }
  if (v && voice==2)
  { "";
  }
  int j, m, counter;
 //- finding the polynomial giving number and degrees of secondary invariants -
  poly p=1;
  for (j=1;j<=n;j=j+1)                 // calculating the denominator of the
  { p=p*(1-var(1)^deg(P[j]));          // Hilbert series of the ring generated
  }                                    // by the primary invariants -
  matrix s[1][2]=M[1,1]*p,M[1,2];      // s is used for canceling
  s=matrix(syz(ideal(s)));
  p=s[2,1];                            // the polynomial telling us where to
                                       // search for secondary invariants
  map slead=br,ideal(0);
  p=1/slead(p)*p;                      // smallest term of p needs to be 1
  if (v)
  { "  Polynomial telling us where to look for secondary invariants:";
    "   "+string(p);
    "";
  }
  matrix dimmat=coeffs(p,var(1));      // 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
  if (v)
  { "  In degree 0 we have: 1";
    "";
  }
 //-------------------------- initializing variables --------------------------
  intmat PP;
  poly pp;
  int k;
  intvec deg_vec;
  ideal sP=std(ideal(P));
  ideal TEST,B,IS;
  ideal S=1;                           // 1 is the first secondary invariant -
 //--------------------- generating secondary invariants ----------------------
  for (i=2;i<=m;i=i+1)                 // going through dimmat -
  { if (int(dimmat[i,1])<>0)           // when it is == 0 we need to find 0
    {                                  // elements in the current degree (i-1)
      if (v)
      { "  Searching in degree "+string(i-1)+", we need to find "+string(int(dimmat[i,1]))+" invariant(s)...";
      }
      TEST=sP;
      counter=0;                       // we'll count up to degvec[i]
      if (IS[1]<>0)
      { PP=power_products(deg_vec,i-1); // finding power products of irreducible
      }                                // secondary invariants
      if (size(ideal(PP))<>0)
      { for (j=1;j<=ncols(PP);j=j+1)   // going through all the power products
        { pp=1;
          for (k=1;k<=nrows(PP);k=k+1)
          { pp=pp*IS[1,k]^PP[k,j];
          }
          if (reduce(pp,TEST)<>0)
          { S=S,pp;
            counter=counter+1;
            if (v)
            { "  We find: "+string(pp);
            }
            if (int(dimmat[i,1])<>counter)
            { TEST=std(TEST+ideal(NF(pp,TEST))); // should be replaced by next
                                                 // line soon
            // TEST=std(TEST,NF(pp,TEST));
            }
            else
            { break;
            }
          }
        }
      }
      if (int(dimmat[i,1])<>counter)
      { B=sort_of_invariant_basis(sP,REY,i-1,int(dimmat[i,1])*6); // B contains
                                       // images of kbase(sP,i-1) under the
                                       // Reynolds operator that are linearly
                                       // independent and that don't reduce to
                                       // 0 modulo sP -
        if (counter==0 && ncols(B)==int(dimmat[i,1])) // then we can take all of
        { S=S,B;                       // B
          IS=IS+B;
          if (deg_vec[1]==0)
          { deg_vec=i-1;
            if (v)
            { "  We find: "+string(B[1]);
            }
            for (j=2;j<=int(dimmat[i,1]);j=j+1)
            { deg_vec=deg_vec,i-1;
              if (v)
              { "  We find: "+string(B[j]);
              }
            }
          }
          else
          { for (j=1;j<=int(dimmat[i,1]);j=j+1)
            { deg_vec=deg_vec,i-1;
              if (v)
              { "  We find: "+string(B[j]);
              }
            }
          }
        }
        else
        { j=0;                         // j goes through all of B -
          while (int(dimmat[i,1])<>counter) // need to find dimmat[i,1]
          {                            // invariants that are linearly
                                       // independent modulo TEST
            j=j+1;
            if (reduce(B[j],TEST)<>0)  // B[j] should be added
            { S=S,B[j];
              IS=IS+ideal(B[j]);
              if (deg_vec[1]==0)
              { deg_vec[1]=i-1;
              }
              else
              { deg_vec=deg_vec,i-1;
              }
              counter=counter+1;
              if (v)
              { "  We find: "+string(B[j]);
              }
              if (int(dimmat[i,1])<>counter)
              { TEST=std(TEST+ideal(NF(B[j],TEST))); // should be replaced by
                                                     // next line
              // TEST=std(TEST,NF(B[j],TEST));
              }
            }
          }
        }
      }
      if (v)
      { "";
      }
    }
  }
  if (v)
  { "  We're done!";
    "";
  }
  return(matrix(S),matrix(IS));
}
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;
         list L=primary_invariants(A);
         matrix S,IS=secondary_char0(L[1..3]);
         print(S);
         print(IS);
}

proc secondary_charp (matrix P, matrix REY, string ring_name, list #)
"USAGE:   secondary_charp(P,REY,ringname[,v]);
         P: a 1xn <matrix> with primary invariants, REY: a gxn <matrix>
         representing the Reynolds operator, ringname: a <string> giving the
         name of a ring of characteristic 0 where the Molien series is stored,
         v: an optional <int>
ASSUME:  n is the number of variables of the basering, g the size of the group,
         REY is the 1st return value of group_reynolds(), reynolds_molien() or
         the second one of primary_invariants(), `ringname` is ring of
         characteristic 0 that has been created by molien() or reynolds_molien()
         or primary_invariants()
RETURN:  secondary invariants of the invariant ring (type <matrix>) and
         irreducible secondary invariants (type <matrix>)
DISPLAY: information if v does not equal 0
EXAMPLE: example secondary_charp; shows an example
THEORY:  The secondary invariants are calculated by finding a basis (in terms of
         monomials) of the basering modulo the primary invariants, mapping those
         to invariants with the Reynolds operator and using these images or
         their power products such that they are linearly independent modulo the
         primary invariants (see paper \"Some Algorithms in Invariant Theory of
         Finite Groups\" by Kemper and Steel (1997)).
"
{ def br=basering;
  degBound=0;
 //---------------- checking input and setting verbose mode -------------------
  if (char(br)==0)
  { "ERROR:   secondary_charp should only be used with rings of characteristic p>0.";
    return();
  }
  int i;
  if (size(#)>0)
  { if (typeof(#[size(#)])=="int")
    { int v=#[size(#)];
    }
    else
    { int v=0;
    }
  }
  else
  { int v=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 (ncols(P)<>n)
  { "ERROR:   The first parameter ought to be the matrix of the primary";
    "         invariants."
    return();
  }
  if (ncols(REY)<>n)
  { "ERROR:   The second parameter ought to be the Reynolds operator."
    return();
  }
  if (typeof(`ring_name`)<>"ring")
  { "ERROR:   The <string> should give the name of the ring where the Molien."
    "         series is stored.";
    return();
  }
  if (v && voice==2)
  { "";
  }
  int j, m, counter, d;
  intvec deg_dim_vec;
 //- finding the polynomial giving number and degrees of secondary invariants -
  for (j=1;j<=n;j=j+1)
  { deg_dim_vec[j]=deg(P[j]);
  }
  setring `ring_name`;
  poly p=1;
  for (j=1;j<=n;j=j+1)                 // calculating the denominator of the
  { p=p*(1-var(1)^deg_dim_vec[j]);     // Hilbert series of the ring generated
  }                                    // by the primary invariants -
  matrix s[1][2]=M[1,1]*p,M[1,2];      // s is used for canceling
  s=matrix(syz(ideal(s)));
  p=s[2,1];                            // the polynomial telling us where to
                                       // search for secondary invariants
  map slead=basering,ideal(0);
  p=1/slead(p)*p;                      // smallest term of p needs to be 1
  if (v)
  { "  Polynomial telling us where to look for secondary invariants:";
    "   "+string(p);
    "";
  }
  matrix dimmat=coeffs(p,var(1));      // 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
  deg_dim_vec=1;
  for (j=2;j<=m;j=j+1)
  { deg_dim_vec=deg_dim_vec,int(dimmat[j,1]);
  }
  if (v)
  { "  In degree 0 we have: 1";
    "";
  }
 //------------------------ initializing variables ----------------------------
  setring br;
  intmat PP;
  poly pp;
  int k;
  intvec deg_vec;
  ideal sP=std(ideal(P));
  ideal TEST,B,IS;
  ideal S=1;                           // 1 is the first secondary invariant
 //------------------- generating secondary invariants ------------------------
  for (i=2;i<=m;i=i+1)                 // going through deg_dim_vec -
  { if (deg_dim_vec[i]<>0)             // when it is == 0 we need to find 0
    {                                  // elements in the current degree (i-1)
      if (v)
      { "  Searching in degree "+string(i-1)+", we need to find "+string(deg_dim_vec[i])+" invariant(s)...";
      }
      TEST=sP;
      counter=0;                       // we'll count up to degvec[i]
      if (IS[1]<>0)
      { PP=power_products(deg_vec,i-1); // generating power products of
      }                                // irreducible secondary invariants
      if (size(ideal(PP))<>0)
      { for (j=1;j<=ncols(PP);j=j+1)   // going through all of those
        { pp=1;
          for (k=1;k<=nrows(PP);k=k+1)
          { pp=pp*IS[1,k]^PP[k,j];
          }
          if (reduce(pp,TEST)<>0)
          { S=S,pp;
            counter=counter+1;
            if (v)
            { "  We find: "+string(pp);
            }
            if (deg_dim_vec[i]<>counter)
            { TEST=std(TEST+ideal(NF(pp,TEST))); // should be soon replaced by
                                                 // next line
            // TEST=std(TEST,NF(pp,TEST));
            }
            else
            { break;
            }
          }
        }
      }
      if (deg_dim_vec[i]<>counter)
      { B=sort_of_invariant_basis(sP,REY,i-1,deg_dim_vec[i]*6); // B contains
                                       // images of kbase(sP,i-1) under the
                                       // Reynolds operator that are linearly
                                       // independent and that don't reduce to
                                       // 0 modulo sP -
        if (counter==0 && ncols(B)==deg_dim_vec[i]) // then we can add all of B
        { S=S,B;
          IS=IS+B;
          if (deg_vec[1]==0)
          { deg_vec=i-1;
            if (v)
            { "  We find: "+string(B[1]);
            }
            for (j=2;j<=deg_dim_vec[i];j=j+1)
            { deg_vec=deg_vec,i-1;
              if (v)
              { "  We find: "+string(B[j]);
              }
            }
          }
          else
          { for (j=1;j<=deg_dim_vec[i];j=j+1)
            { deg_vec=deg_vec,i-1;
              if (v)
              { "  We find: "+string(B[j]);
              }
            }
          }
        }
        else
        { j=0;                         // j goes through all of B -
          while (deg_dim_vec[i]<>counter) // need to find deg_dim_vec[i]
          {                            // invariants that are linearly
                                       // independent modulo TEST
            j=j+1;
            if (reduce(B[j],TEST)<>0)   // B[j] should be added
            { S=S,B[j];
              IS=IS+ideal(B[j]);
              if (deg_vec[1]==0)
              { deg_vec[1]=i-1;
              }
              else
              { deg_vec=deg_vec,i-1;
              }
              counter=counter+1;
              if (v)
              { "  We find: "+string(B[j]);
              }
              if (deg_dim_vec[i]<>counter)
              { TEST=std(TEST+ideal(NF(B[j],TEST))); // should be soon replaced
                                                     // by next line
              // TEST=std(TEST,NF(B[j],TEST));
              }
            }
          }
        }
      }
      if (v)
      { "";
      }
    }
  }
  if (v)
  { "  We're done!";
    "";
  }
  if (ring_name=="aksldfalkdsflkj")
  { kill `ring_name`;
  }
  return(matrix(S),matrix(IS));
}
example
{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed into";
  "         characteristic 3)";
  echo=2;
         ring R=3,(x,y,z),dp;
         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
         list L=primary_invariants(A);
         matrix S,IS=secondary_charp(L[1..size(L)]);
         print(S);
         print(IS);
}

proc secondary_no_molien (matrix P, matrix REY, list #)
"USAGE:   secondary_no_molien(P,REY[,deg_vec,v]);
         P: a 1xn <matrix> with primary invariants, REY: a gxn <matrix>
         representing the Reynolds operator, deg_vec: an optional <intvec>
         listing some degrees where no non-trivial homogeneous invariants can be
         found, v: an optional <int>
ASSUME:  n is the number of variables of the basering, g the size of the group,
         REY is the 1st return value of group_reynolds(), reynolds_molien() or
         the second one of primary_invariants(), deg_vec is the second return
         value of primary_char0_no_molien(), primary_charp_no_molien(),
         primary_char0_no_molien_random() or primary_charp_no_molien_random()
RETURN:  secondary invariants of the invariant ring (type <matrix>)
DISPLAY: information if v does not equal 0
EXAMPLE: example secondary_no_molien; shows an example
THEORY:  The secondary invariants are calculated by finding a basis (in terms of
         monomials) of the basering modulo the primary invariants, mapping those
         to invariants with the Reynolds operator and using these images as
         candidates for secondary invariants.
"
{ int i;
  degBound=0;
 //------------------ checking input and setting verbose ----------------------
  if (size(#)==1 or size(#)==2)
  { if (typeof(#[size(#)])=="int")
    { if (size(#)==2)
      { if (typeof(#[size(#)-1])=="intvec")
        { intvec deg_vec=#[size(#)-1];
        }
        else
        { "ERROR:   the third parameter should be an <intvec>";
          return();
        }
      }
      int v=#[size(#)];
    }
    else
    { if (size(#)==1)
      { if (typeof(#[size(#)])=="intvec")
        { intvec deg_vec=#[size(#)];
          int v=0;
        }
        else
        { "ERROR:   the third parameter should be an <intvec>";
          return();
        }
      }
      else
      { "ERROR:   wrong list of parameters";
        return();
      }
    }
  }
  else
  { if (size(#)>2)
    { "ERROR:   there are too many parameters";
      return();
    }
    int v=0;
  }
  int n=nvars(basering);               // 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 (ncols(P)<>n)
  { "ERROR:   The first parameter ought to be the matrix of the primary";
    "         invariants."
    return();
  }
  if (ncols(REY)<>n)
  { "ERROR:   The second parameter ought to be the Reynolds operator."
    return();
  }
  if (v && voice==2)
  { "";
  }
  int j, m, d;
  int max=1;
  for (j=1;j<=n;j=j+1)
  { max=max*deg(P[j]);
  }
  max=max/nrows(REY);
  if (v)
  { "  We need to find "+string(max)+" secondary invariants.";
    "";
    "  In degree 0 we have: 1";
    "";
  }
 //------------------------- initializing variables ---------------------------
  ideal sP=std(ideal(P));
  ideal B, TEST;
  ideal S=1;                           // 1 is the first secondary invariant
  int counter=1;
  i=0;
  if (defined(deg_vec)<>voice)
  { intvec deg_vec;
  }
  int k=1;
 //--------------------- generating secondary invariants ----------------------
  while (counter<>max)
  { i=i+1;
    if (deg_vec[k]<>i)
    { if (v)
      { "  Searching in degree "+string(i)+"...";
      }
      B=sort_of_invariant_basis(sP,REY,i,max); // B contains images of
                                       // kbase(sP,i) under the Reynolds
                                       // operator that are linearly independent
                                       // and that don't reduce to 0 modulo sP
      TEST=sP;
      for (j=1;j<=ncols(B);j=j+1)
      {                                // that are linearly independent modulo
                                       // TEST
        if (reduce(B[j],TEST)<>0)      // B[j] should be added
        { S=S,B[j];
          counter=counter+1;
          if (v)
          { "  We find: "+string(B[j]);
          }
          if (counter==max)
          { break;
          }
          else
          { if (j<>ncols(B))
            { TEST=std(TEST+ideal(NF(B[j],TEST))); // should soon be replaced by
                                                   // next line
            // TEST=std(TEST,NF(B[j],TEST));
            }
          }
        }
      }
    }
    else
    { if (size(deg_vec)==k)
      { k=1;
      }
      else
      { k=k+1;
      }
    }
  }
  if (v)
  { "";
  }
  if (v)
  { "  We're done!";
    "";
  }
  return(matrix(S));
}
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;
         list L=primary_invariants(A,intvec(1,1,0));
         matrix S=secondary_no_molien(L[1..3]);
         print(S);
}

proc secondary_with_irreducible_ones_no_molien (matrix P, matrix REY, list #)
"USAGE:   secondary_with_irreducible_ones_no_molien(P,REY[,v]);
         P: a 1xn <matrix> with primary invariants, REY: a gxn <matrix>
         representing the Reynolds operator, v: an optional <int>
ASSUME:  n is the number of variables of the basering, g the size of the group,
         REY is the 1st return value of group_reynolds(), reynolds_molien() or
         the second one of primary_invariants()
RETURN:  secondary invariants of the invariant ring (type <matrix>) and
         irreducible secondary invariants (type <matrix>)
DISPLAY: information if v does not equal 0
EXAMPLE: example secondary_with_irreducible_ones_no_molien; shows an example
THEORY:  The secondary invariants are calculated by finding a basis (in terms of
         monomials) of the basering modulo the primary invariants, mapping those
         to invariants with the Reynolds operator and using these images or
         their power products such that they are linearly independent modulo the
         primary invariants (see paper \"Some Algorithms in Invariant Theory of
         Finite Groups\" by Kemper and Steel (1997)).
"
{ int i;
  degBound=0;
 //--------------------- checking input and setting verbose mode --------------
  if (size(#)==1 or size(#)==2)
  { if (typeof(#[size(#)])=="int")
    { if (size(#)==2)
      { if (typeof(#[size(#)-1])=="intvec")
        { intvec deg_vec=#[size(#)-1];
        }
        else
        { "ERROR:   the third parameter should be an <intvec>";
          return();
        }
      }
      int v=#[size(#)];
    }
    else
    { if (size(#)==1)
      { if (typeof(#[size(#)])=="intvec")
        { intvec deg_vec=#[size(#)];
          int v=0;
        }
        else
        { "ERROR:   the third parameter should be an <intvec>";
          return();
        }
      }
      else
      { "ERROR:   wrong list of parameters";
        return();
      }
    }
  }
  else
  { if (size(#)>2)
    { "ERROR:   there are too many parameters";
      return();
    }
    int v=0;
  }
  int n=nvars(basering);               // 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 (ncols(P)<>n)
  { "ERROR:   The first parameter ought to be the matrix of the primary";
    "         invariants."
    return();
  }
  if (ncols(REY)<>n)
  { "ERROR:   The second parameter ought to be the Reynolds operator."
    return();
  }
  if (v && voice==2)
  { "";
  }
  int j, m, d;
  int max=1;
  for (j=1;j<=n;j=j+1)
  { max=max*deg(P[j]);
  }
  max=max/nrows(REY);
  if (v)
  { "  We need to find "+string(max)+" secondary invariants.";
    "";
    "  In degree 0 we have: 1";
    "";
  }
 //------------------------ initializing variables ----------------------------
  intmat PP;
  poly pp;
  int k;
  intvec irreducible_deg_vec;
  ideal sP=std(ideal(P));
  ideal B,TEST,IS;
  ideal S=1;                           // 1 is the first secondary invariant
  int counter=1;
  i=0;
  if (defined(deg_vec)<>voice)
  { intvec deg_vec;
  }
  int l=1;
 //------------------- generating secondary invariants ------------------------
  while (counter<>max)
  { i=i+1;
    if (deg_vec[l]<>i)
    { if (v)
      { "  Searching in degree "+string(i)+"...";
      }
      TEST=sP;
      if (IS[1]<>0)
      { PP=power_products(irreducible_deg_vec,i);  // generating all power
      }                                // products of irreducible secondary
                                       // invariants
      if (size(ideal(PP))<>0)
      { for (j=1;j<=ncols(PP);j=j+1)   // going through all those power products
        { pp=1;
          for (k=1;k<=nrows(PP);k=k+1)
          { pp=pp*IS[1,k]^PP[k,j];
          }
          if (reduce(pp,TEST)<>0)
          { S=S,pp;
            counter=counter+1;
            if (v)
            { "  We find: "+string(pp);
            }
            if (counter<>max)
            { TEST=std(TEST+ideal(NF(pp,TEST))); // should soon be replaced by
                                                 // next line
            // TEST=std(TEST,NF(pp,TEST));
            }
            else
            { break;
            }
          }
        }
      }
      if (max<>counter)
      { B=sort_of_invariant_basis(sP,REY,i,max); // B contains images of
                                       // kbase(sP,i) under the Reynolds
                                       // operator that are linearly independent
                                       // and that don't reduce to 0 modulo sP
        for (j=1;j<=ncols(B);j=j+1)
        { if (reduce(B[j],TEST)<>0)    // B[j] should be added
          { S=S,B[j];
            IS=IS+ideal(B[j]);
            if (irreducible_deg_vec[1]==0)
            { irreducible_deg_vec[1]=i;
            }
            else
            { irreducible_deg_vec=irreducible_deg_vec,i;
            }
            counter=counter+1;
            if (v)
            { "  We find: "+string(B[j]);
            }
            if (counter==max)
            { break;
            }
            else
            { if (j<>ncols(B))
              { TEST=std(TEST+ideal(NF(B[j],TEST))); // should soon be replaced
                                                     // by next line
              // TEST=std(TEST,NF(B[j],TEST));
              }
            }
          }
        }
      }
    }
    else
    { if (size(deg_vec)==l)
      { l=1;
      }
      else
      { l=l+1;
      }
    }
  }
  if (v)
  { "";
  }
  if (v)
  { "  We're done!";
    "";
  }
  return(matrix(S),matrix(IS));
}
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;
         list L=primary_invariants(A,intvec(1,1,0));
         matrix S,IS=secondary_with_irreducible_ones_no_molien(L[1..2]);
         print(S);
         print(IS);
}

proc secondary_not_cohen_macaulay (matrix P, list #)
"USAGE:   secondary_not_cohen_macaulay(P,G1,G2,...[,v]);
         P: a 1xn <matrix> with primary invariants, G1,G2,...: nxn <matrices>
         generating a finite matrix group, v: an optional <int>
ASSUME:  n is the number of variables of the basering
RETURN:  secondary invariants of the invariant ring (type <matrix>)
DISPLAY: information if v does not equal 0
EXAMPLE: example secondary_not_cohen_macaulay; shows an example
THEORY:  The secondary invariants are generated following \"Generating Invariant
         Rings of Finite Groups over Arbitrary Fields\" by Kemper (1996, to
         appear in JSC).
"
{ int i, j;
  degBound=0;
  def br=basering;
  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 (size(#)>0)                       // checking input and setting verbose
  { if (typeof(#[size(#)])=="int")
    { int gen_num=size(#)-1;
      if (gen_num==0)
      { "ERROR:   There are no generators of the finite matrix group given.";
        return();
      }
      int v=#[size(#)];
      for (i=1;i<=gen_num;i=i+1)
      { if (typeof(#[i])<>"matrix")
        { "ERROR:   These parameters should be generators of the finite matrix group.";
          return();
        }
        if ((n<>nrows(#[i])) or (n<>ncols(#[i])))
        { "ERROR:   matrices need to be square and of the same dimensions";
          return();
        }
      }
    }
    else
    { int v=0;
      int gen_num=size(#);
      for (i=1;i<=gen_num;i=i+1)
      { if (typeof(#[i])<>"matrix")
        { "ERROR:   These parameters should be generators of the finite matrix group.";
          return();
        }
        if ((n<>nrows(#[i])) or (n<>ncols(#[i])))
        { "ERROR:   matrices need to be square and of the same dimensions";
          return();
        }
      }
    }
  }
  else
  { "ERROR:   There are no generators of the finite matrix group given.";
    return();
  }
  if (ncols(P)<>n)
  { "ERROR:   The first parameter ought to be the matrix of the primary";
    "         invariants."
    return();
  }
  if (v && voice==2)
  { "";
  }
  ring alskdfalkdsj=0,x,dp;
  matrix M[1][2]=1,(1-x)^n;            // we look at our primary invariants as
  export alskdfalkdsj;
  export M;
  setring br;                          // such of the subgroup that only
  matrix REY=matrix(maxideal(1));      // 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 - M[1,1]/M[1,2] 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
  if (v)
  { "  The procedure secondary_charp() is called to calculate secondary invariants";
    "  of the invariant ring of the trivial group with respect to the primary";
    "  invariants found previously.";
    "";
  }
  matrix trivialS=secondary_charp(P,REY,"alskdfalkdsj",v);
  kill alskdfalkdsj;
  // now we have those secondary invariants
  int k=ncols(trivialS);               // k is the number of the secondary
                                       // invariants, we just calculated
  if (v)
  { "  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)[gen_num][k];             // M(1) will contain a module
  ideal B;
  for (i=1;i<=gen_num;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(trivialS)-trivialS;            // 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)=nres(M(1),2)[2];
  int 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 P=f(P);                       // 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)-P[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;
  ideal substitute=maxideal(1),ideal(P),1;
  f=R,substitute;                      // replacing y(1..n) by primary
                                       // invariants -
  M(2)=f(M);                           // M(2) is the new module
  m=ncols(M(2));
  matrix S[1][m];
  S=matrix(trivialS)*matrix(M(2));     // S now contains the secondary
                                       // invariants
  for (i=1; i<=m;i=i+1)
  { S[1,i]=S[1,i]/leadcoef(S[1,i]); // making elements nice
  }
  S=sort(ideal(S))[1];
  if (v)
  { "  These are the secondary invariants: ";
    for (i=1;i<=m;i=i+1)
    { "   "+string(S[1,i]);
    }
    "";
    "  We're done!";
    "";
  }
  if ((v or (voice==2)) && (m>1))
  { "  WARNING: The invariant ring might not have a Hironaka decomposition";
    "           if the characteristic of the coefficient field divides the";
    "           group order.";
  }
  return(S);
}
example
{ echo=2;
           ring R=2,(x,y,z),dp;
           matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
           list L=primary_invariants(A);
           matrix S=secondary_not_cohen_macaulay(L[1],A);
           print(S);
}

proc invariant_ring (list #)
"USAGE:   invariant_ring(G1,G2,...[,flags]);
         G1,G2,...: <matrices> generating a finite matrix group, flags: an
         optional <intvec> with three entries: 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 Molien series should not
         be computed, if it equals -1 characteristic 0 is simulated, i.e. the
         Molien series is computed as if the base field were characteristic 0
         (the user must choose a field of large prime characteristic, e.g.
         32003) and if the first one is anything else, it means that the
         characteristic of the base field divides the group order (i.e. it will
         not even be attempted to compute the Reynolds operator or Molien
         series), the second component should give the size of intervals
         between canceling common factors in the expansion of the Molien series,
         0 (the default) means only once after generating all terms, in prime
         characteristic also a negative number can be given to indicate that
         common factors should always be canceled when the expansion is simple
         (the root of the extension field does not occur among the coefficients)
RETURN:  primary and secondary invariants (both of type <matrix>) generating the
         invariant ring with respect to the matrix group generated by the
         matrices in the input and irreducible secondary invariants (type
         <matrix>) if the Molien series was available
DISPLAY: information about the various stages of the program if the third flag
         does not equal 0
EXAMPLE: example invariant_ring; shows an example
THEORY:  Bases of homogeneous invariants are generated successively and those
         are chosen as primary invariants that lower the dimension of the ideal
         generated by the previously found invariants (see paper \"Generating a
         Noetherian Normalization of the Invariant Ring of a Finite Group\" by
         Decker, Heydtmann, Schreyer (1997) to appear in JSC). In the
         non-modular case secondary invariants are calculated by finding a
         basis (in terms of monomials) of the basering modulo the primary
         invariants, mapping to invariants with the Reynolds operator and using
         those or their power products such that they are linearly independent
         modulo the primary invariants (see paper \"Some Algorithms in Invariant
         Theory of Finite Groups\" by Kemper and Steel (1997)). In the modular
         case they are generated according to \"Generating Invariant Rings of
         Finite Groups over Arbitrary Fields\" by Kemper (1996, to appear in
         JSC).
"
{ if (size(#)==0)
  { "ERROR:   There are no generators given.";
    return();
  }
  int ch=char(basering);               // the algorithms depend very much on the
                                       // characteristic of the ground field -
  int n=nvars(basering);               // 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
  int gen_num;
  int mol_flag, v;
 //------------------- checking input and setting flags -----------------------
  if (typeof(#[size(#)])=="intvec")
  { if (size(#[size(#)])<>3)
    { "ERROR:   The <intvec> should have three entries.";
      return();
    }
    gen_num=size(#)-1;
    mol_flag=#[size(#)][1];
    if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0)))
    { "ERROR:   the second component of <intvec> should be >=0";
      return();
    }
    int interval=#[size(#)][2];
    v=#[size(#)][3];
  }
  else
  { gen_num=size(#);
    mol_flag=0;
    int interval=0;
    v=0;
  }
 //----------------------------------------------------------------------------
  if (mol_flag==0)                     // calculation Molien series will be
  { if (ch==0)                         // attempted -
    { matrix REY,M=reynolds_molien(#[1..gen_num],intvec(0,interval,v)); // one
                                       // will contain Reynolds operator and the
                                       // other enumerator and denominator of
                                       // Molien series
      matrix P=primary_char0(REY,M,v);
      matrix S,IS=secondary_char0(P,REY,M,v);
      return(P,S,IS);
    }
    else
    { list L=group_reynolds(#[1..gen_num],v);
      if (L[1]<>0)                     // testing whether we are in the modular
      { string newring="aksldfalkdsflkj"; // case
        if (minpoly==0)
        { if (v)
          { "  We are dealing with the non-modular case.";
          }
          molien(L[3..size(L)],newring,L[2],intvec(0,interval,v));
          matrix P=primary_charp(L[1],newring,v);
          matrix S,IS=secondary_charp(P,L[1],newring,v);
          if (defined(aksldfalkdsflkj)==2)
          { kill aksldfalkdsflkj;
          }
          return(P,S,IS);
        }
        else
        { if (v)
          { "  Since it is impossible for this programme to calculate the Molien
 series for";
            "  invariant rings over extension fields of prime characteristic, we
 have to";
            "  continue without it.";
            "";

          }
          list l=primary_charp_no_molien(L[1],v);
          if (size(l)==2)
          { matrix S=secondary_no_molien(l[1],L[1],l[2],v);
          }
          else
          { matix S=secondary_no_molien(l[1],L[1],v);
          }
          return(l[1],S);
        }
      }
      else                             // the modular case
      { if (v)
        { "  There is also no Molien series or Reynolds operator, we can make use of...";
          "";
          "  We can start looking for primary invariants...";
          "";
        }
        matrix P=primary_charp_without(#[1..gen_num],v);
        matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
        return(P,S);
      }
    }
  }
  if (mol_flag==1)                     // the user wants no calculation of the
  { list L=group_reynolds(#[1..gen_num],v); // Molien series
    if (ch==0)
    { list l=primary_char0_no_molien(L[1],v);
      if (size(l)==2)
      { matix S=secondary_no_molien(l[1],L[1],l[2],v);
      }
      else
      { matix S=secondary_no_molien(l[1],L[1],v);
      }
      return(l[1],S);
    }
    else
    { if (L[1]<>0)                     // testing whether we are in the modular
      { list l=primary_charp_no_molien(L[1],v); // case
        if (size(l)==2)
        { matix S=secondary_no_molien(l[1],L[1],l[2],v);
        }
        else
        { matix S=secondary_no_molien(l[1],L[1],v);
        }
        return(l[1],S);
      }
      else                             // the modular case
      { if (v)
        { "  We can start looking for primary invariants...";
          "";
        }
        matrix P=primary_charp_without(#[1..gen_num],v);
        matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
        return(L[1],S);
      }
    }
  }
  if (mol_flag==-1)
  { if (ch==0)
    { "ERROR:   Characteristic 0 can only be simulated in characteristic p>>0.
";
      return();
    }
    list L=group_reynolds(#[1..gen_num],v);
    string newring="aksldfalkdsflkj";
    molien(L[2..size(L)],newring,intvec(1,interval,v));
    matrix P=primary_charp(L[1],newring,v);
    matrix S,IS=secondary_charp(P,L[1],newring,v);
    kill aksldfalkdsflkj;
    return(P,S,IS);
  }
  else                                 // the user specified that the
  { if (ch==0)                         // characteristic divides the group order
    { "ERROR:   The characteristic cannot divide the group order when it is 0.
";
      return();
    }
    if (v)
    { "";
    }
    matrix P=primary_charp_without(#[1..gen_num],v);
    matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
    return(L[1],S);
  }
}
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 P,S,IS=invariant_ring(A);
         print(P);
         print(S);
         print(IS);
}

proc invariant_ring_random (list #)
"USAGE:   invariant_ring_random(G1,G2,...,r[,flags]);
         G1,G2,...: <matrices> generating a finite matrix group, r: an <int>
         where -|r| to |r| is the range of coefficients of random
         combinations of bases elements that serve as primary invariants,
         flags: an optional <intvec> with three entries: 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 Molien series
         should not be computed, if it equals -1 characteristic 0 is simulated,
         i.e. the Molien series is computed as if the base field were
         characteristic 0 (the user must choose a field of large prime
         characteristic, e.g.  32003) and if the first one is anything else, it
         means that the characteristic of the base field divides the group order
         (i.e. it will not even be attempted to compute the Reynolds operator or
         Molien series), the second component should give the size of intervals
         between canceling common factors in the expansion of the Molien series,
         0 (the default) means only once after generating all terms, in prime
         characteristic also a negative number can be given to indicate that
         common factors should always be canceled when the expansion is simple
         (the root of the extension field does not occur among the coefficients)
RETURN:  primary and secondary invariants (both of type <matrix>) generating the
         invariant ring with respect to the matrix group generated by the
         matrices in the input and irreducible secondary invariants (type
         <matrix>) if the Molien series was available
DISPLAY: information about the various stages of the program if the third flag
         does not equal 0
EXAMPLE: example invariant_ring_random; shows an example
THEORY:  is the same as for invariant_ring except that random combinations of
         basis elements are chosen as candidates for primary invariants and
         hopefully they lower the dimension of the previously found primary
         invariants by the right amount.
"
{ if (size(#)<2)
  { "ERROR:   There are too few parameters.";
    return();
  }
  int ch=char(basering);               // the algorithms depend very much on the
                                       // characteristic of the ground field
  int n=nvars(basering);               // 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
  int gen_num;
  int mol_flag, v;
 //------------------- checking input and setting flags -----------------------
  if (typeof(#[size(#)])=="intvec" && typeof(#[size(#)-1])=="int")
  { if (size(#[size(#)])<>3)
    { "ERROR:   <intvec> should have three entries.";
      return();
    }
    gen_num=size(#)-2;
    mol_flag=#[size(#)][1];
    if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0)))
    { "ERROR:   the second component of <intvec> should be >=0";
      return();
    }
    int interval=#[size(#)][2];
    v=#[size(#)][3];
    int max=#[size(#)-1];
    if (gen_num==0)
    { "ERROR:   There are no generators of a finite matrix group given.";
      return();
    }
  }
  else
  { if (typeof(#[size(#)])=="int")
    { gen_num=size(#)-1;
      mol_flag=0;
      int interval=0;
      v=0;
      int max=#[size(#)];
    }
   else
    { "ERROR:   If the two last parameters are not <int> and <intvec>, the last";
      "         parameter should be an <int>.";
      return();
    }
  }
  for (int i=1;i<=gen_num;i=i+1)
  { if (typeof(#[i])=="matrix")
    { if (nrows(#[i])<>n or ncols(#[i])<>n)
      { "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 first parameters should be a list of matrices";
      return();
    }
  }
 //----------------------------------------------------------------------------
  if (mol_flag==0)
  { if (ch==0)
    { matrix REY,M=reynolds_molien(#[1..gen_num],intvec(0,interval,v)); // one
                                       // will contain Reynolds operator and the
                                       // other enumerator and denominator of
                                       // Molien series
      matrix P=primary_char0_random(REY,M,max,v);
      matrix S,IS=secondary_char0(P,REY,M,v);
      return(P,S,IS);
    }
    else
    { list L=group_reynolds(#[1..gen_num],v);
      if (L[1]<>0)                     // testing whether we are in the modular
      { string newring="aksldfalkdsflkj"; // case
        if (minpoly==0)
        { if (v)
          { "  We are dealing with the non-modular case.";
          }
          molien(L[3..size(L)],newring,L[2],intvec(0,interval,v));
          matrix P=primary_charp_random(L[1],newring,max,v);
          matrix S,IS=secondary_charp(P,L[1],newring,v);
          if (voice==2)
          { kill aksldfalkdsflkj;
          }
          return(P,S,IS);
        }
        else
        { if (v)
          { "  Since it is impossible for this programme to calculate the Molien
 series for";
            "  invariant rings over extension fields of prime characteristic, we
 have to";
            "  continue without it.";
            "";

          }
          list l=primary_charp_no_molien_random(L[1],max,v);
          if (size(l)==2)
          { matix S=secondary_no_molien(l[1],L[1],l[2],v);
          }
          else
          { matix S=secondary_no_molien(l[1],L[1],v);
          }
          return(l[1],S);
        }
      }
      else                             // the modular case
      { if (v)
        { "  There is also no Molien series, we can make use of...";
          "";
          "  We can start looking for primary invariants...";
          "";
        }
        matrix P=primary_charp_without_random(#[1..gen_num],max,v);
        matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
        return(L[1],S);
      }
    }
  }
  if (mol_flag==1)                     // the user wants no calculation of the
  { list L=group_reynolds(#[1..gen_num],v); // Molien series
    if (ch==0)
    { list l=primary_char0_no_molien_random(L[1],max,v);
      if (size(l)==2)
      { matix S=secondary_no_molien(l[1],L[1],l[2],v);
      }
      else
      { matix S=secondary_no_molien(l[1],L[1],v);
      }
      return(l[1],S);
    }
    else
    { if (L[1]<>0)                     // testing whether we are in the modular
      { list l=primary_charp_no_molien_random(L[1],max,v); // case
        if (size(l)==2)
        { matix S=secondary_no_molien(l[1],L[1],l[2],v);
        }
        else
        { matix S=secondary_no_molien(l[1],L[1],v);
        }
        return(l[1],S);
      }
      else                             // the modular case
      { if (v)
        { "  We can start looking for primary invariants...";
          "";
        }
        matrix P=primary_charp_without_random(#[1..gen_num],max,v);
        matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
        return(L[1],S);
      }
    }
  }
  if (mol_flag==-1)
  { if (ch==0)
    { "ERROR:   Characteristic 0 can only be simulated in characteristic p>>0.
";
      return();
    }
    list L=group_reynolds(#[1..gen_num],v);
    string newring="aksldfalkdsflkj";
    molien(L[2..size(L)],newring,intvec(1,v));
    matrix P=primary_charp_random(L[1],newring,max,v);
    matrix S,IS=secondary_charp(P,L[1],newring,v);
    kill aksldfalkdsflkj;
    return(P,S,IS);
  }
  else                                 // the user specified that the
  { if (ch==0)                         // characteristic divides the group order
    { "ERROR:   The characteristic cannot divide the group order when it is 0.
";
      return();
    }
    if (v)
    { "";
    }
    matrix P=primary_charp_without_random(#[1..gen_num],max,v);
    matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
    return(L[1],S);
  }
}
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 P,S,IS=invariant_ring_random(A,1);
         print(P);
         print(S);
         print(IS);
}

proc algebra_containment (poly p, matrix A)
"USAGE:   algebra_containment(p,A);
         p: arbitrary <poly>, A: a 1xm <matrix> giving generators of a
         subalgebra of the basering
RETURN:  1 (TRUE) (type <int>) if p is contained in the subalgebra
         0 (FALSE) (type <int>) if <poly> is not contained
DISPLAY: a representation of p in terms of algebra generators A[1,i]=y(i) if p
         is contained in the subalgebra
EXAMPLE: example algebra_containment; shows an example
THEORY:  The ideal of algebraic relations of the algebra generators f1,...,fm
         given by A is computed introducing new variables y(i) and the product
         order: x^a*y^b > y^d*y^e if x^a > x^d or else if y^b > y^e. p reduces
         to a polynomial only in the y(i) <=> p is contained in the subring
         generated by the polynomials in A.
"
{ degBound=0;
  if (nrows(A)==1)
  { def br=basering;
    int n=nvars(br);
    int m=ncols(A);
    string mp=string(minpoly);
    execute "ring R=("+charstr(br)+"),(x(1..n),y(1..m)),(dp(n),dp(m));";
    execute "minpoly=number("+mp+");";
    ideal vars=x(1..n);
    map emb=br,vars;
    ideal A=ideal(emb(A));
    ideal check=emb(p);
    for (int i=1;i<=m;i=i+1)
    { A[i]=A[i]-y(i);
    }
    A=std(A);
    check[1]=reduce(check[1],A);
    A=elim(check,1,n);
    if (A[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 A[1][7]=x2+y2,z2,x4+y4,1,x2z-1y2z,xyz,x3y-1xy3;
         poly p1=x10z3-x8y2z3+2x6y4z3-2x4y6z3+x2y8z3-y10z3+x6z4+3x4y2z4+3x2y4z4+y6z4;
         algebra_containment(p1,A);
         poly p2=z;
         algebra_containment(p2,A);
}

proc module_containment(poly p, matrix P, matrix S)
"USAGE:   module_containment(p,P,S);
         p: arbitrary <poly>, P: a 1xn <matrix> giving generators of an algebra,
         S: a 1xt <matrix> giving generators of a module over the algebra
         generated by P
ASSUME:  n is the number of variables in the basering and the generators in P
         are algebraically independent
RETURNS: 1 (TRUE) (type <int>) if p is contained in the ring
         0 (FALSE) type <int>) if p is not contained
DISPLAY: the representation of p in terms of algebra generators P[1,i]=z(i) and
         module generators S[1,j]=y(j) if p is contained in the module
EXAMPLE: example module_containment; shows an example
THEORY:  The ideal of algebraic relations of all the generators p1,...,pn,
         s1,...,st given by P and S is computed introducing new variables y(j),
         z(i) and the product order: x^a*y^b*z^c > x^d*y^e*z^f if x^a > x^d
         with respect to the lp ordering or else if z^c > z^f with respect to
         the dp ordering or else if y^b > y^e with respect to the lp ordering
         again. p reduces to a polynomial only in the y(j) and z(i) linear in
         the z(i)) <=> p is contained in the module.
"
{ def br=basering;
  degBound=0;
  int n=nvars(br);
  if (ncols(P)==n and nrows(P)==1 and nrows(S)==1)
  { int m=ncols(S);
    string mp=string(minpoly);
    execute "ring R=("+charstr(br)+"),(x(1..n),y(1..m),z(1..n)),(lp(n),dp(m),lp(n));";
    execute "minpoly=number("+mp+");";
    ideal vars=x(1..n);
    map emb=br,vars;
    matrix P=emb(P);
    matrix S=emb(S);
    ideal check=emb(p);
    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[m+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:   the first <matrix> must have the same number of columns as the";
    "         basering and both <matrices> 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=x10z3-x8y2z3+2x6y4z3-2x4y6z3+x2y8z3-y10z3+x6z4+3x4y2z4+3x2y4z4+y6z4;
         module_containment(p1,P,S);
         poly p2=z;
         module_containment(p2,P,S);
}

proc orbit_variety (matrix F,string newring)
"USAGE:   orbit_variety(F,s);
         F: a 1xm <matrix> defing an invariant ring, s: a <string> giving 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 `s`)
EXAMPLE: example orbit_variety; shows an example
THEORY:  The ideal of algebraic relations of the invariant ring generators is
         calculated, then the variables of the original ring are eliminated and
         the polynomials that are left over define the orbit variety
"
{ if (newring=="")
  { "ERROR:   the second parameter 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_variety(F,newring);
         print(G);
         basering;
}

proc relative_orbit_variety(ideal I,matrix F,string newring)
"USAGE:   relative_orbit_variety(I,F,s);
         I: an <ideal> invariant under the action of a group, F: a 1xm
         <matrix> defining the invariant ring of this group, s: a <string>
         giving 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 I in the new ring (named s)
EXAMPLE: example relative_orbit_variety; shows an example
THEORY:  A Groebner basis of the ideal of algebraic relations of the invariant
         ring generators is calculated, then one of the basis elements plus the
         ideal generators. The variables of the original ring are eliminated and
         the polynomials that are left over define thecrelative orbit variety
         with respect to I.
"
{ if (newring=="")
  { "ERROR:   the third parameter may not be empty a <string>";
    return();
  }
  degBound=0;
  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";
         relative_orbit_variety(I,F,newring);
         print(G);
         basering;
}

proc image_of_variety(ideal I,matrix F)
"USAGE:   image_of_variety(I,F);
         I: an arbitray <ideal>, F: a 1xm <matrix> defining an invariant ring
         of a some matrix group
RETURN:  the <ideal> defining the image under that group of the variety defined
         by I
EXAMPLE: example image_of_variety; shows an example
THEORY:  relative_orbit_variety(I,F,s) is called and the newly introduced
         variables in the output are replaced by the generators of the
         invariant ring. This ideal in the original variables defines the image
         of the variety defined by I
"
{ if (nrows(F)==1)
  { def br=basering;
    int n=nvars(br);
    string newring="E";
    relative_orbit_variety(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(image_of_variety(I,F));
}