Changeset e52b9f in git

Timestamp:

Jan 22, 2007, 2:08:25 PM (17 years ago)

Author:

Hans Schönemann <hannes@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

8223cd635042d178325b09664e9e3c94114b9a5c

Parents:

05a31dd1bac7c8886acb526ecf5503742a948378

Message:

*king: new secondary_charp, irred_secondary_char0, secondary_char0


git-svn-id: file:///usr/local/Singular/svn/trunk@9742 2c84dea3-7e68-4137-9b89-c4e89433aadc

File:

• Singular/LIB/finvar.lib (modified) (90 diffs)

Singular/LIB/finvar.lib

@@ -1 +1 @@
 ///////////////////////////////////////////////////////////////////////////////
-version="$Id: finvar.lib,v 1.51 2006-08-02 15:40:47 Singular Exp $"
+version="$Id: finvar.lib,v 1.52 2007-01-22 13:08:25 Singular Exp $"
 category="Invariant theory";
 info="
@@ -40 +40 @@
  irred_secondary_char0()           irreducible secondary invariants in char 0
  secondary_charp()                 secondary invariants in char p
+ secondary_charp_degbound()        finds all (irred.) secondary invariants in char p
+                                   up to a given degree bound, without using Molien
+                                   series or Reynolds operator
  secondary_no_molien()             secondary invariants, without Molien series
+                                   but with Reynolds operator
  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
@@ -54 +58 @@
 // 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
+// int_number_map                  integers 1..q are mapped to q field elements
 // search                          searches a number of p.i., char 0
 // p_search                        searches a number of p.i., char p
