Changeset 7bcc87 in git
 Timestamp:
 Feb 21, 2007, 2:10:11 PM (16 years ago)
 Branches:
 (u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
 Children:
 119c19a74782ca2e52ffe0f1040a7185f0c028d5
 Parents:
 7a7d813adf208ec7a9f3cb07b00a9dab84ea015e
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/finvar.lib
r7a7d813 r7bcc87 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: finvar.lib,v 1.6 1 20070215 22:47:28king Exp $"2 version="$Id: finvar.lib,v 1.62 20070221 13:10:11 king Exp $" 3 3 category="Invariant theory"; 4 4 info=" … … 16 16 primary_invariants() primary invariants (p.i.) 17 17 primary_invariants_random() primary invariants, randomized alg. 18 invariant_algebra_reynolds() minimal (algebra) generating set for the invariant 19 ring of a finite matrix group, in the nonmodular case 20 invariant_algebra_perm() minimal (algebra) generating set for the invariant 21 ring or a permutation group, in the nonmodular case 18 22 19 23 AUXILIARY PROCEDURES: … … 51 55 relative_orbit_variety() ideal of a relative orbit variety (old version) 52 56 image_of_variety() ideal of the image of a variety 57 orbit_sums orbit sums of a set of monomials under the action of 58 a permutation group 53 59 "; 54 60 /////////////////////////////////////////////////////////////////////////////// … … 64 70 // p_search_random searches a # of p.i., char p, randomized 65 71 // concat_intmat concatenates two integer matrices 72 // GetGroup disjoint cycle presentation of perm. group 73 // > permuation matrices 66 74 /////////////////////////////////////////////////////////////////////////////// 67 75 … … 4428 4436 4429 4437 proc secondary_char0 (matrix P, matrix REY, matrix M, list #) 4430 " 4431 USAGE: secondary_char0(P,REY,M[,v][,\"old\"]); 4438 "USAGE: secondary_char0(P,REY,M[,v][,\"old\"]); 4432 4439 @* P: a 1xn <matrix> with homogeneous primary invariants, where 4433 4440 n is the number of variables of the basering; … … 4883 4890 4884 4891 proc irred_secondary_char0 (matrix P, matrix REY, matrix M, list #) 4885 " 4886 USAGE: irred_secondary_char0(P,REY,M[,v][,\"PP\"]); 4892 "USAGE: irred_secondary_char0(P,REY,M[,v][,\"PP\"]); 4887 4893 @* P: a 1xn <matrix> with homogeneous primary invariants, where 4888 4894 n is the number of variables of the basering; … … 7777 7783 7778 7784 proc relative_orbit_variety(ideal I,matrix F,string newring) 7779 " 7780 USAGE: relative_orbit_variety(I,F,s); 7785 "USAGE: relative_orbit_variety(I,F,s); 7781 7786 I: an <ideal> invariant under the action of a group, 7782 7787 @* F: a 1xm <matrix> defining the invariant ring of this group, … … 7906 7911 /////////////////////////////////////////////////////////////////////////////// 7907 7912 7913 static proc GetMatrix(list GEN) 7914 { matrix M[nvars(basering)][nvars(basering)]; 7915 int i,j; 7916 for (i=1;i<=size(GEN);i++) 7917 { for (j=1;j<size(GEN[i]);j++) 7918 { M[GEN[i][j+1],GEN[i][j]] = 1; 7919 } 7920 M[GEN[i][1],GEN[i][size(GEN[i])]] = 1; 7921 } 7922 return(M); 7923 } 7924 7925 /////////////////////////////////////////////////////////////////////////////// 7926 // MINIMAL GENERATING SETS 7927 /////////////////////////////////////////////////////////////////////////////// 7928 7929 7930 /////////////////////////////////////////////////////////////////////////////// 7931 // Input: Disjoint cycle presentation of a permutation group 7932 // Output: Permuation matrices representing that group 7933 /////////////////////////////////////////////////////////////////////////////// 7934 proc GetGroup(list GRP) 7935 { list L; 7936 int i,j,k; 7937 intvec fix; 7938 for (i=1;i<=size(GRP);i++) 7939 { L[i] = GetMatrix(GRP[i]); 7940 // Fixelemente des Zyklus 7941 fix = intvec(0); 7942 fix[nvars(basering)]=0; 7943 for (j=1;j<=size(GRP[i]);j++) 7944 { for (k=1;k<=size(GRP[i][j]);k++) 7945 { fix[GRP[i][j][k]]=1; 7946 } 7947 } 7948 for (j=1;j<=nvars(basering);j++) 7949 { if (fix[j]==0) 7950 { L[i][j,j]=1; 7951 } 7952 } 7953 } 7954 return(L); 7955 } 7956 7957 /////////////////////////////////////////////////////////////////////////////// 7958 7959 proc orbit_sums (ideal mon, list M) 7960 "USAGE: orbit_sums(mon, Per); 7961 mon is of type <ideal> and is formed by monomials of degree d. 7962 Per is a list of ideals providing a permutation of the ring variables. 7963 ASSUME: mon contains at least the minimal monomial in each orbit sum of degree d 7964 under the permutation group given by Per. 7965 RETURN: an <ideal> whose generators are the orbit sums of degree d under the permutation 7966 group given by Per. 7967 SEE ALSO: invariant_algebra_perm 7968 KEYWORDS: orbit sum 7969 EXAMPLE: example orbit_sums; shows an example 7970 " 7971 { 7972 // checking input  7973 mon=simplify(mon,2); // remove 0 7974 if (not homog(sum(mon))) 7975 { "ERROR: Monomials in the first parameter must be of the same degree."; 7976 return(); 7977 } 7978 int k,k2; 7979 for (k=1;k<=ncols(mon);k++) 7980 { if (leadmonom(mon[k])!=mon[k]) 7981 { "ERROR: The first parameter must be formed by monomials."; 7982 return(); 7983 } 7984 } 7985 for (k2=1;k2<=size(M);k2++) 7986 { if (typeof(M[k2])!="ideal") 7987 { "ERROR: The second parameter must be a list of ideals providing a "; 7988 " permutation of the ring variables."; 7989 return(); 7990 } 7991 if ((ncols(M[k2])!=nvars(basering)) 7992 or (ncols(M[k2])!=ncols(simplify(M[k2],6))) 7993 or (M[k2][1]==0)) 7994 { "ERROR: The second parameter must be a list of ideals providing a "; 7995 " permutation of the ring variables."; 7996 return(); 7997 } 7998 for (k=1;k<=ncols(M[k2]);k++) 7999 { if ((leadmonom(M[k2][k])!=M[k2][k]) or (deg(M[k2][k])!=1)) 8000 { "ERROR: The second parameter must be a list of ideals providing a "; 8001 " permutation of the ring variables."; 8002 return(); 8003 } 8004 } 8005 } 8006 // preparing steps  8007 int dgb = degBound; 8008 degBound = 0; 8009 def br = basering; 8010 ideal tmpI; 8011 int g; 8012 // rings that allow fast permutation of variables with >fetch< and >imap< 8013 for (g=1;g<=size(M);g++) 8014 { list tmpL(g); 8015 execute("ring tmpR("+string(g)+") = "+charstr(br)+", ("+string(M[g])+"),("+ordstr(br)+")"); 8016 ideal tmpI; poly OrbitS; list L(g); 8017 setring(br); 8018 } 8019 // Remove nonminimal monomials. 8020 // Strategy: If m1, m2 are monomials related by a group generator 8021 // then remember only the smaller of the two. 8022 // This is where we use the assumption on mon! 8023 for (g=1;g<=size(M);g++) 8024 { setring tmpR(g); 8025 tmpI=fetch(br,mon); 8026 setring br; 8027 tmpI=imap(tmpR(g),tmpI); 8028 for (k=1;k<=ncols(mon);k++) 8029 { if ( mon[k]!=0) 8030 { if (tmpI[k] < mon[k]) 8031 { 8032 mon[k]=0; 8033 } 8034 else 8035 { if (tmpI[k]>mon[k]) 8036 { 8037 mon=NF(mon,std(tmpI[k])); 8038 } 8039 } 8040 } 8041 } 8042 } 8043 mon=simplify(mon,2); 8044 8045 // main part: form orbits 8046 list L,newL; // contains partial orbits 8047 for (k=1;k<=ncols(mon);k++) 8048 { L[k]=ideal(mon[k]); 8049 attrib(L[k],"isSB",1); 8050 } 8051 newL=L; 8052 poly tmp; 8053 int GenOK, j,i; 8054 poly OrbitS=mon[1]; 8055 mon[1]=0; 8056 int countOrbits; 8057 ideal AllOrbits; 8058 while (1) // apply all generators to the monomials that have been found 8059 // in the preceding step 8060 { GenOK = 0; 8061 for (g=1;g<=size(M);g++) // Generator g is applied to partial orbits 8062 { setring tmpR(g); 8063 L(g)=fetch(br,newL); 8064 setring br; 8065 tmpL(g)=imap(tmpR(g),L(g)); 8066 } 8067 for (k=1;k<=size(L);k++) 8068 { if (L[k][1]!=0) // if the orbit was not fused or completed before 8069 { newL[k] = tmpL(1)[k]; 8070 for (g=2;g<=size(M);g++) 8071 { newL[k] = newL[k]+tmpL(g)[k]; 8072 } 8073 newL[k] = simplify(NF(newL[k],L[k]),2); 8074 attrib(newL[k],"isSB",1); 8075 if (newL[k][1] != 0) // i.e., the Orbit is not complete under all generators yet 8076 { GenOK=1; 8077 L[k] = L[k],newL[k]; 8078 attrib(L[k],"isSB",1); 8079 // fuse orbits: 8080 for (k2=k+1;k2<=size(L);k2++) 8081 { if (size(NF(newL[k2],L[k]))!=size(newL[k2])) 8082 { L[k] = L[k]+L[k2]; 8083 attrib(L[k],"isSB",1); 8084 newL[k] = newL[k]+newL[k2]; 8085 attrib(newL[k],"isSB",1); 8086 L[k2]=ideal(0); 8087 newL[k2]=ideal(0); 8088 } 8089 } // fuse orbits 8090 } // the orbit was not complete yet 8091 else 8092 { countOrbits++; 8093 AllOrbits[countOrbits]=sum(L[k]); 8094 L[k]=ideal(0); 8095 newL[k]=ideal(0); 8096 } // the orbit was complete 8097 } // if orbit was not yet fused/completed before 8098 } // loop through partial orbits 8099 if (GenOK == 0) // i.e., all orbits are complete, so we can break the iteration 8100 { break; 8101 } 8102 } // iterate 8103 degBound = dgb; 8104 return(AllOrbits); 8105 } 8106 example 8107 { "EXAMPLE: orbits sums in degree 3 under the action of the alternating group;"; 8108 " it doesn't matter that we are in the modular case"; echo=2; 8109 ring R=3,(a,b,c,d),dp; 8110 list Per=ideal(b,c,a,d),ideal(a,c,d,b); 8111 ideal orb=orbit_sums(kbase(std(0),2),Per); 8112 orb; 8113 } 8114 8115 /////////////////////////////////////////////////////////////////////////////// 8116 8117 proc invariant_algebra_perm (list #) 8118 "USAGE: invariant_algebra_perm(GEN[,v]); 8119 @* GEN: a list of generators of a permutation group. It is given in disjoint cycle 8120 form, where trivial cycles can be omitted; e.g., the generator (1,2)(3,4)(5) is 8121 given by <list(list(1,2),list(3,4))>. 8122 @* v: an optional <int> 8123 RETURN: A minimal homogeneous generating set of the invariant ring of the group presented 8124 by GEN, type <matrix> 8125 ASSUME: We are in the nonmodular case, i.e., the characteristic of the basering 8126 does not divide the group order. Note that the function does not verify whether 8127 the assumption holds or not 8128 DISPLAY: Information on the progress of computations if v does not equal 0 8129 THEORY: We do an incremental search in increasing degree d. Generators of the invariant 8130 ring are found among the orbit sums of degree d. The generators are chosen by 8131 Groebner basis techniques. 8132 SEE ALSO: invariant_algebra_reynolds 8133 KEYWORDS: invariant ring minimal generating set permutation group 8134 EXAMPLE: example invariant_algebra_perm; shows an example 8135 " 8136 { 8137 // checking input and setting verbose mode  8138 if (size(#)==0) 8139 { "ERROR: There are no parameters given."; 8140 return(); 8141 } 8142 if (typeof(#[size(#)])=="int") 8143 { int v=#[size(#)]; 8144 list GEN=#[1..size(#)1]; 8145 } 8146 else 8147 { int v=0; 8148 list GEN=#; 8149 } 8150 int i,j,k; 8151 for (i=1;i<=size(GEN);i++) 8152 { if (typeof(GEN[i])!="list") 8153 { "ERROR: The permutation group must be given by generators in"; 8154 " disjoint cycle presentation"; 8155 return(); 8156 } 8157 for (j=1;j<=size(GEN[i]);j++) 8158 { if (typeof(GEN[i][j])!="list") 8159 { "ERROR: The permutation group must be given by generators in"; 8160 " disjoint cycle presentation"; 8161 return(); 8162 } 8163 for (k=1;k<=size(GEN[i][j]);k++) 8164 { if (typeof(GEN[i][j][k])!="int") 8165 { "ERROR: The permutation group must be given by generators in"; 8166 " disjoint cycle presentation"; 8167 return(); 8168 } 8169 } 8170 } 8171 } 8172 def TST=sort(sum(sum(GEN)))[1]; 8173 if (TST[1]<=0) 8174 { "ERROR: The permutation group must be given by generators in"; 8175 " disjoint cycle presentation"; 8176 return(); 8177 } 8178 for (i=1;i<size(TST);i++) 8179 { if ((TST[i]==TST[i+1]) or (TST[i+1]>nvars(basering))) 8180 { "ERROR: The permutation group must be given by generators in"; 8181 " disjoint cycle presentation"; 8182 return(); 8183 } 8184 } 8185 8186 // starting computation  8187 list GRP = GetGroup(GEN); 8188 def br=basering; 8189 int dgb = degBound; 8190 degBound = 0; 8191 int counter, totalcount, OrbSize, 8192 NextCand, sGoffset, fin,MaxD,LastGDeg; 8193 ideal OrbSums, RedSums, NewReductor; 8194 poly helpP,lmp; 8195 8196 8197 // Transform the matrices generating the group into ring maps 8198 int maxG = size(GRP); 8199 8200 matrix vars=transpose(matrix(maxideal(1))); 8201 8202 ideal tmpI; 8203 int g; 8204 list M; // M provides maps corresponding to the group generators: 8205 for (i=1;i<=maxG;i++) 8206 { tmpI = ideal(GRP[i]*vars); 8207 M[i] = tmpI; 8208 } 8209 8210 ideal G; // G will contain homog. alg.generators 8211 poly Inv; 8212 ideal mon; 8213 8214 ideal sG = std(G); 8215 int d=1; 8216 while ((fin==0) or (d<=MaxD)) 8217 { if (v) 8218 { "Searching generators in degree",d; 8219 } 8220 counter = 0; // counts generators in degree d 8221 OrbSums = orbit_sums(kbase(sG,d),M); 8222 RedSums = reduce(OrbSums,sG); 8223 sGoffset=ncols(sG); 8224 OrbSize = ncols(OrbSums); 8225 if (v) 8226 { " We have "+string(OrbSize)+" orbit sums of degree "+string(d); 8227 } 8228 // loop through orbit_sum, and select generators among them 8229 for (i=1;i<=OrbSize;i++) 8230 { helpP = RedSums[i]; 8231 if (helpP<>0) 8232 { counter++; // found new inv. algebra generator! 8233 totalcount++; 8234 if (v) 8235 { " We found generator number "+string(totalcount)+ 8236 " in degree "+string(d); 8237 } 8238 G[totalcount] = OrbSums[i]; 8239 sG[sGoffset+counter] = helpP; 8240 attrib(sG,"isSB",1); 8241 NewReductor[1]=helpP; 8242 attrib(NewReductor,"isSB",1); 8243 RedSums = reduce(RedSums,NewReductor); 8244 LastGDeg=d; 8245 } // found new generator 8246 } // loop through orbit sums 8247 if ((MaxD>0) and (d==MaxD)) 8248 { if (v) 8249 { "We went beyond the degree bound, so we are done!"; 8250 } 8251 break; 8252 } 8253 d++; 8254 // Prepare computations of the next degree 8255 if (LastGDeg==d2) 8256 { degBound=0; 8257 if (v) {"Presumably all generators have been found";} 8258 sG = groebner(sG); 8259 fin = 0; 8260 } 8261 if (LastGDeg==d1) 8262 { degBound=d; 8263 if (v) { "Computing Groebner basis up to the new degree",d; } 8264 sG = groebner(sG); 8265 attrib(sG,"isSB",1); 8266 fin = 0; 8267 } 8268 if ((fin==0) and (dim(sG)==0)) 8269 { MaxD = maxdeg1(kbase(sG)); 8270 if (v) 8271 { if (MaxD>=d) {"We found the degree bound",MaxD; } 8272 else {"We found the degree bound",d1;} 8273 } 8274 if (MaxD<d) 8275 { if (v) {"We went beyond the degree bound, so, we are done!";} 8276 break; 8277 } 8278 fin = 1; 8279 } // computing degree bound 8280 } // for increasing degrees 8281 kill sG; 8282 degBound=dgb; 8283 return(matrix (G)); 8284 } 8285 example 8286 { "EXAMPLE: A nontransitive action of S_2 on four variables"; echo=2; 8287 ring R=0,(a,b,c,d),dp; 8288 def GEN=list(list(list(1,3),list(2,4))); 8289 matrix G = invariant_algebra_perm(GEN,1); 8290 G; 8291 } 8292 8293 /////////////////////////////////////////////////////////////////////////////// 8294 8295 proc invariant_algebra_reynolds (matrix REY, list #) 8296 "USAGE: invariant_algebra_reynolds(REY[,v]); 8297 @* REY: a gxn <matrix> representing the Reynolds operator of a finite matrix group, 8298 where g ist the group order and n is the number of variables of the basering; 8299 @* v: an optional <int> 8300 RETURN: A minimal homogeneous generating set of the invariant ring, type <matrix> 8301 ASSUME: We are in the nonmodular case, i.e., the characteristic of the basering 8302 does not divide the group order; 8303 REY is the 1st return value of group_reynolds(), reynolds_molien() or 8304 the second one of primary_invariants() 8305 DISPLAY: Information on the progress of computations if v does not equal 0 8306 THEORY: We do an incremental search in increasing degree d. Generators of the invariant 8307 ring are found among the Reynolds images of monomials of degree d. The generators are 8308 chosen by Groebner basis techniques. 8309 SEE ALSO: invariant_algebra_perm 8310 KEYWORDS: invariant ring minimal generating set matrix group 8311 EXAMPLE: example invariant_algebra_reynolds; shows an example 8312 " 8313 { int MonStep = 10; 8314 int v=1; 8315 def br=basering; 8316 int dgb = degBound; 8317 degBound = 0; 8318 int block, counter, totalcount, monsize,i,j, 8319 sGoffset, fin,MaxD,LastGDeg; 8320 poly helpP,lmp; 8321 block = 1; 8322 ideal mon, MON, Inv, G; // G will contain homog. alg.generators 8323 ideal RedInv, NewG, NewReductor; 8324 int NewGcounter; 8325 intvec dimvec; 8326 8327 ideal sG = std(G); 8328 int d=1; 8329 dimvec[d] = dim(sG); 8330 while ((fin==0) or (d<=MaxD)) 8331 { NewG = ideal(); 8332 NewGcounter=0; 8333 mon = kbase(sG,d); 8334 sGoffset=ncols(sG); 8335 monsize = ncols(mon); 8336 if (mon[1]<>0) 8337 { if (v) { " We have "+string(monsize)+ 8338 " relevant monomials in degree "+string(d);} 8339 block=0; // Loops through the monomials 8340 while (block<monsize) 8341 { if ((block mod MonStep) == 0) 8342 { if ((block+MonStep) <= monsize) 8343 { ideal tmp = normalize(evaluate_reynolds(REY,ideal(mon[block+1..block+MonStep]))); 8344 Inv[block+1..block+MonStep] = tmp[1..MonStep]; 8345 tmp = reduce(tmp,sG); 8346 RedInv[block+1..block+MonStep] = tmp[1..MonStep]; 8347 kill tmp; 8348 } 8349 else 8350 { ideal tmp = normalize(evaluate_reynolds(REY,ideal(mon[block+1..monsize]))); 8351 Inv[block+1..monsize] = tmp[1..monsizeblock]; 8352 tmp = reduce(tmp,sG); 8353 RedInv[block+1..monsize] = tmp[1..monsizeblock]; 8354 kill tmp; 8355 } 8356 } 8357 block++; 8358 helpP = RedInv[block]; 8359 if (helpP<>0) 8360 { counter++; // found new inv. algebra generator! 8361 if (v) { " We found generator number "+string(counter)+ 8362 " in degree "+string(d); } 8363 G[counter] = Inv[block]; 8364 sG[sGoffset+counter] = helpP; 8365 attrib(sG,"isSB",1); 8366 NewGcounter++; 8367 NewG[NewGcounter] = helpP; 8368 mon[block]=0; 8369 Inv[block]=0; 8370 RedInv[block]=0; 8371 NewReductor[1]=helpP; 8372 attrib(NewReductor,"isSB",1); 8373 RedInv = reduce(RedInv,NewReductor); 8374 LastGDeg=d; 8375 } 8376 } // loop through monomials 8377 } // if mon[1]<>0 8378 if ((MaxD>0) and (d==MaxD)) 8379 { if (v) 8380 { "We went beyond the degree bound, so we are done!"; 8381 } 8382 break; 8383 } 8384 d++; 8385 // Prepare computations of the next degree 8386 if (LastGDeg==d2) 8387 { degBound=0; 8388 if (v) {"Presumably all generators have been found";} 8389 sG = groebner(sG); 8390 fin = 0; 8391 } 8392 if (LastGDeg==d1) 8393 { degBound=d; 8394 if (v) { "Computing Groebner basis up to the new degree",d; } 8395 sG = groebner(sG); 8396 attrib(sG,"isSB",1); 8397 fin = 0; 8398 } 8399 if ((fin==0) and (dim(sG)==0)) 8400 { MaxD = maxdeg1(kbase(sG)); 8401 if (v) 8402 { if (MaxD>=d) {"We found the degree bound",MaxD; } 8403 else {"We found the degree bound",d1;} 8404 } 8405 if (MaxD<d) 8406 { if (v) {"We went beyond the degree bound, so, we are done!";} 8407 break; 8408 } 8409 fin = 1; 8410 } // computing degree bound 8411 } // for increasing degrees 8412 kill sG; 8413 degBound=dgb; 8414 return(matrix (G)); 8415 } 8416 example 8417 { "EXAMPLE: A nontransitive action of S_2 on four variables"; echo=2; 8418 ring R=0,(a,b,c,d),dp; 8419 matrix A[4][4]= 8420 0,0,1,0, 8421 0,0,0,1, 8422 1,0,0,0, 8423 0,1,0,0; 8424 list L = group_reynolds(A); 8425 matrix G = invariant_algebra_reynolds(L[1],1); 8426 G; 8427 } 8428 8429 ///////////////////////////////////////////////////////////////////////////////
Note: See TracChangeset
for help on using the changeset viewer.