source: git/Singular/LIB/finvar.lib @ 731e67e

///////////////////////////////////////////////////////////////////////////////
version="$Id: finvar.lib,v 1.50 2006-07-18 15:48:13 Singular Exp $"
category="Invariant theory";
info="
LIBRARY:  finvar.lib    Invariant Rings of Finite Groups
AUTHOR: Agnes E. Heydtmann, email: agnes@math.uni-sb.de;
        Simon A. King, email: king@mfo.de
OVERVIEW:
 A library for computing polynomial invariants of finite matrix groups and
 generators of related varieties. The algorithms are based on B. Sturmfels,
 G. Kemper and W. Decker et al..

MAIN PROCEDURES:
 invariant_ring()                  generators of the invariant ring (i.r.)
 invariant_ring_random()           generators of the i.r., randomized alg.
 primary_invariants()              primary invariants (p.i.)
 primary_invariants_random()       primary invariants, randomized alg.

AUXILIARY PROCEDURES:
 cyclotomic()                      cyclotomic polynomial
 group_reynolds()                  finite group and Reynolds operator (R.o.)
 molien()                          Molien series (M.s.)
 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 of a degree
 invariant_basis_reynolds()        as invariant_basis(), with R.o.
 primary_char0()                   primary invariants (p.i.) in char 0
 primary_charp()                   primary invariants in char p
 primary_char0_no_molien()         p.i., char 0, without Molien series
 primary_charp_no_molien()         p.i., char p, without Molien series
 primary_charp_without()           p.i., char p, without R.o. or Molien series
 primary_char0_random()            primary invariants in char 0, randomized
 primary_charp_random()            primary invariants in char p, randomized
 primary_char0_no_molien_random()  p.i., char 0, without M.s., randomized
 primary_charp_no_molien_random()  p.i., char p, without M.s., randomized
 primary_charp_without_random()    p.i., char p, without R.o. or M.s., random.
 power_products()                  exponents for power products
 secondary_char0()                 secondary invariants (s.i.) in char 0
 irred_secondary_char0()           irreducible secondary invariants in char 0
 secondary_charp()                 secondary invariants in char p
 secondary_no_molien()             secondary invariants, without Molien series
 secondary_and_irreducibles_no_molien() s.i. & irreducible s.i., without M.s.
 secondary_not_cohen_macaulay()    s.i. when invariant ring not Cohen-Macaulay
 orbit_variety()                   ideal of the orbit variety
 rel_orbit_variety()               ideal of a relative orbit variety (new version)
 relative_orbit_variety()          ideal of a relative orbit variety (old version)
 image_of_variety()                ideal of the image of a variety
";
///////////////////////////////////////////////////////////////////////////////
// perhaps useful procedures (no help provided):
// unique()                        is a matrix among other matrices?
// exponent()                      gives the exponent of a number
// sort_of_invariant_basis()       lin. ind. invariants of a degree mod p.i.
// next_vector                     lists all of Z^n with first nonzero entry 1
// int_number_map                  integers 1..q are maped to q field elements
// search                          searches a number of p.i., char 0
// p_search                        searches a number of p.i., char p
// search_random                   searches a # of p.i., char 0, randomized
// p_search_random                 searches a # of p.i., char p, randomized
// concat_intmat                   concatenates two integer matrices
///////////////////////////////////////////////////////////////////////////////

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

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

proc cyclotomic (int i)
"USAGE:   cyclotomic(i); i integer > 0
RETURNS: the i-th cyclotomic polynomial (type <poly>) as one in the first ring
         variable
THEORY:  x^i-1 is divided by the j-th cyclotomic polynomial where j takes on
         the value of proper divisors of i