@@ -90 +94 @@
 "
 { if (i<=0)
-  { "ERROR:   the input should be > 0.";
-    return();
+  { ERROR("ERROR:   the input should be > 0.");
   }
   poly v1=var(1);
@@ -138 +141 @@
 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 is non-cyclic 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
@@ -148 +151 @@
          (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 are found in one run. Additionally each time a new
+         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
@@ -161 +164 @@
   if (typeof(#[size(#)])=="int")
   { if (size(#)==1)
-    { "ERROR:   there are no matrices given among the parameters";
-      return();
+    { ERROR("ERROR:   there are no matrices given among the parameters");
     }
     int v=#[size(#)];
@@ -172 +174 @@
   }
   if (typeof(#[1])<>"matrix")
-  { "ERROR:   The parameters must be a list of matrices and maybe an <int>";
-    return();
+  { ERROR("ERROR:   The parameters must be a list of matrices and maybe an <int>");
   }
   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();
+  { ERROR("ERROR:   the number of variables of the basering needs to be the same"+newline+
+    "         as the dimension of the matrices");
   }
   if (n<>ncols(#[1]))
-  { "ERROR:   matrices need to be square and of the same dimensions";
-    return();
+  { ERROR("ERROR:   matrices need to be square and of the same dimensions");
   }
   matrix vars=matrix(maxideal(1));     // creating an nx1-matrix containing the
@@ -209 +208 @@
                                        // procedure
     if (not(typeof(#[j])=="matrix"))
-    { "ERROR:   The parameters must be a list of matrices and maybe an <int>";
-      return();
+    { ERROR("ERROR:   The parameters must be a list of matrices and maybe an <int>");
     }
     if ((n!=nrows(#[j])) or (n!=ncols(#[j])))
-    { "ERROR:   matrices need to be square and of the same dimensions";
-       return();
+    { ERROR("ERROR:   matrices need to be square and of the same dimensions");
     }
     if (unique(G(1..i),#[j]))
@@ -358 +355 @@
          elements generated by group_reynolds(), lcm is the second return value
          of group_reynolds()
-RETURN:  in case of characteristic 0 a 1x2 <matrix> giving numerator and
+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
@@ -379 +376 @@
     { 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();
+      { ERROR("ERROR:   the second component of <intvec> should be >=0");
       }
       int interval=#[size(#)][2];
@@ -386 +382 @@
     }
     else
-    { "ERROR:   <intvec> should have three components";
-      return();
+    { ERROR("ERROR:   <intvec> should have three components");
     }
     if (ch<>0)
@@ -393 +388 @@
       { 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();
+        { ERROR("ERROR:   In characteristic p>0 a <string> must be given for the name of a new"+newline+
+          "         ring where the Molien series can be stored");
         }
         else
         { if (#[size(#)-2]=="")
-          { "ERROR:   <string> may not be empty";
-            return();
+          { ERROR("ERROR:   <string> may not be empty");
           }
           string newring=#[size(#)-2];
@@ -408 +401 @@
       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";
+        { ERROR("ERROR:   In characteristic p>0 a <string> must be given for the name of a new"+newline+
+          "         ring where the Molien series can be stored");
           return();
         }
         else
         { if (#[size(#)-1]=="")
-          { "ERROR:   <string> may not be empty";
-            return();
+          { ERROR("ERROR:   <string> may not be empty");
           }
           string newring=#[size(#)-1];
@@ -435 +427 @@
       { 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();
+        { ERROR("ERROR:   in characteristic p>0 a <string> must be given for the name of a new"+newline+
+          "         ring where the Molien series can be stored");
         }
         else
         { if (#[size(#)-1]=="")
-          { "ERROR:   <string> may not be empty";
-            return();
+          { ERROR("ERROR:   <string> may not be empty");
           }
           string newring=#[size(#)-1];
@@ -450 +440 @@
       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();
+        { ERROR("ERROR:   in characteristic p>0 a <string> must be given for the name of a new"+newline+
+          "         ring where the Molien series can be stored");
         }
         else
         { if (#[size(#)]=="")
-          { "ERROR:   <string> may not be empty";
-            return();
+          { ERROR("ERROR:   <string> may not be empty");
           }
           string newring=#[size(#)];
@@ -471 +459 @@
   if (ch<>0)
   { if ((g/r)*r<>g)
-   { "ERROR:   <int> should divide the group order."
-      return();
+   { ERROR("ERROR:   <int> should divide the group order.");
     }
   }
@@ -507 +494 @@
  //----------------------------------------------------------------------------
   if (not(typeof(#[1])=="matrix"))
-  { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
-    return();
+  { ERROR("ERROR:   the parameters must be a list of matrices and maybe an <intvec>");
   }
   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();
+  { ERROR("ERROR:   the number of variables of the basering needs to be the same"+newline+
+    "         as the dimension of the square matrices");
   }
   if (v && voice<>2)
@@ -540 +525 @@
     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();
+      { ERROR("ERROR:   the parameters must be a list of matrices and maybe an <intvec>");
       }
       if ((n<>nrows(#[j])) or (n<>ncols(#[j])))
-      { "ERROR:   matrices need to be square and of the same dimensions";
-         return();
+      { ERROR("ERROR:   matrices need to be square and of the same dimensions");
       }
       p=det(I-v1*#[j]);                // denominator of new term -
@@ -637 +620 @@
       { setring br;
         if (not(typeof(#[i])=="matrix"))
-        { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
-          return();
+        { ERROR("ERROR:   the parameters must be a list of matrices and maybe an <intvec>");
         }
         if ((n<>nrows(#[i])) or (n<>ncols(#[i])))
-        { "ERROR:   matrices need to be square and of the same dimensions";
-           return();
+        { ERROR("ERROR:   matrices need to be square and of the same dimensions");
         }
         setring bre;
@@ -798 +779 @@
     { setring br;
       if (not(typeof(#[i])=="matrix"))
-      { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
-        return();
+      { ERROR("ERROR:   the parameters must be a list of matrices and maybe an <intvec>");
       }
       if ((n<>nrows(#[i])) or (n<>ncols(#[i])))
-      { "ERROR:   matrices need to be square and of the same dimensions";
-         return();
+      { ERROR("ERROR:   matrices need to be square and of the same dimensions");
       }
       string stM(i)=string(#[i]);
@@ -908 +887 @@
          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 term 1/det(1-xE)
+         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 numerator and denominator of the expanded
+         returned matrix gives enumerator and denominator of the expanded
          version where common factors have been canceled.
 EXAMPLE: example reynolds_molien; shows an example
@@ -925 +904 @@
     { 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();
+      { ERROR("ERROR:   the second component of the <intvec> should be >=0");
       }
       int interval=#[size(#)][2];
@@ -932 +910 @@
     }
     else
-    { "ERROR:   <intvec> should have three components";
-      return();
+    { ERROR("ERROR:   <intvec> should have three components");
     }
     if (ch<>0)
     { if (typeof(#[size(#)-1])<>"string")
-      { "ERROR:   in characteristic p a <string> must be given for the name";
+      { ERROR("ERROR:   in characteristic p a <string> must be given for the name");
         "         of a new ring where the Molien series can be stored";
         return();
@@ -943 +920 @@
       else
       { if (#[size(#)-1]=="")
-        { "ERROR:   <string> may not be empty";
+        { ERROR("ERROR:   <string> may not be empty");
           return();
         }
@@ -960 +937 @@
     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();
+      { ERROR("ERROR:   in characteristic p a <string> must be given for the name"+newline+
+        "         of a new ring where the Molien series can be stored");
       }
       else
       { if (#[size(#)]=="")
-        { "ERROR:   <string> may not be empty";
-          return();
+        { ERROR("ERROR:   <string> may not be empty");
         }
         string newring=#[size(#)];
@@ -1018 +993 @@
   if (ch==0)
   { if (typeof(#[1])<>"matrix")
-    { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
-      return();
+    { ERROR("ERROR:   the parameters must be a list of matrices and maybe an <intvec>");
     }
     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();
+    { ERROR("ERROR:   the number of variables of the basering needs to be the same"+newline+
+      "         as the dimension of the matrices");
     }
     if (n<>ncols(#[1]))
-    { "ERROR:   matrices need to be square and of the same dimensions";
-      return();
+    { ERROR("ERROR:   matrices need to be square and of the same dimensions");
     }
     matrix vars=matrix(maxideal(1));   // creating an nx1-matrix containing the
@@ -1055 +1027 @@
                                        // procedure
       if (not(typeof(#[j])=="matrix"))
-      { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
-        return();
+      { ERROR("ERROR:   the parameters must be a list of matrices and maybe an <intvec>");
       }
       if ((n!=nrows(#[j])) or (n!=ncols(#[j])))
-      { "ERROR:   matrices need to be square and of the same dimensions";
-         return();
+      { ERROR("ERROR:   matrices need to be square and of the same dimensions");
       }
       if (unique(G(1..i),#[j]))
@@ -1651 +1621 @@
            print(invariant_basis(2,A));
 }
+///////////////////////////////////////////////////////////////////////////////
+
+stat proc invariant_basis_modS (int g,ideal sS, list #)
+"USAGE:   invariant_basis(g,sP,G1,G2,...);
+         g: an <int> indicating of which degree (>0) the homogeneous basis
+         shoud be,
+         sS: normal forms of subsequently found sec. inv. in degree g,
+         G1,G2,...: <matrices> generating a finite matrix group
+RETURNS: 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.
+"
+{ if (g<=0)
+  { "ERROR:   the first parameter should be > 0";
+    return();
+  }
+  def br=basering;
+  ideal mon=sort(kbase(sS,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);
+}
+
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -2064 +2120 @@
 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 program if v does not
+DISPLAY: information about the various stages of the programme if v does not
          equal 0
 RETURN:  primary invariants (type <matrix>) of the invariant ring
@@ -2163 +2219 @@
       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
+      }
+      else                             // i.e. we can take all of B
       { for(j=i+1;j<=i+dif;j++)
         { CI=CI+ideal(var(j)^d);
@@ -2207 +2263 @@
          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 program if v does not
+DISPLAY: information about the various stages of the programme if v does not
          equal 0
 RETURN:  primary invariants (type <matrix>) of the invariant ring
@@ -2249 +2305 @@
   }
   if (typeof(`ring_name`)<>"ring")
-  { "ERROR:   Second parameter ought to the name of a ring where the Molien";
-    "         is stored.";
+  { "ERROR:   Second parameter ought to be ring where ";
+    "         the Molien series is stored.";
     return();
   }
@@ -2360 +2416 @@
          <int>
 ASSUME:  REY is the first return value of group_reynolds or reynolds_molien
-DISPLAY: information about the various stages of the program if v does not
+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
@@ -2503 +2559 @@
          <int>
 ASSUME:  REY is the first return value of group_reynolds or reynolds_molien
-DISPLAY: information about the various stages of the program if v does not
+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
@@ -2644 +2700 @@
          G1,G2,...: <matrices> generating a finite matrix group, v: an optional
          <int>
-DISPLAY: information about the various stages of the program if v does not
+DISPLAY: information about the various stages of the programme if v does not
          equal 0
 RETURN:  primary invariants (type <matrix>) of the invariant ring
@@ -2658 +2714 @@
  //--------------------- checking input and setting verbose mode --------------
   if (char(basering)==0)
-  { "ERROR:   primary_charp_without should only be used with rings of";
+  { "ERROR:   primary_charp_without should not be used with rings of";
     "         characteristic 0.";
     return();
@@ -2704 +2760 @@
                                        // space of invariants of degree d,
                                        // newdim: dimension the ideal generated
-                                       // the primary invariants plus basis
+                                       // by the primary invariants plus basis
                                        // elements, dif=n-i-newdim, i.e. the
                                        // number of new primary invairants that
@@ -2784 +2840 @@
          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 program attempts to compute the Molien series and
-         Reynolds operator, if it equals 1, the program is told that the
+         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
+         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
@@ -2797 +2853 @@
          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 program if the third
+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
@@ -2892 +2948 @@
         else
         { if (v)
-          { "  Since it is impossible for this program to calculate the Molien series for";
+          { "  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.";
@@ -3199 +3255 @@
 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 program if v does not
+DISPLAY: information about the various stages of the programme if v does not
          equal 0
 RETURN:  primary invariants (type <matrix>) of the invariant ring
@@ -3346 +3402 @@
          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 program if v does not
+DISPLAY: information about the various stages of the programme if v does not
          equal 0
 RETURN:  primary invariants (type <matrix>) of the invariant ring
@@ -3499 +3555 @@
          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 program if v does not
+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
@@ -3646 +3702 @@
          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 program if v does not
+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
@@ -3792 +3848 @@
          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 program if v does not
+DISPLAY: information about the various stages of the programme if v does not
          equal 0
 RETURN:  primary invariants (type <matrix>) of the invariant ring
@@ -3943 +3999 @@
          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 program
+         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 program is told that the Molien series should not be
+         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
+         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
@@ -3957 +4013 @@
          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 program if the third
+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
@@ -4060 +4116 @@
         else
         { if (v)
-          { "  Since it is impossible for this program to calculate the Molien series for";
+          { "  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.";
@@ -4424 +4480 @@
 proc secondary_char0 (matrix P, matrix REY, matrix M, list #)
 "
-USAGE:    secondary_char0(P,REY,M[,v][,\"AH\"]);
+USAGE:    secondary_char0(P,REY,M[,v][,\"old\"]);
 @*          P: a 1xn <matrix> with homogeneous primary invariants, where
                n is the number of variables of the basering;
@@ -4432 +4488 @@
                series;
 @*          v: an optional <int>;
-@*          \"AH\": if this string occurs as (optional) parameter, then an
+@*          \"old\": 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
+ASSUME:   The characteristic of basering is zero;
+          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()
@@ -4450 +4507 @@
           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.
+          Steel (1997)) using Groebner basis techniques. 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
+          Specifically, the output of secondary_char0(P,REY,M,\"old\") 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.
+@*        There are three internal parameters \"pieces\", \"MonStep\" and \"IrrSwitch\".
+          The default values of the parameters should be fine in most cases. However,
+          in some cases, different values may provide a better balance of memory
+          consumption (smaller values) and speed (bigger values).
 SEE ALSO: irred_secondary_char0;
 EXAMPLE:  example secondary_char0; shows an example
 "
-{ def br=basering;
+{
+ //----------  Internal parameters, whose choice might improve the performance -----
+  int pieces = 15;   // For generating reducible secondaries, blocks of #pieces# secondaries
+                     // are formed.
+                     // If this parameter is small, the memory consumption will decrease.
+  int MonStep = 15;  // The Reynolds operator is applied to blocks of #MonStep# monomials.
+                     // If this parameter is small, the memory consumption will decrease.
+  int IrrSwitch = 7; // Up to degree #IrrSwitch#-1, we use a method for generating
+                     // irreducible secondaries that tends to produce a smaller output
+                     // (which is good for subsequent computations). However, this method
+                     // needs to much memory in high degrees, and so we use a sparser method
+                     // from degree #IrrSwitch# on.
+
+  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();
+  { "ERROR:   secondary_char0 can only be used with rings of characteristic 0.";
+    "         Try secondary_charp";
+      return();
   }
   int i;
@@ -4484 +4556 @@
   if (size(#)>0)
   { if (typeof(#[size(#)])=="string")
-    { if (#[size(#)]=="AH")
+    { if (#[size(#)]=="old")
       { if (typeof(#[1])=="int")
         { matrix S,IS = old_secondary_char0(P,REY,M,#[1]);
@@ -4495 +4567 @@
       }
       else
-      { "ERROR:   If the last optional parameter is a string, it should be \"AH\".";
+      { "ERROR:   If the last optional parameter is a string, it should be \"old\".";
         return(matrix(ideal()),matrix(ideal()));
       }
@@ -4561 +4633 @@
   option(redSB);
   ideal sP = groebner(ideal(P));
+                           // This is the only explicit Groebner basis computation!
   ideal Reductor,SaveRed;  // sP union Reductor is a Groebner basis up to degree i-1
   int SizeSave;
@@ -4576 +4649 @@
   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
+  ideal mon,B,IS;              // IS will contain all irr. sec. inv.
 
   for (i=1;i<m;i++)
@@ -4587 +4657 @@
       attrib(SSort[i][k],"size",0);
     }
-   }
+  }
 //--------------------- generating secondary invariants ----------------------
   for (i=2;i<=m;i++)                   // going through dimmat -
@@ -4642 +4712 @@
                 }
               }
-              ProdCand = Mul1*Mul2;
+//              if (i<10)
+                ProdCand = simplify(Mul1*Mul2,4);
+//              else { ProdCand = Mul1*Mul2;}
               ReducedCandidates = reduce(ProdCand,sP);
                   // sP union SaveRed union Reductor is a homogeneous Groebner basis
@@ -4654 +4726 @@
                              // invariants and previously computed secondary invariants.
                              // Otherwise ProdCand[ii] can be taken as secondary invariant.
-          if (size(Indicator)<>0)
+              if (size(Indicator)<>0)
               { for (ii=1;ii<=ncols(ProdCand);ii++)     // going through all the power products
-                { helpP = reduce(Indicator[ii],Reductor);
+                { helpP = Indicator[ii];
                   if (helpP <> 0)
                   { counter++;
@@ -4667 +4739 @@
                     }
                     if (int(dimmat[i,1])<>counter)
-                    { helpint++;
-                      Reductor[helpint] = helpP;
+                    { Reductor = ideal(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
@@ -4676 +4747 @@
                         // 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);
-                      }
+                      Indicator=reduce(Indicator,Reductor);
+                      SizeSave++;
+                      SaveRed[SizeSave] = helpP;
+                      attrib(SaveRed, "isSB",1);
                     }
                     else
@@ -4694 +4758 @@
                 }
               }
-              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);
+//              attrib(SaveRed, "isSB",1);
             }
           }
@@ -4712 +4769 @@
         { "There are irreducible secondary invariants in degree ", i-1," !!";
         }
-        B=sort_of_invariant_basis(sP,REY,i-1,int(dimmat[i,1])*6);
+        if (i>IrrSwitch)  // use sparse algorithm
+        { mon = kbase(sP,i-1);
+//          mon = ideal(mon[ncols(mon)..1]);
+          if (counter==0 && ncols(mon)==int(dimmat[i,1])) // then we can use all of mon
+          { B=normalize(evaluate_reynolds(REY,mon));
+            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 mon -
+            // Compare the comments on the computation of reducible sec. inv.!
+            while (int(dimmat[i,1])<>counter)
+            { if ((j mod MonStep) == 0)
+              { if ((j+MonStep) <= ncols(mon))
+                { B[j+1..j+MonStep] = normalize(evaluate_reynolds(REY,ideal(mon[j+1..j+MonStep])));
+                  ReducedCandidates[j+1..j+MonStep] = reduce(ideal(B[j+1..j+MonStep]),sP);
+                  Indicator[j+1..j+MonStep] = reduce(ideal(ReducedCandidates[j+1..j+MonStep]),SaveRed);
+                }
+                else
+                { B[j+1..ncols(mon)] = normalize(evaluate_reynolds(REY,ideal(mon[j+1..ncols(mon)])));
+                  ReducedCandidates[j+1..ncols(mon)] = reduce(ideal(B[j+1..ncols(mon)]),sP);
+                  Indicator[j+1..ncols(mon)] = reduce(ideal(ReducedCandidates[j+1..ncols(mon)]),SaveRed);
+                }
+              }
+              j++;
+              helpP = Indicator[j];
+              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]);
+                }
+                Reductor = ideal(helpP);
+                attrib(Reductor, "isSB",1);
+                Indicator=reduce(Indicator,Reductor);
+                SizeSave++;
+                SaveRed[SizeSave] = helpP;
+                attrib(SaveRed, "isSB",1);
+              }
+              B[j]=0;
+              ReducedCandidates[j]=0;
+              Indicator[j]=0;
+            }
+          }
+        }  // if i>IrrSwitch
+        else  // use fast algorithm
+        { 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;
+                                       // independent
+          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
-          { 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];
+          { 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 = Indicator[j];
+              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]);
+                }
+                Reductor = ideal(helpP);
+                attrib(Reductor, "isSB",1);
+                Indicator=reduce(Indicator,Reductor);
+                SizeSave++;
+                SaveRed[SizeSave] = helpP;
+                attrib(SaveRed, "isSB",1);
               }
-              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);
-              }
+              B[j]=0;
+              ReducedCandidates[j]=0;
+              Indicator[j]=0;
             }
           }
-        }
-      }
+        } // i<=IrrSwitch
+      } // Computation of irreducible secondaries
       if (v)
       { "";
       }
-    }
-  }
+    } // if (int(dimmat[i,1])<>0)
+  }  // for i
   if (v)
   { "  We're done!";
@@ -4784 +4901 @@
   degBound = dgb;
 
-// Prepare return:
+  // Prepare return:
   int TotalNumber;
   for (k=1;k<=m;k++)
@@ -4806 +4923 @@
 }
 
-
 example
 { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2;
@@ -4816 +4932 @@
          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 #)
@@ -4841 +4948 @@
 RETURN:   Irreducible homogeneous secondary invariants of the invariant ring
           (type <matrix>)
-ASSUME:   REY is the 1st return value of group_reynolds(), reynolds_molien() or
+ASSUME:   We are in the non-modular case, i.e., the characteristic of the basering
+          does not divide the group order;
+          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
@@ -4852 +4961 @@
           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
+          Steel (1997)) using Groebner basis techniques. 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.
+@*        There are three internal parameters \"pieces\", \"MonStep\" and \"IrrSwitch\".
+          The default values of the parameters should be fine in most cases. However,
+          in some cases, different values may provide a better balance of memory
+          consumption (smaller values) and speed (bigger values).
 SEE ALSO: secondary_char0
 KEYWORDS: irreducible secondary invariant
 EXAMPLE:  example irred_secondary_char0; shows an example
 "
-{ def br=basering;
-
+{ //----------  Internal parameters, whose choice might improve the performance -----
+  int pieces = 15;   // For generating reducible secondaries, blocks of #pieces# secondaries
+                     // are formed.
+                     // If this parameter is small, the memory consumption will decrease.
+  int MonStep = 15;  // The Reynolds operator is applied to blocks of #MonStep# monomials.
+                     // If this parameter is small, the memory consumption will decrease.
+  int IrrSwitch = 7; // Up to degree #IrrSwitch#-1, we use a method for generating
+                     // irreducible secondaries that tends to produce a smaller output
+                     // (which is good for subsequent computations). However, this method
+                     // needs to much memory in high degrees, and so we use a sparser method
+                     // from degree #IrrSwitch# on.
+
+  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();
+  { "ERROR:   irred_secondary_char0 can only be used with rings of characteristic 0.";
+    "         Try irred_secondary_charp";
+      return();
   }
   int i;
@@ -4958 +5080 @@
   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
+  ideal mon,B,IS;              // IS will contain all irr. sec. inv.
 
   for (i=1;i<m;i++)
@@ -5027 +5146 @@
                 }
               }
-              ProdCand = Mul1*Mul2;
+//              if (i<10)
+                ProdCand = simplify(Mul1*Mul2,4);
+//              else { ProdCand = Mul1*Mul2;}
               ReducedCandidates = reduce(ProdCand,sP);
                   // sP union SaveRed union Reductor is a homogeneous Groebner basis
@@ -5042 +5163 @@
               if (size(Indicator)<>0)
               { for (ii=1;ii<=ncols(ProdCand);ii++)     // going through all the power products
-                { helpP = reduce(Indicator[ii],Reductor);
+                { helpP = Indicator[ii];
                   if (helpP <> 0)
                   { counter++;
@@ -5054 +5175 @@
                     }
                     if (int(dimmat[i,1])<>counter)
-                    { helpint++;
-                      Reductor[helpint] = helpP;
+                    { Reductor = ideal(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
@@ -5063 +5183 @@
                         // 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);
-                      }
+                      Indicator=reduce(Indicator,Reductor);
+                      SizeSave++;
+                      SaveRed[SizeSave] = helpP;
+                      attrib(SaveRed, "isSB",1);
+                      attrib(Reductor, "isSB",1);
                     }
                     else
@@ -5081 +5195 @@
                 }
               }
-              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);
@@ -5099 +5207 @@
         { "There are irreducible secondary invariants in degree ", i-1," !!";
         }
-        B=sort_of_invariant_basis(sP,REY,i-1,int(dimmat[i,1])*6);
+        if (i>IrrSwitch)  // use sparse algorithm
+        { mon = kbase(sP,i-1);
+//          mon = ideal(mon[ncols(mon)..1]);
+          if (counter==0 && ncols(mon)==int(dimmat[i,1])) // then we can use all of mon
+          { B=normalize(evaluate_reynolds(REY,mon));
+            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] + NF(B,sP);
+            attrib(RedSSort[i-1][i-1],"size", saveAttr);
+            if (typeof(RedISSort[i-1]) <> "none")
+            { RedISSort[i-1] = RedISSort[i-1] + NF(B,sP);
+            }
+            else
+            { RedISSort[i-1] = NF(B,sP);
+            }
+            if (v) {"  We found: "; print(B);}
+          }
+          else
+          { irrcounter=0;
+            j=0;                         // j goes through all of mon -
+          // Compare the comments on the computation of reducible sec. inv.!
+            while (int(dimmat[i,1])<>counter)
+            { if ((j mod MonStep) == 0)
+              { if ((j+MonStep) <= ncols(mon))
+                { B[j+1..j+MonStep] = normalize(evaluate_reynolds(REY,ideal(mon[j+1..j+MonStep])));
+                  ReducedCandidates[j+1..j+MonStep] = reduce(ideal(B[j+1..j+MonStep]),sP);
+                  Indicator[j+1..j+MonStep] = reduce(ideal(ReducedCandidates[j+1..j+MonStep]),SaveRed);
+                }
+                else
+                { B[j+1..ncols(mon)] = normalize(evaluate_reynolds(REY,ideal(mon[j+1..ncols(mon)])));
+                  ReducedCandidates[j+1..ncols(mon)] = reduce(ideal(B[j+1..ncols(mon)]),sP);
+                  Indicator[j+1..ncols(mon)] = reduce(ideal(ReducedCandidates[j+1..ncols(mon)]),SaveRed);
+                }
+              }
+              j++;
+              helpP = Indicator[j];
+              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]);
+                }
+                Reductor = ideal(helpP);
+                attrib(Reductor, "isSB",1);
+                Indicator=reduce(Indicator,Reductor);
+                SizeSave++;
+                SaveRed[SizeSave] = helpP;
+                attrib(SaveRed, "isSB",1);
+                attrib(Reductor, "isSB",1);
+              }
+              mon[j]=0;
+              B[j]=0;
+              ReducedCandidates[j]=0;
+              Indicator[j]=0;
+            }
+          }
+        }  // if i>IrrSwitch
+        else  // use fast algorithm
+        { 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
@@ -5105 +5280 @@
                                        // 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;
+          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
-          { 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];
+          { 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 = Indicator[j];
+              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];
                 }
-                SizeSave=SizeSave+size(Reductor);
-                Reductor = ideal(0);
-                helpint = 0;
+                else
+                { RedISSort[i-1] = ideal(ReducedCandidates[j]);
+                }
+                if (v)
+                { "  We found the irreducible sec. inv. "+string(B[j]);
+                }
+                Reductor = ideal(helpP);
+                attrib(Reductor, "isSB",1);
+                Indicator=reduce(Indicator,Reductor);
+                SizeSave++;
+                SaveRed[SizeSave] = helpP;
                 attrib(SaveRed, "isSB",1);
                 attrib(Reductor, "isSB",1);
               }
+              B[j]=0;
+              ReducedCandidates[j]=0;
+              Indicator[j]=0;
             }
           }
-        }
-      }
+        } // i<=IrrSwitch
+      } // Computation of irreducible secondaries
       if (v)
       { "";
@@ -5189 +5361 @@
 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...
+list L = primary_invariants(A1,A2);
 matrix IS = irred_secondary_char0(L[1],L[2],L[3],0);
 IS;
@@ -5197 +5369 @@
 ///////////////////////////////////////////////////////////////////////////////
 
-
-proc secondary_charp (matrix P, matrix REY, string ring_name, list #)
+proc old_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>
@@ -5421 +5592 @@
 ///////////////////////////////////////////////////////////////////////////////
 
+
+proc secondary_charp (matrix P, matrix REY, string ring_name, list #)
+"USAGE:   secondary_charp(P,REY,ringname[,v][,\"old\"]);
+@*          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;
+@*          ringname: a <string> giving the name of a ring of characteristic 0
+               containing a 1x2 <matrix> M giving numerator and denominator of the Molien
+               series;
+@*          v: an optional <int>;
+@*          \"old\": if this string occurs as (optional) parameter, then an
+               old version of secondary_char0 is used (for downward
+               compatibility)
+ASSUME:   The characteristic of basering is not zero;
+          REY is the 1st return value of group_reynolds(), reynolds_molien() or
+          the second one of primary_invariants();
+@*        `ringname` is the name of a ring of characteristic 0 that has been created
+          by molien() or reynolds_molien() or primary_invariants()
+RETURN:   secondary invariants of the invariant ring (type <matrix>) and
+          irreducible secondary invariants (type <matrix>)
+DISPLAY:  information if v does not equal 0
+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)) using Groebner basis techniques. The size of this set
+          can be read off from the Molien series.
+EXAMPLE:  example secondary_charp; shows an example
+"
+{
+ //----------  Internal parameters, whose choice might improve the performance -----
+  int pieces = 15;   // For generating reducible secondaries, blocks of #pieces# secondaries
+                     // are formed.
+                     // If this parameter is small, the memory consumption will decrease.
+  int MonStep = 15;  // The Reynolds operator is applied to blocks of #MonStep# monomials.
+                     // If this parameter is small, the memory consumption will decrease.
+  int IrrSwitch = 7; // Up to degree #IrrSwitch#-1, we use a method for generating
+                     // irreducible secondaries that tends to produce a smaller output
+                     // (which is good for subsequent computations). However, this method
+                     // needs to much memory in high degrees, and so we use a sparser method
+                     // from degree #IrrSwitch# on.
+
+  def br=basering;
+
+ //---------------- 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;
+  }
+  if (size(#)>0)
+  { if (typeof(#[size(#)])=="string")
+    { if (#[size(#)]=="old")
+      { if (typeof(#[1])=="int")
+        { matrix S,IS = old_secondary_charp(P,REY,ring_name,#[1]);
+          return(S,IS);
+        }
+        else
+        { matrix S,IS = old_secondary_charp(P,REY,ring_name);
+          return(S,IS);
+        }
+      }
+      else
+      { "ERROR:   If the last optional parameter is a string, it should be \"old\".";
+        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 (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, irrcounter, 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
+
+  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;
+       }
+     }
+  }
+  deg_dim_vec=1;
+  for (j=2;j<=m;j++)
+  { deg_dim_vec=deg_dim_vec,int(dimmat[j,1]); // i.e., deg_dim_vec[i] = dimmat[i-1,1],
+                                              // which is the number of secondaries in degree i-1
+  }
+  if (v)
+  { "  In degree 0 we have: 1";
+    "";
+  }
+ //------------------------ initializing variables ----------------------------
+  setring br;
+  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));
+                           // This is the only explicit Groebner basis computation!
+  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 mon,B,IS;              // IS will contain all irr. sec. inv.
+
+  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 deg_dim_vec -
+  { if (deg_dim_vec[i]<>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(deg_dim_vec[i])+
+        " invariant(s)...";
+        "  Looking for Power Products...";
+      }
+      counter = 0;                       // we'll count up to deg_dim_vec[i]
+      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 (deg_dim_vec[i]==counter)
+        { break;
+        }
+        if (typeof(ISSort[k])<>"none")
+        { Mul1 = ISSort[k];
+        }
+        else
+        { Mul1 = ideal(0);
+        }
+        if ((deg_dim_vec[i-k]>0) && (Mul1[1]<>0))
+        { for (minD=k;minD<i-k;minD++)
+          { if (deg_dim_vec[i]==counter)
+            { break;
+            }
+            for (k2=1;k2 <= ((attrib(SSort[i-k-1][minD],"size")-1) div pieces)+1; k2++)
+            { if (deg_dim_vec[i]==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];
+                }
+              }
+//              if (i<10)
+                ProdCand = simplify(Mul1*Mul2,4);
+//              else { 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 = Indicator[ii];
+                  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", deg_dim_vec[i]," secondaries in degree",i-1;
+                    }
+                    if (deg_dim_vec[i]<>counter)
+                    { Reductor = ideal(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.
+                      Indicator=reduce(Indicator,Reductor);
+                      SizeSave++;
+                      SaveRed[SizeSave] = helpP;
+                      attrib(SaveRed, "isSB",1);
+                    }
+                    else
+                    { break;
+                    }
+                  }
+                }
+              }
+//              attrib(SaveRed, "isSB",1);
+            }
+          }
+        }
+      }
+
+      // The remaining sec. inv. are irreducible!
+      if (deg_dim_vec[i]<>counter)   // need more than all the power products
+      { if (v)
+        { "There are irreducible secondary invariants in degree ", i-1," !!";
+        }
+        if (i>IrrSwitch)  // use sparse algorithm
+        { mon = kbase(sP,i-1);
+//          mon = ideal(mon[ncols(mon)..1]);
+          if (counter==0 && ncols(mon)==deg_dim_vec[i]) // then we can use all of mon
+          { B=normalize(evaluate_reynolds(REY,mon));
+            IS=IS+B;
+            saveAttr = attrib(SSort[i-1][i-1],"size")+deg_dim_vec[i];
+            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 mon -
+            // Compare the comments on the computation of reducible sec. inv.!
+            while (deg_dim_vec[i]<>counter)
+            { if ((j mod MonStep) == 0)
+              { if ((j+MonStep) <= ncols(mon))
+                { B[j+1..j+MonStep] = normalize(evaluate_reynolds(REY,ideal(mon[j+1..j+MonStep])));
+                  ReducedCandidates[j+1..j+MonStep] = reduce(ideal(B[j+1..j+MonStep]),sP);
+                  Indicator[j+1..j+MonStep] = reduce(ideal(ReducedCandidates[j+1..j+MonStep]),SaveRed);
+                }
+                else
+                { B[j+1..ncols(mon)] = normalize(evaluate_reynolds(REY,ideal(mon[j+1..ncols(mon)])));
+                  ReducedCandidates[j+1..ncols(mon)] = reduce(ideal(B[j+1..ncols(mon)]),sP);
+                  Indicator[j+1..ncols(mon)] = reduce(ideal(ReducedCandidates[j+1..ncols(mon)]),SaveRed);
+                }
+              }
+              j++;
+              helpP = Indicator[j];
+              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]);
+                }
+                Reductor = ideal(helpP);
+                attrib(Reductor, "isSB",1);
+                Indicator=reduce(Indicator,Reductor);
+                SizeSave++;
+                SaveRed[SizeSave] = helpP;
+                attrib(SaveRed, "isSB",1);
+              }
+              B[j]=0;
+              ReducedCandidates[j]=0;
+              Indicator[j]=0;
+            }
+          }
+        }  // if i>IrrSwitch
+        else  // use fast algorithm
+        { 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
+          if (counter==0 && ncols(B)==deg_dim_vec[i]) // then we can take all of B
+          { IS=IS+B;
+            saveAttr = attrib(SSort[i-1][i-1],"size")+deg_dim_vec[i];
+            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 (deg_dim_vec[i]<>counter)
+            { j++;
+              helpP = Indicator[j];
+              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]);
+                }
+                Reductor = ideal(helpP);
+                attrib(Reductor, "isSB",1);
+                Indicator=reduce(Indicator,Reductor);
+                SizeSave++;
+                SaveRed[SizeSave] = helpP;
+                attrib(SaveRed, "isSB",1);
+              }
+              B[j]=0;
+              ReducedCandidates[j]=0;
+              Indicator[j]=0;
+            }
+          }
+        } // i<=IrrSwitch
+      } // Computation of irreducible secondaries
+      if (v)
+      { "";
+      }
+    } // if (deg_dim_vec[i]<>0)
+  }  // for i
+  if (v)
+  { "  We're done!";
+    "";
+  }
+  degBound = dgb;
+
+  // Prepare return:
+  int TotalNumber;
+  for (k=1;k<=m;k++)
+  { TotalNumber = TotalNumber + deg_dim_vec[k];
+  }
+  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();
+     }
+    }
+  }
+  if (ring_name=="aksldfalkdsflkj")
+  { kill `ring_name`;
+  }
+  return(matrix(S),matrix(compress(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)],1);
+         print(S);
+         print(IS);
+}
+///////////////////////////////////////////////////////////////////////////////
+
 proc secondary_no_molien (matrix P, matrix REY, list #)
 "USAGE:   secondary_no_molien(P,REY[,deg_vec,v]);
@@ -5437 +6066 @@
          monomials) of the basering modulo primary invariants, mapping those to
          invariants with the Reynolds operator and using these images as
-         candidates for secondary invariants.
+         candidates for secondary invariants. We have the Reynolds operator, hence,
+         we are in the non-modular case. Therefore, the invariant ring is Cohen-Macaulay,
+         hence the number of secondary invariants is the product of the degrees of
+         primary invariants divided by the group order.
 EXAMPLE: example secondary_no_molien; shows an example
 "
@@ -5879 +6511 @@
     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);
+//  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
@@ -5933 +6565 @@
     "           group order.";
   }
-  return(S);
+  return(matrix(S));
 }
 example
@@ -5954 +6586 @@
          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
+         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
@@ -6051 +6683 @@
         else
         { if (v)
-          { "  Since it is impossible for this program to calculate the Molien
+          { "  Since it is impossible for this programme to calculate the Molien
  series for";
             "  invariant rings over extension fields of prime characteristic, we
@@ -6171 +6803 @@
          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,
+         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
@@ -6286 +6918 @@
         else
         { if (v)
-          { "  Since it is impossible for this program to calculate the Molien
+          { "  Since it is impossible for this programme to calculate the Molien
  series for";
             "  invariant rings over extension fields of prime characteristic, we
@@ -6403 +7035 @@
 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.
+         the polynomials that are left over define the orbit variety
 EXAMPLE: example orbit_variety; shows an example
 "
@@ -6620 +7252 @@
 @*       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.
+         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
@@ -6706 +7339 @@
 @*       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.
+         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.
+@*       of the variety defined by I
 EXAMPLE: example image_of_variety; shows an example
 "
@@ -6743 +7376 @@
 
 
-
-
+proc secondary_charp_degbound (matrix P, int dmax, list #)
+"USAGE:   secondary_charp_degbound(P,G1,...,Gk[,v]);
+@*          P: a 1xn <matrix> with homogeneous primary invariants, where
+               n is the number of variables of the basering;
+@*          G1,...,Gk matrices generating a matrix group for which P are homogeneous
+            primary invariants.
+@*          v: an optional <int>;
+ASSUME:   The characteristic of basering is not zero
+RETURN:   secondary invariants of the invariant ring (type <matrix>) and
+          irreducible secondary invariants (type <matrix>) of degree at most dmax.
+DISPLAY:  information if v does not equal 0
+THEORY:   A basis of invariants of degrees d<dmax+1 is computed using invariant_basis.
+          Among the basis elements and 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)) using Groebner basis techniques.
+@*        If dmax is at least the Noether number of the invariant ring, then primary
+          invariants and the irreducible secondaries found by secondary_charp_degbound
+          are algebra generators of the invariant ring. This may work
+          faster than an application of secondary_not_cohen_macaulay.
+SEE ALSO: secondary_charp, secondary_char0, secondary_not_cohen_macaulay
+EXAMPLE:  example secondary_charp_molien; shows an example
+"
+{
+ //----------  Internal parameters, whose choice might improve the performance -----
+  int pieces = 15;   // For generating reducible secondaries, blocks of #pieces# secondaries
+                     // are formed.
+                     // If this parameter is small, the memory consumption will decrease.
+
+  def br=basering;
+
+ //---------------- checking input and setting verbose mode -------------------
+  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 (char(br)==0)
+  { "ERROR:   secondary_charp_degbound should only be used with rings of characteristic p>0.";
+    return();
+  }
+  int i;
+  if (size(#)==0)
+  { "ERROR: There are no generators given.";
+    return();
+  }
+  if (typeof(#[size(#)])=="int")
+  { int v=#[size(#)];
+  }
+  else
+  { int v=0;
+  }
+  list GEN;
+  for (i=1;i<size(#);i++)
+  { if ((typeof(#[i])=="matrix") and
+        (ncols(#[i])==n) and (nrows(#[i])==n))
+    { GEN[i] = #[i];
+    }
+    else
+    { "ERROR: The generators should be square matrices of appropriate size.";
+       return();
+    }
+  }
+  if (typeof(#[size(#)])<>"int")
+  { if ((typeof(#[size(#)])=="matrix") and
+        (ncols(#[size(#)])==n) and (nrows(#[size(#)])==n))
+    { GEN[i] = #[size(#)];
+    }
+    else
+    { "ERROR: The last optional parameter should be an integer.";
+       return();
+    }
+  }
+
+  if (ncols(P)<>n)
+  { "ERROR:   The first parameter ought to be the matrix of the primary";
+    "         invariants."
+    return();
+  }
+  if (v && voice==2)
+  { "";
+  }
+  if (v)
+  { "  In degree 0 we have: 1";
+    "";
+  }
+ //------------------------ initializing variables ----------------------------
+  setring br;
+  ideal ProdCand;          // contains products of secondary invariants,
+                           // i.e., candidates for reducible sec. inv.
+  ideal Mul1,Mul2;
+
+  int j, counter, irrcounter, d;
+  int m = dmax+1;
+  int dgb = degBound;
+  degBound = 0;
+//  degBound = dmax;
+  option(redSB);
+  ideal sP = groebner(ideal(P));
+                           // This is the only explicit Groebner basis computation!
+  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 mon,B,IS;              // IS will contain all irr. sec. inv.
+
+  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 the invariants of degree i-1 -
+    {                                  // starting with reducibles
+      if (v)
+      { "  Looking for Power Products...";
+      }
+      counter = 0;                      // counts the number of secondaries
+      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 (typeof(ISSort[k])<>"none")
+        { Mul1 = ISSort[k];
+        }
+        else
+        { Mul1 = ideal(0);
+        }
+        if (Mul1[1]<>0)
+        { for (minD=k;minD<i-k;minD++)
+          { for (k2=1;k2 <= ((attrib(SSort[i-k-1][minD],"size")-1) div pieces)+1; k2++)
+            { 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];
+                }
+              }
+//              if (i<10)
+                ProdCand = simplify(Mul1*Mul2,4);
+//              else { 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 = Indicator[ii];
+                  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, " secondaries in degree",i-1;
+                    }
+                    Reductor = ideal(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.
+                    Indicator=reduce(Indicator,Reductor);
+                    SizeSave++;
+                    SaveRed[SizeSave] = helpP;
+                    attrib(SaveRed, "isSB",1);
+                  }
+                }
+              }
+//              attrib(SaveRed, "isSB",1);
+            }
+          }
+        }
+      }
+
+      // If there are more sec. inv., then these are irreducible!
+      if (v)
+      { "Searching irreducible secondary invariants in degree ", i-1;
+      }
+      B=invariant_basis_modS(i-1,SaveRed,GEN); // B contains enough invariant polynomials of degree d
+      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 (size(Indicator)>0)
+      { j++;
+        helpP = Indicator[j];
+        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]);
+          }
+          Reductor = ideal(helpP);
+          attrib(Reductor, "isSB",1);
+          Indicator=reduce(Indicator,Reductor);
+          SizeSave++;
+          SaveRed[SizeSave] = helpP;
+          attrib(SaveRed, "isSB",1);
+        }
+        B[j]=0;
+        ReducedCandidates[j]=0;
+        Indicator[j]=0;
+      }
+      if (v)
+      { "";
+      }
+    } // for i
+  if (v)
+  { "  We're done!";
+    "";
+  }
+  degBound = dgb;
+
+  // Prepare return:
+  int TotalNumber;
+  for (k=1;k<m;k++)
+  { for (k2=1;k2<=k;k2++)
+    { if (typeof(attrib(SSort[k][k2],"size"))=="int")
+      { TotalNumber = TotalNumber + attrib(SSort[k][k2],"size");
+      }
+    }
+  }
+  matrix S[1][TotalNumber+1];
+  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(matrix(S),matrix(compress(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;
+         matrix P=primary_charp_without(A);
+         matrix S,IS=secondary_charp_degbound(P,4,A,1);
+         print(S);
+         print(IS);
+}
+