/////////////////////////////////////////////////////////////////////////////// version="$Id: finvar.lib,v 1.80 2009-01-07 16:11:36 Singular Exp $" category="Invariant theory"; info=" LIBRARY: finvar.lib Invariant Rings of Finite Groups AUTHOR: Agnes E. Heydtmann, email: agnes@math.uni-sb.de; Simon A. King, email: king@mathematik.uni-jena.de OVERVIEW: A library for computing polynomial invariants of finite matrix groups and generators of related varieties. The algorithms are based on B. Sturmfels, G. Kemper, S. King and W. Decker et al.. MAIN PROCEDURES: invariant_ring() generators of the invariant ring (i.r.) invariant_ring_random() generators of the i.r., randomized alg. primary_invariants() primary invariants (p.i.) primary_invariants_random() primary invariants, randomized alg. invariant_algebra_reynolds() minimal generating set for the invariant ring of a finite matrix group, in the non-modular case invariant_algebra_perm() minimal generating set for the invariant ring or a permutation group, in the non-modular case AUXILIARY PROCEDURES: cyclotomic() cyclotomic polynomial group_reynolds() finite group and Reynolds operator (R.o.) molien() Molien series (M.s.) reynolds_molien() Reynolds operator and Molien series partial_molien() partial expansion of Molien series evaluate_reynolds() image under the Reynolds operator invariant_basis() basis of homogeneous invariants of a degree invariant_basis_reynolds() as invariant_basis(), with R.o. primary_char0() primary invariants (p.i.) in char 0 primary_charp() primary invariants in char p primary_char0_no_molien() p.i., char 0, without Molien series primary_charp_no_molien() p.i., char p, without Molien series primary_charp_without() p.i., char p, without R.o. or Molien series primary_char0_random() primary invariants in char 0, randomized primary_charp_random() primary invariants in char p, randomized primary_char0_no_molien_random() p.i., char 0, without M.s., randomized primary_charp_no_molien_random() p.i., char p, without M.s., randomized primary_charp_without_random() p.i., char p, without R.o. or M.s., random. power_products() exponents for power products secondary_char0() secondary invariants (s.i.) in char 0 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 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) relative_orbit_variety() ideal of a relative orbit variety (old version) image_of_variety() ideal of the image of a variety orbit_sums orbit sums of a set of monomials under the action of a permutation group "; /////////////////////////////////////////////////////////////////////////////// // perhaps useful procedures (no help provided): // unique() is a matrix among other matrices? // exponent() gives the exponent of a number // 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 mapped to q field elements // search searches a number of p.i., char 0 // p_search searches a number of p.i., char p // search_random searches a # of p.i., char 0, randomized // p_search_random searches a # of p.i., char p, randomized // concat_intmat concatenates two integer matrices // GetGroup gets disjoint cycle presentation of perm. group, // yields permuation matrices /////////////////////////////////////////////////////////////////////////////// LIB "matrix.lib"; LIB "elim.lib"; LIB "general.lib"; LIB "algebra.lib"; /////////////////////////////////////////////////////////////////////////////// // Checks whether the last parameter, being a matrix, is among the previous // parameters, also being matrices /////////////////////////////////////////////////////////////////////////////// proc unique (list #) { int s=size(#); def m=#[s]; for (int i=1;i 0 RETURNS: the i-th cyclotomic polynomial (type ) as one in the first ring variable THEORY: x^i-1 is divided by the j-th cyclotomic polynomial where j takes on the value of proper divisors of i EXAMPLE: example cyclotomic; shows an example " { if (i<=0) { "ERROR: the input should be > 0."; return(); } poly v1=var(1); if (i==1) { return(v1-1); // 1-st cyclotomic polynomial } poly min=v1^i-1; matrix s[1][2]; min=min/(v1-1); // dividing by the 1-st cyclotomic // polynomial int j=2; int n; poly c; int flag=1; while(2*j<=i) // there are no proper divisors of i { if ((i%j)==0) // greater than i/2 { if (flag==1) { n=j; // n stores the first proper divisor of } // i > 1 flag=0; c=cyclotomic(j); // recursive computation s=min,c; s=matrix(syz(ideal(s))); // dividing min=s[2,1]; } if (n*j==i) // the earliest possible point to break { break; } j++; } min=min/leadcoef(min); // making sure that the leading return(min); // coefficient is 1 } example { "EXAMPLE:"; echo=2; ring R=0,(x,y,z),dp; print(cyclotomic(25)); } proc group_reynolds (list #) "USAGE: group_reynolds(G1,G2,...[,v]); G1,G2,...: nxn generating a finite matrix group, v: an optional ASSUME: n is the number of variables of the basering, g the number of group elements RETURN: a , the first list element will be a gxn representing the Reynolds operator if we are in the non-modular case; if the characteristic is >0, minpoly==0 and the finite group non-cyclic the second list element is an giving the lowest common multiple of the matrix group elements' order (used in molien); in general all other list elements are nxn listing all elements of the finite group DISPLAY: information if v does not equal 0 THEORY: The entire matrix group is generated by getting all left products of generators with the new elements from the last run through the loop (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 found in one run. Additionally each time a new group element is found the corresponding ring mapping of which the Reynolds operator is made up is generated. They are stored in the rows of the first return value. EXAMPLE: example group_reynolds; shows an example " { int ch=char(basering); // the existance of the Reynolds operator // is dependent on the characteristic of // the base field int gen_num; // number of generators //------------------------ making sure the input is okay --------------------- if (typeof(#[size(#)])=="int") { if (size(#)==1) { "ERROR: there are no matrices given among the parameters"; return(); } int v=#[size(#)]; gen_num=size(#)-1; } else // last parameter is not { int v=0; // no information is default gen_num=size(#); } if (typeof(#[1])<>"matrix") { "ERROR: The parameters must be a list of matrices and maybe an "; return(); } 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(); } if (n<>ncols(#[1])) { "ERROR: matrices need to be square and of the same dimensions"; return(); } matrix vars=matrix(maxideal(1)); // creating an nx1-matrix containing the vars=transpose(vars); // variables of the ring - matrix REY=#[1]*vars; // calculating the first ring mapping - // REY will contain the Reynolds // operator - matrix G(1)=#[1]; // G(k) are elements of the group - if (ch<>0 && minpoly==0 && gen_num<>1) // finding out of which order the { matrix I=diag(1,n); // group element is matrix TEST=G(1); int o1=1; int o2; while (TEST<>I) { TEST=TEST*G(1); o1++; } } int i=1; // -------------- doubles among the generators should be avoided ------------- for (int j=2;j<=gen_num;j++) // this loop adds the parameters to the { // group, leaving out doubles and // checking whether the parameters are // compatible with the task of the // procedure if (not(typeof(#[j])=="matrix")) { "ERROR: The parameters must be a list of matrices and maybe an "; return(); } if ((n!=nrows(#[j])) or (n!=ncols(#[j]))) { "ERROR: matrices need to be square and of the same dimensions"; return(); } if (unique(G(1..i),#[j])) { i++; matrix G(i)=#[j]; if (ch<>0 && minpoly==0) // finding out of which order the group { TEST=G(i); // element is o2=1; while (TEST<>I) { TEST=TEST*G(i); o2++; } o1=o1*o2/gcd(o1,o2); // lowest common multiple of the element } // orders - REY=concat(REY,#[j]*vars); // adding ring homomorphisms to REY } } int g=i; // G(1)..G(i) are generators without // doubles - g generally is the number // of elements in the group so far - j=i; // j is the number of new elements that // we use as factors int k, m, l; if (v) { ""; " Generating the entire matrix group and the Reynolds operator..."; ""; } // -------------- main loop that finds all the group elements ---------------- while (1) { l=0; // l is the number of products we get in // one going for (m=g-j+1;m<=g;m++) { for (k=1;k<=i;k++) { l=l+1; matrix P(l)=G(k)*G(m); // possible new element } } j=0; for (k=1;k<=l;k++) { if (unique(G(1..g),P(k))) { j++; // a new factor for next run g++; matrix G(g)=P(k); // a new group element - if (ch<>0 && minpoly==0 && i<>1) // finding out of which order the { TEST=G(g); //group element is o2=1; while (TEST<>I) { TEST=TEST*G(g); o2++; } o1=o1*o2/gcd(o1,o2); // lowest common multiple of the element } // orders - REY=concat(REY,P(k)*vars); // adding new mapping to REY if (v) { " Group element ",g," has been found."; } } kill P(k); } if (j==0) // when we didn't add any new elements { break; // in one run through the while loop } // we are done } if (v) { if (g<=i) { " There are only ",g," group elements."; } ""; } REY=transpose(REY); // when we evaluate the Reynolds operator // later on, we actually want 1xn // matrices if (ch<>0) { if ((g%ch)==0) { if (voice==2) { "WARNING: The characteristic of the coefficient field divides the group order."; " Proceed without the Reynolds operator!"; } else { if (v) { " The characteristic of the base field divides the group order."; " We have to continue without Reynolds operator..."; ""; } } kill REY; matrix REY[1][1]=0; return(list(REY,G(1..g))); } if (minpoly==0) { if (i>1) { return(list(REY,o1,G(1..g))); } return(list(REY,G(1..g))); } } if (v) { " Done generating the group and the Reynolds operator."; ""; } return(list(REY,G(1..g))); } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2; ring R=0,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; list L=group_reynolds(A); print(L[1]); print(L[2..size(L)]); } /////////////////////////////////////////////////////////////////////////////// // Returns i such that root^i==n, i.e. it heavily relies on the right input. /////////////////////////////////////////////////////////////////////////////// proc exponent(number n, number root) { int i=0; while((n/root^i)<>1) { i++; } return(i); } /////////////////////////////////////////////////////////////////////////////// proc molien (list #) "USAGE: molien(G1,G2,...[,ringname,lcm,flags]); G1,G2,...: nxn , all elements of a finite matrix group, ringname: a giving a name for a new ring of characteristic 0 for the Molien series in case of prime characteristic, lcm: an giving the lowest common multiple of the elements' orders in case of prime characteristic, minpoly==0 and a non-cyclic group, flags: an optional with three components: if the first element is not equal to 0 characteristic 0 is simulated, i.e. the Molien series is computed as if the base field were characteristic 0 (the user must choose a field of large prime characteristic, e.g. 32003), the second component should give the size of intervals between canceling common factors in the expansion of the Molien series, 0 (the default) means only once after generating all terms, in prime characteristic also a negative number can be given to indicate that common factors should always be canceled when the expansion is simple (the root of the extension field does not occur among the coefficients) ASSUME: n is the number of variables of the basering, G1,G2... are the group elements generated by group_reynolds(), lcm is the second return value of group_reynolds() RETURN: in case of characteristic 0 a 1x2 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 Molien series (named M) is stored DISPLAY: information if the third component of flags does not equal 0 THEORY: In characteristic 0 the terms 1/det(1-xE) for all group elements of the Molien series are computed in a straight forward way. In prime characteristic a Brauer lift is involved. The returned matrix gives enumerator and denominator of the expanded version where common factors have been canceled. EXAMPLE: example molien; shows an example " { def br=basering; // the Molien series depends on the int ch=char(br); // characteristic of the coefficient // field - int g; // size of the group //---------------------- making sure the input is okay ----------------------- if (typeof(#[size(#)])=="intvec") { if (size(#[size(#)])==3) { int mol_flag=#[size(#)][1]; if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0))) { "ERROR: the second component of should be >=0" return(); } int interval=#[size(#)][2]; int v=#[size(#)][3]; } else { "ERROR: should have three components"; return(); } if (ch<>0) { if (typeof(#[size(#)-1])=="int") { int r=#[size(#)-1]; if (typeof(#[size(#)-2])<>"string") { "ERROR: In characteristic p>0 a must be given for the name of a new"; " ring where the Molien series can be stored"; return(); } else { if (#[size(#)-2]=="") { "ERROR: may not be empty"; return(); } string newring=#[size(#)-2]; g=size(#)-3; } } else { if (typeof(#[size(#)-1])<>"string") { "ERROR: In characteristic p>0 a must be given for the name of a new"; " ring where the Molien series can be stored"; return(); } else { if (#[size(#)-1]=="") { "ERROR: may not be empty"; return(); } string newring=#[size(#)-1]; g=size(#)-2; int r=g; } } } else // then ist not needed { g=size(#)-1; } } else // last parameter is not { int v=0; // no information is default int mol_flag=0; // computing of Molien series is default int interval=0; if (ch<>0) { if (typeof(#[size(#)])=="int") { int r=#[size(#)]; if (typeof(#[size(#)-1])<>"string") { "ERROR: in characteristic p>0 a must be given for the name of a new"; " ring where the Molien series can be stored"; return(); } else { if (#[size(#)-1]=="") { "ERROR: may not be empty"; return(); } string newring=#[size(#)-1]; g=size(#)-2; } } else { if (typeof(#[size(#)])<>"string") { "ERROR: in characteristic p>0 a must be given for the name of a new"; " ring where the Molien series can be stored"; return(); } else { if (#[size(#)]=="") { "ERROR: may not be empty"; return(); } string newring=#[size(#)]; g=size(#)-1; int r=g; } } } else { g=size(#); } } if (ch<>0) { if ((g/r)*r<>g) { "ERROR: should divide the group order." return(); } } if (ch<>0) { if ((g%ch)==0) { if (voice==2) { "WARNING: The characteristic of the coefficient field divides the group"; " order. Proceed without the Molien series!"; } else { if (v) { " The characteristic of the base field divides the group order."; " We have to continue without Molien series..."; ""; } } } if (minpoly<>0 && mol_flag==0) { if (voice==2) { "WARNING: It is impossible for this program to calculate the Molien series"; " for finite groups over extension fields of prime characteristic."; } else { if (v) { " Since it is impossible for this program to calculate the Molien series for"; " invariant rings over extension fields of prime characteristic, we have to"; " continue without it."; ""; } } return(); } } //---------------------------------------------------------------------------- if (not(typeof(#[1])=="matrix")) { "ERROR: the parameters must be a list of matrices and maybe an "; return(); } 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(); } if (v && voice<>2) { ""; " Generating the Molien series..."; ""; } if (v && voice==2) { ""; } //------------- calculating Molien series in characteristic 0 ---------------- if (ch==0) // when ch==0 we can calculate the Molien { matrix I=diag(1,n); // series in any case - poly v1=maxideal(1)[1]; // the Molien series will be in terms of // the first variable of the current // ring - matrix M[1][2]; // M will contain the Molien series - M[1,1]=0; // M[1,1] will be the numerator - M[1,2]=1; // M[1,2] will be the denominator - matrix s; // will help us canceling in the // fraction poly p; // will contain the denominator of the // new term of the Molien series //------------ computing 1/det(1+xE) for all E in the group ------------------ for (int j=1;j<=g;j++) { if (not(typeof(#[j])=="matrix")) { "ERROR: the parameters must be a list of matrices and maybe an "; return(); } if ((n<>nrows(#[j])) or (n<>ncols(#[j]))) { "ERROR: matrices need to be square and of the same dimensions"; return(); } p=det(I-v1*#[j]); // denominator of new term - M[1,1]=M[1,1]*p+M[1,2]; // expanding M[1,1]/M[1,2] + 1/p M[1,2]=M[1,2]*p; if (interval<>0) // canceling common terms of denominator { if ((j/interval)*interval==j or j==g) // and enumerator - { s=matrix(syz(ideal(M))); // once gcd() is faster than syz() these M[1,1]=-s[2,1]; // three lines should be replaced by the M[1,2]=s[1,1]; // following three // p=gcd(M[1,1],M[1,2]); // M[1,1]=M[1,1]/p; // M[1,2]=M[1,2]/p; } } if (v) { " Term ",j," of the Molien series has been computed."; } } if (interval==0) // canceling common terms of denominator { // and enumerator - s=matrix(syz(ideal(M))); // once gcd() is faster than syz() these M[1,1]=-s[2,1]; // three lines should be replaced by the M[1,2]=s[1,1]; // following three // p=gcd(M[1,1],M[1,2]); // M[1,1]=M[1,1]/p; // M[1,2]=M[1,2]/p; } map slead=br,ideal(0); s=slead(M); M[1,1]=1/s[1,1]*M[1,1]; // numerator and denominator have to have M[1,2]=1/s[1,2]*M[1,2]; // a constant term of 1 if (v) { ""; " We are done calculating the Molien series."; ""; } return(M); } //---- calculating Molien series in prime characteristic with Brauer lift ---- if (ch<>0 && mol_flag==0) { if (g<>1) { matrix G(1..g)=#[1..g]; if (interval<0) { string Mstring; } //------ preparing everything for Brauer lifts into characteristic 0 --------- ring Q=0,x,dp; // we want to extend our ring as well as // the ring of rational numbers Q to // contain r-th primitive roots of unity // in order to factor characteristic // polynomials of group elements into // linear factors and lift eigenvalues to // characteristic 0 - poly minq=cyclotomic(r); // minq now contains the size-of-group-th // cyclotomic polynomial of Q, it is // irreducible there ring `newring`=(0,e),x,dp; map f=Q,ideal(e); minpoly=number(f(minq)); // e is now a r-th primitive root of // unity - kill Q, f; // no longer needed - poly p=1; // used to build the denominator of the // new term in the Molien series matrix s[1][2]; // used for canceling - matrix M[1][2]=0,1; // will contain Molien series - ring v1br=char(br),x,dp; // we calculate the r-th cyclotomic poly minp=cyclotomic(r); // polynomial of the base field and pick minp=factorize(minp)[1][2]; // an irreducible factor of it - if (deg(minp)==1) // in this case the base field contains { ring bre=char(br),x,dp; // r-th roots of unity already map f1=v1br,ideal(0); number e=-number((f1(minp))); // e is a r-th primitive root of unity } else { ring bre=(char(br),e),x,dp; map f1=v1br,ideal(e); minpoly=number(f1(minp)); // e is a r-th primitive root of unity } map f2=br,ideal(0); // we need f2 to map our group elements // to this new extension field bre matrix xI=diag(x,n); poly p; // used for the characteristic polynomial // to factor - list L; // will contain the linear factors of the ideal F; // characteristic polynomial of the group intvec C; // elements and their powers int i, j, k; // -------------- finding all the terms of the Molien series ----------------- for (i=1;i<=g;i++) { setring br; if (not(typeof(#[i])=="matrix")) { "ERROR: the parameters must be a list of matrices and maybe an "; return(); } if ((n<>nrows(#[i])) or (n<>ncols(#[i]))) { "ERROR: matrices need to be square and of the same dimensions"; return(); } setring bre; p=det(xI-f2(G(i))); // characteristic polynomial of G(i) L=factorize(p); F=L[1]; C=L[2]; for (j=2;j<=ncols(F);j++) { F[j]=-1*(F[j]-x); // F[j] is now an eigenvalue of G(i), // it is a power of a primitive r-th root // of unity - k=exponent(number(F[j]),e); // F[j]==e^k setring `newring`; p=p*(1-x*(e^k))^C[j]; // building the denominator of the new setring bre; // term } // ----------- // k=0; // while(kg) { Mstring=string(M); for (j=1;j<=size(Mstring);j++) { if (Mstring[j]=="e") { interval=0; break; } } } if (interval<>0) { s=matrix(syz(ideal(M))); // once gcd() is faster than syz() M[1,1]=-s[2,1]; // these three lines should be M[1,2]=s[1,1]; // replaced by the following three // p=gcd(M[1,1],M[1,2]); // M[1,1]=M[1,1]/p; // M[1,2]=M[1,2]/p; } else { interval=-1; } } else { if (interval<>0) // canceling common terms of denominator { if ((i/interval)*interval==i or i==g) // and enumerator { s=matrix(syz(ideal(M))); // once gcd() is faster than syz() M[1,1]=-s[2,1]; // these three lines should be M[1,2]=s[1,1]; // replaced by the following three // p=gcd(M[1,1],M[1,2]); // M[1,1]=M[1,1]/p; // M[1,2]=M[1,2]/p; } } } p=1; setring bre; if (v) { " Term ",i," of the Molien series has been computed."; } } if (v) { ""; } setring `newring`; if (interval==0) // canceling common terms of denominator { // and enumerator - s=matrix(syz(ideal(M))); // once gcd() is faster than syz() these M[1,1]=-s[2,1]; // three lines should be replaced by the M[1,2]=s[1,1]; // following three // p=gcd(M[1,1],M[1,2]); // M[1,1]=M[1,1]/p; // M[1,2]=M[1,2]/p; } map slead=`newring`,ideal(0); s=slead(M); // forcing the constant term of numerator M[1,1]=1/s[1,1]*M[1,1]; // and denominator to be 1 M[1,2]=1/s[1,2]*M[1,2]; kill slead; kill s; kill p; } else // if the group only contains an identity { ring `newring`=0,x,dp; // element, it is very easy to calculate matrix M[1][2]=1,(1-x)^n; // the Molien series } exportto(Top,`newring`); // we keep the ring where we computed the export M; // Molien series in such that we can setring br; // keep it if (v) { " We are done calculating the Molien series."; ""; } } else // i.e. char<>0 and mol_flag<>0, the user { // has specified that we are dealing with // a ring of large characteristic which // can be treated like a ring of // characteristic 0; we'll avoid the // Brauer lifts //----------------------- simulating characteristic 0 ------------------------ string chst=charstr(br); for (int i=1;i<=size(chst);i++) { if (chst[i]==",") { break; } } //----------------- generating ring of characteristic 0 ---------------------- if (minpoly==0) { if (i>size(chst)) { execute("ring "+newring+"=0,("+varstr(br)+"),("+ordstr(br)+")"); } else { chst=chst[i..size(chst)]; execute ("ring "+newring+"=(0"+chst+"),("+varstr(br)+"),("+ordstr(br)+")"); } } else { string minp=string(minpoly); minp=minp[2..size(minp)-1]; chst=chst[i..size(chst)]; execute("ring "+newring+"=(0"+chst+"),("+varstr(br)+"),("+ordstr(br)+")"); execute("minpoly="+minp); } matrix I=diag(1,n); poly v1=maxideal(1)[1]; // the Molien series will be in terms of // the first variable of the current // ring - matrix M[1][2]; // M will contain the Molien series - M[1,1]=0; // M[1,1] will be the numerator - M[1,2]=1; // M[1,2] will be the denominator - matrix s; // will help us canceling in the // fraction poly p; // will contain the denominator of the // new term of the Molien series int j; string links, rechts; //----------------- finding all terms of the Molien series ------------------- for (i=1;i<=g;i++) { setring br; if (not(typeof(#[i])=="matrix")) { "ERROR: the parameters must be a list of matrices and maybe an "; return(); } if ((n<>nrows(#[i])) or (n<>ncols(#[i]))) { "ERROR: matrices need to be square and of the same dimensions"; return(); } string stM(i)=string(#[i]); for (j=1;j<=size(stM(i));j++) { if (stM(i)[j]==" ") { links=stM(i)[1..j-1]; rechts=stM(i)[j+1..size(stM(i))]; stM(i)=links+rechts; } } setring `newring`; execute("matrix G(i)["+string(n)+"]["+string(n)+"]="+stM(i)); p=det(I-v1*G(i)); // denominator of new term - M[1,1]=M[1,1]*p+M[1,2]; // expanding M[1,1]/M[1,2] + 1/p M[1,2]=M[1,2]*p; if (interval<>0) // canceling common terms of denominator { if ((i/interval)*interval==i or i==g) // and enumerator { s=matrix(syz(ideal(M))); // once gcd() is faster than syz() these M[1,1]=-s[2,1]; // three lines should be replaced by the M[1,2]=s[1,1]; // following three // p=gcd(M[1,1],M[1,2]); // M[1,1]=M[1,1]/p; // M[1,2]=M[1,2]/p; } } if (v) { " Term ",i," of the Molien series has been computed."; } } if (interval==0) // canceling common terms of denominator { // and enumerator - s=matrix(syz(ideal(M))); // once gcd() is faster than syz() these M[1,1]=-s[2,1]; // three lines should be replaced by the M[1,2]=s[1,1]; // following three // p=gcd(M[1,1],M[1,2]); // M[1,1]=M[1,1]/p; // M[1,2]=M[1,2]/p; } map slead=`newring`,ideal(0); s=slead(M); M[1,1]=1/s[1,1]*M[1,1]; // numerator and denominator have to have M[1,2]=1/s[1,2]*M[1,2]; // a constant term of 1 if (v) { ""; " We are done calculating the Molien series."; ""; } kill G(1..g), s, slead, p, v1, I; export `newring`; // we keep the ring where we computed the export M; // the Molien series such that we can setring br; // keep it } } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; " note the case of prime characteristic"; echo=2; ring R=0,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; list L=group_reynolds(A); matrix M=molien(L[2..size(L)]); print(M); ring S=3,(x,y,z),dp; string newring="alksdfjlaskdjf"; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; list L=group_reynolds(A); molien(L[2..size(L)],newring); setring alksdfjlaskdjf; print(M); setring S; kill alksdfjlaskdjf; } /////////////////////////////////////////////////////////////////////////////// proc reynolds_molien (list #) "USAGE: reynolds_molien(G1,G2,...[,ringname,flags]); G1,G2,...: nxn generating a finite matrix group, ringname: a giving a name for a new ring of characteristic 0 for the Molien series in case of prime characteristic, flags: an optional with three components: if the first element is not equal to 0 characteristic 0 is simulated, i.e. the Molien series is computed as if the base field were characteristic 0 (the user must choose a field of large prime characteristic, e.g. 32003) the second component should give the size of intervals between canceling common factors in the expansion of the Molien series, 0 (the default) means only once after generating all terms, in prime characteristic also a negative number can be given to indicate that common factors should always be canceled when the expansion is simple (the root of the extension field does not occur among the coefficients) ASSUME: n is the number of variables of the basering, G1,G2... are the group elements generated by group_reynolds(), g is the size of the group RETURN: a gxn representing the Reynolds operator is the first return value and in case of characteristic 0 a 1x2 giving enumerator and denominator of Molien series is the second one; in case of prime characteristic a ring with the name `ringname` of characteristic 0 is created where the same Molien series (named M) is stored DISPLAY: information if the third component of flags does not equal 0 THEORY: The entire matrix group is generated by getting all left products of the generators with new elements from the last run through the loop (or the generators themselves during the first run). All the ones that have been generated before are thrown out and the program terminates when are no new elements found in one run. Additionally each time a new group element is found the corresponding ring mapping of which the Reynolds operator is made up is generated. They are stored in the rows of the first return value. In characteristic 0 the terms 1/det(1-xE) is computed whenever a new element E is found. In prime characteristic a Brauer lift is involved and the terms are only computed after the entire matrix group is generated (to avoid the modular case). The returned matrix gives enumerator and denominator of the expanded version where common factors have been canceled. EXAMPLE: example reynolds_molien; shows an example " { def br=basering; // the Molien series depends on the int ch=char(br); // characteristic of the coefficient // field int gen_num; //------------------- making sure the input is okay -------------------------- if (typeof(#[size(#)])=="intvec") { if (size(#[size(#)])==3) { int mol_flag=#[size(#)][1]; if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0))) { "ERROR: the second component of the should be >=0"; return(); } int interval=#[size(#)][2]; int v=#[size(#)][3]; } else { "ERROR: should have three components"; return(); } if (ch<>0) { if (typeof(#[size(#)-1])<>"string") { "ERROR: in characteristic p a must be given for the name"; " of a new ring where the Molien series can be stored"; return(); } else { if (#[size(#)-1]=="") { "ERROR: may not be empty"; return(); } string newring=#[size(#)-1]; gen_num=size(#)-2; } } else // then ist not needed { gen_num=size(#)-1; } } else // last parameter is not { int v=0; // no information is default int interval; int mol_flag=0; // computing of Molien series is default if (ch<>0) { if (typeof(#[size(#)])<>"string") { "ERROR: in characteristic p a must be given for the name"; " of a new ring where the Molien series can be stored"; return(); } else { if (#[size(#)]=="") { "ERROR: may not be empty"; return(); } string newring=#[size(#)]; gen_num=size(#)-1; } } else { gen_num=size(#); } } // ----------------- computing the terms with Brauer lift -------------------- if (ch<>0 && mol_flag==0) { list L=group_reynolds(#[1..gen_num],v); if (L[1]==0) { if (voice==2) { "WARNING: The characteristic of the coefficient field divides the group order."; " Proceed without the Reynolds operator or the Molien series!"; return(); } if (v) { " The characteristic of the base field divides the group order."; " We have to continue without Reynolds operator or the Molien series..."; return(); } } if (minpoly<>0) { if (voice==2) { "WARNING: It is impossible for this program to calculate the Molien series"; " for finite groups over extension fields of prime characteristic."; return(L[1]); } else { if (v) { " Since it is impossible for this program to calculate the Molien series for"; " invariant rings over extension fields of prime characteristic, we have to"; " continue without it."; return(L[1]); } } } if (typeof(L[2])=="int") { molien(L[3..size(L)],newring,L[2],intvec(mol_flag,interval,v)); } else { molien(L[2..size(L)],newring,intvec(mol_flag,interval,v)); } return(L[1]); } //----------- computing Molien series in the straight forward way ------------ if (ch==0) { if (typeof(#[1])<>"matrix") { "ERROR: the parameters must be a list of matrices and maybe an "; return(); } 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(); } if (n<>ncols(#[1])) { "ERROR: matrices need to be square and of the same dimensions"; return(); } matrix vars=matrix(maxideal(1)); // creating an nx1-matrix containing the vars=transpose(vars); // variables of the ring - matrix A(1)=#[1]*vars; // calculating the first ring mapping - // A(1) will contain the Reynolds // operator - poly v1=vars[1,1]; // the Molien series will be in terms of // the first variable of the current // ring matrix I=diag(1,n); matrix A(2)[1][2]; // A(2) will contain the Molien series - A(2)[1,1]=1; // A(2)[1,1] will be the numerator matrix G(1)=#[1]; // G(k) are elements of the group - A(2)[1,2]=det(I-v1*(G(1))); // A(2)[1,2] will be the denominator - matrix s; // will help us canceling in the // fraction poly p; // will contain the denominator of the // new term of the Molien series int i=1; for (int j=2;j<=gen_num;j++) // this loop adds the parameters to the { // group, leaving out doubles and // checking whether the parameters are // compatible with the task of the // procedure if (not(typeof(#[j])=="matrix")) { "ERROR: the parameters must be a list of matrices and maybe an "; return(); } if ((n!=nrows(#[j])) or (n!=ncols(#[j]))) { "ERROR: matrices need to be square and of the same dimensions"; return(); } if (unique(G(1..i),#[j])) { i++; matrix G(i)=#[j]; A(1)=concat(A(1),#[j]*vars); // adding ring homomorphisms to A(1) - p=det(I-v1*#[j]); // denominator of new term - A(2)[1,1]=A(2)[1,1]*p+A(2)[1,2]; // expanding A(2)[1,1]/A(2)[1,2] +1/p A(2)[1,2]=A(2)[1,2]*p; if (interval<>0) // canceling common terms of denominator { if ((i/interval)*interval==i) // and enumerator { s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz() these A(2)[1,1]=-s[2,1]; // three lines should be replaced by the A(2)[1,2]=s[1,1]; // following three // p=gcd(A(2)[1,1],A(2)[1,2]); // A(2)[1,1]=A(2)[1,1]/p; // A(2)[1,2]=A(2)[1,2]/p; } } } } int g=i; // G(1)..G(i) are generators without // doubles - g generally is the number // of elements in the group so far - j=i; // j is the number of new elements that // we use as factors int k, m, l; if (v) { ""; " Generating the entire matrix group. Whenever a new group element is found,"; " the corresponding ring homomorphism of the Reynolds operator and the"; " corresponding term of the Molien series is generated."; ""; } //----------- computing 1/det(I-xE) whenever a new element E is found -------- while (1) { l=0; // l is the number of products we get in // one going for (m=g-j+1;m<=g;m=m+1) { for (k=1;k<=i;k++) { l++; matrix P(l)=G(k)*G(m); // possible new element } } j=0; for (k=1;k<=l;k++) { if (unique(G(1..g),P(k))) { j++; // a new factor for next run g++; matrix G(g)=P(k); // a new group element - A(1)=concat(A(1),P(k)*vars); // adding new mapping to A(1) p=det(I-v1*P(k)); // denominator of new term A(2)[1,1]=A(2)[1,1]*p+A(2)[1,2]; A(2)[1,2]=A(2)[1,2]*p; // expanding A(2)[1,1]/A(2)[1,2] + 1/p - if (interval<>0) // canceling common terms of denominator { if ((g/interval)*interval==g) // and enumerator { s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz() A(2)[1,1]=-s[2,1]; // these three lines should be replaced A(2)[1,2]=s[1,1]; // by the following three // p=gcd(A(2)[1,1],A(2)[1,2]); // A(2)[1,1]=A(2)[1,1]/p; // A(2)[1,2]=A(2)[1,2]/p; } } if (v) { " Group element ",g," has been found."; } } kill P(k); } if (j==0) // when we didn't add any new elements { break; // in one run through the while loop } // we are done } if (v) { if (g<=i) { " There are only ",g," group elements."; } ""; } A(1)=transpose(A(1)); // when we evaluate the Reynolds operator // later on, we actually want 1xn // matrices if (interval==0) // canceling common terms of denominator { // and enumerator - s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz() A(2)[1,1]=-s[2,1]; // these three lines should be replaced A(2)[1,2]=s[1,1]; // by the following three // p=gcd(A(2)[1,1],A(2)[1,2]); // A(2)[1,1]=A(2)[1,1]/p; // A(2)[1,2]=A(2)[1,2]/p; } if (interval<>0) // canceling common terms of denominator { if ((g/interval)*interval<>g) // and enumerator { s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz() A(2)[1,1]=-s[2,1]; // these three lines should be replaced A(2)[1,2]=s[1,1]; // by the following three // p=gcd(A(2)[1,1],A(2)[1,2]); // A(2)[1,1]=A(2)[1,1]/p; // A(2)[1,2]=A(2)[1,2]/p; } } map slead=br,ideal(0); s=slead(A(2)); A(2)[1,1]=1/s[1,1]*A(2)[1,1]; // numerator and denominator have to have A(2)[1,2]=1/s[1,2]*A(2)[1,2]; // a constant term of 1 if (v) { " Now we are done calculating Molien series and Reynolds operator."; ""; } return(A(1..2)); } //------------------------ simulating characteristic 0 ----------------------- else // if ch<>0 and mol_flag<>0 { if (typeof(#[1])<>"matrix") { "ERROR: the parameters must be a list of matrices and maybe an "; return(); } 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(); } if (n<>ncols(#[1])) { "ERROR: matrices need to be square and of the same dimensions"; return(); } matrix vars=matrix(maxideal(1)); // creating an nx1-matrix containing the vars=transpose(vars); // variables of the ring - matrix A(1)=#[1]*vars; // calculating the first ring mapping - // A(1) will contain the Reynolds // operator string chst=charstr(br); for (int i=1;i<=size(chst);i++) { if (chst[i]==",") { break; } } if (minpoly==0) { if (i>size(chst)) { execute("ring "+newring+"=0,("+varstr(br)+"),("+ordstr(br)+")"); } else { chst=chst[i..size(chst)]; execute ("ring "+newring+"=(0"+chst+"),("+varstr(br)+"),("+ordstr(br)+")"); } } else { string minp=string(minpoly); minp=minp[2..size(minp)-1]; chst=chst[i..size(chst)]; execute("ring "+newring+"=(0"+chst+"),("+varstr(br)+"),("+ordstr(br)+")"); execute("minpoly="+minp); } poly v1=var(1); // the Molien series will be in terms of // the first variable of the current // ring matrix I=diag(1,n); int o; setring br; matrix G(1)=#[1]; string links, rechts; string stM(1)=string(#[1]); for (o=1;o<=size(stM(1));o++) { if (stM(1)[o]==" ") { links=stM(1)[1..o-1]; rechts=stM(1)[o+1..size(stM(1))]; stM(1)=links+rechts; } } setring `newring`; execute("matrix G(1)["+string(n)+"]["+string(n)+"]="+stM(1)); matrix A(2)[1][2]; // A(2) will contain the Molien series - A(2)[1,1]=1; // A(2)[1,1] will be the numerator A(2)[1,2]=det(I-v1*(G(1))); // A(2)[1,2] will be the denominator - matrix s; // will help us canceling in the // fraction poly p; // will contain the denominator of the // new term of the Molien series i=1; for (int j=2;j<=gen_num;j++) // this loop adds the parameters to the { // group, leaving out doubles and // checking whether the parameters are // compatible with the task of the // procedure setring br; if (not(typeof(#[j])=="matrix")) { "ERROR: the parameters must be a list of matrices and maybe an "; return(); } if ((n!=nrows(#[j])) or (n!=ncols(#[j]))) { "ERROR: matrices need to be square and of the same dimensions"; return(); } if (unique(G(1..i),#[j])) { i++; matrix G(i)=#[j]; A(1)=concat(A(1),G(i)*vars); // adding ring homomorphisms to A(1) string stM(i)=string(G(i)); for (o=1;o<=size(stM(i));o++) { if (stM(i)[o]==" ") { links=stM(i)[1..o-1]; rechts=stM(i)[o+1..size(stM(i))]; stM(i)=links+rechts; } } setring `newring`; execute("matrix G(i)["+string(n)+"]["+string(n)+"]="+stM(i)); p=det(I-v1*G(i)); // denominator of new term - A(2)[1,1]=A(2)[1,1]*p+A(2)[1,2]; // expanding A(2)[1,1]/A(2)[1,2] +1/p A(2)[1,2]=A(2)[1,2]*p; if (interval<>0) // canceling common terms of denominator { if ((i/interval)*interval==i) // and enumerator { s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz() these A(2)[1,1]=-s[2,1]; // three lines should be replaced by the A(2)[1,2]=s[1,1]; // following three // p=gcd(A(2)[1,1],A(2)[1,2]); // A(2)[1,1]=A(2)[1,1]/p; // A(2)[1,2]=A(2)[1,2]/p; } } setring br; } } int g=i; // G(1)..G(i) are generators without // doubles - g generally is the number // of elements in the group so far - j=i; // j is the number of new elements that // we use as factors int k, m, l; if (v) { ""; " Generating the entire matrix group. Whenever a new group element is found,"; " the corresponding ring homomorphism of the Reynolds operator and the"; " corresponding term of the Molien series is generated."; ""; } // taking all elements in a ring of characteristic 0 and computing the terms // of the Molien series there while (1) { l=0; // l is the number of products we get in // one going for (m=g-j+1;m<=g;m++) { for (k=1;k<=i;k++) { l++; matrix P(l)=G(k)*G(m); // possible new element } } j=0; for (k=1;k<=l;k++) { if (unique(G(1..g),P(k))) { j++; // a new factor for next run g++; matrix G(g)=P(k); // a new group element - A(1)=concat(A(1),G(g)*vars); // adding new mapping to A(1) string stM(g)=string(G(g)); for (o=1;o<=size(stM(g));o++) { if (stM(g)[o]==" ") { links=stM(g)[1..o-1]; rechts=stM(g)[o+1..size(stM(g))]; stM(g)=links+rechts; } } setring `newring`; execute("matrix G(g)["+string(n)+"]["+string(n)+"]="+stM(g)); p=det(I-v1*G(g)); // denominator of new term A(2)[1,1]=A(2)[1,1]*p+A(2)[1,2]; A(2)[1,2]=A(2)[1,2]*p; // expanding A(2)[1,1]/A(2)[1,2] + 1/p - if (interval<>0) // canceling common terms of denominator { if ((g/interval)*interval==g) // and enumerator { s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz() A(2)[1,1]=-s[2,1]; // these three lines should be replaced A(2)[1,2]=s[1,1]; // by the following three // p=gcd(A(2)[1,1],A(2)[1,2]); // A(2)[1,1]=A(2)[1,1]/p; // A(2)[1,2]=A(2)[1,2]/p; } } if (v) { " Group element ",g," has been found."; } setring br; } kill P(k); } if (j==0) // when we didn't add any new elements { break; // in one run through the while loop } // we are done } if (v) { if (g<=i) { " There are only ",g," group elements."; } ""; } A(1)=transpose(A(1)); // when we evaluate the Reynolds operator // later on, we actually want 1xn // matrices setring `newring`; if (interval==0) // canceling common terms of denominator { // and enumerator - s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz() A(2)[1,1]=-s[2,1]; // these three lines should be replaced A(2)[1,2]=s[1,1]; // by the following three // p=gcd(A(2)[1,1],A(2)[1,2]); // A(2)[1,1]=A(2)[1,1]/p; // A(2)[1,2]=A(2)[1,2]/p; } if (interval<>0) // canceling common terms of denominator { if ((g/interval)*interval<>g) // and enumerator { s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz() A(2)[1,1]=-s[2,1]; // these three lines should be replaced A(2)[1,2]=s[1,1]; // by the following three // p=gcd(A(2)[1,1],A(2)[1,2]); // A(2)[1,1]=A(2)[1,1]/p; // A(2)[1,2]=A(2)[1,2]/p; } } map slead=`newring`,ideal(0); s=slead(A(2)); A(2)[1,1]=1/s[1,1]*A(2)[1,1]; // numerator and denominator have to have A(2)[1,2]=1/s[1,2]*A(2)[1,2]; // a constant term of 1 if (v) { " Now we are done calculating Molien series and Reynolds operator."; ""; } matrix M=A(2); kill G(1..g), s, slead, p, v1, I, A(2); export `newring`; // we keep the ring where we computed the export M; // the Molien series such that we can setring br; // keep it return(A(1)); } } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; " note the case of prime characteristic"; echo=2; ring R=0,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; matrix REY,M=reynolds_molien(A); print(REY); print(M); ring S=3,(x,y,z),dp; string newring="Qadjoint"; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; matrix REY=reynolds_molien(A,newring); print(REY); setring Qadjoint; print(M); setring S; kill Qadjoint; } /////////////////////////////////////////////////////////////////////////////// proc partial_molien (matrix M, int n, list #) "USAGE: partial_molien(M,n[,p]); M: a 1x2 , n: an indicating number of terms in the expansion, p: an optional ASSUME: M is the return value of molien or the second return value of reynolds_molien, p ought to be the second return value of a previous run of partial_molien and avoids recalculating known terms RETURN: n terms (type ) of the partial expansion of the Molien series (first n if there is no third parameter given, otherwise the next n terms depending on a previous calculation) and an intermediate result (type ) of the calculation to be used as third parameter in a next run of partial_molien THEORY: The following calculation is implemented: @format (1+a1x+a2x^2+...+anx^n)/(1+b1x+b2x^2+...+bmx^m)=(1+(a1-b1)x+... (1+b1x+b2x^2+...+bmx^m) ----------------------- (a1-b1)x+(a2-b2)x^2+... (a1-b1)x+b1(a1-b1)x^2+... @end format EXAMPLE: example partial_molien; shows an example " { poly A(2); // A(2) will contain the return value of // the intermediate result if (char(basering)<>0) { "ERROR: you have to change to a basering of characteristic 0, one in"; " which the Molien series is defined"; } if (ncols(M)==2 && nrows(M)==1 && n>0 && size(#)<2) { def br=basering; // keeping track of the old ring map slead=br,ideal(0); matrix s=slead(M); if (s[1,1]<>1 || s[1,2]<>1) { "ERROR: the constant terms of enumerator and denominator are not 1"; return(); } if (size(#)==0) { A(2)=M[1,1]; // if a third parameter is not given, the // intermediate result from the last run // corresponds to the numerator - we need } // its smallest term else { if (typeof(#[1])=="poly") { A(2)=#[1]; // if a third term is given we 'start' } // with its smallest term else { "ERROR: as third parameter expected"; return(); } } poly A(1)=M[1,2]; // denominator of Molien series (for now) string mp=string(minpoly); execute("ring R=("+charstr(br)+"),("+varstr(br)+"),ds;"); execute("minpoly=number("+mp+");"); poly A(1)=0; // A(1) will contain the sum of n terms - poly min; // min will be our smallest term - poly A(2)=fetch(br,A(2)); // fetching A(2) and M[1,2] into R poly den=fetch(br,A(1)); for (int i=1; i<=n; i++) // getting n terms and adding them up { min=lead(A(2)); A(1)=A(1)+min; A(2)=A(2)-min*den; } setring br; // moving A(1) and A(2) back in the A(1)=fetch(R,A(1)); // actual ring for output A(2)=fetch(R,A(2)); return(A(1..2)); } else { "ERROR: the first parameter has to be a 1x2-matrix, i.e. the matrix"; " returned by the procedure 'reynolds_molien', the second one"; " should be > 0 and there should be no more than 3 parameters;" return(); } } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2; ring R=0,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; matrix REY,M=reynolds_molien(A); poly p(1..2); p(1..2)=partial_molien(M,5); p(1); p(1..2)=partial_molien(M,5,p(2)); p(1); } /////////////////////////////////////////////////////////////////////////////// proc evaluate_reynolds (matrix REY, ideal I) "USAGE: evaluate_reynolds(REY,I); REY: a representing the Reynolds operator, I: an arbitrary ASSUME: REY is the first return value of group_reynolds() or reynolds_molien() RETURNS: image of the polynomials defining I under the Reynolds operator (type ) NOTE: the characteristic of the coefficient field of the polynomial ring should not divide the order of the finite matrix group THEORY: REY has been constructed in such a way that each row serves as a ring mapping of which the Reynolds operator is made up. EXAMPLE: example evaluate_reynolds; shows an example " { def br=basering; int n=nvars(br); if (ncols(REY)==n) { int m=nrows(REY); // we need m to 'cut' the ring // homomorphisms 'out' of REY and to // divide by the group order in the end int num_poly=ncols(I); matrix MI=matrix(I); matrix MiI[1][num_poly]; map pREY; matrix rowREY[1][n]; for (int i=1;i<=m;i++) { rowREY=REY[i,1..n]; pREY=br,ideal(rowREY); // f is now the i-th ring homomorphism MiI=pREY(MI)+MiI; } MiI=(1/number(m))*MiI; return(ideal(MiI)); } else { "ERROR: the number of columns in the should be the same as the"; " number of variables in the basering; in fact it should be first"; " return value of group_reynolds() or reynolds_molien()."; return(); } } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2; ring R=0,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; list L=group_reynolds(A); ideal I=x2,y2,z2; print(evaluate_reynolds(L[1],I)); } /////////////////////////////////////////////////////////////////////////////// proc invariant_basis (int g,list #) "USAGE: invariant_basis(g,G1,G2,...); g: an indicating of which degree (>0) the homogeneous basis shoud be, G1,G2,...: generating a finite matrix group RETURNS: the basis (type ) of the space of 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. EXAMPLE: example invariant_basis; shows an example " { if (g<=0) { "ERROR: the first parameter should be > 0"; return(); } def br=basering; ideal vars = maxideal(1); ideal mon=sort(maxideal(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); int DGB = degBound; //---------------------- 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(); } } //---------------------------------------------------------------------------- ring tmpR = char(br),(p(1..m)),dp; def T = br+tmpR; setring T; // 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 ideal vars = fetch(br,vars); 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); degBound=1; 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 } } degBound=DGB; return(I); } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2; ring R=0,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; print(invariant_basis(2,A)); } /////////////////////////////////////////////////////////////////////////////// proc invariant_basis_reynolds (matrix REY,int d,list #) "USAGE: invariant_basis_reynolds(REY,d[,flags]); REY: a representing the Reynolds operator, d: an indicating of which degree (>0) the homogeneous basis shoud be, flags: an optional with two entries: its first component gives the dimension of the space (default <0 meaning unknown) and its second component is used as the number of polynomials that should be mapped to invariants during one call of evaluate_reynolds if the dimension of the space is unknown or the number such that number x dimension polynomials are mapped to invariants during one call of evaluate_reynolds ASSUME: REY is the first return value of group_reynolds() or reynolds_molien() and flags[1] given by partial_molien RETURN: the basis (type ) of the space of invariants of degree d THEORY: Monomials of degree d are mapped to invariants with the Reynolds operator. A linearly independent set is generated with the help of minbase. EXAMPLE: example invariant_basis_reynolds; shows an example " { //---------------------- checking that the input is ok ----------------------- if (d<=0) { " ERROR: the second parameter should be > 0"; return(); } if (size(#)>1) { " ERROR: there should be at most three parameters"; return(); } if (size(#)==1) { if (typeof(#[1])<>"intvec") { " ERROR: the third parameter should be of type "; return(); } if (size(#[1])<>2) { " ERROR: there should be two components in "; return(); } else { int cd=#[1][1]; int step_fac=#[1][2]; } if (step_fac<=0) { " ERROR: the second component of should be > 0"; return(); } if (cd==0) { return(ideal(0)); } } else { int step_fac=1; int cd=-1; } if (ncols(REY)<>nvars(basering)) { "ERROR: the number of columns in the should be the same as the"; " number of variables in the basering; in fact it should be first"; " return value of group_reynolds() or reynolds_molien()."; return(); } //---------------------------------------------------------------------------- ideal mon=sort(maxideal(d))[1]; int DEGB = degBound; degBound=d; int j=ncols(mon); mon=sort(mon,intvec(j..1))[1]; ideal B; // will contain the basis if (cd<0) { if (step_fac>j) // all of mon will be mapped to { B=evaluate_reynolds(REY,mon); // invariants at once B=minbase(B); degBound=DEGB; return(B); } } else { if (step_fac*cd>j) // all of mon will be mapped to { B=evaluate_reynolds(REY,mon); // invariants at once B=minbase(B); degBound=DEGB; return(B); } } int i,k; int upper_bound=0; int lower_bound=0; ideal part_mon; // a part of mon of size step_fac*cd while (1) { lower_bound=upper_bound+1; if (cd<0) { upper_bound=upper_bound+step_fac; } else { upper_bound=upper_bound+step_fac*cd; } if (upper_bound>j) { upper_bound=j; } part_mon=mon[lower_bound..upper_bound]; B=minbase(B+evaluate_reynolds(REY,part_mon)); if ((ncols(B)==cd and B[1]<>0) or upper_bound==j) { degBound=DEGB; return(B); } } } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2; ring R=0,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; intvec flags=0,1,0; matrix REY,M=reynolds_molien(A,flags); flags=8,6; print(invariant_basis_reynolds(REY,6,flags)); } /////////////////////////////////////////////////////////////////////////////// // This procedure generates linearly independent invariant polynomials of // degree d that do not reduce to 0 modulo the primary invariants. It does this // by applying the Reynolds operator to the monomials returned by kbase(sP,d). // The result is used when computing secondary invariants. /////////////////////////////////////////////////////////////////////////////// proc sort_of_invariant_basis (ideal sP,matrix REY,int d,int step_fac) { ideal mon=kbase(sP,d); int DEGB=degBound; degBound=d; int j=ncols(mon); int i; mon=sort(mon,intvec(j..1))[1]; ideal B; // will contain the "sort of basis" if (step_fac>j) { B=compress(evaluate_reynolds(REY,mon)); for (i=1;i<=ncols(B);i++) // those are taken our that are o mod sP { if (reduce(B[i],sP)==0) { B[i]=0; } } B=minbase(B); // here are the linearly independent ones degBound=DEGB; return(B); } int upper_bound=0; int lower_bound=0; ideal part_mon; // parts of mon while (1) { lower_bound=upper_bound+1; upper_bound=upper_bound+step_fac; if (upper_bound>j) { upper_bound=j; } part_mon=mon[lower_bound..upper_bound]; part_mon=compress(evaluate_reynolds(REY,part_mon)); for (i=1;i<=ncols(part_mon);i++) { if (reduce(part_mon[i],sP)==0) { part_mon[i]=0; } } B=minbase(B+part_mon); // here are the linearly independent ones if (upper_bound==j) { degBound=DEGB; return(B); } } } /////////////////////////////////////////////////////////////////////////////// // Procedure returning the succeeding vector after vec. It is used to list // all the vectors of Z^n with first nonzero entry 1. They are listed by // increasing sum of the absolute value of their entries. /////////////////////////////////////////////////////////////////////////////// proc next_vector(intmat vec) { int n=ncols(vec); // p: >0, n: <0, p0: >=0, n0: <=0 for (int i=1;i<=n;i++) // finding out which is the first { if (vec[1,i]<>0) // component <>0 { break; } } intmat new[1][n]; if (i>n) // 0,...,0 --> 1,0....,0 { new[1,1]=1; return(new); } if (i==n) // 0,...,1 --> 1,1,0,...,0 { new[1,1..2]=1,1; return(new); } if (i==n-1) { if (vec[1,n]==0) // 0,...,0,1,0 --> 0,...,0,1 { new[1,n]=1; return(new); } if (vec[1,n]>0) // 0,..,0,1,p --> 0,...,0,1,-p { new[1,1..n]=vec[1,1..n-1],-vec[1,n]; return(new); } new[1,1..2]=1,1-vec[1,n]; // 0,..,0,1,n --> 1,1-n,0,..,0 return(new); } if (i>1) { intmat temp[1][n-i+1]=vec[1,i..n]; // 0,...,0,1,*,...,* --> 1,*,...,* temp=next_vector(temp); new[1,i..n]=temp[1,1..n-i+1]; return(new); } // case left: 1,*,...,* for (i=2;i<=n;i++) { if (vec[1,i]>0) // make first positive negative and { vec[1,i]=-vec[1,i]; // return return(vec); } else { vec[1,i]=-vec[1,i]; // make all negatives before positives } // positive } for (i=2;i<=n-1;i++) // case: 1,p,...,p after 1,n,...,n { if (vec[1,i]>0) { vec[1,2]=vec[1,i]-1; // shuffleing things around... if (i>2) // same sum of absolute values of entries { vec[1,i]=0; } vec[1,i+1]=vec[1,i+1]+1; return(vec); } } // case left: 1,0,...,0 --> 1,1,0,...,0 new[1,2..3]=1,vec[1,n]; // and: 1,0,...,0,1 --> 0,1,1,0,...,0 return(new); } /////////////////////////////////////////////////////////////////////////////// // Maps integers to elements of the base field. It is only called if the base // field is of prime characteristic. If the base field has q elements // (depending on minpoly) 1..q is mapped to those q elements. /////////////////////////////////////////////////////////////////////////////// proc int_number_map (int i) { int p=char(basering); if (minpoly==0) // if no minpoly is given, we have p { i=i%p; // elements in the field return(number(i)); } int d=pardeg(minpoly); if (i<0) { int bool=1; i=(-1)*i; } i=i%p^d; // base field has p^d elements - number a=par(1); // a is the root of the minpoly - we have number out=0; // to construct a linear combination of int j=1; // a^k int k; while (1) { if (i0. /////////////////////////////////////////////////////////////////////////////// proc p_search (int n,int d,ideal B,int cd,ideal P,ideal sP,int i,int dif,int dB,ideal CI) { def br=basering; degBound=0; matrix vec(1)[1][cd]; // starting with 0-vector - intmat new[1][cd]; // the coefficients for the next // combination - matrix pnew[1][cd]; // new needs to be mapped into br - int counter=1; // counts the vectors int j; int p=char(br); if (minpoly<>0) { int ext_deg=pardeg(minpoly); // field has p^d elements } else { int ext_deg=1; // field has p^d elements } poly test_poly; // the linear combination to test int test_dim; ring R=0,x,dp; // just to calculate next variable // bound - number bound=(number(p)^(ext_deg*cd)-1)/(number(p)^ext_deg-1)+1; // this is how many linearly independent // vectors of size cd exist having // entries in the base field of br setring br; intvec h; // Hilbert series int k=i+1; if (ncols(B)bound) // otherwise, we are done { setring br; new=next_vector(new); for (j=1;j<=cd;j++) { pnew[1,j]=int_number_map(new[1,j]); // mapping an integer into br } if (unique(vec(1..counter),pnew)) //checking whether we tried pnew before { counter++; matrix vec(counter)=pnew; // keeping track of the ones we tried - test_poly=(vec(counter)*tB)[1,1]; // linear combination - // degBound=dB; TEST=sP+ideal(test_poly); attrib(TEST,"isSB",1); test_dim=dim(TEST); // degBound=0; if (n-test_dim==k) // the dimension has been lowered by one { sP=TEST; setring R; break; } // degBound=dB; //TEST=std(sP+ideal(test_poly)); // should soon to be replaced by next // line TEST=std(sP,test_poly); // or, better: // TEST=std(sP,test_poly,h); // Hilbert driven std-calculation test_dim=dim(TEST); // degBound=0; if (n-test_dim==k) // the dimension has been lowered by one { sP=TEST; setring R; break; } } setring R; } if (number(counter)<=bound) { setring br; P[k]=test_poly; // test_poly ist added to primary } // invariants else { setring br; CI=CI[1..size(CI)-1]; return(P,sP,CI,dB-d+1); } k++; } return(P,sP,CI,dB); } /////////////////////////////////////////////////////////////////////////////// proc primary_char0 (matrix REY,matrix M,list #) "USAGE: primary_char0(REY,M[,v]); REY: a representing the Reynolds operator, M: a 1x2 representing the Molien series, v: an optional 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 programme if v does not equal 0 RETURN: primary invariants (type ) of the invariant ring THEORY: Bases of homogeneous invariants are generated successively and those are chosen as primary invariants that lower the dimension of the ideal generated by the previously found invariants (see paper \"Generating a Noetherian Normalization of the Invariant Ring of a Finite Group\" by Decker, Heydtmann, Schreyer (1998)). EXAMPLE: example primary_char0; shows an example " { degBound=0; if (char(basering)<>0) { "ERROR: primary_char0 should only be used with rings of characteristic 0."; return(); } //----------------- checking input and setting verbose mode ------------------ if (size(#)>1) { "ERROR: primary_char0 can only have three parameters."; return(); } if (size(#)==1) { if (typeof(#[1])<>"int") { "ERROR: The third parameter should be of type ."; return(); } else { int v=#[1]; } } else { 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(REY)<>n) { "ERROR: First parameter ought to be the Reynolds operator." return(); } if (ncols(M)<>2 or nrows(M)<>1) { "ERROR: Second parameter ought to be the Molien series." return(); } //---------------------------------------------------------------------------- if (v && voice<>2) { " We can start looking for primary invariants..."; ""; } if (v && voice==2) { ""; } //------------------------- initializing variables --------------------------- int dB; poly p(1..2); // p(1) will be used for single terms of // the partial expansion, p(2) to store p(1..2)=partial_molien(M,1); // the intermediate result - poly v1=var(1); // we need v1 to split off coefficients // in the partial expansion of M (which // is in terms of the first variable) - int j,d,cd,newdim,dif; // d: current degree, cd: dimension of // space of invariants of degree d, // newdim: dimension the ideal generated // the primary invariants plus basis // elements, dif=n-i-newdim, i.e. the // number of new primary invairants that // should be added in this degree - ideal P,Pplus,sPplus,CI,B; // P: will contain primary invariants, // Pplus: P+B, CI: a complete // intersection with the same Hilbert // function as P ideal sP=groebner(P); dB=1; // used as degree bound int i=0; //-------------- loop that searches for primary invariants ------------------ while(1) // repeat until n primary invariants are { // found - p(1..2)=partial_molien(M,1,p(2)); // next term of the partial expansion - d=deg(p(1)); // degree where we'll search - cd=int(coef(p(1),v1)[2,1]); // dimension of the homogeneous space of // inviarants of degree d if (v) { " Computing primary invariants in degree ",d,":"; } B=invariant_basis_reynolds(REY,d,intvec(cd,6)); // basis of invariants of // degree d if (B[1]<>0) { Pplus=P+B; sPplus=groebner(Pplus); newdim=dim(sPplus); dif=n-i-newdim; } else { dif=0; } if (dif<>0) // we have to find dif new primary { // invariants if (cd<>dif) { P,sP,CI,dB=search(n,d,B,cd,P,sP,i,dif,dB,CI); // searching for dif invariants } else // i.e. we can take all of B { for(j=i+1;j<=i+dif;j++) { CI=CI+ideal(var(j)^d); } dB=dB+dif*(d-1); P=Pplus; sP=sPplus; } if (v) { for (j=1;j<=dif;j++) { " We find: "+string(P[i+j]); } } i=i+dif; if (i==n) // found all primary invariants { if (v) { ""; " We found all primary invariants."; ""; } return(matrix(P)); } } // done with degree d } } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2; ring R=0,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; matrix REY,M=reynolds_molien(A); matrix P=primary_char0(REY,M); print(P); } /////////////////////////////////////////////////////////////////////////////// proc primary_charp (matrix REY,string ring_name,list #) "USAGE: primary_charp(REY,ringname[,v]); REY: a representing the Reynolds operator, ringname: a giving the name of a ring where the Molien series is stored, v: an optional ASSUME: REY is the first return value of group_reynolds or reynolds_molien and 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 programme if v does not equal 0 RETURN: primary invariants (type ) of the invariant ring THEORY: Bases of homogeneous invariants are generated successively and those are chosen as primary invariants that lower the dimension of the ideal generated by the previously found invariants (see paper \"Generating a Noetherian Normalization of the Invariant Ring of a Finite Group\" by Decker, Heydtmann, Schreyer (1998)). EXAMPLE: example primary_charp; shows an example " { degBound=0; // ---------------- checking input and setting verbose mode ------------------- if (char(basering)==0) { "ERROR: primary_charp should only be used with rings of characteristic p>0."; return(); } if (size(#)>1) { "ERROR: primary_charp can only have three parameters."; return(); } if (size(#)==1) { if (typeof(#[1])<>"int") { "ERROR: The third parameter should be of type ."; return(); } else { int v=#[1]; } } else { int v=0; } def br=basering; 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(REY)<>n) { "ERROR: First parameter ought to be the Reynolds operator." return(); } if (typeof(`ring_name`)<>"ring") { "ERROR: Second parameter ought to be ring where "; " the Molien series is stored."; return(); } //---------------------------------------------------------------------------- if (v && voice<>2) { " We can start looking for primary invariants..."; ""; } if (v && voice==2) { ""; } //----------------------- initializing variables ----------------------------- int dB; setring `ring_name`; // the Molien series is stores here - poly p(1..2); // p(1) will be used for single terms of // the partial expansion, p(2) to store p(1..2)=partial_molien(M,1); // the intermediate result - poly v1=var(1); // we need v1 to split off coefficients // in the partial expansion of M (which // is in terms of the first variable) setring br; int j,d,cd,newdim,dif; // d: current degree, cd: dimension of // space of invariants of degree d, // newdim: dimension the ideal generated // the primary invariants plus basis // elements, dif=n-i-newdim, i.e. the // number of new primary invairants that // should be added in this degree - ideal P,Pplus,sPplus,CI,B; // P: will contain primary invariants, // Pplus: P+B, CI: a complete // intersection with the same Hilbert // function as P ideal sP=groebner(P); dB=1; // used as degree bound int i=0; //---------------- loop that searches for primary invariants ----------------- while(1) // repeat until n primary invariants are { // found setring `ring_name`; p(1..2)=partial_molien(M,1,p(2)); // next term of the partial expansion - d=deg(p(1)); // degree where we'll search - cd=int(coef(p(1),v1)[2,1]); // dimension of the homogeneous space of // inviarants of degree d setring br; if (v) { " Computing primary invariants in degree ",d,":"; } B=invariant_basis_reynolds(REY,d,intvec(cd,6)); // basis of invariants of // degree d if (ncols(B)0) { Pplus=P+B; sPplus=groebner(Pplus); newdim=dim(sPplus); dif=n-i-newdim; } else { dif=0; } if (dif<>0) // we have to find dif new primary { // invariants if (cd<>dif) { P,sP,CI,dB=p_search(n,d,B,cd,P,sP,i,dif,dB,CI); } else // i.e. we can take all of B { for(j=i+1;j>i+dif;j++) { CI=CI+ideal(var(j)^d); } dB=dB+dif*(d-1); P=Pplus; sP=sPplus; } if (v) { for (j=1;j<=size(P)-i;j++) { " We find: "+string(P[i+j]); } } i=size(P); if (i==n) // found all primary invariants { if (v) { ""; " We found all primary invariants."; ""; } return(matrix(P)); } } // done with degree d } } 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=group_reynolds(A); string newring="alskdfj"; molien(L[2..size(L)],newring); matrix P=primary_charp(L[1],newring); kill `newring`; print(P); } /////////////////////////////////////////////////////////////////////////////// proc primary_char0_no_molien (matrix REY, list #) "USAGE: primary_char0_no_molien(REY[,v]); REY: a representing the Reynolds operator, v: an optional ASSUME: REY is the first return value of group_reynolds or reynolds_molien DISPLAY: information about the various stages of the programme if v does not equal 0 RETURN: primary invariants (type ) of the invariant ring and an listing some of the degrees where no non-trivial homogeneous invariants are to be found THEORY: Bases of homogeneous invariants are generated successively and those are chosen as primary invariants that lower the dimension of the ideal generated by the previously found invariants (see paper \"Generating a Noetherian Normalization of the Invariant Ring of a Finite Group\" by Decker, Heydtmann, Schreyer (1998)). EXAMPLE: example primary_char0_no_molien; shows an example " { degBound=0; //-------------- checking input and setting verbose mode --------------------- if (char(basering)<>0) { "ERROR: primary_char0_no_molien should only be used with rings of"; " characteristic 0."; return(); } if (size(#)>1) { "ERROR: primary_char0_no_molien can only have two parameters."; return(); } if (size(#)==1) { if (typeof(#[1])<>"int") { "ERROR: The second parameter should be of type ."; return(); } else { int v=#[1]; } } else { 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(REY)<>n) { "ERROR: First parameter ought to be the Reynolds operator." return(); } //---------------------------------------------------------------------------- if (v && voice<>2) { " We can start looking for primary invariants..."; ""; } if (v && voice==2) { ""; } //----------------------- initializing variables ----------------------------- int dB; int j,d,cd,newdim,dif; // d: current degree, cd: dimension of // space of invariants of degree d, // newdim: dimension the ideal generated // the primary invariants plus basis // elements, dif=n-i-newdim, i.e. the // number of new primary invairants that // should be added in this degree - ideal P,Pplus,CI,B; // P: will contain primary invariants, // Pplus: P+B, CI: a complete // intersection with the same Hilbert // function as P ideal sP=groebner(P); dB=1; // used as degree bound - d=0; // initializing int i=0; intvec deg_vector; //------------------ loop that searches for primary invariants --------------- while(1) // repeat until n primary invariants are { // found - d++; // degree where we'll search if (v) { " Computing primary invariants in degree ",d,":"; } B=invariant_basis_reynolds(REY,d,intvec(-1,6)); // basis of invariants of // degree d if (B[1]<>0) { Pplus=P+B; newdim=dim(groebner(Pplus)); dif=n-i-newdim; } else { dif=0; deg_vector=deg_vector,d; } if (dif<>0) // we have to find dif new primary { // invariants cd=size(B); if (cd<>dif) { P,sP,CI,dB=search(n,d,B,cd,P,sP,i,dif,dB,CI); } else // i.e. we can take all of B { for(j=i+1;j<=i+dif;j++) { CI=CI+ideal(var(j)^d); } dB=dB+dif*(d-1); P=Pplus; sP=groebner(P); } if (v) { for (j=1;j<=dif;j++) { " We find: "+string(P[i+j]); } } i=i+dif; if (i==n) // found all primary invariants { if (v) { ""; " We found all primary invariants."; ""; } if (deg_vector==0) { return(matrix(P)); } else { return(matrix(P),compress(deg_vector)); } } } // done with degree d else { if (v) { " None here..."; } } } } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2; ring R=0,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; list L=group_reynolds(A); list l=primary_char0_no_molien(L[1]); print(l[1]); } /////////////////////////////////////////////////////////////////////////////// proc primary_charp_no_molien (matrix REY, list #) "USAGE: primary_charp_no_molien(REY[,v]); REY: a representing the Reynolds operator, v: an optional ASSUME: REY is the first return value of group_reynolds or reynolds_molien DISPLAY: information about the various stages of the programme if v does not equal 0 RETURN: primary invariants (type ) of the invariant ring and an listing some of the degrees where no non-trivial homogeneous invariants are to be found THEORY: Bases of homogeneous invariants are generated successively and those are chosen as primary invariants that lower the dimension of the ideal generated by the previously found invariants (see paper \"Generating a Noetherian Normalization of the Invariant Ring of a Finite Group\" by Decker, Heydtmann, Schreyer (1998)). EXAMPLE: example primary_charp_no_molien; shows an example " { degBound=0; //----------------- checking input and setting verbose mode ------------------ if (char(basering)==0) { "ERROR: primary_charp_no_molien should only be used with rings of"; " characteristic p>0."; return(); } if (size(#)>1) { "ERROR: primary_charp_no_molien can only have two parameters."; return(); } if (size(#)==1) { if (typeof(#[1])<>"int") { "ERROR: The second parameter should be of type ."; return(); } else { int v=#[1]; } } else { 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(REY)<>n) { "ERROR: First parameter ought to be the Reynolds operator." return(); } //---------------------------------------------------------------------------- if (v && voice<>2) { " We can start looking for primary invariants..."; ""; } if (v && voice==2) { ""; } //-------------------- initializing variables -------------------------------- int dB; int j,d,cd,newdim,dif; // d: current degree, cd: dimension of // space of invariants of degree d, // newdim: dimension the ideal generated // the primary invariants plus basis // elements, dif=n-i-newdim, i.e. the // number of new primary invairants that // should be added in this degree - ideal P,Pplus,sPplus,CI,B; // P: will contain primary invariants, // Pplus: P+B, CI: a complete // intersection with the same Hilbert // function as P ideal sP=groebner(P); dB=1; // used as degree bound - d=0; // initializing int i=0; intvec deg_vector; //------------------ loop that searches for primary invariants --------------- while(1) // repeat until n primary invariants are { // found - d++; // degree where we'll search if (v) { " Computing primary invariants in degree ",d,":"; } B=invariant_basis_reynolds(REY,d,intvec(-1,6)); // basis of invariants of // degree d if (B[1]<>0) { Pplus=P+B; sPplus=groebner(Pplus); newdim=dim(sPplus); dif=n-i-newdim; } else { dif=0; deg_vector=deg_vector,d; } if (dif<>0) // we have to find dif new primary { // invariants cd=size(B); if (cd<>dif) { P,sP,CI,dB=p_search(n,d,B,cd,P,sP,i,dif,dB,CI); } else // i.e. we can take all of B { for(j=i+1;j<=i+dif;j++) { CI=CI+ideal(var(j)^d); } dB=dB+dif*(d-1); P=Pplus; sP=sPplus; } if (v) { for (j=1;j<=size(P)-i;j++) { " We find: "+string(P[i+j]); } } i=size(P); if (i==n) // found all primary invariants { if (v) { ""; " We found all primary invariants."; ""; } if (deg_vector==0) { return(matrix(P)); } else { return(matrix(P),compress(deg_vector)); } } } // done with degree d else { if (v) { " None here..."; } } } } 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=group_reynolds(A); list l=primary_charp_no_molien(L[1]); print(l[1]); } /////////////////////////////////////////////////////////////////////////////// proc primary_charp_without (list #) "USAGE: primary_charp_without(G1,G2,...[,v]); G1,G2,...: generating a finite matrix group, v: an optional DISPLAY: information about the various stages of the programme if v does not equal 0 RETURN: primary invariants (type ) of the invariant ring THEORY: Bases of homogeneous invariants are generated successively and those are chosen as primary invariants that lower the dimension of the ideal generated by the previously found invariants (see paper \"Generating a Noetherian Normalization of the Invariant Ring of a Finite Group\" by Decker, Heydtmann, Schreyer (1998)). No Reynolds operator or Molien series is used. EXAMPLE: example primary_charp_without; shows an example " { degBound=0; //--------------------- checking input and setting verbose mode -------------- if (char(basering)==0) { "ERROR: primary_charp_without should not be used with rings of"; " characteristic 0."; return(); } if (size(#)==0) { "ERROR: There are no parameters."; return(); } if (typeof(#[size(#)])=="int") { int v=#[size(#)]; int gen_num=size(#)-1; if (gen_num==0) { "ERROR: There are no generators of a finite matrix group given."; return(); } } else { int v=0; int gen_num=size(#); } 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 for (int i=1;i<=gen_num;i++) { if (typeof(#[i])=="matrix") { if (nrows(#[i])<>n or ncols(#[i])<>n) { "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 first parameters should be a list of matrices"; return(); } } //---------------------------------------------------------------------------- if (v && voice==2) { ""; } //---------------------------- initializing variables ------------------------ int dB; int j,d,cd,newdim,dif; // d: current degree, cd: dimension of // space of invariants of degree d, // newdim: dimension the ideal generated // by the primary invariants plus basis // elements, dif=n-i-newdim, i.e. the // number of new primary invairants that // should be added in this degree - ideal P,Pplus,sPplus,CI,B; // P: will contain primary invariants, // Pplus: P+B, CI: a complete // intersection with the same Hilbert // function as P ideal sP=groebner(P); dB=1; // used as degree bound - d=0; // initializing i=0; intvec deg_vector; //-------------------- loop that searches for primary invariants ------------- while(1) // repeat until n primary invariants are { // found - d++; // degree where we'll search if (v) { " Computing primary invariants in degree ",d,":"; } B=invariant_basis(d,#[1..gen_num]); // basis of invariants of degree d if (B[1]<>0) { Pplus=P+B; sPplus=groebner(Pplus); newdim=dim(sPplus); dif=n-i-newdim; } else { dif=0; deg_vector=deg_vector,d; } if (dif<>0) // we have to find dif new primary { // invariants cd=size(B); if (cd<>dif) { P,sP,CI,dB=p_search(n,d,B,cd,P,sP,i,dif,dB,CI); } else // i.e. we can take all of B { for(j=i+1;j<=i+dif;j++) { CI=CI+ideal(var(j)^d); } dB=dB+dif*(d-1); P=Pplus; sP=sPplus; } if (v) { for (j=1;j<=size(P)-i;j++) { " We find: "+string(P[i+j]); } } i=size(P); if (i==n) // found all primary invariants { if (v) { ""; " We found all primary invariants."; ""; } return(matrix(P)); } } // done with degree d else { if (v) { " None here..."; } } } } example { "EXAMPLE:"; echo=2; ring R=2,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; matrix P=primary_charp_without(A); print(P); } /////////////////////////////////////////////////////////////////////////////// proc primary_invariants (list #) "USAGE: primary_invariants(G1,G2,...[,flags]); G1,G2,...: generating a finite matrix group, flags: an optional with three 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 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 characteristic of the base field divides the group order, the second component should give the size of intervals between canceling common factors in the expansion of the Molien series, 0 (the default) means only once after generating all terms, in prime characteristic also a negative number can be given to indicate that common factors should always be canceled when the expansion is simple (the root of the extension field occurs not among the coefficients) DISPLAY: information about the various stages of the programme if the third flag does not equal 0 RETURN: primary invariants (type ) of the invariant ring and if computable Reynolds operator (type ) and Molien series (type ) or ring name (type string) where the Molien series can be found in the char p case; if the first flag is 1 and we are in the non-modular case then an is returned giving some of the degrees where no non-trivial homogeneous invariants can be found THEORY: Bases of homogeneous invariants are generated successively and those are chosen as primary invariants that lower the dimension of the ideal generated by the previously found invariants (see paper \"Generating a Noetherian Normalization of the Invariant Ring of a Finite Group\" by Decker, Heydtmann, Schreyer (1998)). EXAMPLE: example primary_invariants; shows an example " { // ----------------- checking input and setting flags ------------------------ if (size(#)==0) { "ERROR: There are no parameters."; return(); } int ch=char(basering); // the algorithms depend very much on the // characteristic of the ground field 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 int gen_num; int mol_flag,v; if (typeof(#[size(#)])=="intvec") { if (size(#[size(#)])<>3) { "ERROR: should have three entries."; return(); } gen_num=size(#)-1; mol_flag=#[size(#)][1]; if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag==-1))) { "ERROR: the second component of should be >=0"; return(); } int interval=#[size(#)][2]; v=#[size(#)][3]; if (gen_num==0) { "ERROR: There are no generators of a finite matrix group given."; return(); } } else { gen_num=size(#); mol_flag=0; int interval=0; v=0; } for (int i=1;i<=gen_num;i++) { if (typeof(#[i])=="matrix") { if (nrows(#[i])<>n or ncols(#[i])<>n) { "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 first parameters should be a list of matrices"; return(); } } //---------------------------------------------------------------------------- if (mol_flag==0) { if (ch==0) { matrix REY,M=reynolds_molien(#[1..gen_num],intvec(mol_flag,interval,v)); // one will contain Reynolds operator and // the other enumerator and denominator // of Molien series matrix P=primary_char0(REY,M,v); return(P,REY,M); } else { list L=group_reynolds(#[1..gen_num],v); if (L[1]<>0) // testing whether we are in the modular { string newring="aksldfalkdsflkj"; // case if (minpoly==0) { if (v) { " We are dealing with the non-modular case."; } if (typeof(L[2])=="int") { molien(L[3..size(L)],newring,L[2],intvec(mol_flag,interval,v)); } else { molien(L[2..size(L)],newring,intvec(mol_flag,interval,v)); } matrix P=primary_charp(L[1],newring,v); return(P,L[1],newring); } else { if (v) { " 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."; ""; } list l=primary_charp_no_molien(L[1],v); if (size(l)==2) { return(l[1],L[1],l[2]); } else { return(l[1],L[1]); } } } else // the modular case { if (v) { " There is also no Molien series, we can make use of..."; ""; " We can start looking for primary invariants..."; ""; } return(primary_charp_without(#[1..gen_num],v)); } } } if (mol_flag==1) // the user wants no calculation of the { list L=group_reynolds(#[1..gen_num],v); // Molien series if (ch==0) { list l=primary_char0_no_molien(L[1],v); if (size(l)==2) { return(l[1],L[1],l[2]); } else { return(l[1],L[1]); } } else { if (L[1]<>0) // testing whether we are in the modular { list l=primary_charp_no_molien(L[1],v); // case if (size(l)==2) { return(l[1],L[1],l[2]); } else { return(l[1],L[1]); } } else // the modular case { if (v) { " We can start looking for primary invariants..."; ""; } return(primary_charp_without(#[1..gen_num],v)); } } } if (mol_flag==-1) { if (ch==0) { "ERROR: Characteristic 0 can only be simulated in characteristic p>>0."; return(); } list L=group_reynolds(#[1..gen_num],v); string newring="aksldfalkdsflkj"; molien(L[2..size(L)],newring,intvec(1,interval,v)); matrix P=primary_charp(L[1],newring,v); return(P,L[1],newring); } else // the user specified that the { if (ch==0) // characteristic divides the group order { "ERROR: The characteristic cannot divide the group order when it is 0."; return(); } if (v) { ""; } return(primary_charp_without(#[1..gen_num],v)); } } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2; ring R=0,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; list L=primary_invariants(A); print(L[1]); } /////////////////////////////////////////////////////////////////////////////// // This procedure finds dif primary invariants in degree d. It returns all // primary invariants found so far. The coefficients lie in a field of // characteristic 0. /////////////////////////////////////////////////////////////////////////////// proc search_random (int n,int d,ideal B,int cd,ideal P,int i,int dif,int dB,ideal CI,int max) { string answer; degBound=0; int j,k,test_dim,flag; matrix test_matrix[1][dif]; // the linear combination to test intvec h; // Hilbert series for (j=i+1;j<=i+dif;j++) { CI=CI+ideal(var(j)^d); // homogeneous polynomial of the same // degree as the one we're looking for // is added } ideal TEST; // h=hilb(std(CI),1); dB=dB+dif*(d-1); // used as degBound while (1) { test_matrix=matrix(B)*random(max,cd,dif); // degBound=dB; TEST=P+ideal(test_matrix); attrib(TEST,"isSB",1); test_dim=dim(TEST); // degBound=0; if (n-test_dim==i+dif) { break; } // degBound=dB; test_dim=dim(groebner(TEST)); // test_dim=dim(std(TEST,h)); // Hilbert driven std-calculation // degBound=0; if (n-test_dim==i+dif) { break; } else { "HELP: The ",dif," random combination(s) of the ",cd," basis elements with"; " coefficients in the range from -",max," to ",max," did not lower the"; " dimension by ",dif,". You can abort, try again or give a new range:"; answer=""; while (answer<>"n " && answer<>"y ") { " Do you want to abort (y/n)?"; answer=read(""); } if (answer=="y ") { flag=1; break; } answer=""; while (answer<>"n " && answer<>"y ") { " Do you want to try again (y/n)?"; answer=read(""); } if (answer=="n ") { flag=1; while (flag) { " Give a new > ",max," that bounds the range of coefficients:"; answer=read(""); for (j=1;j<=size(answer)-1;j++) { for (k=0;k<=9;k++) { if (answer[j]==string(k)) { break; } } if (k>9) { flag=1; break; } flag=0; } if (not(flag)) { execute("test_dim="+string(answer[1..size(answer)])); if (test_dim<=max) { flag=1; } else { max=test_dim; } } } } } } if (not(flag)) { P[(i+1)..(i+dif)]=test_matrix[1,1..dif]; } return(P,CI,dB); } /////////////////////////////////////////////////////////////////////////////// // This procedure finds at most dif primary invariants in degree d. It returns // all primary invariants found so far. The coefficients lie in the field of // characteristic p>0. /////////////////////////////////////////////////////////////////////////////// proc p_search_random (int n,int d,ideal B,int cd,ideal P,int i,int dif,int dB,ideal CI,int max) { string answer; degBound=0; int j,k,test_dim,flag; matrix test_matrix[1][dif]; // the linear combination to test intvec h; // Hilbert series ideal TEST; while (dif>0) { for (j=i+1;j<=i+dif;j++) { CI=CI+ideal(var(j)^d); // homogeneous polynomial of the same // degree as the one we're looking for // is added } // h=hilb(std(CI),1); dB=dB+dif*(d-1); // used as degBound test_matrix=matrix(B)*random(max,cd,dif); // degBound=dB; TEST=P+ideal(test_matrix); attrib(TEST,"isSB",1); test_dim=dim(TEST); // degBound=0; if (n-test_dim==i+dif) { break; } // degBound=dB; test_dim=dim(groebner(TEST)); // test_dim=dim(std(TEST,h)); // Hilbert driven std-calculation // degBound=0; if (n-test_dim==i+dif) { break; } else { "HELP: The ",dif," random combination(s) of the ",cd," basis elements with"; " coefficients in the range from -",max," to ",max," did not lower the"; " dimension by ",dif,". You can abort, try again, lower the number of"; " combinations searched for by 1 or give a larger coefficient range:"; answer=""; while (answer<>"n " && answer<>"y ") { " Do you want to abort (y/n)?"; answer=read(""); } if (answer=="y ") { flag=1; break; } answer=""; while (answer<>"n " && answer<>"y ") { " Do you want to try again (y/n)?"; answer=read(""); } if (answer=="n ") { answer=""; while (answer<>"n " && answer<>"y ") { " Do you want to lower the number of combinations by 1 (y/n)?"; answer=read(""); } if (answer=="y ") { dif=dif-1; } else { flag=1; while (flag) { " Give a new > ",max," that bounds the range of coefficients:"; answer=read(""); for (j=1;j<=size(answer)-1;j++) { for (k=0;k<=9;k++) { if (answer[j]==string(k)) { break; } } if (k>9) { flag=1; break; } flag=0; } if (not(flag)) { execute("test_dim="+string(answer[1..size(answer)])); if (test_dim<=max) { flag=1; } else { max=test_dim; } } } } } } CI=CI[1..i]; dB=dB-dif*(d-1); } if (dif && not(flag)) { P[(i+1)..(i+dif)]=test_matrix[1,1..dif]; } if (dif && flag) { P[n+1]=0; } return(P,CI,dB); } /////////////////////////////////////////////////////////////////////////////// proc primary_char0_random (matrix REY,matrix M,int max,list #) "USAGE: primary_char0_random(REY,M,r[,v]); REY: a representing the Reynolds operator, M: a 1x2 representing the Molien series, r: an where -|r| to |r| is the range of coefficients of the random combinations of bases elements, v: an optional 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 programme if v does not equal 0 RETURN: primary invariants (type ) of the invariant ring THEORY: Bases of homogeneous invariants are generated successively and random linear combinations are chosen as primary invariants that lower the dimension of the ideal generated by the previously found invariants (see \"Generating a Noetherian Normalization of the Invariant Ring of a Finite Group\" by Decker, Heydtmann, Schreyer (1998)). EXAMPLE: example primary_char0_random; shows an example " { degBound=0; if (char(basering)<>0) { "ERROR: primary_char0_random should only be used with rings of"; " characteristic 0."; return(); } //----------------- checking input and setting verbose mode ------------------ if (size(#)>1) { "ERROR: primary_char0_random can only have four parameters."; return(); } if (size(#)==1) { if (typeof(#[1])<>"int") { "ERROR: The fourth parameter should be of type ."; return(); } else { int v=#[1]; } } else { 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(REY)<>n) { "ERROR: First parameter ought to be the Reynolds operator." return(); } if (ncols(M)<>2 or nrows(M)<>1) { "ERROR: Second parameter ought to be the Molien series." return(); } //---------------------------------------------------------------------------- if (v && voice<>2) { " We can start looking for primary invariants..."; ""; } if (v && voice==2) { ""; } //------------------------- initializing variables --------------------------- int dB; poly p(1..2); // p(1) will be used for single terms of // the partial expansion, p(2) to store p(1..2)=partial_molien(M,1); // the intermediate result - poly v1=var(1); // we need v1 to split off coefficients // in the partial expansion of M (which // is in terms of the first variable) - int j,d,cd,newdim,dif; // d: current degree, cd: dimension of // space of invariants of degree d, // newdim: dimension the ideal generated // the primary invariants plus basis // elements, dif=n-i-newdim, i.e. the // number of new primary invairants that // should be added in this degree - ideal P,Pplus,CI,B; // P: will contain primary invariants, // Pplus: P+B,CI: a complete // intersection with the same Hilbert // function as P - dB=1; // used as degree bound int i=0; //-------------- loop that searches for primary invariants ------------------ while(1) // repeat until n primary invariants are { // found - p(1..2)=partial_molien(M,1,p(2)); // next term of the partial expansion - d=deg(p(1)); // degree where we'll search - cd=int(coef(p(1),v1)[2,1]); // dimension of the homogeneous space of // inviarants of degree d if (v) { " Computing primary invariants in degree ",d,":"; } B=invariant_basis_reynolds(REY,d,intvec(cd,6)); // basis of invariants of // degree d if (B[1]<>0) { Pplus=P+B; newdim=dim(groebner(Pplus)); dif=n-i-newdim; } else { dif=0; } if (dif<>0) // we have to find dif new primary { // invariants if (cd<>dif) { P,CI,dB=search_random(n,d,B,cd,P,i,dif,dB,CI,max); // searching for } // dif invariants - else // i.e. we can take all of B { for(j=i+1;j>i+dif;j++) { CI=CI+ideal(var(j)^d); } dB=dB+dif*(d-1); P=Pplus; } if (ncols(P)==i) { "WARNING: The return value is not a set of primary invariants, but"; " polynomials qualifying as the first ",i," primary invariants."; return(matrix(P)); } if (v) { for (j=1;j<=dif;j++) { " We find: "+string(P[i+j]); } } i=i+dif; if (i==n) // found all primary invariants { if (v) { ""; " We found all primary invariants."; ""; } return(matrix(P)); } } // done with degree d } } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2; ring R=0,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; matrix REY,M=reynolds_molien(A); matrix P=primary_char0_random(REY,M,1); print(P); } /////////////////////////////////////////////////////////////////////////////// proc primary_charp_random (matrix REY,string ring_name,int max,list #) "USAGE: primary_charp_random(REY,ringname,r[,v]); REY: a representing the Reynolds operator, ringname: a giving the name of a ring where the Molien series is stored, r: an where -|r| to |r| is the range of coefficients of the random combinations of bases elements, v: an optional ASSUME: REY is the first return value of group_reynolds or reynolds_molien and 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 programme if v does not equal 0 RETURN: primary invariants (type ) of the invariant ring THEORY: Bases of homogeneous invariants are generated successively and random linear combinations are chosen as primary invariants that lower the dimension of the ideal generated by the previously found invariants (see \"Generating a Noetherian Normalization of the Invariant Ring of a Finite Group\" by Decker, Heydtmann, Schreyer (1998)). EXAMPLE: example primary_charp_random; shows an example " { degBound=0; // ---------------- checking input and setting verbose mode ------------------ if (char(basering)==0) { "ERROR: primary_charp_random should only be used with rings of"; " characteristic p>0."; return(); } if (size(#)>1) { "ERROR: primary_charp_random can only have four parameters."; return(); } if (size(#)==1) { if (typeof(#[1])<>"int") { "ERROR: The fourth parameter should be of type ."; return(); } else { int v=#[1]; } } else { int v=0; } def br=basering; 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(REY)<>n) { "ERROR: First parameter ought to be the Reynolds operator." return(); } if (typeof(`ring_name`)<>"ring") { "ERROR: Second parameter ought to the name of a ring where the Molien"; " is stored."; return(); } //---------------------------------------------------------------------------- if (v && voice<>2) { " We can start looking for primary invariants..."; ""; } if (v && voice==2) { ""; } //----------------------- initializing variables ----------------------------- int dB; setring `ring_name`; // the Molien series is stores here - poly p(1..2); // p(1) will be used for single terms of // the partial expansion, p(2) to store p(1..2)=partial_molien(M,1); // the intermediate result - poly v1=var(1); // we need v1 to split off coefficients // in the partial expansion of M (which // is in terms of the first variable) setring br; int j,d,cd,newdim,dif; // d: current degree, cd: dimension of // space of invariants of degree d, // newdim: dimension the ideal generated // the primary invariants plus basis // elements, dif=n-i-newdim, i.e. the // number of new primary invairants that // should be added in this degree - ideal P,Pplus,CI,B; // P: will contain primary invariants, // Pplus: P+B, CI: a complete // intersection with the same Hilbert // function as P - dB=1; // used as degree bound int i=0; //---------------- loop that searches for primary invariants ----------------- while(1) // repeat until n primary invariants are { // found setring `ring_name`; p(1..2)=partial_molien(M,1,p(2)); // next term of the partial expansion - d=deg(p(1)); // degree where we'll search - cd=int(coef(p(1),v1)[2,1]); // dimension of the homogeneous space of // inviarants of degree d setring br; if (v) { " Computing primary invariants in degree ",d,":"; } B=invariant_basis_reynolds(REY,d,intvec(cd,6)); // basis of invariants of // degree d if (B[1]<>0) { Pplus=P+B; newdim=dim(groebner(Pplus)); dif=n-i-newdim; } else { dif=0; } if (dif<>0) // we have to find dif new primary { // invariants if (cd<>dif) { P,CI,dB=p_search_random(n,d,B,cd,P,i,dif,dB,CI,max); } else // i.e. we can take all of B { for(j=i+1;j>i+dif;j++) { CI=CI+ideal(var(j)^d); } dB=dB+dif*(d-1); P=Pplus; } if (ncols(P)==n+1) { "WARNING: The first return value is not a set of primary invariants,"; " but polynomials qualifying as the first ",i," primary invariants."; return(matrix(P)); } if (v) { for (j=1;j<=size(P)-i;j++) { " We find: "+string(P[i+j]); } } i=size(P); if (i==n) // found all primary invariants { if (v) { ""; " We found all primary invariants."; ""; } return(matrix(P)); } } // done with degree d } } 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=group_reynolds(A); string newring="alskdfj"; molien(L[2..size(L)],newring); matrix P=primary_charp_random(L[1],newring,1); kill `newring`; print(P); } /////////////////////////////////////////////////////////////////////////////// proc primary_char0_no_molien_random (matrix REY, int max, list #) "USAGE: primary_char0_no_molien_random(REY,r[,v]); REY: a representing the Reynolds operator, r: an where -|r| to |r| is the range of coefficients of the random combinations of bases elements, v: an optional ASSUME: REY is the first return value of group_reynolds or reynolds_molien DISPLAY: information about the various stages of the programme if v does not equal 0 RETURN: primary invariants (type ) of the invariant ring and an listing some of the degrees where no non-trivial homogeneous invariants are to be found THEORY: Bases of homogeneous invariants are generated successively and random linear combinations are chosen as primary invariants that lower the dimension of the ideal generated by the previously found invariants (see \"Generating a Noetherian Normalization of the Invariant Ring of a Finite Group\" by Decker, Heydtmann, Schreyer (1998)). EXAMPLE: example primary_char0_no_molien_random; shows an example " { degBound=0; //-------------- checking input and setting verbose mode --------------------- if (char(basering)<>0) { "ERROR: primary_char0_no_molien_random should only be used with rings of"; " characteristic 0."; return(); } if (size(#)>1) { "ERROR: primary_char0_no_molien_random can only have three parameters."; return(); } if (size(#)==1) { if (typeof(#[1])<>"int") { "ERROR: The third parameter should be of type ."; return(); } else { int v=#[1]; } } else { 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(REY)<>n) { "ERROR: First parameter ought to be the Reynolds operator." return(); } //---------------------------------------------------------------------------- if (v && voice<>2) { " We can start looking for primary invariants..."; ""; } if (v && voice==2) { ""; } //----------------------- initializing variables ----------------------------- int dB; int j,d,cd,newdim,dif; // d: current degree, cd: dimension of // space of invariants of degree d, // newdim: dimension the ideal generated // the primary invariants plus basis // elements, dif=n-i-newdim, i.e. the // number of new primary invairants that // should be added in this degree - ideal P,Pplus,CI,B; // P: will contain primary invariants, // Pplus: P+B, CI: a complete // intersection with the same Hilbert // function as P - dB=1; // used as degree bound - d=0; // initializing int i=0; intvec deg_vector; //------------------ loop that searches for primary invariants --------------- while(1) // repeat until n primary invariants are { // found - d++; // degree where we'll search if (v) { " Computing primary invariants in degree ",d,":"; } B=invariant_basis_reynolds(REY,d,intvec(-1,6)); // basis of invariants of // degree d if (B[1]<>0) { Pplus=P+B; newdim=dim(groebner(Pplus)); dif=n-i-newdim; } else { dif=0; deg_vector=deg_vector,d; } if (dif<>0) // we have to find dif new primary { // invariants cd=size(B); if (cd<>dif) { P,CI,dB=search_random(n,d,B,cd,P,i,dif,dB,CI,max); } else // i.e. we can take all of B { for(j=i+1;j<=i+dif;j++) { CI=CI+ideal(var(j)^d); } dB=dB+dif*(d-1); P=Pplus; } if (ncols(P)==i) { "WARNING: The first return value is not a set of primary invariants,"; " but polynomials qualifying as the first ",i," primary invariants."; return(matrix(P)); } if (v) { for (j=1;j<=dif;j++) { " We find: "+string(P[i+j]); } } i=i+dif; if (i==n) // found all primary invariants { if (v) { ""; " We found all primary invariants."; ""; } if (deg_vector==0) { return(matrix(P)); } else { return(matrix(P),compress(deg_vector)); } } } // done with degree d else { if (v) { " None here..."; } } } } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2; ring R=0,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; list L=group_reynolds(A); list l=primary_char0_no_molien_random(L[1],1); print(l[1]); } /////////////////////////////////////////////////////////////////////////////// proc primary_charp_no_molien_random (matrix REY, int max, list #) "USAGE: primary_charp_no_molien_random(REY,r[,v]); REY: a representing the Reynolds operator, r: an where -|r| to |r| is the range of coefficients of the random combinations of bases elements, v: an optional ASSUME: REY is the first return value of group_reynolds or reynolds_molien DISPLAY: information about the various stages of the programme if v does not equal 0 RETURN: primary invariants (type ) of the invariant ring and an listing some of the degrees where no non-trivial homogeneous invariants are to be found THEORY: Bases of homogeneous invariants are generated successively and random linear combinations are chosen as primary invariants that lower the dimension of the ideal generated by the previously found invariants (see \"Generating a Noetherian Normalization of the Invariant Ring of a Finite Group\" by Decker, Heydtmann, Schreyer (1998)). EXAMPLE: example primary_charp_no_molien_random; shows an example " { degBound=0; //----------------- checking input and setting verbose mode ------------------ if (char(basering)==0) { "ERROR: primary_charp_no_molien_random should only be used with rings of"; " characteristic p>0."; return(); } if (size(#)>1) { "ERROR: primary_charp_no_molien_random can only have three parameters."; return(); } if (size(#)==1) { if (typeof(#[1])<>"int") { "ERROR: The third parameter should be of type ."; return(); } else { int v=#[1]; } } else { 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(REY)<>n) { "ERROR: First parameter ought to be the Reynolds operator." return(); } //---------------------------------------------------------------------------- if (v && voice<>2) { " We can start looking for primary invariants..."; ""; } if (v && voice==2) { ""; } //-------------------- initializing variables -------------------------------- int dB; int j,d,cd,newdim,dif; // d: current degree, cd: dimension of // space of invariants of degree d, // newdim: dimension the ideal generated // the primary invariants plus basis // elements, dif=n-i-newdim, i.e. the // number of new primary invairants that // should be added in this degree - ideal P,Pplus,CI,B; // P: will contain primary invariants, // Pplus: P+B, CI: a complete // intersection with the same Hilbert // function as P - dB=1; // used as degree bound - d=0; // initializing int i=0; intvec deg_vector; //------------------ loop that searches for primary invariants --------------- while(1) // repeat until n primary invariants are { // found - d++; // degree where we'll search if (v) { " Computing primary invariants in degree ",d,":"; } B=invariant_basis_reynolds(REY,d,intvec(-1,6)); // basis of invariants of // degree d if (B[1]<>0) { Pplus=P+B; newdim=dim(groebner(Pplus)); dif=n-i-newdim; } else { dif=0; deg_vector=deg_vector,d; } if (dif<>0) // we have to find dif new primary { // invariants cd=size(B); if (cd<>dif) { P,CI,dB=p_search_random(n,d,B,cd,P,i,dif,dB,CI,max); } else // i.e. we can take all of B { for(j=i+1;j<=i+dif;j++) { CI=CI+ideal(var(j)^d); } dB=dB+dif*(d-1); P=Pplus; } if (ncols(P)==n+1) { "WARNING: The first return value is not a set of primary invariants,"; " but polynomials qualifying as the first ",i," primary invariants."; return(matrix(P)); } if (v) { for (j=1;j<=size(P)-i;j++) { " We find: "+string(P[i+j]); } } i=size(P); if (i==n) // found all primary invariants { if (v) { ""; " We found all primary invariants."; ""; } if (deg_vector==0) { return(matrix(P)); } else { return(matrix(P),compress(deg_vector)); } } } // done with degree d else { if (v) { " None here..."; } } } } 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=group_reynolds(A); list l=primary_charp_no_molien_random(L[1],1); print(l[1]); } /////////////////////////////////////////////////////////////////////////////// proc primary_charp_without_random (list #) "USAGE: primary_charp_without_random(G1,G2,...,r[,v]); G1,G2,...: generating a finite matrix group, r: an where -|r| to |r| is the range of coefficients of the random combinations of bases elements, v: an optional DISPLAY: information about the various stages of the programme if v does not equal 0 RETURN: primary invariants (type ) of the invariant ring THEORY: Bases of homogeneous invariants are generated successively and random linear combinations are chosen as primary invariants that lower the dimension of the ideal generated by the previously found invariants (see \"Generating a Noetherian Normalization of the Invariant Ring of a Finite Group\" by Decker, Heydtmann, Schreyer (1998)). No Reynolds operator or Molien series is used. EXAMPLE: example primary_charp_without_random; shows an example " { degBound=0; //--------------------- checking input and setting verbose mode -------------- if (char(basering)==0) { "ERROR: primary_charp_without_random should only be used with rings of"; " characteristic 0."; return(); } if (size(#)<2) { "ERROR: There are too few parameters."; return(); } if (typeof(#[size(#)])=="int" && typeof(#[size(#)-1])=="int") { int v=#[size(#)]; int max=#[size(#)-1]; int gen_num=size(#)-2; if (gen_num==0) { "ERROR: There are no generators of a finite matrix group given."; return(); } } else { if (typeof(#[size(#)])=="int") { int max=#[size(#)]; int v=0; int gen_num=size(#)-1; } else { "ERROR: The last parameter should be an ."; return(); } } 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 for (int i=1;i<=gen_num;i++) { if (typeof(#[i])=="matrix") { if (nrows(#[i])<>n or ncols(#[i])<>n) { "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 first parameters should be a list of matrices"; return(); } } //---------------------------------------------------------------------------- if (v && voice==2) { ""; } //---------------------------- initializing variables ------------------------ int dB; int j,d,cd,newdim,dif; // d: current degree, cd: dimension of // space of invariants of degree d, // newdim: dimension the ideal generated // the primary invariants plus basis // elements, dif=n-i-newdim, i.e. the // number of new primary invairants that // should be added in this degree - ideal P,Pplus,CI,B; // P: will contain primary invariants, // Pplus: P+B, CI: a complete // intersection with the same Hilbert // function as P - dB=1; // used as degree bound - d=0; // initializing i=0; intvec deg_vector; //-------------------- loop that searches for primary invariants ------------- while(1) // repeat until n primary invariants are { // found - d++; // degree where we'll search if (v) { " Computing primary invariants in degree ",d,":"; } B=invariant_basis(d,#[1..gen_num]); // basis of invariants of degree d if (B[1]<>0) { Pplus=P+B; newdim=dim(groebner(Pplus)); dif=n-i-newdim; } else { dif=0; deg_vector=deg_vector,d; } if (dif<>0) // we have to find dif new primary { // invariants cd=size(B); if (cd<>dif) { P,CI,dB=p_search_random(n,d,B,cd,P,i,dif,dB,CI,max); } else // i.e. we can take all of B { for(j=i+1;j<=i+dif;j++) { CI=CI+ideal(var(j)^d); } dB=dB+dif*(d-1); P=Pplus; } if (ncols(P)==n+1) { "WARNING: The first return value is not a set of primary invariants,"; " but polynomials qualifying as the first ",i," primary invariants."; return(matrix(P)); } if (v) { for (j=1;j<=size(P)-i;j++) { " We find: "+string(P[i+j]); } } i=size(P); if (i==n) // found all primary invariants { if (v) { ""; " We found all primary invariants."; ""; } return(matrix(P)); } } // done with degree d else { if (v) { " None here..."; } } } } example { "EXAMPLE:"; echo=2; ring R=2,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; matrix P=primary_charp_without_random(A,1); print(P); } /////////////////////////////////////////////////////////////////////////////// proc primary_invariants_random (list #) "USAGE: primary_invariants_random(G1,G2,...,r[,flags]); G1,G2,...: generating a finite matrix group, r: an where -|r| to |r| is the range of coefficients of the random combinations of bases elements, flags: an optional with three 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 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 characteristic of the base field divides the group order, the second component should give the size of intervals between canceling common factors in the expansion of the Molien series, 0 (the default) means only once after generating all terms, in prime characteristic also a negative number can be given to indicate that common factors should always be canceled when the expansion is simple (the root of the extension field does not occur among the coefficients) DISPLAY: information about the various stages of the programme if the third flag does not equal 0 RETURN: primary invariants (type ) of the invariant ring and if computable Reynolds operator (type ) and Molien series (type ), if the first flag is 1 and we are in the non-modular case then an is returned giving some of the degrees where no non-trivial homogeneous invariants can be found THEORY: Bases of homogeneous invariants are generated successively and random linear combinations are chosen as primary invariants that lower the dimension of the ideal generated by the previously found invariants (see \"Generating a Noetherian Normalization of the Invariant Ring of a Finite Group\" by Decker, Heydtmann, Schreyer (1998)). EXAMPLE: example primary_invariants_random; shows an example " { // ----------------- checking input and setting flags ------------------------ if (size(#)<2) { "ERROR: There are too few parameters."; return(); } int ch=char(basering); // the algorithms depend very much on the // characteristic of the ground field 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 int gen_num; int mol_flag,v; if (typeof(#[size(#)])=="intvec" && typeof(#[size(#)-1])=="int") { if (size(#[size(#)])<>3) { "ERROR: should have three entries."; return(); } gen_num=size(#)-2; mol_flag=#[size(#)][1]; if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0))) { "ERROR: the second component of should be >=0"; return(); } int interval=#[size(#)][2]; v=#[size(#)][3]; int max=#[size(#)-1]; if (gen_num==0) { "ERROR: There are no generators of a finite matrix group given."; return(); } } else { if (typeof(#[size(#)])=="int") { gen_num=size(#)-1; mol_flag=0; int interval=0; v=0; int max=#[size(#)]; } else { "ERROR: If the two last parameters are not and , the last"; " parameter should be an ."; return(); } } for (int i=1;i<=gen_num;i++) { if (typeof(#[i])=="matrix") { if (nrows(#[i])<>n or ncols(#[i])<>n) { "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 first parameters should be a list of matrices"; return(); } } //---------------------------------------------------------------------------- if (mol_flag==0) { if (ch==0) { matrix REY,M=reynolds_molien(#[1..gen_num],intvec(0,interval,v)); // one will contain Reynolds operator and // the other enumerator and denominator // of Molien series matrix P=primary_char0_random(REY,M,max,v); return(P,REY,M); } else { list L=group_reynolds(#[1..gen_num],v); if (L[1]<>0) // testing whether we are in the modular { string newring="aksldfalkdsflkj"; // case if (minpoly==0) { if (v) { " We are dealing with the non-modular case."; } if (typeof(L[2])=="int") { molien(L[3..size(L)],newring,L[2],intvec(0,interval,v)); } else { molien(L[2..size(L)],newring,intvec(0,interval,v)); } matrix P=primary_charp_random(L[1],newring,max,v); return(P,L[1],newring); } else { if (v) { " 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."; ""; } list l=primary_charp_no_molien_random(L[1],max,v); if (size(l)==2) { return(l[1],L[1],l[2]); } else { return(l[1],L[1]); } } } else // the modular case { if (v) { " There is also no Molien series, we can make use of..."; ""; " We can start looking for primary invariants..."; ""; } return(primary_charp_without_random(#[1..gen_num],max,v)); } } } if (mol_flag==1) // the user wants no calculation of the { list L=group_reynolds(#[1..gen_num],v); // Molien series if (ch==0) { list l=primary_char0_no_molien_random(L[1],max,v); if (size(l)==2) { return(l[1],L[1],l[2]); } else { return(l[1],L[1]); } } else { if (L[1]<>0) // testing whether we are in the modular { list l=primary_charp_no_molien_random(L[1],max,v); // case if (size(l)==2) { return(l[1],L[1],l[2]); } else { return(l[1],L[1]); } } else // the modular case { if (v) { " We can start looking for primary invariants..."; ""; } return(primary_charp_without_random(#[1..gen_num],max,v)); } } } if (mol_flag==-1) { if (ch==0) { "ERROR: Characteristic 0 can only be simulated in characteristic p>>0."; return(); } list L=group_reynolds(#[1..gen_num],v); string newring="aksldfalkdsflkj"; if (typeof(L[2])=="int") { molien(L[3..size(L)],newring,L[2],intvec(0,interval,v)); } else { molien(L[2..size(L)],newring,intvec(0,interval,v)); } matrix P=primary_charp_random(L[1],newring,max,v); return(P,L[1],newring); } else // the user specified that the { if (ch==0) // characteristic divides the group order { "ERROR: The characteristic cannot divide the group order when it is 0."; return(); } if (v) { ""; } return(primary_charp_without_random(#[1..gen_num],max,v)); } } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2; ring R=0,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; list L=primary_invariants_random(A,1); print(L[1]); } /////////////////////////////////////////////////////////////////////////////// proc concat_intmat(intmat A,intmat B) { int n=nrows(A); int m1=ncols(A); int m2=ncols(B); intmat C[n][m1+m2]; C[1..n,1..m1]=A[1..n,1..m1]; C[1..n,m1+1..m1+m2]=B[1..n,1..m2]; return(C); } /////////////////////////////////////////////////////////////////////////////// proc power_products(intvec deg_vec,int d) "USAGE: power_products(dv,d); dv: an giving the degrees of homogeneous polynomials, d: the degree of the desired power products RETURN: a size(dv)*m where each column ought to be interpreted as containing the exponents of the corresponding polynomials. The product of the powers is then homogeneous of degree d. EXAMPLE: example power_products; shows an example " { ring R=0,x,dp; if (d<=0) { "ERROR: The may not be <= 0"; return(); } int d_neu,j,nc; int s=size(deg_vec); intmat PP[s][1]; intmat TEST[s][1]; for (int i=1;i<=s;i++) { if (i<0) { "ERROR: The entries of may not be <= 0"; return(); } d_neu=d-deg_vec[i]; if (d_neu>0) { intmat PPd_neu=power_products(intvec(deg_vec[i..s]),d_neu); if (size(ideal(PPd_neu))<>0) { nc=ncols(PPd_neu); intmat PPd_neu_gross[s][nc]; PPd_neu_gross[i..s,1..nc]=PPd_neu[1..s-i+1,1..nc]; for (j=1;j<=nc;j++) { PPd_neu_gross[i,j]=PPd_neu_gross[i,j]+1; } PP=concat_intmat(PP,PPd_neu_gross); kill PPd_neu_gross; } kill PPd_neu; } if (d_neu==0) { intmat PPd_neu[s][1]; PPd_neu[i,1]=1; PP=concat_intmat(PP,PPd_neu); kill PPd_neu; } } if (matrix(PP)<>matrix(TEST)) { PP=compress(PP); } return(PP); } example { "EXAMPLE:"; echo=2; intvec dv=5,5,5,10,10; print(power_products(dv,10)); print(power_products(dv,7)); } /////////////////////////////////////////////////////////////////////////////// static proc old_secondary_char0 (matrix P, matrix REY, matrix M, list #) "USAGE: secondary_char0(P,REY,M[,v]); P: a 1xn with primary invariants, REY: a gxn representing the Reynolds operator, M: a 1x2 giving numerator and denominator of the Molien series, v: an optional 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(), M the return value of molien() or the second one of reynolds_molien() or the third one of primary_invariants() RETURN: secondary invariants of the invariant ring (type ) and irreducible secondary invariants (type ) 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 and using these images or their power products such that they are linearly independent modulo the primary invariants (see paper \"Some Algorithms in Invariant Theory of Finite Groups\" by Kemper and Steel (1997)). " { def br=basering; degBound=0; //----------------- checking input and setting verbose mode ------------------ if (char(br)<>0) { "ERROR: secondary_char0 should only be used with rings of characteristic 0."; return(); } int i; if (size(#)>0) { if (typeof(#[size(#)])=="int") { int v=#[size(#)]; } else { int v=0; } } else { int v=0; } 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 (ncols(M)<>2 or nrows(M)<>1) { "ERROR: The third parameter ought to be the Molien series." return(); } if (v && voice==2) { ""; } int j, m, counter; //- finding the polynomial giving number and degrees of secondary invariants - poly p=1; for (j=1;j<=n;j++) // calculating the denominator of the { p=p*(1-var(1)^deg(P[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=br,ideal(0); p=1/slead(p)*p; // smallest term of p needs to be 1 if (v) { " Polynomial telling us where to look for secondary invariants:"; " "+string(p); ""; } 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) { "In degree 0 we have: 1"; ""; } //-------------------------- initializing variables -------------------------- intmat PP; poly pp; int k; intvec deg_vec; ideal sP=groebner(ideal(P)); ideal TEST,B,IS; ideal S=1; // 1 is the first secondary invariant - //--------------------- 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 0 { // elements in the current degree (i-1) if (v) { "Searching in degree ",i-1,", we need to find ",int(dimmat[i,1])," invariant(s)..."; } TEST=sP; counter=0; // we'll count up to degvec[i] if (IS[1]<>0) { PP=power_products(deg_vec,i-1); // finding power products of irreducible } // secondary invariants if (size(ideal(PP))<>0) { for (j=1;j<=ncols(PP);j++) // going through all the power products { pp=1; for (k=1;k<=nrows(PP);k++) { pp=pp*IS[1,k]^PP[k,j]; } if (reduce(pp,TEST)<>0) { S=S,pp; counter++; if (v) { " We find: "+string(pp); } if (int(dimmat[i,1])<>counter) { // TEST=std(TEST+ideal(NF(pp,TEST))); // should be replaced by next // line soon TEST=std(TEST,pp); } else { break; } } } } if (int(dimmat[i,1])<>counter) { 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 { S=S,B; // B IS=IS+B; if (deg_vec[1]==0) { deg_vec=i-1; if (v) { " We find: "+string(B[1]); } for (j=2;j<=int(dimmat[i,1]);j++) { deg_vec=deg_vec,i-1; if (v) { " We find: "+string(B[j]); } } } else { for (j=1;j<=int(dimmat[i,1]);j++) { deg_vec=deg_vec,i-1; if (v) { " We find: "+string(B[j]); } } } } else { j=0; // j goes through all of B - while (int(dimmat[i,1])<>counter) // need to find dimmat[i,1] { // invariants that are linearly // independent modulo TEST j++; if (reduce(B[j],TEST)<>0) // B[j] should be added { S=S,B[j]; IS=IS+ideal(B[j]); if (deg_vec[1]==0) { deg_vec[1]=i-1; } else { deg_vec=deg_vec,i-1; } counter++; if (v) { " We find: "+string(B[j]); } if (int(dimmat[i,1])<>counter) { //TEST=std(TEST+ideal(NF(B[j],TEST))); // should be replaced by // next line TEST=std(TEST,B[j]); } } } } } if (v) { ""; } } } if (v) { " We're done!"; ""; } return(matrix(S),matrix(IS)); } /////////////////////////////////////////////////////////////////////////////// proc secondary_char0 (matrix P, matrix REY, matrix M, list #) "USAGE: secondary_char0(P,REY,M[,v][,\"old\"]); @* P: a 1xn with homogeneous primary invariants, where n is the number of variables of the basering; @* REY: a gxn representing the Reynolds operator, where g the size of the corresponding group; @* M: a 1x2 giving numerator and denominator of the Molien series; @* v: an optional ; @* \"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 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() or the second one of reynolds_molien() or the third one of primary_invariants() RETURN: Homogeneous secondary invariants and irreducible secondary invariants of the invariant ring (both type ) DISPLAY: Information on the progress of the computations if v is an integer different from 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 secondary invariants using Groebner basis techniques (see S. King: Fast Computation of Secondary Invariants). 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,\"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 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 " { //---------- Internal parameters, whose choice might improve the performance ----- 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. def br=basering; //----------------- checking input and setting verbose mode ------------------ if (char(br)<>0) { "ERROR: secondary_char0 can only be used with rings of characteristic 0."; " Try secondary_charp"; return(); } int i; if (size(#)>0) { if (typeof(#[1])=="int") { int v=#[1]; } 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_char0(P,REY,M,#[1]); return(S,IS); } else { matrix S,IS = old_secondary_char0(P,REY,M); 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 (ncols(M)<>2 or nrows(M)<>1) { "ERROR: The third parameter ought to be the Molien series."; return(); } if (v && voice==2) { ""; } int j, m, counter, irrcounter; //- finding the polynomial giving number and degrees of secondary invariants - poly p=1; for (j=1;j<=n;j++) // calculating the denominator of the { p=p*(1-var(1)^deg(P[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=br,ideal(0); p=1/slead(p)*p; // smallest [constant] 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; } } } if (v) { "In degree 0 we have: 1"; ""; } //-------------------------- initializing variables -------------------------- 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 ideal sP_Reductor; // will contain sP union Reductor. int sPOffset = ncols(sP); 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 int NrIS; poly helpP; ideal helpI; ideal Indicator; // will tell us which candidates for sec. inv. we can choose 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;i0) // when it is == 0 we need to find no { // elements in the current degree (i-1) if (v) { "Searching in degree ",i-1,", we need to find ",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); 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"none") { Mul1 = ISSort[k]; } else { Mul1 = ideal(0); } if ((int(dimmat[i-k,1])>0) && (Mul1[1]<>0)) { for (minD=k;minD=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]; } } ProdCand = simplify(Mul1*Mul2,4); // sP union SaveRed union Reductor is a homogeneous Groebner basis // up to degree i-1. 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 { 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", int(dimmat[i,1]),"secondaries in degree",i-1; } if (int(dimmat[i,1])<>counter) { Reductor = ideal(helpP); attrib(Reductor, "isSB",1); 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 { break; } } } } } } } } NrIS = int(dimmat[i,1])-counter; // The remaining NrIS sec. inv. are irreducible! if (NrIS>0) // need more than all the power products { if (v) { " There are ",NrIS,"irreducible secondary invariants in degree ", i-1; } mon = kbase(sP_Reductor,i-1); ii = ncols(mon); if ((i>IrrSwitch) or (ii>200)) // use sparse algorithm { if (ii+counter==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")+ii; 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 all",NrIS,"irreducibles in degree ",i-1;} } else { irrcounter=0; j=0; // j goes through all of mon - // Compare the comments on the computation of reducible sec. inv.! while (counter <> int(dimmat[i,1])) { if ((j mod MonStep) == 0) { if ((j+MonStep) <= ii) { ideal tmp = normalize(evaluate_reynolds(REY,ideal(mon[j+1..j+MonStep]))); B[j+1..j+MonStep] = tmp[1..MonStep]; tmp = reduce(tmp,sP_Reductor); Indicator[j+1..j+MonStep] = tmp[1..MonStep]; kill tmp; } else { ideal tmp = normalize(evaluate_reynolds(REY,ideal(mon[j+1..ii]))); B[j+1..ii] = tmp[1..ii-j]; tmp = reduce(tmp,sP_Reductor); Indicator[j+1..ii] = tmp[1..ii-j]; kill tmp; } } 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", irrcounter,"of", NrIS ,"irr. sec. inv. in degree ",i-1; } Reductor = ideal(helpP); attrib(Reductor, "isSB",1); Indicator=reduce(Indicator,Reductor); SizeSave++; SaveRed[SizeSave] = helpP; attrib(SaveRed, "isSB",1); sP_Reductor[sPOffset+SizeSave] = helpP; attrib(sP_Reductor, "isSB",1); } B[j]=0; Indicator[j]=0; } } } // if i>IrrSwitch else // use fast algorithm { 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 // Reynolds operator that are linearly // 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 all",NrIS,"irreducibles in degree ",i-1;} } else { irrcounter=0; j=0; // j goes through all of B - // Compare the comments on the computation of reducible sec. inv.! Indicator = reduce(B,sP_Reductor); 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", irrcounter,"of", NrIS,"irr. sec. inv. in degree ",i-1; } Reductor = ideal(helpP); attrib(Reductor, "isSB",1); Indicator=reduce(Indicator,Reductor); SizeSave++; SaveRed[SizeSave] = helpP; attrib(SaveRed, "isSB",1); sP_Reductor[sPOffset+SizeSave] = helpP; attrib(sP_Reductor, "isSB",1); } B[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!"; ""; } degBound = dgb; // Prepare return: int TotalNumber; for (k=1;k<=m;k++) { TotalNumber = TotalNumber + int(dimmat[k,1]); } matrix S[1][TotalNumber]; S[1,1]=1; j=1; for (k=1;k with homogeneous primary invariants, where n is the number of variables of the basering; @* REY: a gxn representing the Reynolds operator, where g the size of the corresponding group; @* M: a 1x2 giving numerator and denominator of the Molien series; @* v: an optional ; @* \"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 ) 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 reynolds_molien() or the third one of primary_invariants() DISPLAY: Information on the progress of computations 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 secondary invariants using Groebner basis techniques (see S. King: Fast Computation of Secondary Invariants). 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, 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 " { //---------- Internal parameters, whose choice might improve the performance ----- 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, 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; //----------------- checking input and setting verbose mode ------------------ if (char(br)<>0) { "ERROR: irred_secondary_char0 can only be used with rings of characteristic 0."; " Try irred_secondary_charp"; return(); } int i; if (size(#)>0) { if (typeof(#[1])=="int") { int v=#[1]; } else { int v=0; } } 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 // 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 (ncols(M)<>2 or nrows(M)<>1) { "ERROR: The third parameter ought to be the Molien series." return(); } if (v && voice==2) { ""; } int j, m, counter, irrcounter; //- finding the polynomial giving number and degrees of secondary invariants - poly p=1; for (j=1;j<=n;j++) // calculating the denominator of the { p=p*(1-var(1)^deg(P[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=br,ideal(0); p=1/slead(p)*p; // smallest [constant] 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) { "There are"; 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; } } } if (v) { "In degree 0 we have: 1"; ""; } //-------------------------- initializing variables -------------------------- 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)); 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 NoPP; // It NoPP==1, we must not use Power Products ideal Reductor,SaveRed; // sP union Reductor is a Groebner basis up to degree i-1 int SizeSave; list RedSSort; // sec. Inv. reduced by sP, sorted first by degree and then // by the minimal degree of a non-constant invariant factor list RedISSort; // irr. sec. Inv. reduced by sP, 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;i0) // when it is == 0 we need to find no { // elements in the current degree (i-1) if (v) { "Searching in degree ",i-1,". There are ", int(dimmat[i,1])," secondary invariant(s)..."; } counter = 0; // we'll count up to dimmat[i,1] Reductor = ideal(0); helpint = 0; SaveRed = Reductor; SizeSave = 0; attrib(Reductor,"isSB",1); attrib(SaveRed,"isSB",1); // Case A: We use PP if ( NoPP==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. // 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". // Finally, as we are not interested in the reducible sec. inv., we will only // work with their reduction modulo sP --- this allows to detect a secondary invariant // without to actually compute it! for (k=1;k"none") { Mul1 = RedISSort[k]; } else { Mul1 = ideal(0); } if ((int(dimmat[i-k,1])>0) && (Mul1[1]<>0)) { for (minD=i-k-1;minD>0;minD--) { if (int(dimmat[i,1])==counter) { break; } for (k2=1;k2 <= ((attrib(RedSSort[i-k-1][minD],"size")-1) div pieces)+1; k2++) { if (int(dimmat[i,1])==counter) { break; } Mul2=ideal(0); if (attrib(RedSSort[i-k-1][minD],"size")>=k2*pieces) { for (k3=1;k3<=pieces;k3++) { Mul2[k3] = RedSSort[i-k-1][minD][((k2-1)*pieces)+k3]; } } else { for (k3=1;k3<=(attrib(RedSSort[i-k-1][minD],"size") mod pieces);k3++) { Mul2[k3] = RedSSort[i-k-1][minD][((k2-1)*pieces)+k3]; } } ProdCand = simplify(Mul1*Mul2,4); 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 ReducedCandidates[ii] is the reduction of an invariant // that we can take 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(RedSSort[i-1][k],"size")+1; RedSSort[i-1][k][saveAttr] = ReducedCandidates[ii]; // By construction, this is the reduction of a reducible s.i. // of degree i-1. attrib(RedSSort[i-1][k],"size",saveAttr); if (v) { " We found reducible sec. inv. number ",counter; } if (int(dimmat[i,1])<>counter) { Reductor = ideal(helpP); attrib(Reductor, "isSB",1); 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 { break; } } // new reducible sec. inv. found } // loop through ProdCand } // if there is some reducible sec. inv. attrib(SaveRed, "isSB",1); attrib(Reductor, "isSB",1); } } } } TotalReductor = sP_Reductor; } // if NoPP==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. } // dealing with reducible sec. inv. // The remaining sec. inv. are irreducible! // Reason: If we use PP, TotalReductor detects the <