EXAMPLE: example cyclotomic; shows an example
"
{ 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++;
  }
  min=min/leadcoef(min);               // making sure that the leading
  return(min);                         // coefficient is 1
}
example
{ "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' order (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
THEORY:  The entire matrix group is generated by getting all left products of
         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 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.
EXAMPLE: example group_reynolds; shows an example
"
{ 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
  { matrix I=diag(1,n);                  // group element is
    matrix TEST=G(1);
    int o1=1;
    int o2;
    while (TEST<>I)
    { TEST=TEST*G(1);
      o1++;
    }
  }
  int i=1;
 // -------------- doubles among the generators should be avoided -------------
  for (int j=2;j<=gen_num;j++)         // 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++;
      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++;
        }
        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++)
    { for (k=1;k<=i;k++)
      { l=l+1;
        matrix P(l)=G(k)*G(m);         // possible new element
      }
    }
    j=0;
    for (k=1;k<=l;k++)
    { if (unique(G(1..g),P(k)))
      { j++;                           // a new factor for next run
        g++;
        matrix G(g)=P(k);              // a new group element -
        if (ch<>0 && minpoly==0 && i<>1) // finding out of which order the
        { TEST=G(g);                     //group element is
          o2=1;
          while (TEST<>I)
          { TEST=TEST*G(g);
            o2++;
          }
          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++;
   }
   return(i);
}
///////////////////////////////////////////////////////////////////////////////

proc molien (list #)
"USAGE:   molien(G1,G2,...[,ringname,lcm,flags]);
         G1,G2,...: nxn <matrices>, all elements of 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
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.
EXAMPLE: example molien; shows an example
"
{ 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++)
    { 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++)
      { 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++)
        { 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++)
            { 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
    }
    exportto(Top,`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++)
    { 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++)
    { 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++)
      { 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
THEORY:  The entire matrix group is generated by getting all left products of
         the generators with 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 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.
EXAMPLE: example reynolds_molien; shows an example
"
{ 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++)       // 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++;
        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 corresponding 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++)
        { l++;
          matrix P(l)=G(k)*G(m);       // possible new element
        }
      }
      j=0;
      for (k=1;k<=l;k++)
      { if (unique(G(1..g),P(k)))
        { j++;                         // a new factor for next run
          g++;
          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++)
    { 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++)
    { 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++)       // 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++;
        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++)
        { 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 corresponding 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++)
      { for (k=1;k<=i;k++)
        { l++;
          matrix P(l)=G(k)*G(m);       // possible new element
        }
      }
      j=0;
      for (k=1;k<=l;k++)
      { if (unique(G(1..g),P(k)))
        { j++;                         // a new factor for next run
          g++;
          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++)
          { 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:
@format
(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+...
@end format
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++)          // 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
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.
EXAMPLE: example evaluate_reynolds; shows an example
"
{ 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++)
    { 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
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.
EXAMPLE: example invariant_basis; shows an example
"
{ 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++)
  { 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--)
  { 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=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++)
    { temp=P/mon[j]-Pnew/mon[j];
      for (k=1;k<=m;k++)
      { 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[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
THEORY:  Monomials of degree d are mapped to invariants with the Reynolds
         operator. A linearly independent set is generated with the help of
         minbase.
EXAMPLE: example invariant_basis_reynolds; shows an example
"
{
 //---------------------- 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];
  int DEGB = degBound;
  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=DEGB;
      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=DEGB;
      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=DEGB;
      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);
  int DEGB=degBound;
  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++)          // 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=DEGB;
    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++)
    { 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=DEGB;
      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++)               // 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++)
  { 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++)                // 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++)
      { 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++;
  }
}

///////////////////////////////////////////////////////////////////////////////
// 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);            // or, better:
      //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++;                               // 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;
  if (ncols(B)<cd) { B[cd]=0; }
  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++)
      { 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++;
        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);        // or, better:
        // 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++;
  }
  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
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 (1998)).
EXAMPLE: example primary_char0; shows an example
"
{ 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=groebner(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=groebner(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++)
        { CI=CI+ideal(var(j)^d);
        }
        dB=dB+dif*(d-1);
        P=Pplus;
        sP=sPplus;
      }
      if (v)
      { for (j=1;j<=dif;j++)
        { "  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
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 (1998)).
EXAMPLE: example primary_charp; shows an example
"
{ 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=groebner(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 (ncols(B)<cd)
    {
      " warning: expected ",cd," invars, found ",ncols(B);
    }
    if (B[1]<>0)
    { Pplus=P+B;
      sPplus=groebner(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++)
        { 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++)
        { "  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 char 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
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 (1998)).
EXAMPLE: example primary_char0_no_molien; shows an example
"
{ 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=groebner(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++;                               // 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(groebner(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++)
        { CI=CI+ideal(var(j)^d);
        }
        dB=dB+dif*(d-1);
        P=Pplus;
        sP=groebner(P);
      }
      if (v)
      { for (j=1;j<=dif;j++)
        { "  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
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 (1998)).
EXAMPLE: example primary_charp_no_molien; shows an example
"
{ 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=groebner(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++;                               // 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=groebner(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++)
        { 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++)
        { "  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 char 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
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 (1998)). No Reynolds
         operator or Molien series is used.
EXAMPLE: example primary_charp_without; shows an example
"
{ 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++)
  { 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=groebner(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++;                               // 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=groebner(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++)
        { 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++)
        { "  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
{ "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 occurs not 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>) or ring name (type string) where the Molien series
         can be found in the char p case; 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
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 (1998)).
EXAMPLE: example primary_invariants; shows an example
"
{
 // ----------------- 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++)
  { 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++)
  { 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(groebner(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++)
          { for (k=0;k<=9;k++)
            { 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++)
    { 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(groebner(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++)
            { for (k=0;k<=9;k++)
              { 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
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 \"Generating a Noetherian Normalization of the Invariant Ring of
         a Finite Group\" by Decker, Heydtmann, Schreyer (1998)).
EXAMPLE: example primary_char0_random; shows an example
"
{ 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(groebner(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++)
        { 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++)
        { "  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
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 \"Generating a Noetherian Normalization of the Invariant Ring of
         a Finite Group\" by Decker, Heydtmann, Schreyer (1998)).
EXAMPLE: example primary_charp_random; shows an example
"
{ 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(groebner(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++)
        { 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++)
        { "  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 char 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
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 \"Generating a Noetherian Normalization of the Invariant Ring of
         a Finite Group\" by Decker, Heydtmann, Schreyer (1998)).
EXAMPLE: example primary_char0_no_molien_random; shows an example
"
{ 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++;                               // 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(groebner(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++)
        { 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++)
        { "  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
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 \"Generating a Noetherian Normalization of the Invariant Ring of
         a Finite Group\" by Decker, Heydtmann, Schreyer (1998)).
EXAMPLE: example primary_charp_no_molien_random; shows an example
"
{ 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++;                               // 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(groebner(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++)
        { 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++)
        { "  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 char 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
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 \"Generating a Noetherian Normalization of the Invariant Ring of
         a Finite Group\" by Decker, Heydtmann, Schreyer (1998)). No Reynolds
         operator or Molien series is used.
EXAMPLE: example primary_charp_without_random; shows an example
"
{ 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++)
  { 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++;                               // 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(groebner(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++)
        { 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++)
        { "  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
{ "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
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 \"Generating a Noetherian Normalization of the Invariant Ring of
         a Finite Group\" by Decker, Heydtmann, Schreyer (1998)).
EXAMPLE: example primary_invariants_random; shows an example
"
{
 // ----------------- 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++)
  { 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.";
          }
          if (typeof(L[2])=="int")
          { molien(L[3..size(L)],newring,L[2],intvec(0,interval,v));
          }
          else
          { 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";
    if (typeof(L[2])=="int")
    { molien(L[3..size(L)],newring,L[2],intvec(0,interval,v));
    }
    else
    { 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                                 // 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; shows 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++)
  { 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++)
        { 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
{ "EXAMPLE:"; echo=2;
         intvec dv=5,5,5,10,10;
         print(power_products(dv,10));
         print(power_products(dv,7));
}
///////////////////////////////////////////////////////////////////////////////

static proc old_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 numerator 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
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++)                   // 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=groebner(ideal(P));
  ideal TEST,B,IS;
  ideal S=1;                           // 1 is the first secondary invariant -
 //--------------------- generating secondary invariants ----------------------
  for (i=2;i<=m;i++)                   // 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++)     // going through all the power products
        { pp=1;
          for (k=1;k<=nrows(PP);k++)
          { pp=pp*IS[1,k]^PP[k,j];
          }
          if (reduce(pp,TEST)<>0)
          { S=S,pp;
            counter++;
            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,pp);
            }
            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++)
            { deg_vec=deg_vec,i-1;
              if (v)
              { "  We find: "+string(B[j]);
              }
            }
          }
          else
          { for (j=1;j<=int(dimmat[i,1]);j++)
            { 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++;
            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++;
              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,B[j]);
              }
            }
          }
        }
      }
      if (v)
      { "";
      }
    }
  }
  if (v)
  { "  We're done!";
    "";
  }
  return(matrix(S),matrix(IS));
}

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


proc secondary_char0 (matrix P, matrix REY, matrix M, list #)
"
USAGE:    secondary_char0(P,REY,M[,v][,\"AH\"]);
@*          P: a 1xn <matrix> with homogeneous primary invariants, where
               n is the number of variables of the basering;
@*          REY: a gxn <matrix> representing the Reynolds operator, where
               g the size of the corresponding group;
@*          M: a 1x2 <matrix> giving numerator and denominator of the Molien
               series;
@*          v: an optional <int>;
@*          \"AH\": if this string occurs as (optional) parameter, then an
               old version of secondary_char0 is used (for downward
               compatibility)
ASSUME:   REY is the 1st return value of group_reynolds(), reynolds_molien() or
          the second one of primary_invariants();
@*        M is the return value of molien()
          or the second one of reynolds_molien() or the third one of
          primary_invariants()
RETURN:   Homogeneous secondary invariants and irreducible secondary
          invariants of the invariant ring (both type <matrix>)
DISPLAY:  Information on the progress of the computations if v is an integer
          different from 0.
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. Among these images
          or their power products we pick a maximal subset that is linearly
          independent modulo the primary invariants (see paper \"Some
          Algorithms in Invariant Theory of Finite Groups\" by Kemper and
          Steel (1997)). The size of this set can be read off from the
          Molien series.
NOTE:     Secondary invariants are not uniquely determined by the given data.
          Specifically, the output of secondary_char0(P,REY,M,"AH") will
          differ from the output of secondary_char0(P,REY,M). However, the
          ideal generated by the irreducible homogeneous
          secondary invariants will be the same in both cases.
@*        There is an internal parameter \"pieces\" that, by default, equals 15.
          Setting \"pieces\" to a smaller value might decrease the memory consumption,
          but increase the running time.
SEE ALSO: irred_secondary_char0;
EXAMPLE:  example secondary_char0; shows an example
"
{ def br=basering;

 //----------------- 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(#[1])=="int")
    { int v=#[1];
    }
    else
    { int v=0;
    }
  }
  else
  { int v=0;
  }
  if (size(#)>0)
  { if (typeof(#[size(#)])=="string")
    { if (#[size(#)]=="AH")
      { if (typeof(#[1])=="int")
        { matrix S,IS = old_secondary_char0(P,REY,M,#[1]);
          return(S,IS);
        }
        else
        { matrix S,IS = old_secondary_char0(P,REY,M);
          return(S,IS);
        }
      }
      else
      { "ERROR:   If the last optional parameter is a string, it should be \"AH\".";
        return(matrix(ideal()),matrix(ideal()));
      }
    }
  }
  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, irrcounter;
 //- finding the polynomial giving number and degrees of secondary invariants -
  poly p=1;
  for (j=1;j<=n;j++)                   // 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 [constant] term of
                                       // p needs to be 1
  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)
  { "We need to find";
    for (j=1;j<=m;j++)
     { if (int(dimmat[j,1])<>1)
       { int(dimmat[j,1]), "secondary invariants in degree",j-1;
       }
       else
       { "1 secondary invariant in degree",j-1;
       }
     }
  }
  if (v)
  { "  In degree 0 we have: 1";
    "";
  }
 //-------------------------- initializing variables --------------------------
  ideal ProdCand;          // contains products of secondary invariants,
                           // i.e., candidates for reducible sec. inv.
  ideal Mul1,Mul2;

  int dgb=degBound;
  degBound = 0;
  option(redSB);
  ideal sP = groebner(ideal(P));
  ideal Reductor,SaveRed;  // sP union Reductor is a Groebner basis up to degree i-1
  int SizeSave;

  list SSort;         // sec. inv. first sorted by degree and then sorted by the
                      // minimal degree of a non-constant invariant factor.
  list ISSort;        // irr. sec. inv. sorted by degree

  poly helpP;
  ideal helpI;
  ideal Indicator;        // will tell us which candidates for sec. inv. we can choose
  ideal ReducedCandidates;
  int helpint;
  int k,k2,k3,minD;
  int ii;
  int saveAttr;
  ideal B,IS;              // IS will contain all irr. sec. inv.

  int pieces = 15;   // If this parameter is small, the memory consumption will decrease.
                     // In some cases, a careful choice will speed up computations

  for (i=1;i<m;i++)
  { SSort[i] = list();
    for (k=1;k<=i;k++)
    { SSort[i][k]=ideal();
      attrib(SSort[i][k],"size",0);
    }
   }
//--------------------- generating secondary invariants ----------------------
  for (i=2;i<=m;i++)                   // going through dimmat -
  { if (int(dimmat[i,1])<>0)           // when it is == 0 we need to find no
    {                                  // 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)...";
        "  Looking for Power Products...";
      }
      counter = 0;                       // we'll count up to dimmat[i,1]
      Reductor = ideal(0);
      helpint = 0;
      SaveRed = Reductor;
      SizeSave = 0;
      attrib(Reductor,"isSB",1);
      attrib(SaveRed,"isSB",1);

// We start searching for reducible secondary invariants in degree i-1, i.e., those
// that are power products of irreducible secondary invariants.
// It suffices to restrict the search at products of one _irreducible_ sec. inv. (Mul1)
// with some sec. inv. (Mul2).
// Moreover, we avoid to consider power products twice since we take a product
// into account only if the minimal degree of a non-constant invariant factor in "Mul2" is not
// smaller than the degree of "Mul1".
      for (k=1;k<i-1;k++)
      { if (int(dimmat[i,1])==counter)
        { break;
        }
        if (typeof(ISSort[k])<>"none")
        { Mul1 = ISSort[k];
        }
        else
        { Mul1 = ideal(0);
        }
        if ((int(dimmat[i-k,1])>0) && (Mul1[1]<>0))
        { for (minD=k;minD<i-k;minD++)
          { if (int(dimmat[i,1])==counter)
            { break;
            }
            for (k2=1;k2 <= ((attrib(SSort[i-k-1][minD],"size")-1) div pieces)+1; k2++)
            { if (int(dimmat[i,1])==counter)
              { break;
              }
              Mul2=ideal(0);
              if (attrib(SSort[i-k-1][minD],"size")>=k2*pieces)
              { for (k3=1;k3<=pieces;k3++)
                { Mul2[k3] = SSort[i-k-1][minD][((k2-1)*pieces)+k3];
                }
              }
              else
              { for (k3=1;k3<=(attrib(SSort[i-k-1][minD],"size") mod pieces);k3++)
                { Mul2[k3] = SSort[i-k-1][minD][((k2-1)*pieces)+k3];
                }
              }
              ProdCand = Mul1*Mul2;
              ReducedCandidates = reduce(ProdCand,sP);
                  // sP union SaveRed union Reductor is a homogeneous Groebner basis
                  // up to degree i-1.
                  // We first reduce by sP (which is fixed, so we can do it once for all),
                  // then by SaveRed resp. by Reductor (which is modified during
                  // the computations).
              Indicator = reduce(ReducedCandidates,SaveRed);
                             // If Indicator[ii]==0 then ReducedCandidates it the reduction
                             // of an invariant that is in the algebra generated by primary
                             // invariants and previously computed secondary invariants.
                             // Otherwise ProdCand[ii] can be taken as secondary invariant.
          if (size(Indicator)<>0)
              { for (ii=1;ii<=ncols(ProdCand);ii++)     // going through all the power products
                { helpP = reduce(Indicator[ii],Reductor);
                  if (helpP <> 0)
                  { counter++;
                    saveAttr = attrib(SSort[i-1][k],"size")+1;
                    SSort[i-1][k][saveAttr] = ProdCand[ii];
                        // By construction, this is a _reducible_ s.i.
                    attrib(SSort[i-1][k],"size",saveAttr);
                    if (v)
                    { "We found",counter, " of", int(dimmat[i,1])," secondaries in degree",i-1;
                    }
                    if (int(dimmat[i,1])<>counter)
                    { helpint++;
                      Reductor[helpint] = helpP;
                        // Lemma: If G is a homogeneous Groebner basis up to degree i-1 and p is a
                        // homogeneous polynomial of degree i-1 then G union NF(p,G) is
                        // a homogeneous Groebner basis up to degree i-1.
                      attrib(Reductor, "isSB",1);
                        // if Reductor becomes too large, we reduce the whole of Indicator by
                        // it, save it in SaveRed, and work with a smaller Reductor. This turns
                        // out to save a little time.
                      if (ncols(Reductor)>100)
                      { Indicator=reduce(Indicator,Reductor);
                        for (k3=1;k3<=ncols(Reductor);k3++)
                        { SaveRed[SizeSave+k3] = Reductor[k3];
                        }
                        SizeSave=SizeSave+size(Reductor);
                        Reductor = ideal(0);
                        helpint = 0;
                        attrib(SaveRed, "isSB",1);
                        attrib(Reductor, "isSB",1);
                      }
                    }
                    else
                    { break;
                    }
                  }
                }
              }
              for (k3=1;k3<=size(Reductor);k3++)
              { SaveRed[SizeSave+k3] = Reductor[k3];
              }
              SizeSave=SizeSave+size(Reductor);
              Reductor = ideal(0);
              helpint = 0;
              attrib(SaveRed, "isSB",1);
              attrib(Reductor, "isSB",1);
            }
          }
        }
      }

      // The remaining sec. inv. are irreducible!
      if (int(dimmat[i,1])<>counter)   // need more than all the power products
      { if (v)
        { "There are irreducible secondary invariants in degree ", i-1," !!";
        }
        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 B
        { IS=IS+B;
          saveAttr = attrib(SSort[i-1][i-1],"size")+int(dimmat[i,1]);
          SSort[i-1][i-1] = SSort[i-1][i-1] + B;
          attrib(SSort[i-1][i-1],"size", saveAttr);
          if (typeof(ISSort[i-1]) <> "none")
          { ISSort[i-1] = ISSort[i-1] + B;
          }
          else
          { ISSort[i-1] = B;
          }
          if (v) {"  We found: "; print(B);}
        }
        else
        { irrcounter=0;
          j=0;                         // j goes through all of B -
        // Compare the comments on the computation of reducible sec. inv.!
          ReducedCandidates = reduce(B,sP);
          Indicator = reduce(ReducedCandidates,SaveRed);
          while (int(dimmat[i,1])<>counter)
          { j++;
            helpP = reduce(Indicator[j],Reductor);
            if (helpP <>0)                  // B[j] should be added
            { counter++; irrcounter++;
              IS=IS,B[j];
              saveAttr = attrib(SSort[i-1][i-1],"size")+1;
              SSort[i-1][i-1][saveAttr] = B[j];
              attrib(SSort[i-1][i-1],"size",saveAttr);
              if (typeof(ISSort[i-1]) <> "none")
              { ISSort[i-1][irrcounter] = B[j];
              }
              else
              { ISSort[i-1] = ideal(B[j]);
              }
              if (v)
              { "  We found the irreducible sec. inv. "+string(B[j]);
              }
              helpint++;
              Reductor[helpint] = helpP;
              attrib(Reductor, "isSB",1);
              if (ncols(Reductor)>100)
              { Indicator=reduce(Indicator,Reductor);
                for (k3=1;k3<=ncols(Reductor);k3++)
                { SaveRed[SizeSave+k3] = Reductor[k3];
                }
                SizeSave=SizeSave+size(Reductor);
                Reductor = ideal(0);
                helpint = 0;
                attrib(SaveRed, "isSB",1);
                attrib(Reductor, "isSB",1);
              }
            }
          }
        }
      }
      if (v)
      { "";
      }
    }
  }
  if (v)
  { "  We're done!";
    "";
  }
  degBound = dgb;

// Prepare return:
  int TotalNumber;
  for (k=1;k<=m;k++)
  { TotalNumber = TotalNumber + int(dimmat[k,1]);
  }
  matrix S[1][TotalNumber];
  S[1,1]=1;
  j=1;
  for (k=1;k<m;k++)
  { for (k2=1;k2<=k;k2++)
    { if (typeof(attrib(SSort[k][k2],"size"))=="int")
     {for (i=1;i<=attrib(SSort[k][k2],"size");i++)
      { j++;
        S[1,j] = SSort[k][k2][i];
      }
      SSort[k][k2]=ideal();
     }
    }
  }
  return(S,matrix(compress(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],1);
         print(S);
         print(IS);
}
//example
//{ "EXAMPLE: S. King"; echo=2;
//ring r= 0, (V1,V2,V3,V4,V5,V6),dp;
//matrix A1[6][6] = 0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0;
//matrix A2[6][6] = 0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0;
//list L = primary_invariants(A1,A2); // sorry, this takes a while...
//matrix S, IS = secondary_char0(L[1],L[2],L[3],1);
//S;
//IS;
//}
///////////////////////////////////////////////////////////////////////////////


proc irred_secondary_char0 (matrix P, matrix REY, matrix M, list #)
"
USAGE:    irred_secondary_char0(P,REY,M[,v]);
@*          P: a 1xn <matrix> with homogeneous primary invariants, where
               n is the number of variables of the basering;
@*          REY: a gxn <matrix> representing the Reynolds operator, where
               g the size of the corresponding group;
@*          M: a 1x2 <matrix> giving numerator and denominator of the
               Molien series;
@*          v: an optional <int>;
RETURN:   Irreducible homogeneous secondary invariants of the invariant ring
          (type <matrix>)
ASSUME:   REY is the 1st return value of group_reynolds(), reynolds_molien() or
          the second one of primary_invariants();
          M is the return value of molien() or the second one of
          reynolds_molien() or the third one of primary_invariants()
DISPLAY:  Information on the progress of computations if v does not equal 0
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. Among these images
          or their power products we pick a maximal subset that is linearly
          independent modulo the primary invariants (see paper \"Some
          Algorithms in Invariant Theory of Finite Groups\" by Kemper and
          Steel (1997)). The size of this set can be read off from the
          Molien series. Here, only irreducible secondary
          invariants are explicitly computed, which saves time and memory.
NOTE:     There is an internal parameter \"pieces\" that, by default, equals 15.
          Setting \"pieces\" to a smaller value might decrease the memory consumption,
          but increase the running time.
SEE ALSO: secondary_char0
KEYWORDS: irreducible secondary invariant
EXAMPLE:  example irred_secondary_char0; shows an example
"
{ def br=basering;

 //----------------- checking input and setting verbose mode ------------------
  if (char(br)<>0)
  { "ERROR:   irred_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, irrcounter;
 //- finding the polynomial giving number and degrees of secondary invariants -
  poly p=1;
  for (j=1;j<=n;j++)                   // 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 [constant] term of
                                       // p needs to be 1
  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)
  { "We need to find";
    for (j=1;j<=m;j++)
     { if (int(dimmat[j,1])<>1)
       { int(dimmat[j,1]), "secondary invariants in degree",j-1;
       }
       else
       { "1 secondary invariant in degree",j-1;
       }
     }
  }
  if (v)
  { "  In degree 0 we have: 1";
    "";
  }
 //-------------------------- initializing variables --------------------------
  ideal ProdCand;          // contains products of secondary invariants,
                           // i.e., candidates for reducible sec. inv.

  ideal Mul1,Mul2;

  int dgb=degBound;
  degBound = 0;
  option(redSB);
  ideal sP = groebner(ideal(P));
  ideal Reductor,SaveRed;  // sP union Reductor is a Groebner basis up to degree i-1
  int SizeSave;

  list RedSSort;           // sec. Inv. reduced by sP, sorted first by degree and then
                           // by the minimal degree of a non-constant invariant factor
  list RedISSort;          // irr. sec. Inv. reduced by sP, sorted by degree

  poly helpP;
  ideal helpI;
  ideal Indicator;        // will tell us which candidates for sec. inv. we can choose
  ideal ReducedCandidates;
  int helpint;
  int k,k2,k3,minD;
  int ii;
  int saveAttr;
  ideal B,IS;              // IS will contain all irr. sec. inv.

  int pieces = 15;   // If this parameter is small, the memory consumption will decrease.
                     // In some cases, a careful choice will speed up computations

  for (i=1;i<m;i++)
  { RedSSort[i] = list();
    for (k=1;k<=i;k++)
    { RedSSort[i][k]=ideal();
      attrib(RedSSort[i][k],"size",0);
    }
   }
//--------------------- generating secondary invariants ----------------------
  for (i=2;i<=m;i++)                   // going through dimmat -
  { if (int(dimmat[i,1])<>0)           // when it is == 0 we need to find no
    {                                  // 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)...";
        "  Looking for Power Products...";
      }
      counter = 0;                       // we'll count up to dimmat[i,1]
      Reductor = ideal(0);
      helpint = 0;
      SaveRed = Reductor;
      SizeSave = 0;
      attrib(Reductor,"isSB",1);
      attrib(SaveRed,"isSB",1);

// We start searching for reducible secondary invariants in degree i-1, i.e., those
// that are power products of irreducible secondary invariants.
// It suffices to restrict the search at products of one _irreducible_ sec. inv. (Mul1)
// with some sec. inv. (Mul2).
// Moreover, we avoid to consider power products twice since we take a product
// into account only if the minimal degree of a non-constant invariant factor in "Mul2" is not
// smaller than the degree of "Mul1".
// Finally, as we are not interested in the reducible sec. inv., we will only
// work with their reduction modulo sP --- this allows to detect a secondary invariant
// without to actually compute it!
      for (k=1;k<i-1;k++)
      { if (int(dimmat[i,1])==counter)
        { break;
        }
        if (typeof(RedISSort[k])<>"none")
        { Mul1 = RedISSort[k];
        }
        else
        { Mul1 = ideal(0);
        }
        if ((int(dimmat[i-k,1])>0) && (Mul1[1]<>0))
        { for (minD=k;minD<i-k;minD++)
          { if (int(dimmat[i,1])==counter)
            { break;
            }
            for (k2=1;k2 <= ((attrib(RedSSort[i-k-1][minD],"size")-1) div pieces)+1; k2++)
            { if (int(dimmat[i,1])==counter)
              { break;
              }
              Mul2=ideal(0);
              if (attrib(RedSSort[i-k-1][minD],"size")>=k2*pieces)
              { for (k3=1;k3<=pieces;k3++)
                { Mul2[k3] = RedSSort[i-k-1][minD][((k2-1)*pieces)+k3];
                }
              }
              else
              { for (k3=1;k3<=(attrib(RedSSort[i-k-1][minD],"size") mod pieces);k3++)
                { Mul2[k3] = RedSSort[i-k-1][minD][((k2-1)*pieces)+k3];
                }
              }
              ProdCand = Mul1*Mul2;
              ReducedCandidates = reduce(ProdCand,sP);
                  // sP union SaveRed union Reductor is a homogeneous Groebner basis
                  // up to degree i-1.
                  // We first reduce by sP (which is fixed, so we can do it once for all),
                  // then by SaveRed resp. by Reductor (which is modified during
                  // the computations).
              Indicator = reduce(ReducedCandidates,SaveRed);
                             // If Indicator[ii]==0 then ReducedCandidates it the reduction
                             // of an invariant that is in the algebra generated by primary
                             // invariants and previously computed secondary invariants.
                             // Otherwise ReducedCandidates[ii] is the reduction of an invariant
                             // that we can take as secondary invariant.
              if (size(Indicator)<>0)
              { for (ii=1;ii<=ncols(ProdCand);ii++)     // going through all the power products
                { helpP = reduce(Indicator[ii],Reductor);
                  if (helpP <> 0)
                  { counter++;
                    saveAttr = attrib(RedSSort[i-1][k],"size")+1;
                    RedSSort[i-1][k][saveAttr] = ReducedCandidates[ii];
                        // By construction, this is the reduction of a reducible s.i.
                        // of degree i-1.
                    attrib(RedSSort[i-1][k],"size",saveAttr);
                    if (v)
                    { "We found",counter, " of", int(dimmat[i,1])," secondaries in degree",i-1;
                    }
                    if (int(dimmat[i,1])<>counter)
                    { helpint++;
                      Reductor[helpint] = helpP;
                        // Lemma: If G is a homogeneous Groebner basis up to degree i-1 and p is a
                        // homogeneous polynomial of degree i-1 then G union NF(p,G) is
                        // a homogeneous Groebner basis up to degree i-1.
                      attrib(Reductor, "isSB",1);
                        // if Reductor becomes too large, we reduce the whole of Indicator by
                        // it, save it in SaveRed, and work with a smaller Reductor. This turns
                        // out to save a little time.
                      if (ncols(Reductor)>100)
                      { Indicator=reduce(Indicator,Reductor);
                        for (k3=1;k3<=ncols(Reductor);k3++)
                        { SaveRed[SizeSave+k3] = Reductor[k3];
                        }
                        SizeSave=SizeSave+size(Reductor);
                        Reductor = ideal(0);
                        helpint = 0;
                        attrib(SaveRed, "isSB",1);
                        attrib(Reductor, "isSB",1);
                      }
                    }
                    else
                    { break;
                    }
                  }
                }
              }
              for (k3=1;k3<=size(Reductor);k3++)
              { SaveRed[SizeSave+k3] = Reductor[k3];
              }
              SizeSave=SizeSave+size(Reductor);
              Reductor = ideal(0);
              helpint = 0;
              attrib(SaveRed, "isSB",1);
              attrib(Reductor, "isSB",1);
            }
          }
        }
      }

      // The remaining sec. inv. are irreducible!
      if (int(dimmat[i,1])<>counter)   // need more than all the power products
      { if (v)
        { "There are irreducible secondary invariants in degree ", i-1," !!";
        }
        B=sort_of_invariant_basis(sP,REY,i-1,int(dimmat[i,1])*6);
                                       // B is a set of
                                       // images of kbase(sP,i-1) under the
                                       // Reynolds operator that are linearly
                                       // independent and that do not reduce to
                                       // 0 modulo sP -
        if (counter==0 && ncols(B)==int(dimmat[i,1])) // then we can take all of B
        { IS=IS+B;
          saveAttr = attrib(RedSSort[i-1][i-1],"size")+int(dimmat[i,1]);
          RedSSort[i-1][i-1] = RedSSort[i-1][i-1] + B;
          attrib(RedSSort[i-1][i-1],"size", saveAttr);
          if (typeof(RedISSort[i-1]) <> "none")
          { RedISSort[i-1] = RedISSort[i-1] + B;
          }
          else
          { RedISSort[i-1] = B;
          }
          if (v) {"  We found: ";print(B);}
        }
        else
        { irrcounter=0;
          j=0;                         // j goes through all of B -
        // Compare the comments on the computation of reducible sec. inv.!
          ReducedCandidates = reduce(B,sP);
          Indicator = reduce(ReducedCandidates,SaveRed);
          while (int(dimmat[i,1])<>counter)    // need to find dimmat[i,1]
          {                                    // invariants that are linearly independent
            j++;
            helpP = reduce(Indicator[j],Reductor);
            if (helpP <>0)                     // B[j] should be added
            { counter++; irrcounter++;
              IS=IS,B[j];
              saveAttr = attrib(RedSSort[i-1][i-1],"size")+1;
              RedSSort[i-1][i-1][saveAttr] = ReducedCandidates[j];
              attrib(RedSSort[i-1][i-1],"size",saveAttr);
              if (typeof(RedISSort[i-1]) <> "none")
              { RedISSort[i-1][irrcounter] = ReducedCandidates[j];
              }
              else
              { RedISSort[i-1] = ideal(ReducedCandidates[j]);
              }
              if (v)
              { "  We found the irreducible sec. inv. "+string(B[j]);
              }
              helpint++;
              Reductor[helpint] = helpP;
              attrib(Reductor, "isSB",1);
              if (ncols(Reductor)>100)
              { Indicator=reduce(Indicator,Reductor);
                for (k3=1;k3<=ncols(Reductor);k3++)
                { SaveRed[SizeSave+k3] = Reductor[k3];
                }
                SizeSave=SizeSave+size(Reductor);
                Reductor = ideal(0);
                helpint = 0;
                attrib(SaveRed, "isSB",1);
                attrib(Reductor, "isSB",1);
              }
            }
          }
        }
      }
      if (v)
      { "";
      }
    }
  }
  if (v)
  { "  We're done!";
    "";
  }
  degBound = dgb;

  return(matrix(compress(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 IS=irred_secondary_char0(L[1..3],1);
//         print(IS);
//}
example
{ "EXAMPLE: S. King"; echo=2;
ring r= 0, (a,b,c,d,e,f),dp;
matrix A1[6][6] = 0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0;
matrix A2[6][6] = 0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0;
list L = primary_invariants(A1,A2); // sorry, this takes a while...
matrix IS = irred_secondary_char0(L[1],L[2],L[3],0);
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 a ring of
         char 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
THEORY:  Secondary invariants are calculated by finding a basis (in terms of
         monomials) of the basering modulo 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)).
EXAMPLE: example secondary_charp; shows an example
"
{ 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++)
  { deg_dim_vec[j]=deg(P[j]);
  }
  setring `ring_name`;
  poly p=1;
  for (j=1;j<=n;j++)                   // 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++)
  { 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=groebner(ideal(P));
  ideal TEST,B,IS;
  ideal S=1;                           // 1 is the first secondary invariant
 //------------------- generating secondary invariants ------------------------
  for (i=2;i<=m;i++)                   // 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++)     // going through all of those
        { pp=1;
          for (k=1;k<=nrows(PP);k++)
          { pp=pp*IS[1,k]^PP[k,j];
          }
          if (reduce(pp,TEST)<>0)
          { S=S,pp;
            counter++;
            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,pp);
            }
            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++)
            { deg_vec=deg_vec,i-1;
              if (v)
              { "  We find: "+string(B[j]);
              }
            }
          }
          else
          { for (j=1;j<=deg_dim_vec[i];j++)
            { 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++;
            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++;
              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,B[j]);
              }
            }
          }
        }
      }
      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 char 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
THEORY:  Secondary invariants are calculated by finding a basis (in terms of
         monomials) of the basering modulo primary invariants, mapping those to
         invariants with the Reynolds operator and using these images as
         candidates for secondary invariants.
EXAMPLE: example secondary_no_molien; shows an example
"
{ 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++)
  { 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=groebner(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++;
    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++)
      {                                // that are linearly independent modulo
                                       // TEST
        if (reduce(B[j],TEST)<>0)      // B[j] should be added
        { S=S,B[j];
          counter++;
          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,B[j]);
            }
          }
        }
      }
    }
    else
    { if (size(deg_vec)==k)
      { k=1; }
      else
      { k++; }
    }
  }
  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_and_irreducibles_no_molien (matrix P, matrix REY, list #)
"USAGE:   secondary_and_irreducibles_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
THEORY:  Secondary invariants are calculated by finding a basis (in terms of
         monomials) of the basering modulo 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)).
EXAMPLE: example secondary_and_irreducibles_no_molien; shows an example
"
{ 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++)
  { 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=groebner(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++;
    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++)     // going through all those power products
        { pp=1;
          for (k=1;k<=nrows(PP);k++)
          { pp=pp*IS[1,k]^PP[k,j];
          }
          if (reduce(pp,TEST)<>0)
          { S=S,pp;
            counter++;
            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,pp);
            }
            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++)
        { 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++;
            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,B[j]);
              }
            }
          }
        }
      }
    }
    else
    { if (size(deg_vec)==l)
      { l=1;
      }
      else
      { l++;
      }
    }
  }
  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_and_irreducibles_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
THEORY:  Secondary invariants are generated following \"Generating Invariant
         Rings of Finite Groups over Arbitrary Fields\" by Kemper (1996).
EXAMPLE: example secondary_not_cohen_macaulay; shows an example
"
{ 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++)
      { 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++)
      { 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, trivialSI=secondary_charp(P,REY,"alskdfalkdsj",v);
  kill trivialSI;
  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++)
  { 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)
  }
  //intvec save_opts=option(get);
  //option(returnSB,redSB);
  //module M(2)=syz(M(1));               // nres(M(1),2)[2];
  //option(set,save_opts);
  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++)
  { for (j=1;j<=k;j++)
    { 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++)
  { 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++)
    { "   "+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
{ "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 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 occurs not 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
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 \"Generating a
         Noetherian Normalization of the Invariant Ring of a Finite Group\" by
         Decker, Heydtmann, Schreyer (1998)). 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 \"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).
EXAMPLE: example invariant_ring; shows an example
"
{ 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.";
          }
          if (typeof(L[2])=="int")
          { molien(L[3..size(L)],newring,L[2],intvec(0,interval,v));
          }
          else
          { molien(L[2..size(L)],newring,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
          { matrix 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)
      { matrix S=secondary_no_molien(l[1],L[1],l[2],v);
      }
      else
      { matrix 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)
        { matrix S=secondary_no_molien(l[1],L[1],l[2],v);
        }
        else
        { matrix 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";
    if (typeof(L[2])=="int")
    { molien(L[3..size(L)],newring,L[2],intvec(1,interval,v));
    }
    else
    { 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 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,
         then the characteristic of the base field divides the group order
         (i.e. we will not even attempt 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
         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
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.
EXAMPLE: example invariant_ring_random; shows an example
"
{ 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++)
  { 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.";
          }
          if (typeof(L[2])=="int")
          { molien(L[3..size(L)],newring,L[2],intvec(0,interval,v));
          }
          else
          { molien(L[2..size(L)],newring,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)
          { matrix S=secondary_no_molien(l[1],L[1],l[2],v);
          }
          else
          { matrix 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(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_random(L[1],max,v);
      if (size(l)==2)
      { matrix S=secondary_no_molien(l[1],L[1],l[2],v);
      }
      else
      { matrix 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)
        { matrix S=secondary_no_molien(l[1],L[1],l[2],v);
        }
        else
        { matrix 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";
    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_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 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`)
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
EXAMPLE: example orbit_variety; shows an example
"
{ 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[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++)
    { vars[i]=0;
    }
    vars=vars,y(1..m);
    map emb=R,vars;
    ideal G=emb(I);
    kill emb, vars, R;
    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;
}
///////////////////////////////////////////////////////////////////////////////

static proc degvec(ideal I)
"USAGE:  degvec(I);
         I an <ideal>.
RETURN:  the <intvec> of degrees of the generators of I.
"
{ intvec v;
  for (int j = 1;j<=ncols(I);j++)
  { v[j]=deg(I[j]);
  }
  return(v);
}

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

proc rel_orbit_variety(ideal I,matrix F,list #)
"USAGE:    rel_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: optional <string>; if s is present then (for downward
               compatibility) the old procedure <relative_orbit_variety>
               is called, and in this case s gives the name of a new <ring>.
RETURN:   Without optional s, a list L of two rings is returned.
@*        The ring L[1] carries a weighted degree order with variables
            y(1..m), the weight of y(k) equal to the degree of the
            k-th generators F[1,k] of the invariant ring.
            L[1] contains a Groebner basis (type <ideal>, named G) of the
            ideal defining the relative orbit variety with  respect to I.
@*        The ring L[2] has the variables of the basering together with y(1..m)
            and carries a block order: The first block is the order of the
            basering, the second is the weighted degree order occuring in L[1].
            L[2] contains G and a Groebner basis (type <ideal>, named Conv)
            such that if p is any invariant polynomial expressed in the
            variables of the basering then reduce(p,Conv) is a polynomial in
            the new variables y(1..m) such that evaluation at the generators
            of the invariant ring yields p. This can be used to avoid the
            application of <algebra_containment>
            (see @ref{algebra_containment}).
@*        For the case of optional s, see @ref{relative_orbit_variety}.
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 define the relative
          orbit variety with respect to I. The elimination is done by a
          weighted blockorder that has the advantage of dealing with
          quasi-homogeneous ideals.
NOTE:     We provide the ring L[1] for the sake of downward compatibility,
          since it is closer to the ring returned by relative_orbit_variety
          than L[2]. However, L[1] carries a weighted degree order, whereas
          the ring returned by relative_orbit_variety is lexicographically
          ordered.
SEE ALSO: relative_orbit_variety
EXAMPLE:  example rel_orbit_variety; shows an example.
"
{ if (size(#)>0)
  { if (typeof(#[1])=="string")
    { relative_orbit_variety(I,F,#[1]);
      keepring basering;
      return();
    }
    else
    { "ERROR:    the third parameter may either be empty or a <string>.";
      return();
    }
  }
  degBound=0;
  if (nrows(F)==1)
  { option(redSB,noredTail);
    def br=basering;
    int n=nvars(br);
    int m=ncols(F);
    // In the following ring definition, any elimination order would work.
    list rlist = ringlist(br);
    int i;
    for (i=1;i<=m;i++)
    { rlist[2][n+i]="y("+string(i)+")";
    }
    rlist[3][size(rlist[3])+1] = rlist[3][size(rlist[3])];
    rlist[3][size(rlist[3])-1][1] = "wp";
    rlist[3][size(rlist[3])-1][2] = degvec(ideal(F));
    def newring = ring(rlist);

    list smallrlist;
    smallrlist[1]=rlist[1];
    smallrlist[2]=list();
    for (i=1;i<=m;i++)
    { smallrlist[2][i]="y("+string(i)+")";
    }
    smallrlist[3]=list();
    smallrlist[3][1]=list();
    smallrlist[3][1][1] = "wp";
    smallrlist[3][1][2] = degvec(ideal(F));
    smallrlist[3][2] = rlist[3][size(rlist[3])];
    smallrlist[4] = rlist[4];
    def smallring = ring(smallrlist);

    setring(newring);
    ideal J=ideal(imap(br,F));
    ideal I=imap(br,I);
    for (i=1;i<=m;i++)
    { J[i]=J[i]-y(i);
    }
    // We chose the weighted block order since this makes J quasi-homogeneous!
    ideal Conv=groebner(J);
    J=Conv,I;
    J=groebner(J);
    ideal vars;
    vars[n]=0;
    vars=vars,y(1..m);
    map emb=newring,vars;
    ideal G=emb(J);
    J=J-G;
    for (i=1;i<=ncols(G);i++)
    { if (J[i]<>0)
      { G[i]=0;
      }
    }
    G=compress(G);
    vars=ideal();
    for (i=2;i<=n;i++)
    { vars[i]=0;
    }
    vars=vars,y(1..m);
    emb=newring,vars;
    G=compress(emb(G));
    export(G);
    export(Conv);

    setring(smallring);
    ideal G = compress(imap(newring,G));
    export(G);

    list L;
    L[1]=smallring;
    L[2]=newring;

    dbprint( printlevel-voice+3,"
// 'rel_orbit_variety' created a list of two rings.
// If L is the name of that list, you can access
// the first ring by
//   def R = L[1]; setring R;
// (similarly for the second ring)");
    return(L);
  }
  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;
         list L = rel_orbit_variety(I,F);
         def AllR = L[2];
         setring(AllR);
         print(G);
         print(Conv);
         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:  The procedure ends with a new ring named s.
         It contains a Groebner basis
         (type <ideal>, named G) for the ideal defining the
         relative orbit variety with respect to I in the new ring.
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 define the relative orbit variety
         with respect to I.
NOTE:    This procedure is now replaced by rel_orbit_variety
         (see @ref{rel_orbit_variety}), which uses a different elemination
         order that should usually allow faster computations.
SEE ALSO: rel_orbit_variety
EXAMPLE: example relative_orbit_variety; shows an example
"
{ 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++)
    { 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++)
    { 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++)
    { vars[i]=0;
    }
    vars=vars,y(1..m);
    map emb=R,vars;
    ideal G=emb(G);
    kill vars, emb;
    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 some matrix group
RETURN:  The <ideal> defining the image under that group of the variety defined
         by I
THEORY:  rel_orbit_variety(I,F) 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
EXAMPLE: example image_of_variety; shows an example
"
{ if (nrows(F)==1)
  { def br=basering;
    int n=nvars(br);
    list L = rel_orbit_variety(I,F);
    def newring=L[2];
    setring(newring);
    ideal F=imap(br,F);
    for (int i=1;i<=n;i++)
    { F=0,F;
    }
    map emb2=newring,F;
    ideal trafo = compress(emb2(G));
    setring br;
    return(imap(newring,trafo));
  }
  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));
}
///////////////////////////////////////////////////////////////////////////////