Index: Singular/LIB/finvar.lib
===================================================================
--- Singular/LIB/finvar.lib	(revision 2244d82c47e474fe1afdfedde1939f9d8ff39771)
+++ Singular/LIB/finvar.lib	(revision 9f79aa10351ffef3887fd9394c2306929f6564f4)
@@ -1,7 +1,7 @@
-// $Id: finvar.lib,v 1.20 1999-07-30 10:37:56 obachman Exp $
+// $Id: finvar.lib,v 1.21 1999-08-27 14:14:21 Singular Exp $
 // author: Agnes Eileen Heydtmann, email:agnes@math.uni-sb.de
 // last change: 98/11/05
-//////////////////////////////////////////////////////////////////////////////
-version="$Id: finvar.lib,v 1.20 1999-07-30 10:37:56 obachman Exp $"
+///////////////////////////////////////////////////////////////////////////////
+version="$Id: finvar.lib,v 1.21 1999-08-27 14:14:21 Singular Exp $"
 info="
 LIBRARY:  finvar.lib       LIBRARY TO CALCULATE INVARIANT RINGS & MORE
@@ -43,4 +43,6 @@
  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
+ algebra_containment()             query of algebra containment
+ module_containment()              query of module containment
  orbit_variety()                   ideal of the orbit variety
  relative_orbit_variety()          ideal of a relative orbit variety
@@ -64,5 +66,5 @@
 LIB "elim.lib";
 LIB "general.lib";
-LIB "algebra.lib";
+
 ///////////////////////////////////////////////////////////////////////////////
 // Checks whether the last parameter, being a matrix, is among the previous
@@ -77,5 +79,5 @@
   return(1);
 }
-//////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
 
 proc cyclotomic (int i)
@@ -127,5 +129,4 @@
           print(cyclotomic(25));
 }
