Changeset df39c9 in git

Timestamp:

Feb 8, 2007, 4:20:54 PM (17 years ago)

Author:

Simon King <king@…>

Branches:

(u'spielwiese', '4a9821a93ffdc22a6696668bd4f6b8c9de3e6c5f')

Children:

8b2345a5e557bd889ff1e0f5110102e0b4ae3c0c

Parents:

302a5d7a57af9072b9b2e8d50715cd8780fae77e

Message:

New version of finvar.lib.
- New algorithm for irred_secondary_char0.
  In high degrees, irreducible secondary invariants are computed without power
  products, i.e., without generating reducible secondary invariants. In some examples,
  this improves the performance by a factor of more than 1000. Particularly useful if
  the total number of secondary invariants is astronomic (>200000), while the number
  and the degrees of irreducible secondary invariants are relatively small.
- New function irred_secondary_no_molien. Computes irreducible secondary
  invariants using the Reynolds operator (hence, in the non-modular case)
  but without Molien series. In all degrees, irreducible secondary invariants are computed
  without generating reducible secondary invariants.
- mild improvements in secondary_char0, secondary_charp and secondary_and_irreducibles_no_molien.
- invariant_ring and invariant_ring_random now return irreducible secondary
  invariants also in characteristic p>0, in the non-modular case.


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

File:

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

Singular/LIB/finvar.lib

@@ -1 +1 @@
 ///////////////////////////////////////////////////////////////////////////////
-version="$Id: finvar.lib,v 1.58 2007-02-02 13:25:09 Singular Exp $"
+version="$Id: finvar.lib,v 1.59 2007-02-08 15:20:54 king Exp $"
 category="Invariant theory";
 info="
@@ -38 +38 @@
  power_products()                  exponents for power products
  secondary_char0()                 secondary invariants (s.i.) in char 0
- irred_secondary_char0()           irreducible secondary invariants in char 0
- secondary_charp()                 secondary invariants in char p
- secondary_no_molien()             secondary invariants, without Molien series
+ irred_secondary_char0()           irreducible s.i. in char 0
+ secondary_charp()                 s.i. in char p, with Molien series and 
+                                   Reynolds operator
+ secondary_no_molien()             s.i., without Molien series but with 
+                                   Reynolds operator
+ irred_secondary_no_molien()       irreducible s.i., 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
+ secondary_not_cohen_macaulay()    s.i. when the invariant ring is not Cohen-Macaulay
  orbit_variety()                   ideal of the orbit variety
  rel_orbit_variety()               ideal of a relative orbit variety (new version)
@@ -4305 +4308 @@
   m=nrows(dimmat);                     // m-1 is the highest degree
   if (v)
