Changeset e0e0e2 in git for Singular/LIB/finvar.lib

Timestamp:

Dec 14, 2007, 11:12:31 AM (16 years ago)

Author:

Simon King <king@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

4f818c4604d3f8dd7bb40260971c1bdd9f709746

Parents:

d617b278e0d4576d32b4c6f260fd73bae2daa907

Message:

Improvements of various details of secondary_not_cohen_macaulay. New (experimental) option for Hilbert driven computation. Still, magma is *much* faster. This ought to change!


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

File:

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

Singular/LIB/finvar.lib

@@ -1 +1 @@
 ///////////////////////////////////////////////////////////////////////////////
-version="$Id: finvar.lib,v 1.70 2007-12-12 15:02:46 king Exp $"
+version="$Id: finvar.lib,v 1.71 2007-12-14 10:12:31 king Exp $"
 category="Invariant theory";
 info="
@@ -6947 +6947 @@
 
 proc secondary_not_cohen_macaulay (matrix P, list #)
-"USAGE:   secondary_not_cohen_macaulay(P,G1,G2,...[,v]);
+"USAGE:   secondary_not_cohen_macaulay(P,G1,G2,...[,v,h]);
          P: a 1xn <matrix> with primary invariants, G1,G2,...: nxn <matrices>
-         generating a finite matrix group, v: an optional <int>
+         generating a finite matrix group, v,h: optional <int>
 ASSUME:  n is the number of variables of the basering
 RETURN:  secondary invariants of the invariant ring (type <matrix>)
@@ -6955 +6955 @@
 THEORY:  Secondary invariants are generated following \"Generating Invariant
          Rings of Finite Groups over Arbitrary Fields\" by Kemper (1996).
+         Hilbert driven computation is used if h does not equal 0
 EXAMPLE: example secondary_not_cohen_macaulay; shows an example
 "
 { int i, j;
+  int v, useH;
   degBound=0;
   def br=basering;
@@ -6966 +6968 @@
   if (size(#)>0)                       // checking input and setting verbose
   { if (typeof(#[size(#)])=="int")
-    { int gen_num=size(#)-1;
+    { if ((size(#)>1) and (typeof(#[size(#)-1])=="int"))
+      { int gen_num=size(#)-2;
+        v=#[size(#)-1];
+        useH=#[size(#)];
+      }
+      else
+      { int gen_num=size(#)-1;
+        v=#[size(#)];
+      }
       if (gen_num==0)
       { "ERROR:   There are no generators of the finite matrix group given.";
         return();
       }
-      int v=#[size(#)];
       for (i=1;i<=gen_num;i++)
       { if (typeof(#[i])<>"matrix")
@@ -6984 +6993 @@
     }
     else
-    { int v=0;
-      int gen_num=size(#);
+    { int gen_num=size(#);
       for (i=1;i<=gen_num;i++)
       { if (typeof(#[i])<>"matrix")
@@ -7010 +7018 @@
   { "";
   }
-  ring alskdfalkdsj=0,x,dp;
-  matrix M[1][2]=1,(1-x)^n;            
-  export alskdfalkdsj;
-  export M;                            // we look at our primary invariants as
-  setring br;                          // such of the subgroup that only
-  matrix REY=matrix(maxideal(1));      // contains the identity, this means that
+  // ring alskdfalkdsj=0,x,dp;
+  // matrix M[1][2]=1,(1-x)^n;            
+  // export alskdfalkdsj;
+  // export M;                            // we look at our primary invariants as
+  // setring br;                          // such of the subgroup that only
+  // matrix REY=matrix(maxideal(1));      // contains the identity, this means that
                                        // ch does not divide the order anymore,
                                        // this means that we can make use of the
@@ -7025 +7033 @@
                                        // ones from the bigger group
   if (v)
-  { "  The procedure secondary_charp() is called to calculate secondary invariants";
-    "  of the invariant ring of the trivial group with respect to the primary";
-    "  invariants found previously.";
+  { "  First, get secondary invariants for the trivial action, with respect to";
+    "  the primary invariants found previously. By the graded Nakayama lemma,";
+    "  this means to compute the standard monomials for the ideal generated by";
+    "  the primary invariants.";
     "";
   }
-  matrix trivialS, trivialIS=secondary_charp(P,REY,"alskdfalkdsj",v);
-  kill trivialIS;
-  kill alskdfalkdsj;
+  ideal GP = groebner(ideal(P));
+  if (dim(GP)<>0)
+  { "ERROR:   The ideal spanned by primary invariants of a finite group";
+    "         always is of dimension zero.";
+    return();
+  }
+  matrix trivialS = matrix(kbase(GP));// , trivialIS=secondary_charp(P,REY,"alskdfalkdsj",v);
+  // kill trivialIS;
+  // kill alskdfalkdsj;
   // now we have those secondary invariants
   int k=ncols(trivialS);               // k is the number of the secondary
@@ -7057 +7072 @@
 //  module M(2)=syz(M(1));               // nres(M(1),2)[2];
 //  option(set,save_opts);
+  if (v) 
+  { "Syzygies of the \"trivial\" secondaries"; }
   module M(2)=nres(M(1),2)[2];
   int m=ncols(M(2));                   // number of generators of the module
@@ -7063 +7080 @@
   // the algebra A^k where A denote the subalgebra of the usual polynomial
   // ring, generated by the primary invariants
-  string mp=string(minpoly);           // generating a ring where we can do
-                                       // elimination
-  execute("ring R=("+charstr(br)+"),(x(1..n),y(1..n),h),dp;");
-  execute("minpoly=number("+mp+");");
-  // map f=br,maxideal(1);                // canonical mapping
-  matrix M[k][m+k*n];
-  M[1..k,1..m]=matrix(fetch(br,M(2)));        // will contain a module -
+
+  // generate a ring where we can do elimination  
+  // in order to make it weighted homogeneous, we choose appropriate weights
+  if (v)
+  { "By elimination: Intersect with the invariant ring"; }
+  intvec dv = degvec(ideal(P)),1;
+  ring newR = char(br), (y(1..n),h),(a(dv),dp);
+  def R = br + newR;
+  setring R;
   matrix P=fetch(br,P);                       // primary invariants in the new ring
+  ideal ElimP;
   for (i=1;i<=n;i++)
+  { ElimP[i]=y(i)-P[1,i]; }
+  def RHilb = newR+br; // the new variables only appear linearly, so that block order wouldn't hurt
+  setring RHilb;
+  ideal ElimP=imap(R,ElimP);
+  intvec EPvec = hilb(groebner(ElimP),1,degvec(maxideal(1)));
+  setring R;
+  ideal GElimP = std(ElimP,EPvec,degvec(maxideal(1)));
+  attrib(GElimP,"isSB",1);
+  int newN = ncols(GElimP);  
+  matrix M[k][m+k*newN];
+  M[1..k,1..m]=matrix(imap(br,M(2)));        // will contain a module -
+  for (i=1;i<=newN;i++)
   { for (j=1;j<=k;j++)
-    { M[j,m+(i-1)*k+j]=y(i)-P[1,i];
-    }
-  }
-  M=elim(module(M),1,n);               // eliminating x(1..n), std-calculation
-                                       // is done internally -
-  M=homog(module(M),h);                // homogenize for 'minbase'
+    { M[j,m+(i-1)*k+j]=GElimP[i]; // y(i)-P[1,i];
+    }
+  }
+
+  if (not useH)                             
+  { M=elim(module(M),1,n);               // eliminating x(1..n), std-calculation
+                                         // is done internally -
+    M=homog(module(M),h);            // homogenize for 'minbase'
+  }
+  else
+  {
+    // Manually, as long as Hilbert driven eliminations are not supported for modules.
+    // ----- Compute the Hilbert function, using an "easy" ordering ----
+    if (v)
+    { "   Compute Hilbert Function"; }
+    setring RHilb;
+    matrix M=imap(R,M);
+    module GM=groebner(module(M));
+    intvec hv = hilb(GM,1,degvec(maxideal(1)));
+    setring R;
+    module GM=imap(RHilb,GM);
+    kill RHilb; 
+    // ----- Hilbert driven standard basis is hopefully reasonably fast ----
+    if (v)
+    { "   Hilbert driven standard basis computation"; }
+    module MGb = std(module(GM),hv,degvec(maxideal(1)));
+    ideal vars;
+    vars[n]=0;
+    vars=vars,y(1..n),h;
+    map emb=R,vars;
+    module MElim=emb(MGb);
+    MGb = MGb - MElim;
+    for (j=1;j<=ncols(MGb);j++)
+    { if (MGb[j]<>0)
+      { MElim[j]=0;
+      }
+    }
+    M=homog(module(MElim),h);                // homogenize for 'minbase'
+  }
+
   M=minbase(module(M));
-  setring br;
-  ideal substitute=maxideal(1),ideal(P),1;
-  f=R,substitute;                      // replacing y(1..n) by primary
+  ideal substitute;
+  for (i=1;i<=n;i++)
+  { substitute[i]=var(i);
+  }
+  substitute = substitute,ideal(P),1;
+  map f=R,substitute;                      // replacing y(1..n) by primary
                                        // invariants and h by 1 -
-  M(2)=f(M);                           // M(2) is the new module
-  m=ncols(M(2));
+  M=f(M);                              // we can't directly map from R to br
+  setring br;                          // if there is a minpoly.
+  module M(3)=imap(R,M);                     // M(2) is the new module
+  m=ncols(M(3));
   matrix S[1][m];
-  S=matrix(trivialS)*matrix(M(2));     // S now contains the secondary
+  S=matrix(trivialS)*matrix(M(3));     // S now contains the secondary
                                        // invariants
   for (i=1; i<=m;i++)