>-module generated by // all secondaries // If we don't, TotalReductor detects the algebra <

> if (int(dimmat[i,1])<>counter) // need more than all the power products // or work without power products { if (v) { " Looking for irreducible secondary invariants in degree ", i-1; } mon = kbase(TotalReductor,i-1); ii = ncols(mon); if (mon[1]<>0) { if ((i>IrrSwitch) or (ii>200)) // use sparse algorithm { if (NoPP==0 && ii+counter==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")+ii; 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 ",size(B)," irred. sec. inv.";} } else { irrcounter=0; j=0; // j goes through all of mon - // Compare the comments on the computation of reducible sec. inv.! while (j < ii) // we will break the loop if counter==int(dimmat[i,1]) { if ((j mod MonStep) == 0) { if ((j+MonStep) <= ii) { ideal tmp = normalize(evaluate_reynolds(REY,ideal(mon[j+1..j+MonStep]))); B[j+1..j+MonStep] = tmp[1..MonStep]; if (NoPP == 0) // we are still working with PowerProducts, // hence, we need to store NF(sec.inv.,sP) { tmp = reduce(tmp,sP); ReducedCandidates[j+1..j+MonStep] = tmp[1..MonStep]; tmp = reduce(tmp,SaveRed); } else { tmp = reduce(tmp,TotalReductor); } Indicator[j+1..j+MonStep] = tmp[1..MonStep]; kill tmp; } else { ideal tmp = normalize(evaluate_reynolds(REY,ideal(mon[j+1..ii]))); B[j+1..ii] = tmp[1..ii-j]; if (NoPP == 0) { tmp = reduce(tmp,sP); ReducedCandidates[j+1..ii] = tmp[1..ii-j]; tmp = reduce(tmp,SaveRed); } else { tmp = reduce(tmp,TotalReductor); } Indicator[j+1..ii] = tmp[1..ii-j]; kill tmp; } } 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; if (NoPP==0) // we will still be working with Power Products { 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 irred. sec. inv. number ",irrcounter; } Reductor = ideal(helpP); attrib(Reductor, "isSB",1); Indicator=reduce(Indicator,Reductor); SizeSave++; SaveRed[SizeSave] = helpP; attrib(SaveRed, "isSB",1); attrib(Reductor, "isSB",1); LastNewIS = i-1; if (int(dimmat[i,1])==counter) { break; } } mon[j]=0; B[j]=0; ReducedCandidates[j]=0; Indicator[j]=0; } } } // if i>IrrSwitch else // use fast algorithm { B=sort_of_invariant_basis(TotalReductor,REY,i-1,int(dimmat[i,1])*6); // B is a linearly independent set of reynolds // images of monomials that do not occur // as leading monomials of the ideal spanned // by primary and previously found irreducible // secondary invariants. if (NoPP==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]); 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 ",size(B)," irred. sec. inv.";} } 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 = Indicator[j]; if (helpP <>0) // B[j] should be added { NoPP=0; // i.e., TotalReductor needs to be re-computed 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 irred. sec. inv. number ",irrcounter; } Reductor = ideal(helpP); attrib(Reductor, "isSB",1); Indicator=reduce(Indicator,Reductor); SizeSave++; SaveRed[SizeSave] = helpP; attrib(SaveRed, "isSB",1); attrib(Reductor, "isSB",1); LastNewIS = i-1; if (int(dimmat[i,1])==counter) { break; } } B[j]=0; ReducedCandidates[j]=0; Indicator[j]=0; } } } // i<=IrrSwitch } // if mon[1]<>0 } // Computation of irreducible secondaries if (v) { ""; } } } if (v) { " We're done!"; ""; } degBound = dgb; return(matrix(compress(IS))); } example { "EXAMPLE: S. King"; echo=2; ring r= 0, (a,b,c,d,e,f),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); matrix IS = irred_secondary_char0(L[1],L[2],L[3],0); IS; } /////////////////////////////////////////////////////////////////////////////// proc old_secondary_charp (matrix P, matrix REY, string ring_name, list #) "USAGE: secondary_charp(P,REY,ringname[,v]); P: a 1xn with primary invariants, REY: a gxn representing the Reynolds operator, ringname: a giving the name of a ring of characteristic 0 where the Molien series is stored, v: an optional 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(), `ringname` is a ring of char 0 that has been created by molien() or reynolds_molien() or primary_invariants() RETURN: secondary invariants of the invariant ring (type ) and irreducible secondary invariants (type ) DISPLAY: information if v does not equal 0 THEORY: Secondary invariants are calculated by finding a basis (in terms of monomials) of the basering modulo primary invariants, mapping those to invariants with the Reynolds operator and using these images or their power products such that they are linearly independent modulo the primary invariants (see paper \"Some Algorithms in Invariant Theory of Finite Groups\" by Kemper and Steel (1997)). EXAMPLE: example secondary_charp; shows an example " { def br=basering; degBound=0; //---------------- 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; } 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 should give the name of the ring where the Molien." " series is stored."; return(); } if (v && voice==2) { ""; } int j, m, counter, 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 if (v) { " Polynomial telling us where to look for secondary invariants:"; " "+string(p); ""; } 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 deg_dim_vec=1; for (j=2;j<=m;j++) { deg_dim_vec=deg_dim_vec,int(dimmat[j,1]); } if (v) { "In degree 0 we have: 1"; ""; } //------------------------ initializing variables ---------------------------- setring br; intmat PP; poly pp; int k; intvec deg_vec; ideal sP=groebner(ideal(P)); ideal TEST,B,IS; ideal S=1; // 1 is the first secondary invariant //------------------- 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 0 { // elements in the current degree (i-1) if (v) { "Searching in degree ",i-1,", we need to find ",deg_dim_vec[i]," invariant(s)..."; } TEST=sP; counter=0; // we'll count up to degvec[i] if (IS[1]<>0) { PP=power_products(deg_vec,i-1); // generating power products of } // irreducible secondary invariants if (size(ideal(PP))<>0) { for (j=1;j<=ncols(PP);j++) // going through all of those { pp=1; for (k=1;k<=nrows(PP);k++) { pp=pp*IS[1,k]^PP[k,j]; } if (reduce(pp,TEST)<>0) { S=S,pp; counter++; if (v) { " We find: "+string(pp); } if (deg_dim_vec[i]<>counter) { //TEST=std(TEST+ideal(NF(pp,TEST))); // should be soon replaced by // next line TEST=std(TEST,pp); } else { break; } } } } if (deg_dim_vec[i]<>counter) { 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 and that don't reduce to // 0 modulo sP - if (counter==0 && ncols(B)==deg_dim_vec[i]) // then we can add all of B { S=S,B; IS=IS+B; if (deg_vec[1]==0) { deg_vec=i-1; if (v) { " We find: "+string(B[1]); } for (j=2;j<=deg_dim_vec[i];j++) { deg_vec=deg_vec,i-1; if (v) { " We find: "+string(B[j]); } } } else { for (j=1;j<=deg_dim_vec[i];j++) { deg_vec=deg_vec,i-1; if (v) { " We find: "+string(B[j]); } } } } else { j=0; // j goes through all of B - while (deg_dim_vec[i]<>counter) // need to find deg_dim_vec[i] { // invariants that are linearly // independent modulo TEST j++; if (reduce(B[j],TEST)<>0) // B[j] should be added { S=S,B[j]; IS=IS+ideal(B[j]); if (deg_vec[1]==0) { deg_vec[1]=i-1; } else { deg_vec=deg_vec,i-1; } counter++; if (v) { " We find: "+string(B[j]); } if (deg_dim_vec[i]<>counter) { //TEST=std(TEST+ideal(NF(B[j],TEST))); // should be soon replaced // by next line TEST=std(TEST,B[j]); } } } } } if (v) { ""; } } } if (v) { " We're done!"; ""; } if (ring_name=="aksldfalkdsflkj") { kill `ring_name`; } return(matrix(S),matrix(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)]); print(S); print(IS); } /////////////////////////////////////////////////////////////////////////////// proc secondary_charp (matrix P, matrix REY, string ring_name, list #) "USAGE: secondary_charp(P,REY,ringname[,v][,\"old\"]); @* P: a 1xn with homogeneous primary invariants, where n is the number of variables of the basering; @* REY: a gxn representing the Reynolds operator, where g the size of the corresponding group; @* ringname: a giving the name of a ring of characteristic 0 containing a 1x2 M giving numerator and denominator of the Molien series; @* v: an optional ; @* \"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 ) and irreducible secondary invariants (type ) 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 secondary invariants using Groebner basis techniques (see S. King: Fast Computation of Secondary Invariants). 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 = 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. 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 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 ideal sP_Reductor; // will contain sP union Reductor. int sPOffset = ncols(sP); 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 int NrIS; 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;i0) // when it is == 0 we need to find no { // elements in the current degree (i-1) if (v) { "Searching in degree ",i-1,", we need to find ",deg_dim_vec[i], " invariant(s)..."; " Looking for Power Products..."; } counter = 0; // we'll count up to deg_dim_vec[i] Reductor = ideal(0); sP_Reductor = sP; 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"none") { Mul1 = ISSort[k]; } else { Mul1 = ideal(0); } if ((deg_dim_vec[i-k]>0) && (Mul1[1]<>0)) { for (minD=k;minD=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]; } } ProdCand = simplify(Mul1*Mul2,4); 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); 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. Hence, sP_Reductor // is a homog. GB up to degree i-1. SaveRed[SizeSave] = helpP; attrib(SaveRed, "isSB",1); } else { break; } } } } } } } } NrIS = deg_dim_vec[i]-counter; // The remaining NrIS sec. inv. are irreducible! if (NrIS>0) // need more than all the power products { if (v) { " There are ",NrIS,"irreducible secondary invariants in degree ", i-1; } mon = kbase(sP_Reductor,i-1); ii = ncols(mon); if ((i>IrrSwitch) or (ii>200)) // use sparse algorithm { 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 all",NrIS,"irreducibles in degree ",i-1;} } 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) <= ii) { 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); Indicator[j+1..j+MonStep] = tmp[1..MonStep]; kill tmp; } else { 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); Indicator[j+1..ii] = tmp[1..ii-j]; kill tmp; } } 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", irrcounter,"of", NrIS ,"irr. sec. inv. in degree ",i-1; } 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 all",NrIS,"irreducibles in degree ",i-1;} } 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", irrcounter,"of", NrIS ,"irr. sec. inv. in degree ",i-1; } 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 with primary invariants, REY: a gxn representing the Reynolds operator, deg_vec: an optional listing some degrees where no non-trivial homogeneous invariants can be found, v: an optional 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(), deg_vec is the second return value of primary_char0_no_molien(), primary_charp_no_molien(), primary_char0_no_molien_random() or primary_charp_no_molien_random() RETURN: secondary invariants of the invariant ring (type ) DISPLAY: information if v does not equal 0 THEORY: Secondary invariants are calculated by finding a basis (in terms of monomials) of the basering modulo primary invariants, mapping those to invariants with the Reynolds operator and using these images as 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. NOTE: should usually be faster and of more useful functionality. SEE ALSO: secondary_and_irreducibles_no_molien EXAMPLE: example secondary_no_molien; shows an example " { 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 i; 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 "; 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 "; 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(); } } else { "WARNING: In char 0, secondary_char0 should be faster."; } if (v && voice==2) { ""; } int j, m, d; int max=1; for (j=1;j<=n;j++) { max=max*deg(P[j]); } max=max/nrows(REY); if (v) { " We need to find ",max," secondary invariants."; ""; "In degree 0 we have: 1"; ""; } //------------------------- initializing variables --------------------------- int dgb = degBound; degBound=0; ideal sP=groebner(ideal(P)); ideal Reductor,SaveRed; // sP union Reductor is a Groebner basis up to degree i-1 int SizeSave; poly helpP; ideal helpI; ideal Indicator; // will tell us which candidates for sec. inv. we can choose ideal ReducedCandidates; int helpint; int k2,k3,minD; int ii; int saveAttr; ideal mon,B; ideal S=1; // 1 is the first secondary invariant int counter=1; i=0; if (defined(deg_vec)<>voice) { intvec deg_vec; } int k=1; //--------------------- generating secondary invariants ---------------------- while (counter<>max) { i++; Reductor = ideal(0); helpint = 0; SaveRed = Reductor; SizeSave = 0; attrib(Reductor,"isSB",1); attrib(SaveRed,"isSB",1); if (deg_vec[k]<>i) { if (v) { "Searching in degree ",i,"..."; } mon = kbase(sP,i); ii = ncols(mon); if ((i>IrrSwitch) or (ii>200)) // use sparse algorithm { j=0; // j goes through all of mon - // Compare the comments on the computation of reducible sec. inv.! while ((j0) // B[j] should be added { counter++; S[counter] = B[j]; if (v) { " We found sec. inv. number ",counter," in degree ",i; } if (counter == max) {break;} 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,max); // B contains images of // kbase(sP,i) under the Reynolds // operator that are linearly independent // and that don't reduce to 0 modulo sP 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); ii = ncols(Indicator); while ((j0) // B[j] should be added { counter++; S[counter] = B[j]; if (v) { " We found sec. inv. number ",counter," in degree ",i; } if (counter == max) {break;} 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 } else // if deg_vec[k]=i: { if (size(deg_vec)==k) { k=1; } else { k++; } } } if (v) { ""; } if (v) { " We're done!"; ""; } degBound = dgb; return(matrix(S)); } 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 S=secondary_no_molien(L[1..2],intvec(1,2),1); print(S); } /////////////////////////////////////////////////////////////////////////////// proc secondary_and_irreducibles_no_molien (matrix P, matrix REY, list #) "USAGE: secondary_and_irreducibles_no_molien(P,REY[,deg_vec,v]); P: a 1xn with primary invariants, REY: a gxn representing the Reynolds operator, deg_vec: an optional listing some degrees where no non-trivial homogeneous invariants can be found, v: an optional 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: secondary invariants of the invariant ring (type ) and irreducible secondary invariants (type ) DISPLAY: information if v does not equal 0 THEORY: Secondary invariants are calculated by finding a basis (in terms of monomials) of the basering modulo primary invariants, mapping those to invariants with the Reynolds operator. Among these images or their power products we pick secondary invariants using Groebner basis techniques (see S. King: Fast Computation of 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. SEE ALSO: secondary_no_molien EXAMPLE: example secondary_and_irreducibles_no_molien; shows an example " { 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 i; 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 "; 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 "; 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(); } } else { "WARNING: In char 0, secondary_char0 should be faster."; } if (v && voice==2) { ""; } int j, m, d; int max=1; for (j=1;j<=n;j++) { max=max*deg(P[j]); } max=max/nrows(REY); if (v) { " We need to find ",max," secondary invariants."; ""; "In degree 0 we have: 1"; ""; } //-------------------------- initializing variables -------------------------- 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. int counter=1; int irrcounter; i=1; // the degree under consideration is i-1 SSort[i] = list(); for (k=1;k<=i;k++) { SSort[i][k]=ideal(); attrib(SSort[i][k],"size",0); } if (defined(deg_vec)<>voice) { intvec deg_vec; } int l=1; //------------------- generating secondary invariants ------------------------ while (counter<>max) { i++; SSort[i] = list(); for (k=1;k<=i;k++) { SSort[i][k]=ideal(); attrib(SSort[i][k],"size",0); } if (deg_vec[l]<>i-1) { if (v) { "Searching in degree ",i-1; " Looking for Power Products..."; } 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"none") { Mul1 = ISSort[k]; } else { Mul1 = ideal(0); } if (Mul1[1]<>0) { for (minD=k;minD=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]; } } ProdCand = simplify(Mul1*Mul2,4); 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 sec. inv. number",counter, "in degree",i-1; } if (max<>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; } } } } } } } } // The remaining sec. inv. are irreducible! if (max<>counter) { if (v) { " Looking for irreducible secondary invariants in degree ", i-1; } mon = kbase(sP,i-1); ii = ncols(mon); if ((i>IrrSwitch) or (ii>200)) // use sparse algorithm { irrcounter=0; j=0; // j goes through all of mon - // Compare the comments on the computation of reducible sec. inv.! while ((jcounter)) { if ((j mod MonStep) == 0) { if ((j+MonStep) <= ii) { 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); Indicator[j+1..j+MonStep] = tmp[1..MonStep]; kill tmp; } else { 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); Indicator[j+1..ii] = tmp[1..ii-j]; kill tmp; } } 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 irreducible sec. inv. number", irrcounter, "in degree", i-1; } if (counter == max) { break; } 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,max); // B contains // images of kbase(sP,i-1) under the // Reynolds operator that are linearly // independent 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); ii = ncols(Indicator); while ((jcounter)) { 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 irreducible sec. inv. number", irrcounter, "in degree", i-1; } Reductor = ideal(helpP); attrib(Reductor, "isSB",1); Indicator=reduce(Indicator,Reductor); SizeSave++; SaveRed[SizeSave] = helpP; attrib(SaveRed, "isSB",1); } if (max == counter) { break;} B[j]=0; ReducedCandidates[j]=0; Indicator[j]=0; } } // i<=IrrSwitch } // Computation of irreducible secondaries if (v) { ""; } } // if (deg_vec[l]<>i-1) else { if (size(deg_vec)==l) { l=1; } else { l++; } } } // while counter <> max if (v) { ""; } if (v) { " We're done!"; ""; } degBound = dgb; // Prepare return: matrix S[1][max]; S[1,1]=1; j=1; k=1; while (j<>max) { 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(); } } k++; } return(S,matrix(compress(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 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 with primary invariants, REY: a gxn representing the Reynolds operator, deg_vec: an optional listing some degrees where no irreducible secondary invariants can be found, v: an optional 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 ) 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 secondary invariants, using Groebner basis techniques. SEE ALSO: irred_secondary_char0 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 "; 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 "; 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 first parameter ought to be a matrix of primary invariants."); } 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 ",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 ",monsize," candidates for irred. secondaries";} block=0; // Loops through the monomials while (block0) { counter++; // found new irr. sec. inv.! if (v) { " We found irr. sec. inv. number ",counter," in degree ",d; } IS[counter] = Inv[block]; sIS = sIS,helpP; attrib(sIS,"isSB",1); NewIScounter++; NewIS[NewIScounter] = helpP; 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); } /////////////////////////////////////////////////////////////////////////////// proc secondary_not_cohen_macaulay (matrix P, list #) "USAGE: secondary_not_cohen_macaulay(P,G1,G2,...[,v]); P: a 1xn with primary invariants, G1,G2,...: nxn generating a finite matrix group, v: optional ASSUME: n is the number of variables of the basering RETURN: secondary invariants of the invariant ring (type ) DISPLAY: information on the progress of computation if v does not equal 0 THEORY: Secondary invariants are generated following \"Generating Invariant Rings of Finite Groups over Arbitrary Fields\" by Kemper (1996). EXAMPLE: example secondary_not_cohen_macaulay; shows an example " { int i, j; int v; degBound=0; def br=basering; 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 (size(#)>0) // checking input and setting verbose { if (typeof(#[size(#)])=="int") { int gen_num=size(#)-1; v=#[size(#)]; if (gen_num==0) { "ERROR: There are no generators of the finite matrix group given."; return(); } for (i=1;i<=gen_num;i++) { if (typeof(#[i])<>"matrix") { "ERROR: These parameters should be generators of the finite matrix group."; return(); } if ((n<>nrows(#[i])) or (n<>ncols(#[i]))) { "ERROR: matrices need to be square and of the same dimensions"; return(); } } } else { int gen_num=size(#); for (i=1;i<=gen_num;i++) { if (typeof(#[i])<>"matrix") { "ERROR: These parameters should be generators of the finite matrix group."; return(); } if ((n<>nrows(#[i])) or (n<>ncols(#[i]))) { "ERROR: matrices need to be square and of the same dimensions"; return(); } } } } else { "ERROR: There are no generators of the finite matrix group given."; return(); } if (ncols(P)<>n) { "ERROR: The first parameter ought to be the matrix of the primary"; " invariants." return(); } if (v && voice==2) { ""; } // 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 // Molien series again - M[1,1]/M[1,2] is // the Molien series of that group, we // now calculate the secondary invariants // of this subgroup in the usual fashion // where the primary invariants are the // ones from the bigger group if (v) { " 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."; ""; } 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 // invariants, we just calculated if (v) { " We calculate secondary invariants from the ones found for the trivial"; " subgroup."; ""; } map f; // used to let generators act on // secondary invariants with respect to // the trivial group - matrix M(1)[gen_num][k]; // M(1) will contain a module ideal B; for (i=1;i<=gen_num;i++) { B=ideal(matrix(maxideal(1))*transpose(#[i])); // image of the various // variables under the i-th generator - f=br,B; // the corresponding mapping - B=f(trivialS)-trivialS; // these relations should be 0 - 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); 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 // M(2) - // the following steps calculates the intersection of the module M(2) with // the algebra A^k where A denote the subalgebra of the usual polynomial // ring, generated by the primary invariants // 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]=GElimP[i]; // y(i)-P[1,i]; } } M=elim2(module(M),1..n); // eliminating x(1..n), std-calculation // is done internally - //M=std(module(M)); // we have already an elimination ordering //M=nselect(M,1..n); M=homog(module(M),h); // homogenize for 'minbase' M=minbase(module(M)); 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=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(3)); // S now contains the secondary // invariants for (i=1; i<=m;i++) { S[1,i]=S[1,i]/leadcoef(S[1,i]); // making elements nice } S=sort(ideal(S))[1]; if (v) { " These are the secondary invariants: "; for (i=1;i<=m;i++) { " "+string(S[1,i]); } ""; " We're done!"; ""; } if ((v or (voice==2)) && (m>1)) { " WARNING: The invariant ring might not have a Hironaka decomposition"; " if the characteristic of the coefficient field divides the"; " group order."; } return(matrix(S)); } example { "EXAMPLE:"; echo=2; ring R=2,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; list L=primary_invariants(A); matrix S=secondary_not_cohen_macaulay(L[1],A); print(S); } /////////////////////////////////////////////////////////////////////////////// proc invariant_ring (list #) "USAGE: invariant_ring(G1,G2,...[,flags]); G1,G2,...: generating a finite matrix group, flags: an optional with three entries: if the first one equals 0, the program attempts to compute the Molien series and Reynolds operator, if it equals 1, the program is told that the Molien series should not be computed, if it equals -1 characteristic 0 is simulated, i.e. the Molien series is computed as if the base field were characteristic 0 (the user must choose a field of large prime characteristic, e.g. 32003) and if the first one is anything else, it 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 series), the second component should give the size of intervals between canceling common factors in the expansion of Molien series, 0 (the default) means only once after generating all terms, in prime characteristic also a negative number can be given to indicate that common factors should always be canceled when the expansion is simple (the root of the extension field occurs not among the coefficients) RETURN: primary and secondary invariants (both of type ) generating 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 THEORY: Bases of homogeneous invariants are generated successively and those are chosen as primary invariants that lower the dimension of the ideal generated by the previously found invariants (see \"Generating a Noetherian Normalization of the Invariant Ring of a Finite Group\" by Decker, Heydtmann, Schreyer (1998)). In the non-modular case secondary invariants are calculated by finding a basis (in terms of monomials) of the basering modulo the primary invariants, mapping to invariants with the Reynolds operator and using those or their power products such that they are linearly independent modulo the primary invariants (see \"Some Algorithms in Invariant Theory of Finite Groups\" by Kemper and Steel (1997)). In the modular case they are generated according to \"Generating Invariant Rings of Finite Groups over Arbitrary Fields\" by Kemper (1996). EXAMPLE: example invariant_ring; shows an example " { if (size(#)==0) { "ERROR: There are no generators given."; return(); } int ch=char(basering); // the algorithms depend very much on the // characteristic of the ground field - 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 int gen_num; int mol_flag, v; //------------------- checking input and setting flags ----------------------- if (typeof(#[size(#)])=="intvec") { if (size(#[size(#)])<>3) { "ERROR: The should have three entries."; return(); } gen_num=size(#)-1; mol_flag=#[size(#)][1]; if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0))) { "ERROR: the second component of should be >=0"; return(); } int interval=#[size(#)][2]; v=#[size(#)][3]; } else { gen_num=size(#); mol_flag=0; int interval=0; v=0; } //---------------------------------------------------------------------------- if (mol_flag==0) // calculation Molien series will be { if (ch==0) // attempted - { matrix REY,M=reynolds_molien(#[1..gen_num],intvec(0,interval,v)); // one // will contain Reynolds operator and the // other enumerator and denominator of // Molien series matrix P=primary_char0(REY,M,v); matrix S,IS=secondary_char0(P,REY,M,v); return(P,S,IS); } else { list L=group_reynolds(#[1..gen_num],v); if (L[1]<>0) // testing whether we are in the modular { string newring="aksldfalkdsflkj"; // case if (minpoly==0) { if (v) { " We are dealing with the non-modular case."; } if (typeof(L[2])=="int") { molien(L[3..size(L)],newring,L[2],intvec(0,interval,v)); } else { molien(L[2..size(L)],newring,intvec(0,interval,v)); } matrix P=primary_charp(L[1],newring,v); matrix S,IS=secondary_charp(P,L[1],newring,v); if (defined(aksldfalkdsflkj)==2) { kill aksldfalkdsflkj; } return(P,S,IS); } else { if (v) { " 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."; ""; } list l=primary_charp_no_molien(L[1],v); if (size(l)==2) { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v); } else { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v); } return(l[1],S,IS); } } else // the modular case { if (v) { " There is also no Molien series or Reynolds operator, we can make use of..."; ""; " We can start looking for primary invariants..."; ""; } matrix P=primary_charp_without(#[1..gen_num],v); matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v); return(P,S); } } } if (mol_flag==1) // the user wants no calculation of the { list L=group_reynolds(#[1..gen_num],v); // Molien series if (ch==0) { list l=primary_char0_no_molien(L[1],v); if (size(l)==2) { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v); } else { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v); } return(l[1],S,IS); } else { if (L[1]<>0) // testing whether we are in the modular { list l=primary_charp_no_molien(L[1],v); // case if (size(l)==2) { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v); } else { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v); } return(l[1],S,IS); } else // the modular case { if (v) { " We can start looking for primary invariants..."; ""; } matrix P=primary_charp_without(#[1..gen_num],v); matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v); return(P,S); } } } if (mol_flag==-1) { if (ch==0) { "ERROR: Characteristic 0 can only be simulated in characteristic p>>0. "; return(); } list L=group_reynolds(#[1..gen_num],v); string newring="aksldfalkdsflkj"; if (typeof(L[2])=="int") { molien(L[3..size(L)],newring,L[2],intvec(1,interval,v)); } else { molien(L[2..size(L)],newring,intvec(1,interval,v)); } matrix P=primary_charp(L[1],newring,v); matrix S,IS=secondary_charp(P,L[1],newring,v); kill aksldfalkdsflkj; return(P,S,IS); } else // the user specified that the { if (ch==0) // characteristic divides the group order { "ERROR: The characteristic cannot divide the group order when it is 0. "; return(); } if (v) { ""; } matrix P=primary_charp_without(#[1..gen_num],v); matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v); return(L[1],S); } } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2; ring R=0,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; matrix P,S,IS=invariant_ring(A); print(P); print(S); print(IS); } /////////////////////////////////////////////////////////////////////////////// proc invariant_ring_random (list #) "USAGE: invariant_ring_random(G1,G2,...,r[,flags]); G1,G2,...: generating a finite matrix group, r: an where -|r| to |r| is the range of coefficients of random combinations of bases elements that serve as primary invariants, flags: an optional with three entries: if the first equals 0, the program attempts to compute the Molien series and Reynolds operator, if it equals 1, the program is told that the Molien series should not be computed, if it equals -1 characteristic 0 is simulated, i.e. the Molien series is computed as if the base field were characteristic 0 (the user must choose a field of large prime characteristic, e.g. 32003) and if the first one is anything else, then the characteristic of the base field divides the group order (i.e. we will not even attempt to compute the Reynolds operator or Molien series), the second component should give the size of intervals between canceling common factors in the expansion of the Molien series, 0 (the default) means only once after generating all terms, in prime characteristic also a negative number can be given to indicate that common factors should always be canceled when the expansion is simple (the root of the extension field does not occur among the coefficients) RETURN: primary and secondary invariants (both of type ) generating 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 THEORY: is the same as for invariant_ring except that random combinations of basis elements are chosen as candidates for primary invariants and hopefully they lower the dimension of the previously found primary invariants by the right amount. EXAMPLE: example invariant_ring_random; shows an example " { if (size(#)<2) { "ERROR: There are too few parameters."; return(); } int ch=char(basering); // the algorithms depend very much on the // characteristic of the ground field 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 int gen_num; int mol_flag, v; //------------------- checking input and setting flags ----------------------- if (typeof(#[size(#)])=="intvec" && typeof(#[size(#)-1])=="int") { if (size(#[size(#)])<>3) { "ERROR: should have three entries."; return(); } gen_num=size(#)-2; mol_flag=#[size(#)][1]; if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0))) { "ERROR: the second component of should be >=0"; return(); } int interval=#[size(#)][2]; v=#[size(#)][3]; int max=#[size(#)-1]; if (gen_num==0) { "ERROR: There are no generators of a finite matrix group given."; return(); } } else { if (typeof(#[size(#)])=="int") { gen_num=size(#)-1; mol_flag=0; int interval=0; v=0; int max=#[size(#)]; } else { "ERROR: If the two last parameters are not and , the last"; " parameter should be an ."; return(); } } for (int i=1;i<=gen_num;i++) { if (typeof(#[i])=="matrix") { if (nrows(#[i])<>n or ncols(#[i])<>n) { "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 first parameters should be a list of matrices"; return(); } } //---------------------------------------------------------------------------- if (mol_flag==0) { if (ch==0) { matrix REY,M=reynolds_molien(#[1..gen_num],intvec(0,interval,v)); // one // will contain Reynolds operator and the // other enumerator and denominator of // Molien series matrix P=primary_char0_random(REY,M,max,v); matrix S,IS=secondary_char0(P,REY,M,v); return(P,S,IS); } else { list L=group_reynolds(#[1..gen_num],v); if (L[1]<>0) // testing whether we are in the modular { string newring="aksldfalkdsflkj"; // case if (minpoly==0) { if (v) { " We are dealing with the non-modular case."; } if (typeof(L[2])=="int") { molien(L[3..size(L)],newring,L[2],intvec(0,interval,v)); } else { molien(L[2..size(L)],newring,intvec(0,interval,v)); } matrix P=primary_charp_random(L[1],newring,max,v); matrix S,IS=secondary_charp(P,L[1],newring,v); if (voice==2) { kill aksldfalkdsflkj; } return(P,S,IS); } else { if (v) { " 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."; ""; } list l=primary_charp_no_molien_random(L[1],max,v); if (size(l)==2) { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v); } else { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v); } return(l[1],S,IS); } } else // the modular case { if (v) { " There is also no Molien series, we can make use of..."; ""; " We can start looking for primary invariants..."; ""; } matrix P=primary_charp_without_random(#[1..gen_num],max,v); matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v); return(P,S); } } } if (mol_flag==1) // the user wants no calculation of the { list L=group_reynolds(#[1..gen_num],v); // Molien series if (ch==0) { list l=primary_char0_no_molien_random(L[1],max,v); if (size(l)==2) { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v); } else { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v); } return(l[1],S,IS); } else { if (L[1]<>0) // testing whether we are in the modular { list l=primary_charp_no_molien_random(L[1],max,v); // case if (size(l)==2) { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],l[2],v); } else { matrix S,IS=secondary_and_irreducibles_no_molien(l[1],L[1],v); } return(l[1],S,IS); } else // the modular case { if (v) { " We can start looking for primary invariants..."; ""; } matrix P=primary_charp_without_random(#[1..gen_num],max,v); matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v); return(L[1],S); } } } if (mol_flag==-1) { if (ch==0) { "ERROR: Characteristic 0 can only be simulated in characteristic p>>0. "; return(); } list L=group_reynolds(#[1..gen_num],v); string newring="aksldfalkdsflkj"; if (typeof(L[2])=="int") { molien(L[3..size(L)],newring,L[2],intvec(mol_flag,interval,v)); } else { molien(L[2..size(L)],newring,intvec(mol_flag,interval,v)); } matrix P=primary_charp_random(L[1],newring,max,v); matrix S,IS=secondary_charp(P,L[1],newring,v); kill aksldfalkdsflkj; return(P,S,IS); } else // the user specified that the { if (ch==0) // characteristic divides the group order { "ERROR: The characteristic cannot divide the group order when it is 0. "; return(); } if (v) { ""; } matrix P=primary_charp_without_random(#[1..gen_num],max,v); matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v); return(L[1],S); } } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2; ring R=0,(x,y,z),dp; matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; matrix P,S,IS=invariant_ring_random(A,1); print(P); print(S); print(IS); } /////////////////////////////////////////////////////////////////////////////// proc orbit_variety (matrix F,string newring) "USAGE: orbit_variety(F,s); F: a 1xm defing an invariant ring, s: a giving the name for a new ring RETURN: a Groebner basis (type , named G) for the ideal defining the orbit variety (i.e. the syzygy ideal) in the new ring (named `s`) 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 EXAMPLE: example orbit_variety; shows an example " { if (newring=="") { "ERROR: the second parameter may not be an empty "; return(); } if (nrows(F)==1) { def br=basering; int n=nvars(br); int m=ncols(F); string mp=string(minpoly); execute("ring R=("+charstr(br)+"),("+varstr(br)+",y(1..m)),dp;"); execute("minpoly=number("+mp+");"); ideal I=ideal(imap(br,F)); for (int i=1;i<=m;i++) { I[i]=I[i]-y(i); } I=elim(I,1..n); execute("ring "+newring+"=("+charstr(br)+"),(y(1..m)),dp(m);"); execute("minpoly=number("+mp+");"); ideal vars; for (i=2;i<=n;i++) { vars[i]=0; } vars=vars,y(1..m); map emb=R,vars; ideal G=emb(I); kill emb, vars, R; keepring `newring`; return(); } else { "ERROR: the may only have one row"; return(); } } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2; ring R=0,(x,y,z),dp; matrix F[1][7]=x2+y2,z2,x4+y4,1,x2z-1y2z,xyz,x3y-1xy3; string newring="E"; orbit_variety(F,newring); print(G); basering; } /////////////////////////////////////////////////////////////////////////////// static proc degvec(ideal I) "USAGE: degvec(I); I an . RETURN: the of degrees of the generators of I. " { intvec v; for (int j = 1;j<=ncols(I);j++) { v[j]=deg(I[j]); } return(v); } /////////////////////////////////////////////////////////////////////////////// proc rel_orbit_variety(ideal I,matrix F,list #) "USAGE: rel_orbit_variety(I,F[,s]); @* I: an invariant under the action of a group, @* F: a 1xm defining the invariant ring of this group. @* s: optional ; if s is present then (for downward compatibility) the old procedure is called, and in this case s gives the name of a new . RETURN: Without optional string s, a list L of two rings is returned. @* The ring L[1] carries a weighted degree order with variables y(1..m), the weight of y(k) equal to the degree of the k-th generators F[1,k] of the invariant ring. L[1] contains a Groebner basis (type , named G) of the ideal defining the relative orbit variety with respect to I. @* The ring L[2] has the variables of the basering together with y(1..m) and carries a block order: The first block is the order of the basering, the second is the weighted degree order occuring in L[1]. L[2] contains G and a Groebner basis (type , named Conv) such that if p is any invariant polynomial expressed in the variables of the basering then reduce(p,Conv) is a polynomial in the new variables y(1..m) such that evaluation at the generators of the invariant ring yields p. This can be used to avoid the application of (see @ref{algebra_containment}). @* For the case of optional string s, the function is equivalent to @ref{relative_orbit_variety}. 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 ideal generators. The variables of the original ring are eliminated and the polynomials that are left define the relative orbit variety with respect to I. The elimination is done by a weighted blockorder that has the advantage of dealing with quasi-homogeneous ideals. NOTE: We provide the ring L[1] for the sake of downward compatibility, since it is closer to the ring returned by relative_orbit_variety than L[2]. However, L[1] carries a weighted degree order, whereas the ring returned by relative_orbit_variety is lexicographically ordered. SEE ALSO: relative_orbit_variety EXAMPLE: example rel_orbit_variety; shows an example. " { if (size(#)>0) { if (typeof(#[1])=="string") { relative_orbit_variety(I,F,#[1]); keepring basering; return(); } else { "ERROR: the third parameter may either be empty or a ."; return(); } } degBound=0; if (nrows(F)==1) { option(redSB,noredTail); def br=basering; int n=nvars(br); int m=ncols(F); // In the following ring definition, any elimination order would work. list rlist = ringlist(br); int i; for (i=1;i<=m;i++) { rlist[2][n+i]="y("+string(i)+")"; } rlist[3][size(rlist[3])+1] = rlist[3][size(rlist[3])]; rlist[3][size(rlist[3])-1][1] = "wp"; rlist[3][size(rlist[3])-1][2] = degvec(ideal(F)); def newring = ring(rlist); list smallrlist; smallrlist[1]=rlist[1]; smallrlist[2]=list(); for (i=1;i<=m;i++) { smallrlist[2][i]="y("+string(i)+")"; } smallrlist[3]=list(); smallrlist[3][1]=list(); smallrlist[3][1][1] = "wp"; smallrlist[3][1][2] = degvec(ideal(F)); smallrlist[3][2] = rlist[3][size(rlist[3])]; smallrlist[4] = rlist[4]; def smallring = ring(smallrlist); setring(newring); ideal J=ideal(imap(br,F)); ideal I=imap(br,I); for (i=1;i<=m;i++) { J[i]=J[i]-y(i); } // We chose the weighted block order since this makes J quasi-homogeneous! ideal Conv=groebner(J); J=Conv,I; J=groebner(J); ideal vars; vars[n]=0; vars=vars,y(1..m); map emb=newring,vars; ideal G=emb(J); J=J-G; for (i=1;i<=ncols(G);i++) { if (J[i]<>0) { G[i]=0; } } G=compress(G); vars=ideal(); for (i=2;i<=n;i++) { vars[i]=0; } vars=vars,y(1..m); emb=newring,vars; G=compress(emb(G)); export(G); export(Conv); setring(smallring); ideal G = compress(imap(newring,G)); export(G); list L; L[1]=smallring; L[2]=newring; dbprint( printlevel-voice+3," // 'rel_orbit_variety' created a list of two rings. // If L is the name of that list, you can access // the first ring by // def R = L[1]; setring R; // (similarly for the second ring)"); return(L); } else { "ERROR: the may only have one row"; return(); } } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.6.3:"; echo=2; ring R=0,(x,y,z),dp; matrix F[1][3]=x+y+z,xy+xz+yz,xyz; ideal I=x2+y2+z2-1,x2y+y2z+z2x-2x-2y-2z,xy2+yz2+zx2-2x-2y-2z; list L = rel_orbit_variety(I,F); def AllR = L[2]; setring(AllR); print(G); print(Conv); basering; } /////////////////////////////////////////////////////////////////////////////// proc relative_orbit_variety(ideal I,matrix F,string newring) "USAGE: relative_orbit_variety(I,F,s); I: an invariant under the action of a group, @* F: a 1xm defining the invariant ring of this group, @* s: a giving a name for a new ring RETURN: The procedure ends with a new ring named s. It contains a Groebner basis (type , 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 ideal generators. The variables of the original ring are eliminated and the polynomials that are left define the relative orbit variety with respect to I. NOTE: This procedure is now replaced by rel_orbit_variety (see @ref{rel_orbit_variety}), which uses a different elemination order that should usually allow faster computations. SEE ALSO: rel_orbit_variety EXAMPLE: example relative_orbit_variety; shows an example " { if (newring=="") { "ERROR: the third parameter may not be empty a "; return(); } degBound=0; if (nrows(F)==1) { def br=basering; int n=nvars(br); int m=ncols(F); string mp=string(minpoly); execute("ring R=("+charstr(br)+"),("+varstr(br)+",y(1..m)),lp;"); execute("minpoly=number("+mp+");"); ideal J=ideal(imap(br,F)); ideal I=imap(br,I); for (int i=1;i<=m;i++) { J[i]=J[i]-y(i); } J=std(J); J=J,I; option(redSB); J=std(J); ideal vars; //for (i=1;i<=n;i=i+1) //{ vars[i]=0; //} vars[n]=0; vars=vars,y(1..m); map emb=R,vars; ideal G=emb(J); J=J-G; for (i=1;i<=ncols(G);i++) { if (J[i]<>0) { G[i]=0; } } G=compress(G); execute("ring "+newring+"=("+charstr(br)+"),(y(1..m)),lp;"); execute("minpoly=number("+mp+");"); ideal vars; for (i=2;i<=n;i++) { vars[i]=0; } vars=vars,y(1..m); map emb=R,vars; ideal G=emb(G); kill vars, emb; keepring `newring`; return(); } else { "ERROR: the may only have one row"; return(); } } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.6.3:"; echo=2; ring R=0,(x,y,z),dp; matrix F[1][3]=x+y+z,xy+xz+yz,xyz; ideal I=x2+y2+z2-1,x2y+y2z+z2x-2x-2y-2z,xy2+yz2+zx2-2x-2y-2z; string newring="E"; relative_orbit_variety(I,F,newring); print(G); basering; } /////////////////////////////////////////////////////////////////////////////// proc image_of_variety(ideal I,matrix F) "USAGE: image_of_variety(I,F); @* I: an arbitray , @* F: a 1xm defining an invariant ring of some matrix group RETURN: The defining the image under that group of the variety defined 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 EXAMPLE: example image_of_variety; shows an example " { if (nrows(F)==1) { def br=basering; int n=nvars(br); list L = rel_orbit_variety(I,F); def newring=L[2]; setring(newring); ideal F=imap(br,F); for (int i=1;i<=n;i++) { F=0,F; } map emb2=newring,F; ideal trafo = compress(emb2(G)); setring br; return(imap(newring,trafo)); } else { "ERROR: the may only have one row"; return(); } } example { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.6.8:"; echo=2; ring R=0,(x,y,z),dp; matrix F[1][3]=x+y+z,xy+xz+yz,xyz; ideal I=xy; print(image_of_variety(I,F)); } /////////////////////////////////////////////////////////////////////////////// static proc GetMatrix(list GEN) { matrix M[nvars(basering)][nvars(basering)]; int i,j; for (i=1;i<=size(GEN);i++) { for (j=1;j and is formed by monomials of degree d. Per is a list of ideals providing a permutation of the ring variables. ASSUME: mon contains at least the minimal monomial in each orbit sum of degree d under the permutation group given by Per. RETURN: an whose generators are the orbit sums of degree d under the permutation group given by Per. SEE ALSO: invariant_algebra_perm KEYWORDS: orbit sum EXAMPLE: example orbit_sums; shows an example " { //----------------- checking input ------------------ mon=simplify(mon,2); // remove 0 if (not homog(sum(mon))) { "ERROR: Monomials in the first parameter must be of the same degree."; return(); } int k,k2; for (k=1;k<=ncols(mon);k++) { if (leadmonom(mon[k])!=mon[k]) { "ERROR: The first parameter must be formed by monomials."; return(); } } for (k2=1;k2<=size(M);k2++) { if (typeof(M[k2])!="ideal") { "ERROR: The second parameter must be a list of ideals providing a "; " permutation of the ring variables."; return(); } if ((ncols(M[k2])!=nvars(basering)) or (ncols(M[k2])!=ncols(simplify(M[k2],6))) or (M[k2][1]==0)) { "ERROR: The second parameter must be a list of ideals providing a "; " permutation of the ring variables."; return(); } for (k=1;k<=ncols(M[k2]);k++) { if ((leadmonom(M[k2][k])!=M[k2][k]) or (deg(M[k2][k])!=1)) { "ERROR: The second parameter must be a list of ideals providing a "; " permutation of the ring variables."; return(); } } } //--------------- preparing steps ------------------- int dgb = degBound; degBound = 0; def br = basering; ideal tmpI; int g; // rings that allow fast permutation of variables with >fetch< and >imap< for (g=1;g<=size(M);g++) { list tmpL(g); execute("ring tmpR("+string(g)+") = "+charstr(br)+", ("+string(M[g])+"),("+ordstr(br)+")"); ideal tmpI; poly OrbitS; list L(g); setring(br); } // Remove non-minimal monomials. // Strategy: If m1, m2 are monomials related by a group generator // then remember only the smaller of the two. // This is where we use the assumption on mon! for (g=1;g<=size(M);g++) { setring tmpR(g); tmpI=fetch(br,mon); setring br; tmpI=imap(tmpR(g),tmpI); for (k=1;k<=ncols(mon);k++) { if ( mon[k]!=0) { if (tmpI[k] < mon[k]) { mon[k]=0; } else { if (tmpI[k]>mon[k]) { mon=NF(mon,std(tmpI[k])); } } } } } mon=simplify(mon,2); // main part: form orbits list L,newL; // contains partial orbits for (k=1;k<=ncols(mon);k++) { L[k]=ideal(mon[k]); attrib(L[k],"isSB",1); } newL=L; poly tmp; int GenOK, j,i; poly OrbitS=mon[1]; mon[1]=0; int countOrbits; ideal AllOrbits; while (1) // apply all generators to the monomials that have been found // in the preceding step { GenOK = 0; for (g=1;g<=size(M);g++) // Generator g is applied to partial orbits { setring tmpR(g); L(g)=fetch(br,newL); setring br; tmpL(g)=imap(tmpR(g),L(g)); } for (k=1;k<=size(L);k++) { if (L[k][1]!=0) // if the orbit was not fused or completed before { newL[k] = tmpL(1)[k]; for (g=2;g<=size(M);g++) { newL[k] = newL[k]+tmpL(g)[k]; } newL[k] = simplify(NF(newL[k],L[k]),2); attrib(newL[k],"isSB",1); if (newL[k][1] != 0) // i.e., the Orbit is not complete under all generators yet { GenOK=-1; L[k] = L[k],newL[k]; attrib(L[k],"isSB",1); // fuse orbits: for (k2=k+1;k2<=size(L);k2++) { if (size(NF(newL[k2],L[k]))!=size(newL[k2])) { L[k] = L[k]+L[k2]; attrib(L[k],"isSB",1); newL[k] = newL[k]+newL[k2]; attrib(newL[k],"isSB",1); L[k2]=ideal(0); newL[k2]=ideal(0); } } // fuse orbits } // the orbit was not complete yet else { countOrbits++; AllOrbits[countOrbits]=sum(L[k]); L[k]=ideal(0); newL[k]=ideal(0); } // the orbit was complete } // if orbit was not yet fused/completed before } // loop through partial orbits if (GenOK == 0) // i.e., all orbits are complete, so we can break the iteration { break; } } // iterate degBound = dgb; return(AllOrbits); } example { "EXAMPLE: orbit sums in degree 2 under the action of the alternating group;"; " it doesn't matter that we are in the modular case"; echo=2; ring R=3,(a,b,c,d),dp; list Per=ideal(b,c,a,d),ideal(a,c,d,b); ideal orb=orbit_sums(kbase(std(0),2),Per); orb; } /////////////////////////////////////////////////////////////////////////////// proc invariant_algebra_perm (list #) "USAGE: invariant_algebra_perm(GEN[,v]); @* GEN: a list of generators of a permutation group. It is given in disjoint cycle form, where trivial cycles can be omitted; e.g., the generator (1,2)(3,4)(5) is given by . @* v: an optional RETURN: A minimal homogeneous generating set of the invariant ring of the group presented by GEN, type ASSUME: We are in the non-modular case, i.e., the characteristic of the basering does not divide the group order. Note that the function does not verify whether this assumption holds or not DISPLAY: Information on the progress of computations if v does not equal 0 THEORY: We do an incremental search in increasing degree d. Generators of the invariant ring are found among the orbit sums of degree d. The generators are chosen by Groebner basis techniques (see S. King: Minimal generating sets of non-modular invariant rings of finite groups). NOTE: invariant_algebra_perm should not be used in rings with weighted orders. SEE ALSO: invariant_algebra_reynolds KEYWORDS: invariant ring minimal generating set permutation group EXAMPLE: example invariant_algebra_perm; shows an example " { //----------------- checking input and setting verbose mode ------------------ if (size(#)==0) { "ERROR: There are no parameters given."; return(); } if (typeof(#[size(#)])=="int") { int v=#[size(#)]; list GEN=#[1..size(#)-1]; } else { int v=0; list GEN=#; } int i,j,k,k2; list TST; for (i=1;i<=size(GEN);i++) { if (typeof(GEN[i])!="list") { "ERROR: The permutation group must be given by generators in"; " disjoint cycle presentation"; return(); } for (j=1;j<=size(GEN[i]);j++) { if (typeof(GEN[i][j])!="list") { "ERROR: The permutation group must be given by generators in"; " disjoint cycle presentation"; return(); } for (k=1;k<=size(GEN[i][j]);k++) { if (typeof(GEN[i][j][k])!="int") { "ERROR: The permutation group must be given by generators in"; " disjoint cycle presentation"; return(); } } } TST=sort(sum(sum(GEN[i])))[1]; if (TST[1]<=0) { "ERROR: The permutation group must be given by generators in"; " disjoint cycle presentation"; return(); } for (k2=1;k2nvars(basering))) { "ERROR: The permutation group must be given by generators in"; " disjoint cycle presentation"; return(); } } } //-------------- starting computation ------------------------- list GRP = GetGroup(GEN); def br=basering; int dgb = degBound; degBound = 0; int counter, totalcount, OrbSize, NextCand, sGoffset, fin,MaxD,LastGDeg; ideal OrbSums, RedSums, NewReductor; poly helpP,lmp; // Transform the matrices generating the group into ring maps int maxG = size(GRP); matrix vars=transpose(matrix(maxideal(1))); ideal tmpI; int g; list M; // M provides maps corresponding to the group generators: for (i=1;i<=maxG;i++) { tmpI = ideal(GRP[i]*vars); M[i] = tmpI; } ideal G; // G will contain homog. alg.generators poly Inv; ideal mon; ideal sG = std(G); int d=1; while ((fin==0) or (d<=MaxD)) { if (v) { "Searching generators in degree",d; } counter = 0; // counts generators in degree d OrbSums = orbit_sums(kbase(sG,d),M); OrbSize = ncols(OrbSums); if (v) { " We have ",OrbSize," orbit sums of degree ",d; } RedSums = reduce(OrbSums,sG); sGoffset=ncols(sG); // loop through orbit_sum, and select generators among them for (i=1;i<=OrbSize;i++) { helpP = RedSums[i]; if (helpP<>0) { counter++; // found new inv. algebra generator! totalcount++; if (v) { " We found generator number ",totalcount, " in degree ",d; } G[totalcount] = OrbSums[i]; sG[sGoffset+counter] = helpP; attrib(sG,"isSB",1); NewReductor[1]=helpP; attrib(NewReductor,"isSB",1); RedSums = reduce(RedSums,NewReductor); LastGDeg=d; } // found new generator } // loop through orbit sums if ((MaxD>0) and (d==MaxD)) { if (v) { "We went beyond the degree bound, so we are done!"; } break; } d++; // Prepare computations of the next degree if (LastGDeg==d-2) { degBound=0; if (v) {"Presumably all generators have been found";} sG = groebner(sG); fin = 0; } if (LastGDeg==d-1) { degBound=d; if (v) { "Computing Groebner basis up to the new degree",d; } sG = groebner(sG); attrib(sG,"isSB",1); fin = 0; } if ((fin==0) and (dim(sG)==0)) { MaxD = maxdeg1(kbase(sG)); if (v) { if (MaxD>=d) {"We found the degree bound",MaxD; } else {"We found the degree bound",d-1;} } if (MaxD representing the Reynolds operator of a finite matrix group, where g ist the group order and n is the number of variables of the basering; @* v: an optional RETURN: A minimal homogeneous generating set of the invariant ring, type 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() DISPLAY: Information on the progress of computations if v does not equal 0 THEORY: We do an incremental search in increasing degree d. Generators of the invariant ring are found among the Reynolds images of monomials of degree d. The generators are chosen by Groebner basis techniques (see S. King: Minimal generating sets of non-modular invariant rings of finite groups). NOTE: invariant_algebra_reynolds should not be used in rings with weighted orders. SEE ALSO: invariant_algebra_perm KEYWORDS: invariant ring minimal generating set matrix group EXAMPLE: example invariant_algebra_reynolds; shows an example " { int MonStep = 10; int v; if (size(#)) { if (typeof(#[1])=="int") { v=#[1]; } } def br=basering; int dgb = degBound; degBound = 0; int block, counter, totalcount, monsize,i,j, sGoffset, fin,MaxD,LastGDeg; poly helpP,lmp; block = 1; ideal mon, MON, Inv, G; // G will contain homog. alg.generators ideal RedInv, NewG, NewReductor; int NewGcounter; intvec dimvec; ideal sG = std(G); int d=1; dimvec[d] = dim(sG); while ((fin==0) or (d<=MaxD)) { NewG = ideal(); NewGcounter=0; mon = kbase(sG,d); sGoffset=ncols(sG); monsize = ncols(mon); if (mon[1]<>0) { if (v) { " We have ",monsize, " relevant monomials in degree ",d;} block=0; // Loops through the monomials while (block0) { counter++; // found new inv. algebra generator! if (v) { " We found generator number ",counter," in degree ",d; } G[counter] = Inv[block]; sG[sGoffset+counter] = helpP; attrib(sG,"isSB",1); NewGcounter++; NewG[NewGcounter] = helpP; mon[block]=0; Inv[block]=0; RedInv[block]=0; NewReductor[1]=helpP; attrib(NewReductor,"isSB",1); RedInv = reduce(RedInv,NewReductor); LastGDeg=d; } } // loop through monomials } // if mon[1]<>0 if ((MaxD>0) and (d==MaxD)) { if (v) { "We went beyond the degree bound, so we are done!"; } break; } d++; // Prepare computations of the next degree if (LastGDeg==d-2) { degBound=0; if (v) {"Presumably all generators have been found";} sG = groebner(sG); fin = 0; } if (LastGDeg==d-1) { degBound=d; if (v) { "Computing Groebner basis up to the new degree",d; } sG = groebner(sG); attrib(sG,"isSB",1); fin = 0; } if ((fin==0) and (dim(sG)==0)) { MaxD = maxdeg1(kbase(sG)); if (v) { if (MaxD>=d) {"We found the degree bound",MaxD; } else {"We found the degree bound",d-1;} } if (MaxD