-  { "  In degree 0 we have: 1";
+  { "In degree 0 we have: 1";
     "";
   }
@@ -4321 +4324 @@
     {                                  // elements in the current degree (i-1)
       if (v)
-      { "  Searching in degree "+string(i-1)+", we need to find "+string(int(dimmat[i,1]))+" invariant(s)...";
+      { "Searching in degree "+string(i-1)+", we need to find "+string(int(dimmat[i,1]))+" invariant(s)...";
       }
       TEST=sP;
@@ -4468 +4471 @@
 {
  //----------  Internal parameters, whose choice might improve the performance -----
-  int pieces = 15;   // For generating reducible secondaries, blocks of #pieces# secondaries
+  int pieces = 3;   // For generating reducible secondaries, blocks of #pieces# secondaries
                      // are formed.
                      // If this parameter is small, the memory consumption will decrease.
@@ -4566 +4569 @@
   }
   if (v)
-  { "  In degree 0 we have: 1";
+  { "In degree 0 we have: 1";
     "";
   }
@@ -4580 +4583 @@
                            // This is the only explicit Groebner basis computation!
   ideal Reductor,SaveRed;  // sP union Reductor is a Groebner basis up to degree i-1
+ideal sP_Reductor;          // will contain sP union Reductor.
+int sPOffset = ncols(sP);
   int SizeSave;
 
@@ -4585 +4590 @@
                       // minimal degree of a non-constant invariant factor.
   list ISSort;        // irr. sec. inv. sorted by degree
+int NrIS;
 
   poly helpP;
   ideal helpI;
   ideal Indicator;        // will tell us which candidates for sec. inv. we can choose
-  ideal ReducedCandidates;
+//  ideal ReducedCandidates;
   int helpint;
   int k,k2,k3,minD;
@@ -4608 +4614 @@
     {                                  // elements in the current degree (i-1)
       if (v)
-      { "  Searching in degree "+string(i-1)+", we need to find "+string(int(dimmat[i,1]))+
+      { "Searching in degree "+string(i-1)+", we need to find "+string(int(dimmat[i,1]))+
         " invariant(s)...";
         "  Looking for Power Products...";
       }
+sP_Reductor = sP;
       counter = 0;                       // we'll count up to dimmat[i,1]
       Reductor = ideal(0);
@@ -4660 +4667 @@
                 ProdCand = simplify(Mul1*Mul2,4);
 //              else { ProdCand = Mul1*Mul2;}
-              ReducedCandidates = reduce(ProdCand,sP);
+//              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.
+//              Indicator = reduce(ReducedCandidates,SaveRed);
+              Indicator = reduce(ProdCand,sP_Reductor);
+                             // If Indicator[ii]<>0 then 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
@@ -4681 +4683 @@
                     attrib(SSort[i-1][k],"size",saveAttr);
                     if (v)
-                    { "We found",counter, " of", int(dimmat[i,1])," secondaries in degree",i-1;
+                    { "    We found",counter, "of", int(dimmat[i,1]),"secondaries in degree",i-1;
                     }
                     if (int(dimmat[i,1])<>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;
+sP_Reductor[sPOffset+SizeSave] = helpP;
+attrib(sP_Reductor, "isSB",1);
+                        // 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. Hence, sP_Reductor
+                        // is a homog. GB up to degree i-1.
                       attrib(SaveRed, "isSB",1);
+                      attrib(Reductor, "isSB",1);
                     }
                     else
@@ -4709 +4712 @@
       }
 
+NrIS = int(dimmat[i,1])-counter;
       // The remaining sec. inv. are irreducible!
-      if (int(dimmat[i,1])<>counter)   // need more than all the power products
+      if (NrIS>0)   // need more than all the power products
       { if (v)
-        { "There are irreducible secondary invariants in degree ", i-1," !";
-        }
-        mon = kbase(sP,i-1);
+        { "  There are ",NrIS,"irreducible secondary invariants in degree ", i-1;
+        }
+        mon = kbase(sP_Reductor,i-1);
 //        mon = ideal(mon[ncols(mon)..1]);
         ii = ncols(mon);
         if ((i>IrrSwitch) or (ii>200)) // use sparse algorithm
-        { if (counter==0 && ii==int(dimmat[i,1])) // then we can use all of mon
+        { if (ii+counter==int(dimmat[i,1])) // then we can use all of mon
           { B=normalize(evaluate_reynolds(REY,mon));
             IS=IS+B;
@@ -4730 +4734 @@
             { ISSort[i-1] = B;
             }
-            if (v) {"  We found: "; print(B);}
+            if (v) {"    We found all",NrIS,"irreducibles in degree ",i-1;}
           }
           else
@@ -4741 +4745 @@
                 { ideal tmp = normalize(evaluate_reynolds(REY,ideal(mon[j+1..j+MonStep])));
                   B[j+1..j+MonStep] = tmp[1..MonStep];
-                  tmp = reduce(tmp,sP);
-                  ReducedCandidates[j+1..j+MonStep] = tmp[1..MonStep];
-                  tmp = reduce(tmp,SaveRed);
+//                  tmp = reduce(tmp,sP);
+//                  ReducedCandidates[j+1..j+MonStep] = tmp[1..MonStep];
+                  tmp = reduce(tmp,sP_Reductor);
                   Indicator[j+1..j+MonStep] = tmp[1..MonStep];
                   kill tmp;
@@ -4750 +4754 @@
                 { ideal tmp = normalize(evaluate_reynolds(REY,ideal(mon[j+1..ii])));
                   B[j+1..ii] = tmp[1..ii-j];
-                  tmp = reduce(tmp,sP);
-                  ReducedCandidates[j+1..ii] = tmp[1..ii-j];
-                  tmp = reduce(tmp,SaveRed);
+//                  tmp = reduce(tmp,sP);
+//                  ReducedCandidates[j+1..ii] = tmp[1..ii-j];
+                  tmp = reduce(tmp,sP_Reductor);
                   Indicator[j+1..ii] = tmp[1..ii-j];
                   kill tmp;
@@ -4771 +4775 @@
                 { ISSort[i-1] = ideal(B[j]);
                 }
-                if (v)
-                { "  We found the irreducible sec. inv. "+string(B[j]);
+                if (v) 
+                { "    We found", irrcounter,"of", NrIS ,"irr. sec. inv. in degree ",i-1; 
                 }
                 Reductor = ideal(helpP);
@@ -4780 +4784 @@
                 SaveRed[SizeSave] = helpP;
                 attrib(SaveRed, "isSB",1);
+sP_Reductor[sPOffset+SizeSave] = helpP;
+attrib(sP_Reductor, "isSB",1);
               }
               B[j]=0;
-              ReducedCandidates[j]=0;
+//              ReducedCandidates[j]=0;
               Indicator[j]=0;
             }
@@ -4788 +4794 @@
         }  // if i>IrrSwitch
         else  // use fast algorithm
-        { B=sort_of_invariant_basis(sP,REY,i-1,int(dimmat[i,1])*6);
+        { B=sort_of_invariant_basis(sP_Reductor,REY,i-1,int(dimmat[i,1])*6);
                                        // B contains 
                                        // images of kbase(sP,i-1) under the
@@ -4804 +4810 @@
             { ISSort[i-1] = B;
             }
-            if (v) {"  We found: "; print(B);}
+            if (v) {"    We found all",NrIS,"irreducibles in degree ",i-1;}
           }
           else
@@ -4810 +4816 @@
             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);
+//            ReducedCandidates = reduce(B,sP);
+            Indicator = reduce(B,sP_Reductor);
             while (int(dimmat[i,1])<>counter)
             { j++;
@@ -4827 +4833 @@
                 { ISSort[i-1] = ideal(B[j]);
                 }
-                if (v)
-                { "  We found the irreducible sec. inv. "+string(B[j]);
+                if (v) 
+                { "    We found", irrcounter,"of", NrIS,"irr. sec. inv. in degree ",i-1; 
                 }
                 Reductor = ideal(helpP);
@@ -4836 +4842 @@
                 SaveRed[SizeSave] = helpP;
                 attrib(SaveRed, "isSB",1);
+sP_Reductor[sPOffset+SizeSave] = helpP;
+attrib(sP_Reductor, "isSB",1);
               }
               B[j]=0;
-              ReducedCandidates[j]=0;
+//              ReducedCandidates[j]=0;
               Indicator[j]=0;
             }
@@ -4892 +4900 @@
 proc irred_secondary_char0 (matrix P, matrix REY, matrix M, list #)
 "
-USAGE:    irred_secondary_char0(P,REY,M[,v]);
+USAGE:    irred_secondary_char0(P,REY,M[,v][,\"PP\"]);
 @*          P: a 1xn <matrix> with homogeneous primary invariants, where
                n is the number of variables of the basering;
@@ -4900 +4908 @@
                Molien series;
 @*          v: an optional <int>;
+@*          \"PP\": if this string occurs as (optional) parameter, then in
+               all degrees power products of irr. sec. inv. will be computed.
 RETURN:   Irreducible homogeneous secondary invariants of the invariant ring 
           (type <matrix>)
@@ -4915 +4925 @@
           independent modulo the primary invariants using Groebner basis techniques
           (see paper \"Fast Computation of Secondary Invariants\" by S. King).
-          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.
+          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.
+@*        Moreover, if no irr. sec. inv. in degree d-1 have been found and unless the last
+          optional paramter \"PP\" is used, a Groebner basis of primary invariants and 
+          irreducible secondary invariants up to degree d-2 is computed, which allows to 
+          detect irr. sec. inv. in degree d without computing power products. 
 @*        There are three internal parameters \"pieces\", \"MonStep\" and \"IrrSwitch\".
           The default values of the parameters should be fine in most cases. However,
@@ -4927 +4941 @@
 "
 { //----------  Internal parameters, whose choice might improve the performance -----
-  int pieces = 15;   // For generating reducible secondaries, blocks of #pieces# secondaries
+  int pieces = 3;    // 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.
+  int IrrSwitch = 7; // Up to degree #IrrSwitch#-1, or if there are few monomials in kbase(sP,i-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# upwards.
 
   def br=basering;
@@ -4947 +4961 @@
   int i;
   if (size(#)>0)
-  { if (typeof(#[size(#)])=="int")
-    { int v=#[size(#)];
+  { if (typeof(#[1])=="int")
+    { int v=#[1];
     }
     else
@@ -4956 +4970 @@
   else
   { int v=0;
+  }
+
+  int UsePP;
+  if (size(#)>0)
+  { if (typeof(#[size(#)])=="string")
+    { if (#[size(#)]=="PP")
+      { UsePP = 1;
+      }
+      else
+      { "ERROR:   If the last optional parameter is a string, it should be \"PP\".";
+        return(matrix(ideal()),matrix(ideal()));
+      }
+    }
   }
   int n=nvars(br);                // n is the number of variables, as well
@@ -4995 +5022 @@
   m=nrows(dimmat);                     // m-1 is the highest degree
   if (v) 
-  { "We need to find";
+  { "There are";
     for (j=1;j<=m;j++)
      { if (int(dimmat[j,1])<>1)
@@ -5006 +5033 @@
   }
   if (v)
-  { "  In degree 0 we have: 1";
+  { "In degree 0 we have: 1";
     "";
   }
@@ -5019 +5046 @@
   option(redSB);
   ideal sP = groebner(ideal(P));
+ideal TotalReductor = sP;        // will eventually be groebner(ideal(P)+IS)
+ideal sP_Reductor;          // will contain sP union Reductor.
+int sPOffset = ncols(sP);
+int LastNewIS=-10;           // if the Last New Irr. Sec. was found in degree i-2 then
+                           // we compute TotalReductor again
+int TotalSet = 1;         // It TotalSet is 0, we eventually need to re-compute TotalReductor.
   ideal Reductor,SaveRed;  // sP union Reductor is a Groebner basis up to degree i-1
   int SizeSave;
@@ -5045 +5078 @@
 //--------------------- generating secondary invariants ----------------------
   for (i=2;i<=m;i++)                   // going through dimmat -
-  { if (int(dimmat[i,1])<>0)           // when it is == 0 we need to find no
+  { 
+if ((LastNewIS == i-3) and (UsePP==0))
+{ if (v) {"Computing Groebner basis for primaries and previously found irreducibles...";}
+  TotalReductor = groebner(sP+IS); 
+  TotalSet = 1; 
+}
+    if (int(dimmat[i,1])<>0)           // when it is == 0 we need to find no
     {                                  // elements in the current degree (i-1)
       if (v)
-      { "  Searching in degree "+string(i-1)+", we need to find "+
-        string(int(dimmat[i,1]))+" invariant(s)...";
-        "  Looking for Power Products...";
+      { "Searching in degree "+string(i-1)+". There are "+
+        string(int(dimmat[i,1]))+" secondary invariant(s)...";
       }
       counter = 0;                       // we'll count up to dimmat[i,1]
@@ -5060 +5098 @@
       attrib(SaveRed,"isSB",1);
 
+if ( TotalSet==0 )
+{
+      if (v)
+      {
+        "  Looking for Power Products...";
+      }
+sP_Reductor = sP;
 // We start searching for reducible secondary invariants in degree i-1, i.e., those
 // that are power products of irreducible secondary invariants.
@@ -5126 +5171 @@
                     attrib(RedSSort[i-1][k],"size",saveAttr);
                     if (v)
-                    { "We found",counter, " of", int(dimmat[i,1])," secondaries in degree",i-1;
+                    { 
+                      "    We found reducible sec. inv. number "+string(counter);
                     }
                     if (int(dimmat[i,1])<>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;
+sP_Reductor[sPOffset+SizeSave] = helpP;
+attrib(sP_Reductor, "isSB",1);
+                        // 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. Hence, sP_Reductor
+                        // is a homog. GB up to degree i-1.
                       attrib(SaveRed, "isSB",1);
                       attrib(Reductor, "isSB",1);
@@ -5146 +5192 @@
                     { break;
                     }
-                  }
-                }
-              }
+                  } // new reducible sec. inv. found
+                } // loop through ProdCand
+              } // if there is some reducible sec. inv.
               attrib(SaveRed, "isSB",1);
               attrib(Reductor, "isSB",1);
@@ -5155 +5201 @@
         }
       }
+} // if TotalSet==0
+else 
+{ if (v) {"  We don't compute Power Products!"; }
+// Instead of computing Power Products, we can compute a Groebner basis of the ideal
+// generated by the primary and previously found irr. sec. invariants up to degree i-2. An invariant
+// polynomial of degree i-1 belongs to this ideal if and only if it belongs to the sub-algebra
+// generated by primary and irreducible secondary invariants up to degree i-2.
+// Hence, if we find an invariant outside this ideal, it is an irreducible secondary 
+// invariant of degree i-1.
+}
 
       // The remaining sec. inv. are irreducible!
       if (int(dimmat[i,1])<>counter)   // need more than all the power products
+                                       // or work without power products
       { if (v)
-        { "There are irreducible secondary invariants in degree ", i-1," !";
-        }
-        mon = kbase(sP,i-1);
+        { "  Looking for irreducible secondary invariants in degree ", i-1;
+        }
+//        mon = kbase(sP,i-1);
+        if (TotalSet==0)
+        { mon = kbase(sP_Reductor,i-1);
+        }
+        else
+        { mon = kbase(TotalReductor,i-1);
+        }
         ii = ncols(mon);
-        if ((i>IrrSwitch) or (ii>200))  // use sparse algorithm
-        { if (counter==0 && ii==int(dimmat[i,1])) // then we can use all of mon
+if (mon[1]<>0)
+      { if ((i>IrrSwitch) or (ii>200))  // use sparse algorithm
+        { if (TotalSet==0 && ii+counter==int(dimmat[i,1])) // then we can use all of mon
           { B=normalize(evaluate_reynolds(REY,mon));
             IS=IS+B;
@@ -5176 +5240 @@
             { RedISSort[i-1] = NF(B,sP);
             }
-            if (v) {"  We found: "; print(B);}
+            TotalSet=0;  // i.e., TotalReductor needs to be re-computed
+//            if (v) {"  We found: "; print(B);}
+            if (v) {"    We found "+string(size(B))+" irred. sec. inv.";}
           }
           else
@@ -5182 +5248 @@
             j=0;                         // j goes through all of mon -
           // Compare the comments on the computation of reducible sec. inv.!
-            while (int(dimmat[i,1])<>counter) 
+while (j < ii) // we can break the loop if counter==int(dimmat[i,1])
+// while (int(dimmat[i,1])<>counter) 
             { if ((j mod MonStep) == 0) 
               { if ((j+MonStep) <= ii)
@@ -5189 +5256 @@
                   tmp = reduce(tmp,sP);
                   ReducedCandidates[j+1..j+MonStep] = tmp[1..MonStep];
+if (TotalSet == 0)
+{
                   tmp = reduce(tmp,SaveRed);
+}
+else
+{
+                  tmp = reduce(tmp,TotalReductor);
+}
                   Indicator[j+1..j+MonStep] = tmp[1..MonStep];
                   kill tmp;
@@ -5198 +5272 @@
                   tmp = reduce(tmp,sP);
                   ReducedCandidates[j+1..ii] = tmp[1..ii-j];
+if (TotalSet == 0)
+{
                   tmp = reduce(tmp,SaveRed);
+}
+else
+{
+                  tmp = reduce(tmp,TotalReductor);
+}
                   Indicator[j+1..ii] = tmp[1..ii-j];
                   kill tmp;
@@ -5206 +5287 @@
               helpP = Indicator[j];
               if (helpP <>0)                     // B[j] should be added
-              { counter++; irrcounter++;
+              { TotalSet=0;  // i.e., TotalReductor needs to be re-computed
+                counter++; irrcounter++;
                 IS=IS,B[j];
                 saveAttr = attrib(RedSSort[i-1][i-1],"size")+1;
@@ -5218 +5300 @@
                 }
                 if (v)
-                { "  We found the irreducible sec. inv. "+string(B[j]);
+                { //"  We found the irreducible sec. inv. "+string(B[j]);
+                  "    We found irred. sec. inv. number "+string(irrcounter);
                 }
                 Reductor = ideal(helpP);
@@ -5227 +5310 @@
                 attrib(SaveRed, "isSB",1);
                 attrib(Reductor, "isSB",1);
+LastNewIS = i-1;
+if (int(dimmat[i,1])==counter) { break; }
               }
               mon[j]=0;
@@ -5236 +5321 @@
         }  // if i>IrrSwitch
         else  // use fast algorithm
-        { B=sort_of_invariant_basis(sP,REY,i-1,int(dimmat[i,1])*6);
+        { if (TotalSet == 0)
+          { B=sort_of_invariant_basis(sP_Reductor,REY,i-1,int(dimmat[i,1])*6);
+          }
+          else
+          { B=sort_of_invariant_basis(TotalReductor,REY,i-1,int(dimmat[i,1])*6);
+          }
+//        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
@@ -5242 +5333 @@
                                        // 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
+          if (TotalSet==0 && ncols(B)+counter==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]);
@@ -5253 +5344 @@
             { RedISSort[i-1] = B;
             }
-            if (v) {"  We found: ";print(B);}
+//            if (v) {"  We found: ";print(B);}
+            if (v) {"    We found "+string(size(B))+" irred. sec. inv.";}
           }
           else
@@ -5266 +5358 @@
               helpP = Indicator[j];
               if (helpP <>0)                     // B[j] should be added
-              { counter++; irrcounter++;
+              { TotalSet=0;  // i.e., TotalReductor needs to be re-computed
+                counter++; irrcounter++;
                 IS=IS,B[j];
                 saveAttr = attrib(RedSSort[i-1][i-1],"size")+1;
@@ -5278 +5371 @@
                 }
                 if (v)
-                { "  We found the irreducible sec. inv. "+string(B[j]);
+                { //"  We found the irreducible sec. inv. "+string(B[j]);
+                  "    We found irred. sec. inv. number "+string(irrcounter);
                 }
                 Reductor = ideal(helpP);
@@ -5287 +5381 @@
                 attrib(SaveRed, "isSB",1);
                 attrib(Reductor, "isSB",1);
+LastNewIS = i-1;
+if (int(dimmat[i,1])==counter) { break; }
               }
               B[j]=0;
@@ -5294 +5390 @@
           }
         } // i<=IrrSwitch
+} // if mon[1]<>0
       } // Computation of irreducible secondaries
       if (v)
@@ -5424 +5521 @@
   }
   if (v)
-  { "  In degree 0 we have: 1";
+  { "In degree 0 we have: 1";
     "";
   }
@@ -5441 +5538 @@
     {                                  // 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)...";
+      { "Searching in degree "+string(i-1)+", we need to find "+string(deg_dim_vec[i])+" invariant(s)...";
       }
       TEST=sP;
@@ -5587 +5684 @@
 {
  //----------  Internal parameters, whose choice might improve the performance -----
-  int pieces = 15;   // For generating reducible secondaries, blocks of #pieces# secondaries
+  int pieces = 3;   // For generating reducible secondaries, blocks of #pieces# secondaries
                      // are formed.
                      // If this parameter is small, the memory consumption will decrease.
@@ -5695 +5792 @@
   }
   if (v)
-  { "  In degree 0 we have: 1";
+  { "In degree 0 we have: 1";
     "";
   }
@@ -5710 +5807 @@
                            // This is the only explicit Groebner basis computation!
   ideal Reductor,SaveRed;  // sP union Reductor is a Groebner basis up to degree i-1
+ideal sP_Reductor;          // will contain sP union Reductor.
+int sPOffset = ncols(sP);
   int SizeSave;
 
@@ -5715 +5814 @@
                       // minimal degree of a non-constant invariant factor.
   list ISSort;        // irr. sec. inv. sorted by degree
+int NrIS;
 
   poly helpP;
@@ -5738 +5838 @@
     {                                  // elements in the current degree (i-1)
       if (v)
-      { "  Searching in degree "+string(i-1)+", we need to find "+string(deg_dim_vec[i])+
+      { "Searching in degree "+string(i-1)+", we need to find "+string(deg_dim_vec[i])+
         " invariant(s)...";
         "  Looking for Power Products...";
@@ -5744 +5844 @@
       counter = 0;                       // we'll count up to deg_dim_vec[i]
       Reductor = ideal(0);
+sP_Reductor = sP;
       helpint = 0;
       SaveRed = Reductor;
@@ -5811 +5912 @@
                     attrib(SSort[i-1][k],"size",saveAttr);
                     if (v)
-                    { "We found",counter, " of", deg_dim_vec[i]," secondaries in degree",i-1;
+                    { "    We found", counter, "of", deg_dim_vec[i], "secondaries in degree",i-1;
                     }
                     if (deg_dim_vec[i]<>counter)
                     { Reductor = ideal(helpP);
+                      attrib(Reductor, "isSB",1);
+                      Indicator=reduce(Indicator,Reductor);
+                      SizeSave++;
+sP_Reductor[sPOffset+SizeSave] = helpP;
+attrib(sP_Reductor, "isSB",1);
                         // 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++;
+                        // a homogeneous Groebner basis up to degree i-1. Hence, sP_Reductor
+                        // is a homog. GB up to degree i-1.
                       SaveRed[SizeSave] = helpP;
                       attrib(SaveRed, "isSB",1);
@@ -5839 +5940 @@
       }
 
+NrIS = deg_dim_vec[i]-counter;
       // The remaining sec. inv. are irreducible!
-      if (deg_dim_vec[i]<>counter)   // need more than all the power products
+      if (NrIS>0)   // need more than all the power products
       { if (v)
-        { "There are irreducible secondary invariants in degree ", i-1," !";
-        }
-        mon = kbase(sP,i-1);
+        { "  There are ",NrIS,"irreducible secondary invariants in degree ", i-1;
+        }
+        mon = kbase(sP_Reductor,i-1);
         ii = ncols(mon);
 //        mon = ideal(mon[ncols(mon)..1]);
@@ -5860 +5962 @@
             { ISSort[i-1] = B;
             }
-            if (v) {"  We found: "; print(B);}
+            if (v) {"    We found all",NrIS,"irreducibles in degree ",i-1;}
           }
           else
@@ -5901 +6003 @@
                 { ISSort[i-1] = ideal(B[j]);
                 }
-                if (v)
-                { "  We found the irreducible sec. inv. "+string(B[j]);
+                if (v) 
+                { "    We found", irrcounter,"of", NrIS ,"irr. sec. inv. in degree ",i-1; 
                 }
                 Reductor = ideal(helpP);
@@ -5934 +6036 @@
             { ISSort[i-1] = B;
             }
-            if (v) {"  We found: "; print(B);}
+            if (v) {"    We found all",NrIS,"irreducibles in degree ",i-1;}
           }
           else
@@ -5957 +6059 @@
                 { ISSort[i-1] = ideal(B[j]);
                 }
-                if (v)
-                { "  We found the irreducible sec. inv. "+string(B[j]);
+                if (v) 
+                { "    We found", irrcounter,"of", NrIS ,"irr. sec. inv. in degree ",i-1; 
                 }
                 Reductor = ideal(helpP);
@@ -6040 +6142 @@
           hence the number of secondary invariants is the product of the degrees of 
           primary invariants divided by the group order.
+          <secondary_and_irreducible_no_molien> should usually be faster and of more useful
+          functionality.
+SEE ALSO: secondary_and_irreducible_no_molien
 EXAMPLE:  example secondary_no_molien; shows an example
 "
@@ -6123 +6228 @@
   { "  We need to find "+string(max)+" secondary invariants.";
     "";
-    "  In degree 0 we have: 1";
+    "In degree 0 we have: 1";
     "";
   }
@@ -6160 +6265 @@
     if (deg_vec[k]<>i)
     { if (v)
-      { "  Searching in degree "+string(i)+"...";
+      { "Searching in degree "+string(i)+"...";
       }
       mon = kbase(sP,i);
@@ -6194 +6299 @@
             S[counter] = B[j];
             if (v)
-            { "  We found sec. inv. number "+string(counter)+" in degree "+string(i);
+            { "    We found sec. inv. number "+string(counter)+" in degree "+string(i);
             }
             if (counter == max) {break;}
@@ -6226 +6331 @@
             S[counter] = B[j];
             if (v)
-            { "  We found sec. inv. number "+string(counter)+" in degree "+string(i);
+            { "    We found sec. inv. number "+string(counter)+" in degree "+string(i);
             }
             if (counter == max) {break;}
@@ -6292 +6397 @@
 EXAMPLE: example secondary_and_irreducibles_no_molien; shows an example
 "
-{ int pieces = 15;   // For generating reducible secondaries, blocks of #pieces# secondaries
+{ int pieces = 3;    // For generating reducible secondaries, blocks of #pieces# secondaries
                      // are formed.
                      // If this parameter is small, the memory consumption will decrease.
@@ -6376 +6481 @@
   { "  We need to find "+string(max)+" secondary invariants.";
     "";
-    "  In degree 0 we have: 1";
+    "In degree 0 we have: 1";
     "";
   }
@@ -6429 +6534 @@
     if (deg_vec[l]<>i-1)
     { if (v)
-      { "  Searching in degree "+string(i-1);
+      { "Searching in degree "+string(i-1);
         "  Looking for Power Products...";
       }
@@ -6499 +6604 @@
                     attrib(SSort[i-1][k],"size",saveAttr);
                     if (v)
-                    { "We found sec. inv. number",counter, " in degree",i-1;
+                    { "    We found sec. inv. number",counter, "in degree",i-1;
                     }
                     if (max<>counter)
@@ -6530 +6635 @@
       if (max<>counter)
       { if (v)
-        { "Looking for irreducible secondary invariants in degree ", i-1;
+        { "  Looking for irreducible secondary invariants in degree ", i-1;
         }
         mon = kbase(sP,i-1);
@@ -6575 +6680 @@
               }
               if (v)
-              { "  We found the irreducible sec. inv. "+string(B[j]);
+              { "    We found irreducible sec. inv. number", irrcounter, "in degree", i-1;
               }
               if (counter == max) { break; }
@@ -6618 +6723 @@
               }
               if (v)
-              { "  We found the irreducible sec. inv. "+string(B[j]);
+              { "    We found irreducible sec. inv. number", irrcounter, "in degree", i-1;
               }
               Reductor = ideal(helpP);
@@ -6684 +6789 @@
          matrix S,IS=secondary_and_irreducibles_no_molien(L[1..2],intvec(1,2),1);
          print(S);
+         print(IS);
+}
+///////////////////////////////////////////////////////////////////////////////
+
+proc irred_secondary_no_molien (matrix P, matrix REY, list #)
+"USAGE:   irred_secondary_no_molien(P,REY[,deg_vec,v]);
+         P: a 1xn <matrix> with primary invariants, REY: a gxn <matrix>
+         representing the Reynolds operator, deg_vec: an optional <intvec>
+         listing some degrees where no irreducible secondary invariants can
+         be found, v: an optional <int>
+ASSUME:  n is the number of variables of the basering, g the size of the group,
+         REY is the 1st return value of group_reynolds(), reynolds_molien() or
+         the second one of primary_invariants()
+RETURN:  Irreducible secondary invariants of the invariant ring (type <matrix>)
+DISPLAY: information if v does not equal 0
+THEORY:  Irred. secondary invariants are calculated by finding a basis (in terms of
+         monomials) of the basering modulo primary and previously found secondary
+         invariants, mapping those to invariants with the Reynolds operator. Among 
+         these images we pick a maximal subset that is linearly independent modulo 
+         the primary and previously found secondary invariants, using Groebner basis 
+         techniques.
+EXAMPLE: example irred_secondary_no_molien; shows an example
+"
+{ int MonStep = 29;
+
+  def br=basering;
+ //------------------ checking input and setting verbose ----------------------
+  if (size(#)==1 or size(#)==2)
+  { if (typeof(#[size(#)])=="int")
+    { if (size(#)==2)
+      { if (typeof(#[size(#)-1])=="intvec")
+        { intvec deg_vec=#[size(#)-1];
+        }
+        else
+        { "ERROR:   the third parameter should be an <intvec>";
+          return();
+        }
+      }
+      int v=#[size(#)];
+    }
+    else
+    { if (size(#)==1)
+      { if (typeof(#[size(#)])=="intvec")
+        { intvec deg_vec=#[size(#)];
+          int v=0;
+        }
+        else
+        { "ERROR:   the third parameter should be an <intvec>";
+          return();
+        }
+      }
+      else
+      { "ERROR:   wrong list of parameters";
+        return();
+      }
+    }
+  }
+  else
+  { if (size(#)>2)
+    { "ERROR:   there are too many parameters";
+      return();
+    }
+    int v=0;
+  }
+  int n=nvars(basering);               // n is the number of variables, as well
+                                       // as the size of the matrices, as well
+                                       // as the number of primary invariants,
+                                       // we should get
+  if (ncols(P)<>n)
+  { "ERROR:   The first parameter ought to be the matrix of the primary";
+    "         invariants."
+    return();
+  }
+  if (ncols(REY)<>n) 
+  { "ERROR:   The second parameter ought to be the Reynolds operator.";
+    return();
+  }
+  if (char(br)<>0)
+  { if (nrows(REY) mod char(br) == 0)
+    { "ERROR:   We need to be in the non-modular case.";
+      return();
+    }
+  }
+  if (v && voice==2)
+  { "";
+  }
+  if (defined(deg_vec)<>voice)
+  { intvec deg_vec;
+  }
+  int l=1;
+
+ //-------------------------- initializing variables --------------------------
+
+  int dgb = degBound;
+  degBound = 0;
+  int d,block, counter, totalcount, monsize,i,j, fin,hlp;
+  poly helpP,lmp; 
+  block = 1;
+  ideal mon, MON, Inv, IS,RedInv, NewIS, NewReductor;
+  int NewIScounter;
+  ideal sP = slimgb(ideal(P));
+  int Max = vdim(sP);
+  if (Max<0) 
+  { "ERROR:   The second parameter ought to be the Reynolds operator.";
+    return();
+  }
+  ideal sIS = sP;   // sIS will be a Groebner basis up to degree d of P+irred. sec. inv.
+  while (totalcount < Max)
+  { d++;
+    if (deg_vec[l]<>d)
+    { if (v) { "Searching irred. sec. inv. in degree "+string(d); }
+      NewIS = ideal();
+      NewIScounter=0;
+      MON = kbase(sP,d);
+      mon = kbase(sIS,d);
+      monsize=ncols(mon);
+      totalcount = totalcount + ncols(MON);
+      if (mon[1]<>0)
+      {
+        if (v) { "  We have "+string(monsize)+" candidates for irred. secondaries";}
+        block=0; // Loops through the monomials
+        while (block<monsize)
+        { if ((block mod MonStep) == 0) 
+          { if ((block+MonStep) <= monsize)
+            { ideal tmp = normalize(evaluate_reynolds(REY,ideal(mon[block+1..block+MonStep])));
+              Inv[block+1..block+MonStep] = tmp[1..MonStep];
+              tmp = reduce(tmp,sIS);
+              RedInv[block+1..block+MonStep] = tmp[1..MonStep];
+              kill tmp;
+            }
+            else
+            { ideal tmp = normalize(evaluate_reynolds(REY,ideal(mon[block+1..monsize])));
+              Inv[block+1..monsize] = tmp[1..monsize-block];
+              tmp = reduce(tmp,sIS);
+              RedInv[block+1..monsize] = tmp[1..monsize-block];
+              kill tmp;
+            }
+          } 
+          block++;
+          helpP = RedInv[block];
+          if (helpP<>0)
+          { counter++; // found new irr. sec. inv.!
+            if (v) { "    We found irred. sec. inv. number "+string(counter)+" in degree "+string(d); }
+            IS[counter] = Inv[block];
+            sIS = sIS,helpP;
+            attrib(sIS,"isSB",1);
+            NewIScounter++;
+            NewIS[NewIScounter] = helpP;
+// for permutation of variables, the following might yield a speed-up;
+// however, for general group actions it does not work:
+//          for (i=block;i<=monsize;i++)
+//          { if (Inv[block]<>reduce(Inv[block],std(mon[i])))
+//            { mon[i] = 0;
+//            }
+//          }
+
+            mon[block]=0;
+            Inv[block]=0; 
+            RedInv[block]=0;
+            NewReductor[1]=helpP;
+            attrib(NewReductor,"isSB",1);
+            RedInv = reduce(RedInv,NewReductor);
+          }
+        } // loop through monomials
+      } // if mon[1]<>0
+      if (totalcount < Max)
+      { if (NewIScounter==0)
+        { if (fin==0)
+          { degBound = 0;
+            sIS = groebner(sIS);
+            fin=1;  
+          }
+        }
+        else
+        { degBound = d+1;
+          sIS = groebner(sIS+NewIS);  
+          degBound = 0;
+          fin = 0;
+        }
+      }
+      attrib(sIS,"isSB",1);
+    } // if degree d is not excluded by deg_vec 
+    else
+    { if (size(deg_vec)==l)
+      { l=1;
+      }
+      else
+      { l++;
+      }
+    }
+  } // loop through degree
+  degBound = dgb;
+  return(matrix(IS));
+}
+
+example
+{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; 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,intvec(1,1,0));
+         // In that example, there are no secondary invariants
+         // in degree 1 or 2.
+         matrix IS=irred_secondary_no_molien(L[1..2],intvec(1,2),1);
          print(IS);
 }
@@ -6881 +7189 @@
 RETURN:  primary and secondary invariants (both of type <matrix>) generating
          the invariant ring with respect to the matrix group generated by the
-         matrices in the input and irreducible secondary invariants (type
-         <matrix>) if the Molien series was available
+         matrices in the input, and irreducible secondary invariants if we are
+         in the non-modular case.
 DISPLAY: information about the various stages of the program if the third flag
          does not equal 0
@@ -6977 +7285 @@
           list l=primary_charp_no_molien(L[1],v);
           if (size(l)==2)
-          { matrix S=secondary_no_molien(l[1],L[1],l[2],v);
+          { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v);
           }
           else
-          { matrix S=secondary_no_molien(l[1],L[1],v);
+          { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v);
           }
-          return(l[1],S);
+          return(l[1],S,IS);
         }
       }
@@ -7003 +7311 @@
     { list l=primary_char0_no_molien(L[1],v);
       if (size(l)==2)
-      { matrix S=secondary_no_molien(l[1],L[1],l[2],v);
+      { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v);
       }
       else
-      { matrix S=secondary_no_molien(l[1],L[1],v);
-      }
-      return(l[1],S);
+      { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v);
+      }
+      return(l[1],S,IS);
     }
     else
@@ -7014 +7322 @@
       { list l=primary_charp_no_molien(L[1],v); // case
         if (size(l)==2)
-        { matrix S=secondary_no_molien(l[1],L[1],l[2],v);
+        { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v);
         }
         else
-        { matrix S=secondary_no_molien(l[1],L[1],v);
-        }
-        return(l[1],S);
+        { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v);
+        }
+        return(l[1],S,IS);
       }
       else                             // the modular case
@@ -7098 +7406 @@
          among the coefficients)
 RETURN:  primary and secondary invariants (both of type <matrix>) generating
-         invariant ring with respect to the matrix group generated by the
-         matrices in the input and irreducible secondary invariants (type
-         <matrix>) if the Molien series was available
+         the invariant ring with respect to the matrix group generated by the
+         matrices in the input, and irreducible secondary invariants if we are
+         in the non-modular case.
 DISPLAY: information about the various stages of the program if the third flag
          does not equal 0
@@ -7212 +7520 @@
           list l=primary_charp_no_molien_random(L[1],max,v);
           if (size(l)==2)
-          { matrix S=secondary_no_molien(l[1],L[1],l[2],v);
+          { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v);
           }
           else
-          { matrix S=secondary_no_molien(l[1],L[1],v);
+          { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v);
           }
-          return(l[1],S);
+          return(l[1],S,IS);
         }
       }
@@ -7238 +7546 @@
     { list l=primary_char0_no_molien_random(L[1],max,v);
       if (size(l)==2)
-      { matrix S=secondary_no_molien(l[1],L[1],l[2],v);
+      { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v);
       }
       else
-      { matrix S=secondary_no_molien(l[1],L[1],v);
-      }
-      return(l[1],S);
+      { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v);
+      }
+      return(l[1],S,IS);
     }
     else
@@ -7249 +7557 @@
       { list l=primary_charp_no_molien_random(L[1],max,v); // case
         if (size(l)==2)
-        { matrix S=secondary_no_molien(l[1],L[1],l[2],v);
+        { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v);
         }
         else
-        { matrix S=secondary_no_molien(l[1],L[1],v);
-        }
-        return(l[1],S);
+        { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v);
+        }
+        return(l[1],S,IS);
       }
       else                             // the modular case