-//////////////////////////////////////////////////////////////////////////////
 
 proc group_reynolds (list #)
@@ -139,6 +140,7 @@
          characteristic is >0, minpoly==0 and the finite group non-cyclic the
          second list element is an <int> giving the lowest common multiple of
-         the matrix group elements (used in molien); in general all other list
-         elements are nxn <matrices> listing all elements of the finite group
+         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
@@ -146,14 +148,14 @@
          (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 were 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.
+         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
+{ 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")
@@ -170,5 +172,5 @@
   }
   if (typeof(#[1])<>"matrix")
-  { "ERROR:   the parameters must be a list of matrices and maybe an <int>";
+  { "ERROR:   The parameters must be a list of matrices and maybe an <int>";
     return();
   }
@@ -189,5 +191,5 @@
                                        // operator -
   matrix G(1)=#[1];                    // G(k) are elements of the group -
-  if (ch<>0 && minpoly==0 && gen_num<>1) //finding out of which order the group
+  if (ch<>0 && minpoly==0 && gen_num<>1) // finding out of which order the group
   { matrix I=diag(1,n);                // element is
     matrix TEST=G(1);
@@ -207,5 +209,5 @@
                                        // procedure
     if (not(typeof(#[j])=="matrix"))
-    { "ERROR:   the parameters must be a list of matrices and maybe an <int>";
+    { "ERROR:   The parameters must be a list of matrices and maybe an <int>";
       return();
     }
@@ -246,5 +248,5 @@
     for (m=g-j+1;m<=g;m++)
     { for (k=1;k<=i;k++)
-      { l++;
+      { l=l+1;
         matrix P(l)=G(k)*G(m);         // possible new element
       }
@@ -256,5 +258,5 @@
         g++;
         matrix G(g)=P(k);              // a new group element -
-        if (ch<>0 && minpoly==0 && i<>1) //finding out of which order the group
+        if (ch<>0 && minpoly==0 && i<>1) // finding out of which order the group
         { TEST=G(g);                   // element is
           o2=1;
@@ -282,5 +284,5 @@
     "";
   }
-  REY=transpose(REY);                  //when we evaluate the Reynolds operator
+  REY=transpose(REY);                  // when we evaluate the Reynolds operator
                                        // later on, we actually want 1xn
                                        // matrices
@@ -316,5 +318,5 @@
 }
 example
-{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";echo=2;
+{ "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;
@@ -334,23 +336,23 @@
    return(i);
 }
-//////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
 
 proc molien (list #)
 "USAGE:   molien(G1,G2,...[,ringname,lcm,flags]);
-         G1,G2,...: nxn <matrices> generating a finite matrix group, ringname:
-         a <string> giving a name for a new ring of characteristic 0 for the
-         Molien series in case of prime characteristic, lcm: an <int> giving
-         the lowest common multiple of the elements' orders in case of prime
-         characteristic, minpoly==0 and a non-cyclic group, flags: an optional
-         <intvec> with three components: if the first element is not equal to 0
-         characteristic 0 is simulated, i.e. the Molien series is computed as
-         if the base field were characteristic 0 (the user must choose a field
-         of large 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)
+         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
@@ -391,5 +393,5 @@
       { 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";
+        { "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();
@@ -406,5 +408,5 @@
       else
       { if (typeof(#[size(#)-1])<>"string")
-        { "ERROR:   in characteristic p>0 a <string> must be given for the name of a new";
+        { "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();
@@ -670,5 +672,5 @@
 //           { break;
 //           }
-//           k++;
+//           k=k+1;
 //         }
         setring `newring`;
@@ -902,13 +904,13 @@
          (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 were 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.
+         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
 "
@@ -1096,5 +1098,5 @@
     { l=0;                             // l is the number of products we get in
                                        // one going
-      for (m=g-j+1;m<=g;m++)
+      for (m=g-j+1;m<=g;m=m+1)
       { for (k=1;k<=i;k++)
         { l++;
@@ -1242,5 +1244,5 @@
                                        // new term of the Molien series
     i=1;
-    for (int j=2;j<=gen_num;j++)     // this loop adds the parameters to the
+    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
@@ -1434,11 +1436,9 @@
          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
+         (1+a1x+a2x^2+...+anx^n)/(1+b1x+b2x^2+...+bmx^m)=(1+(a1-b1)x+...
+         (1+b1x+b2x^2+...+bmx^m)
+         -----------------------
+            (a1-b1)x+(a2-b2)x^2+...
+            (a1-b1)x+b1(a1-b1)x^2+...
 EXAMPLE: example partial_molien; shows an example
 "
@@ -1480,5 +1480,5 @@
     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
+    for (int i=1; i<=n; i++)           // getting n terms and adding them up
     { min=lead(A(2));
       A(1)=A(1)+min;
@@ -1508,5 +1508,4 @@
          p(1);
 }
-///////////////////////////////////////////////////////////////////////////////
 
 proc evaluate_reynolds (matrix REY, ideal I)
@@ -1557,5 +1556,4 @@
          print(evaluate_reynolds(L[1],I));
 }
-///////////////////////////////////////////////////////////////////////////////
 
 proc invariant_basis (int g,list #)
@@ -1605,5 +1603,5 @@
   ideal mon=imap(br,mon);
   poly P=0;
-  for (i=m;i>=1;i=i-1)
+  for (i=m;i>=1;i--)
   { P=P+p(i)*mon[i];                   // P is the general polynomial
   }
@@ -1757,21 +1755,21 @@
 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));
+           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.
+// 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;
+  int DEGB=degBound;
   degBound=d;
   int j=ncols(mon);
@@ -2064,10 +2062,10 @@
          generated by the previously found invariants (see paper \"Generating a
          Noetherian Normalization of the Invariant Ring of a Finite Group\" by
-         Decker, Heydtmann, Schreyer (1997) to appear in JSC).
+         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.";
+  { "ERROR:   primary_char0 should only be used with rings of characteristic 0.";
     return();
   }
@@ -2207,5 +2205,5 @@
          generated by the previously found invariants (see paper \"Generating a
          Noetherian Normalization of the Invariant Ring of a Finite Group\" by
-         Decker, Heydtmann, Schreyer (1997) to appear in JSC).
+         Decker, Heydtmann, Schreyer (1998)).
 EXAMPLE: example primary_charp; shows an example
 "
@@ -2332,5 +2330,5 @@
 }
 example
-{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed to char 3)";  echo=2;
+{ "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;
@@ -2359,5 +2357,5 @@
          generated by the previously found invariants (see paper \"Generating a
          Noetherian Normalization of the Invariant Ring of a Finite Group\" by
-         Decker, Heydtmann, Schreyer (1997) to appear in JSC).
+         Decker, Heydtmann, Schreyer (1998)).
 EXAMPLE: example primary_char0_no_molien; shows an example
 "
@@ -2479,5 +2477,5 @@
 }
 example
-{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";echo=2;
+{ "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;
@@ -2486,5 +2484,4 @@
          print(l[1]);
 }
-///////////////////////////////////////////////////////////////////////////////
 
 proc primary_charp_no_molien (matrix REY, list #)
@@ -2502,5 +2499,5 @@
          generated by the previously found invariants (see paper \"Generating a
          Noetherian Normalization of the Invariant Ring of a Finite Group\" by
-         Decker, Heydtmann, Schreyer (1997) to appear in JSC).
+         Decker, Heydtmann, Schreyer (1998)).
 EXAMPLE: example primary_charp_no_molien; shows an example
 "
@@ -2623,5 +2620,5 @@
 }
 example
-{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed to char 3)"; echo=2;
+{ "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;
@@ -2643,5 +2640,5 @@
          generated by the previously found invariants (see paper \"Generating a
          Noetherian Normalization of the Invariant Ring of a Finite Group\" by
-         Decker, Heydtmann, Schreyer (1997) to appear in JSC). No Reynolds
+         Decker, Heydtmann, Schreyer (1998)). No Reynolds
          operator or Molien series is used.
 EXAMPLE: example primary_charp_without; shows an example
@@ -2764,5 +2761,5 @@
 }
 example
-{"EXAMPLE:";echo=2;
+{ "EXAMPLE:" echo=2;
          ring R=2,(x,y,z),dp;
          matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
@@ -2793,12 +2790,13 @@
 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
+         <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 (1997) to appear in JSC).
+         Decker, Heydtmann, Schreyer (1998)).
 EXAMPLE: example primary_invariants; shows an example
 "
@@ -3196,7 +3194,6 @@
          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 (1997) to appear in
-         JSC).
+         (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
 "
@@ -3320,5 +3317,5 @@
 }
 example
-{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";echo=2;
+{ "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;
@@ -3344,7 +3341,6 @@
          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 (1997) to appear in
-         JSC).
+         (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
 "
@@ -3474,5 +3470,5 @@
 }
 example
-{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed to char 3)"; echo=2;
+{ "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;
@@ -3501,7 +3497,6 @@
          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 (1997) to appear in
-         JSC).
+         (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
 "
@@ -3626,5 +3621,5 @@
 }
 example
-{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";  echo=2;
+{ "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;
@@ -3649,7 +3644,6 @@
          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 (1997) to appear in
-         JSC).
+         (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
 "
@@ -3774,5 +3768,5 @@
 }
 example
-{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed to char 3)"; echo=2;
+{ "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;
@@ -3794,7 +3788,7 @@
          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 (1997) to appear in
-         JSC). No Reynolds operator or Molien series is used.
+         (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
 "
@@ -3926,5 +3920,5 @@
 }
 example
-{  "EXAMPLE:"; echo=2;
+{ "EXAMPLE:"; echo=2;
          ring R=2,(x,y,z),dp;
          matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
@@ -3963,7 +3957,6 @@
          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 (1997) to appear in
-         JSC).
+         (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
 "
@@ -4046,5 +4039,10 @@
           { "  We are dealing with the non-modular case.";
           }
-          molien(L[2..size(L)],newring,intvec(0,interval,v));
+          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);
@@ -4115,5 +4113,10 @@
     list L=group_reynolds(#[1..gen_num],v);
     string newring="aksldfalkdsflkj";
-    molien(L[2..size(L)],newring,intvec(1,interval,v));
+    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);
@@ -4131,5 +4134,5 @@
 }
 example
-{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";  echo=2;
+{ "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;
@@ -4201,5 +4204,5 @@
 }
 example
-{  "EXAMPLE:"; echo=2;
+{ "EXAMPLE:"; echo=2;
          intvec dv=5,5,5,10,10;
          print(power_products(dv,10));
@@ -4224,7 +4227,7 @@
          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 s. t. they are linearly independent modulo
-         primary invariants (see paper \"Some Algorithms in Invariant Theory of
-         Finite Groups\" by Kemper and Steel (1997)).
+         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_char0; shows an example
 "
@@ -4425,6 +4428,6 @@
          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()
+         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>)
@@ -4433,7 +4436,7 @@
          monomials) of the basering modulo primary invariants, mapping those
          to invariants with the Reynolds operator and using these images or
-         their power products s. t. they are linearly independent modulo the
-         primary invariants (see paper \"Some Algorithms in Invariant Theory of
-         Finite Groups\" by Kemper and Steel (1997)).
+         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
 "
@@ -4629,5 +4632,5 @@
 }
 example
-{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed to char 3)"; echo=2;
+{ "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;
@@ -4653,6 +4656,6 @@
 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
+         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
@@ -4809,7 +4812,7 @@
 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 s.t. they are linearly independent modulo the
+         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)).
@@ -5003,6 +5006,5 @@
 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, to
-         appear in JSC).
+         Rings of Finite Groups over Arbitrary Fields\" by Kemper (1996).
 EXAMPLE: example secondary_not_cohen_macaulay; shows an example
 "
@@ -5080,5 +5082,6 @@
     "";
   }
-  matrix trivialS=secondary_charp(P,REY,"alskdfalkdsj",v);
+  matrix trivialS, trivialSI=secondary_charp(P,REY,"alskdfalkdsj",v);
+  kill trivialSI;
   kill alskdfalkdsj;
   // now we have those secondary invariants
@@ -5096,5 +5099,5 @@
   ideal B;
   for (i=1;i<=gen_num;i++)
-  { B=ideal(matrix(maxideal(1))*transpose(#[i]));   // image of the various
+  { B=ideal(matrix(maxideal(1))*transpose(#[i])); // image of the various
                                        // variables under the i-th generator -
     f=br,B;                            // the corresponding mapping -
@@ -5135,5 +5138,5 @@
                                        // invariants
   for (i=1; i<=m;i++)
-  { S[1,i]=S[1,i]/leadcoef(S[1,i]); // making elements nice
+  { S[1,i]=S[1,i]/leadcoef(S[1,i]);    // making elements nice
   }
   S=sort(ideal(S))[1];
@@ -5155,5 +5158,5 @@
 }
 example
-{  "EXAMPLE:"; echo=2;
+{ "EXAMPLE:"; echo=2;
            ring R=2,(x,y,z),dp;
            matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
@@ -5177,11 +5180,11 @@
          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
+         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
+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
@@ -5190,7 +5193,7 @@
 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
+         generated by the previously found invariants (see \"Generating a
          Noetherian Normalization of the Invariant Ring of a Finite Group\" by
-         Decker, Heydtmann, Schreyer (1997) to appear in JSC). In the
+         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
@@ -5200,6 +5203,5 @@
          Theory of Finite Groups\" by Kemper and Steel (1997)). In the modular
          case they are generated according to \"Generating Invariant Rings of
-         Finite Groups over Arbitrary Fields\" by Kemper (1996, to appear in
-         JSC).
+         Finite Groups over Arbitrary Fields\" by Kemper (1996).
 EXAMPLE: example invariant_ring; shows an example
 "
@@ -5256,5 +5258,10 @@
           { "  We are dealing with the non-modular case.";
           }
-          molien(L[3..size(L)],newring,L[2],intvec(0,interval,v));
+          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);
@@ -5339,5 +5346,10 @@
     list L=group_reynolds(#[1..gen_num],v);
     string newring="aksldfalkdsflkj";
-    molien(L[2..size(L)],newring,intvec(1,interval,v));
+    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);
@@ -5360,6 +5372,5 @@
 }
 example
-{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
-  echo=2;
+{ "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;
@@ -5376,21 +5387,21 @@
          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
+         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, it
+         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 occurs not among the coefficients)
-RETURN:  primary and secondary invariants (both of type matrix) generating the
+         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
@@ -5482,5 +5493,10 @@
           { "  We are dealing with the non-modular case.";
           }
-          molien(L[3..size(L)],newring,L[2],intvec(0,interval,v));
+          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);
@@ -5519,5 +5535,5 @@
         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);
+        return(P,S);
       }
     }
@@ -5565,5 +5581,10 @@
     list L=group_reynolds(#[1..gen_num],v);
     string newring="aksldfalkdsflkj";
-    molien(L[2..size(L)],newring,intvec(1,v));
+    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);
@@ -5593,4 +5614,125 @@
          print(S);
          print(IS);
+}
+///////////////////////////////////////////////////////////////////////////////
+
+proc algebra_containment (poly p, ideal A)
+"USAGE:   algebra_containment(p,A); p <poly>, A <ideal> giving generators, say
+         f1,...,fm, subalgebra of the basering
+RETURN:  1 (TRUE) (type <int>) if p is contained in the subalgebra K[f1,...,fm]
+         0 (FALSE) (type <int>) if p is not contained
+DISPLAY: a representation of p in terms of algebra generators fi=y(i) if p
+         is contained in the subalgebra, provided printlevel>=1
+         (default: printlevel=0)
+THEORY:  The ideal of algebraic relations of the algebra generators f1,...,fm
+         given by A is computed introducing new variables y(i) and the product
+         order with x(i) >> y(i) where x(i) are the variable of the basering.
+         p reduces to a polynomial only in the y(i) <=> p is contained in the
+         subring generated by the polynomials f1,...,fm in A.
+EXAMPLE: example algebra_containment; shows an example
+"
+{ int DEGB = degBound;
+  degBound=0;
+  def br=basering;
+  int n=nvars(br);
+  int m=ncols(A);
+  string mp=string(minpoly);
+  execute "ring R=("+charstr(br)+"),(x(1..n),y(1..m)),(dp(n),dp(m));";
+  execute "minpoly=number("+mp+");";
+  ideal vars=x(1..n);
+  map emb=br,vars;
+  ideal A=ideal(emb(A));
+  poly check=emb(p);
+  for (int i=1;i<=m;i++)
+  { A[i]=A[i]-y(i);
+  }
+  A=std(A);
+  check=reduce(check,A);
+  degBound = DEGB;
+  if( sum(leadexp(check),1..n) == 0 )
+  { dbprint(printlevel-voice+2,"// "+string(check));
+    return(1);
+  }
+  else
+  { return(0);
+  }
+}
+example
+{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2;
+         ring R=0,(x,y,z),dp;
+         matrix A[1][7]=x2+y2,z2,x4+y4,1,x2z-1y2z,xyz,x3y-1xy3;
+         poly p1=x10z3-x8y2z3+2x6y4z3-2x4y6z3+x2y8z3-y10z3+x6z4+3x4y2z4+3x2y4z4+y6z4;
+         printlevel=2;
+         algebra_containment(p1,A);
+         poly p2=z;
+         algebra_containment(p2,A);
+}
+///////////////////////////////////////////////////////////////////////////////
+
+proc module_containment(poly p, matrix P, matrix S)
+"USAGE:  module_containment(p,P,S); p <poly>, P <ideal>, generators of an
+         algebra, S <ideal>, generators of a module over the algebra generated
+         by P
+ASSUME:  ncols(P)=number of variables in the basering and the generators in P
+RETURNS: 1 (TRUE) (type <int>) if p is contained in the module
+         0 (FALSE) type <int>) if p is not contained
+DISPLAY: the representation of p in terms of algebra generators P[i]=z(i) and
+         module generators S[j]=y(j) if p is contained in the module
+THEORY:  The ideal of algebraic relations of all the generators p1,...,pn,
+         s1,...,st given by P and S is computed introducing new variables y(j),
+         z(i) and the product order: x^a*y^b*z^c > x^d*y^e*z^f if x^a > x^d
+         with respect to the lp ordering or else if z^c > z^f with respect to
+         the dp ordering or else if y^b > y^e with respect to the lp ordering
+         again. p reduces to a polynomial only in the y(j) and z(i), linear in
+         the z(i) <=> p is contained in the module.
+EXAMPLE: example module_containment; shows an example
+"
+{ def br=basering;
+  int DEGB=degBound;
+  degBound=0;
+  int n=nvars(br);
+  if (ncols(P)==n)
+  { int m=ncols(S);
+    string mp=string(minpoly);
+    execute "ring R=("+charstr(br)+"),(x(1..n),y(1..m),z(1..n)),(lp(n),dp(m),lp(n));";
+    execute "minpoly=number("+mp+");";
+    ideal vars=x(1..n);
+    map emb=br,vars;
+    matrix P=emb(P);
+    matrix S=emb(S);
+    ideal check=emb(p);
+    ideal I;
+    for (int i=1;i<=m;i++)
+    { I[i]=S[1,i]-y(i);
+    }
+    for (i=1;i<=n;i++)
+    { I[m+i]=P[1,i]-z(i);
+    }
+    I=std(I);
+    check[1]=reduce(check[1],I);
+    I=elim(check,1,n);                 // checking whether all variables from
+    if (I[1]<>0)                       // the former ring have disappeared
+    { dbprint(printlevel-voice+2,"// "+string(check));
+      return(1);
+    }
+    else
+    { return(0);
+    }
+  }
+  else
+   { "ERROR: the first ideal must have nvars(basering) entries";
+    return();
+  }
+}
+example
+{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2;
+         ring R=0,(x,y,z),dp;
+         matrix P[1][3]=x2+y2,z2,x4+y4;
+         matrix S[1][4]=1,x2z-1y2z,xyz,x3y-1xy3;
+         poly p1=x10z3-x8y2z3+2x6y4z3-2x4y6z3+x2y8z3-y10z3+x6z4+3x4y2z4+3x2y4z4+y6z4;
+         printlevel=2;
+         module_containment(p1,P,S);
+         poly p2=z;
+         module_containment(p2,P,S);
 }
 ///////////////////////////////////////////////////////////////////////////////
@@ -5689,5 +5831,5 @@
     J=std(J);
     ideal vars;
-    //for (i=1;i<=n;i++)
+    //for (i=1;i<=n;i=i+1)
     //{ vars[i]=0;
     //}
@@ -5772,3 +5914,2 @@
          print(image_of_variety(I,F));
 }
-///////////////////////////////////////////////////////////////////////////////
