| 144 | | //////////////////////////////////////////////////////////////////////////////// |
| 145 | | // Returns i such that root^i==n, i.e. it heavily relies on the right input. |
| 146 | | //////////////////////////////////////////////////////////////////////////////// |
| 147 | | proc power(number n, number root) |
| 148 | | { int i=0; |
| 149 | | while((n/root^i)<>1) |
| 150 | | { i=i+1; |
| 151 | | } |
| 152 | | return(i); |
| 153 | | } |
| 154 | | |
| 155 | | //////////////////////////////////////////////////////////////////////////////// |
| 156 | | // Generates the Molien series when the characteristic of the base field is p>0 |
| 157 | | // and p does not divide the group order. Input is the entire group and a name |
| 158 | | // for a new ring. |
| 159 | | //////////////////////////////////////////////////////////////////////////////// |
| 160 | | proc p_molien(list #) |
| 161 | | { def br=basering; // keeping track of the base ring since |
| 162 | | int n=nvars(br); // we have to go into an extension of the |
| 163 | | int g=size(#)-2; // basefield - |
| 164 | | matrix G(1..g)=#[1..g]; // rewriting the group elements |
| 165 | | string newring=#[g+1]; |
| 166 | | int flag=#[g+2]; |
| 167 | | if (g<>1) |
| 168 | | { ring Q=0,x,dp; // we want to extend our ring as well as |
| 169 | | // the ring of rational numbers Q to |
| 170 | | // contain g-th primitive roots of unity |
| 171 | | // in order to factor characteristic |
| 172 | | // polynomials of group elements into |
| 173 | | // linear factors and lift eigenvalues to |
| 174 | | // characteristic 0 - |
| 175 | | poly minq=cycle(g); // minq now contains the size-of-group-th |
| 176 | | // cyclotomic polynomial of Q, it is |
| 177 | | // irreducible there |
| 178 | | ring `newring`=(0,e),x,dp; |
| 179 | | map f=Q,ideal(e); |
| 180 | | minpoly=number(f(minq)); // e is now a g-th primitive root of |
| 181 | | // unity - |
| 182 | | kill Q, f; // no longer needed - |
| 183 | | poly p=1; // used to build the denominator of the |
| 184 | | // new term in the Molien series |
| 185 | | matrix s[1][2]; // used for canceling - |
| 186 | | matrix M[1][2]=0,1; // will contain Molien series - |
| 187 | | ring v1br=char(br),x,dp; // we calculate the g-th cyclotomic |
| 188 | | poly minp=cycle(g); // polynomial of the base field and pick |
| 189 | | minp=factorize(minp)[1][2]; // an irreducible factor of it - |
| 190 | | if (deg(minp)==1) // in this case the base field contains |
| 191 | | { ring bre=char(br),x,dp; // g-th roots of unity already |
| 192 | | map f1=v1br,ideal(0); |
| 193 | | number e=-number((f1(minp))); // e is a g-th primitive root of unity |
| 194 | | } |
| 195 | | else |
| 196 | | { ring bre=(char(br),e),x,dp; |
| 197 | | map f1=v1br,ideal(e); |
| 198 | | minpoly=number(f1(minp)); // e is a g-th primitive root of unity |
| 199 | | } |
| 200 | | map f2=br,ideal(0); // we need f2 to map our group elements |
| 201 | | // to this new extension field bre |
| 202 | | matrix I=unitmat(n); |
| 203 | | poly p; // used for the characteristic polynomial |
| 204 | | // to factor - |
| 205 | | list L; // will contain the linear factors of the |
| 206 | | ideal F; // characteristic polynomial of the group |
| 207 | | intvec C; // elements and their powers |
| 208 | | int i, j, k; |
| 209 | | for (i=1;i<=g;i=i+1) |
| 210 | | { p=det(x*I-f2(G(i))); // characteristic polynomial of G(i) |
| 211 | | L=factorize(p); |
| 212 | | F=L[1]; |
| 213 | | C=L[2]; |
| 214 | | for (j=2;j<=ncols(F);j=j+1) |
| 215 | | { F[j]=-1*(F[j]-x); // F[j] is now an eigenvalue of G(i), |
| 216 | | // it is a power of a primitive g-th root |
| 217 | | // of unity - |
| 218 | | k=power(number(F[j]),e); // F[j]==e^k |
| 219 | | setring `newring`; |
| 220 | | p=p*(1-x*(e^k))^C[j]; // building the denominator of the new |
| 221 | | setring bre; // term |
| 222 | | } |
| 223 | | setring `newring`; |
| 224 | | M[1,1]=M[1,1]*p+M[1,2]; // expanding M[1,1]/M[1,2] + 1/p |
| 225 | | M[1,2]=M[1,2]*p; |
| 226 | | p=1; |
| 227 | | s=matrix(syz(ideal(M))); // canceling common terms of denominator |
| 228 | | M[1,1]=-s[2,1]; // and enumerator |
| 229 | | M[1,2]=s[1,1]; |
| 230 | | setring bre; |
| 231 | | if (flag) |
| 232 | | { " Term "+string(i)+" has been computed."; |
| 233 | | } |
| 234 | | } |
| 235 | | if (flag) |
| 236 | | { ""; |
| 237 | | } |
| 238 | | setring `newring`; |
| 239 | | map slead=`newring`,ideal(0); |
| 240 | | s=slead(M); // forcing the constant term of numerator |
| 241 | | M[1,1]=1/s[1,1]*M[1,1]; // and denominator to be 1 |
| 242 | | M[1,2]=1/s[1,2]*M[1,2]; |
| 243 | | kill slead; |
| 244 | | kill s; |
| 245 | | kill p; |
| 246 | | } |
| 247 | | else // if the group only contains an identity |
| 248 | | { ring `newring`=0,x,dp; // element, it is very easy to calculate |
| 249 | | matrix M[1][2]=1,(1-x)^n; // the Molien series |
| 250 | | } |
| 251 | | // keepring `newring`; |
| 252 | | export `newring`; // TTO we keep the ring where we computed |
| 253 | | // the Molien series |
| 254 | | export M; // TTO so that we can keep the Molien |
| 255 | | // series |
| 256 | | setring br; |
| 257 | | } |
| 258 | | |
| 259 | | //////////////////////////////////////////////////////////////////////////////// |
| 260 | | // This procedure calculates all members of a finite matrix group in terms of |
| 261 | | // the given generators. In one run trough the main loop, all left products of |
| 262 | | // the generators with the new elements from the last run through the loop (or |
| 263 | | // the generators themselves in the first run) will be formed. After that the |
| 264 | | // newly generated elements will be added to the group and the loop starts over |
| 265 | | // again unless no elements were added. |
| 266 | | // Additionally, every time a new matrix is added to the group, its |
| 267 | | // corresponding ring mapping in the Reynolds operator and if the |
| 268 | | // characteristic is 0, its corresponding summand of the Molien series is |
| 269 | | // calculated. |
| 270 | | // When the characteristic of the basefield is p>0 such that it does not |
| 271 | | // divide the group order, the Molien series is calculated at the end of the |
| 272 | | // procedure. |
| 273 | | // No matter when the Molien series is calculated, the procedure expands after |
| 274 | | // every step to obtain a rational function. |
| 275 | | // The first result of the procedure is the Reynolds operator, presented in |
| 276 | | // form of a matrix; each row can be transformed into an ideal and from |
| 277 | | // there can be used as a ring homomorphism via the command 'map'. |
| 278 | | // If the characteristic is 0, the second result is a matrix, containing |
| 279 | | // enumerator and denominator (with no common divisor) of the final |
| 280 | | // rational function representing the Molien series. |
| 281 | | // When the characteristic of the basefield is p>0 such that it does not |
| 282 | | // divide the group order, the Molien series is returned in a ring of |
| 283 | | // characteristic 0. It names was specified in the list of parameters. |
| 284 | | //////////////////////////////////////////////////////////////////////////////// |
| 285 | | proc rey_mol (list #) |
| 286 | | USAGE: rey_mol(<generators of a finite matrix group>[,<string>,<int>]); |
| 287 | | if the characteristic of the coefficient field is prime, <string> |
| 288 | | has to contain the name for a new polynomials ring with coefficient |
| 289 | | field of characteristic 0 that stores the Molien series - if <int> is |
| 290 | | not not equal to 0, some information will be printed during the run |
| 291 | | RETURNS: if the characteristic is 0: Reynolds operator (type <matrix>), Molien |
| 292 | | series (type <matrix> with two components, first being the numerator, |
| 293 | | second the denominator) |
| 294 | | if the characteristic is p>0 not dividing the group order: Reynolds |
| 295 | | operator (type <matrix>) - the Molien series will directly be stored |
| 296 | | under the name M (type <matrix>) in the ring `<string>` |
| 297 | | if the characteristic is p>0 dividing the group order: Reynolds |
| 298 | | operator (type <matrix>) |
| 299 | | EXAMPLE: example rey_mol; shows an example |
| 300 | | { def br=basering; // the Molien series depends on the |
| 301 | | int ch=char(br); // characteristic of the coefficient |
| 302 | | int flag; // field - |
| 303 | | if (ch<>0) // making sure the input is 'correct'... |
| 304 | | { if (typeof(#[size(#)])=="string") |
| 305 | | { flag=size(#)-1; |
| 306 | | string newring=#[size(#)]; |
| 307 | | int v=0; // no information is default |
| 308 | | } |
| 309 | | else |
| 310 | | { if (typeof(#[size(#)-1])=="string") |
| 311 | | { flag=size(#)-2; |
| 312 | | string newring=#[size(#)-1]; |
| 313 | | if (typeof(#[size(#)])<>"int") |
| 314 | | { " ERROR: if the second last parameter is <string>, the last must be"; |
| 315 | | " of type <int>"; |
| 316 | | return(); |
| 317 | | } |
| 318 | | int v=#[size(#)]; |
| 319 | | } |
| 320 | | else |
| 321 | | { " ERROR: in characteristic p a <string> must be given for the name"; |
| 322 | | " of a new ring"; |
| 323 | | return(); |
| 324 | | } |
| 325 | | } |
| 326 | | if (newring=="") |
| 327 | | { " ERROR: <string> may not be empty"; |
| | 111 | proc group_reynolds (list #) |
| | 112 | USAGE: group_reynolds(G1,G2,...[,v]); |
| | 113 | G1,G2,...: nxn <matrices> generating a finite matrix group, v: an |
| | 114 | optional <int> |
| | 115 | ASSUME: n is the number of variables of the basering, g the number of group |
| | 116 | elements |
| | 117 | RETURN: a <list>, the first list element will be a gxn <matrix> representing |
| | 118 | the Reynolds operator if we are in the non-modular case; if the |
| | 119 | characteristic is >0, minpoly==0 and the finite group non-cyclic the |
| | 120 | second list element is an <int> giving the lowest common multiple of |
| | 121 | the matrix group elements (used in molien); in general all other list |
| | 122 | elements are nxn <matrices> listing all elements of the finite group |
| | 123 | DISPLAY: information if v does not equal 0 |
| | 124 | EXAMPLE: example group_reynolds; shows an example |
| | 125 | THEORY: The entire matrix group is generated by getting all left products of |
| | 126 | the generators with the new elements from the last run through the loop |
| | 127 | (or the generators themselves during the first run). All the ones that |
| | 128 | have been generated before are thrown out and the program terminates |
| | 129 | when there are no new elements found in one run. Additionally each time |
| | 130 | a new group element is found the corresponding ring mapping of which |
| | 131 | the Reynolds operator is made up is generated. They are stored in the |
| | 132 | rows of the first return value. |
| | 133 | { int ch=char(basering); // the existance of the Reynolds operator |
| | 134 | // is dependent on the characteristic of |
| | 135 | // the base field |
| | 136 | int gen_num; // number of generators |
| | 137 | //------------------------ making sure the input is okay --------------------- |
| | 138 | if (typeof(#[size(#)])=="int") |
| | 139 | { if (size(#)==1) |
| | 140 | { "ERROR: there are no matrices given among the parameters"; |
| 495 | | return(A(1)); |
| 496 | | } |
| 497 | | } |
| 498 | | if (ch<>0) |
| 499 | | { if ((g%ch)==0) |
| | 480 | return(); |
| | 481 | } |
| | 482 | } |
| | 483 | //---------------------------------------------------------------------------- |
| | 484 | if (not(typeof(#[1])=="matrix")) |
| | 485 | { "ERROR: the parameters must be a list of matrices and maybe an <intvec>"; |
| | 486 | return(); |
| | 487 | } |
| | 488 | int n=nrows(#[1]); |
| | 489 | if (n<>nvars(br)) |
| | 490 | { "ERROR: the number of variables of the basering needs to be the same"; |
| | 491 | " as the dimension of the square matrices"; |
| | 492 | return(); |
| | 493 | } |
| | 494 | if (v && voice<>2) |
| | 495 | { ""; |
| | 496 | " Generating the Molien series..."; |
| | 497 | ""; |
| | 498 | } |
| | 499 | if (v && voice==2) |
| | 500 | { ""; |
| | 501 | } |
| | 502 | //------------- calculating Molien series in characteristic 0 ---------------- |
| | 503 | if (ch==0) // when ch==0 we can calculate the Molien |
| | 504 | { matrix I=diag(1,n); // series in any case - |
| | 505 | poly v1=maxideal(1)[1]; // the Molien series will be in terms of |
| | 506 | // the first variable of the current |
| | 507 | // ring - |
| | 508 | matrix M[1][2]; // M will contain the Molien series - |
| | 509 | M[1,1]=0; // M[1,1] will be the numerator - |
| | 510 | M[1,2]=1; // M[1,2] will be the denominator - |
| | 511 | matrix s; // will help us canceling in the |
| | 512 | // fraction |
| | 513 | poly p; // will contain the denominator of the |
| | 514 | // new term of the Molien series |
| | 515 | //------------ computing 1/det(1+xE) for all E in the group ------------------ |
| | 516 | for (int j=1;j<=g;j=j+1) |
| | 517 | { if (not(typeof(#[j])=="matrix")) |
| | 518 | { "ERROR: the parameters must be a list of matrices and maybe an <intvec>"; |
| | 519 | return(); |
| | 520 | } |
| | 521 | if ((n<>nrows(#[j])) or (n<>ncols(#[j]))) |
| | 522 | { "ERROR: matrices need to be square and of the same dimensions"; |
| | 523 | return(); |
| | 524 | } |
| | 525 | p=det(I-v1*#[j]); // denominator of new term - |
| | 526 | M[1,1]=M[1,1]*p+M[1,2]; // expanding M[1,1]/M[1,2] + 1/p |
| | 527 | M[1,2]=M[1,2]*p; |
| | 528 | if (interval<>0) // canceling common terms of denominator |
| | 529 | { if ((j/interval)*interval==j or j==g) // and enumerator - |
| | 530 | { s=matrix(syz(ideal(M))); // once gcd() is faster than syz() these |
| | 531 | M[1,1]=-s[2,1]; // three lines should be replaced by the |
| | 532 | M[1,2]=s[1,1]; // following three |
| | 533 | // p=gcd(M[1,1],M[1,2]); |
| | 534 | // M[1,1]=M[1,1]/p; |
| | 535 | // M[1,2]=M[1,2]/p; |
| | 536 | } |
| | 537 | } |
| | 538 | if (v) |
| | 539 | { " Term "+string(j)+" of the Molien series has been computed."; |
| | 540 | } |
| | 541 | } |
| | 542 | if (interval==0) // canceling common terms of denominator |
| | 543 | { // and enumerator - |
| | 544 | s=matrix(syz(ideal(M))); // once gcd() is faster than syz() these |
| | 545 | M[1,1]=-s[2,1]; // three lines should be replaced by the |
| | 546 | M[1,2]=s[1,1]; // following three |
| | 547 | // p=gcd(M[1,1],M[1,2]); |
| | 548 | // M[1,1]=M[1,1]/p; |
| | 549 | // M[1,2]=M[1,2]/p; |
| | 550 | } |
| | 551 | map slead=br,ideal(0); |
| | 552 | s=slead(M); |
| | 553 | M[1,1]=1/s[1,1]*M[1,1]; // numerator and denominator have to have |
| | 554 | M[1,2]=1/s[1,2]*M[1,2]; // a constant term of 1 |
| | 555 | if (v) |
| | 556 | { ""; |
| | 557 | " We are done calculating the Molien series."; |
| | 558 | ""; |
| | 559 | } |
| | 560 | return(M); |
| | 561 | } |
| | 562 | //---- calculating Molien series in prime characteristic with Brauer lift ---- |
| | 563 | if (ch<>0 && mol_flag==0) |
| | 564 | { if (g<>1) |
| | 565 | { matrix G(1..g)=#[1..g]; |
| | 566 | if (interval<0) |
| | 567 | { string Mstring; |
| | 568 | } |
| | 569 | //------ preparing everything for Brauer lifts into characteristic 0 --------- |
| | 570 | ring Q=0,x,dp; // we want to extend our ring as well as |
| | 571 | // the ring of rational numbers Q to |
| | 572 | // contain r-th primitive roots of unity |
| | 573 | // in order to factor characteristic |
| | 574 | // polynomials of group elements into |
| | 575 | // linear factors and lift eigenvalues to |
| | 576 | // characteristic 0 - |
| | 577 | poly minq=cyclotomic(r); // minq now contains the size-of-group-th |
| | 578 | // cyclotomic polynomial of Q, it is |
| | 579 | // irreducible there |
| | 580 | ring `newring`=(0,e),x,dp; |
| | 581 | map f=Q,ideal(e); |
| | 582 | minpoly=number(f(minq)); // e is now a r-th primitive root of |
| | 583 | // unity - |
| | 584 | kill Q, f; // no longer needed - |
| | 585 | poly p=1; // used to build the denominator of the |
| | 586 | // new term in the Molien series |
| | 587 | matrix s[1][2]; // used for canceling - |
| | 588 | matrix M[1][2]=0,1; // will contain Molien series - |
| | 589 | ring v1br=char(br),x,dp; // we calculate the r-th cyclotomic |
| | 590 | poly minp=cyclotomic(r); // polynomial of the base field and pick |
| | 591 | minp=factorize(minp)[1][2]; // an irreducible factor of it - |
| | 592 | if (deg(minp)==1) // in this case the base field contains |
| | 593 | { ring bre=char(br),x,dp; // r-th roots of unity already |
| | 594 | map f1=v1br,ideal(0); |
| | 595 | number e=-number((f1(minp))); // e is a r-th primitive root of unity |
| | 596 | } |
| | 597 | else |
| | 598 | { ring bre=(char(br),e),x,dp; |
| | 599 | map f1=v1br,ideal(e); |
| | 600 | minpoly=number(f1(minp)); // e is a r-th primitive root of unity |
| | 601 | } |
| | 602 | map f2=br,ideal(0); // we need f2 to map our group elements |
| | 603 | // to this new extension field bre |
| | 604 | matrix xI=diag(x,n); |
| | 605 | poly p; // used for the characteristic polynomial |
| | 606 | // to factor - |
| | 607 | list L; // will contain the linear factors of the |
| | 608 | ideal F; // characteristic polynomial of the group |
| | 609 | intvec C; // elements and their powers |
| | 610 | int i, j, k; |
| | 611 | // -------------- finding all the terms of the Molien series ----------------- |
| | 612 | for (i=1;i<=g;i=i+1) |
| | 613 | { setring br; |
| | 614 | if (not(typeof(#[i])=="matrix")) |
| | 615 | { "ERROR: the parameters must be a list of matrices and maybe an <intvec>"; |
| | 616 | return(); |
| | 617 | } |
| | 618 | if ((n<>nrows(#[i])) or (n<>ncols(#[i]))) |
| | 619 | { "ERROR: matrices need to be square and of the same dimensions"; |
| | 620 | return(); |
| | 621 | } |
| | 622 | setring bre; |
| | 623 | p=det(xI-f2(G(i))); // characteristic polynomial of G(i) |
| | 624 | L=factorize(p); |
| | 625 | F=L[1]; |
| | 626 | C=L[2]; |
| | 627 | for (j=2;j<=ncols(F);j=j+1) |
| | 628 | { F[j]=-1*(F[j]-x); // F[j] is now an eigenvalue of G(i), |
| | 629 | // it is a power of a primitive r-th root |
| | 630 | // of unity - |
| | 631 | k=exponent(number(F[j]),e); // F[j]==e^k |
| | 632 | setring `newring`; |
| | 633 | p=p*(1-x*(e^k))^C[j]; // building the denominator of the new |
| | 634 | setring bre; // term |
| | 635 | } |
| | 636 | // ----------- |
| | 637 | // k=0; |
| | 638 | // while(k<r) |
| | 639 | // { map f3=basering,ideal(e^k); |
| | 640 | // while (f3(p)==0) |
| | 641 | // { p=p/(x-e^k); |
| | 642 | // setring `newring`; |
| | 643 | // p=p*(1-x*(e^k)); // building the denominator of the new |
| | 644 | // setring bre; |
| | 645 | // } |
| | 646 | // kill f3; |
| | 647 | // if (p==1) |
| | 648 | // { break; |
| | 649 | // } |
| | 650 | // k=k+1; |
| | 651 | // } |
| | 652 | setring `newring`; |
| | 653 | M[1,1]=M[1,1]*p+M[1,2]; // expanding M[1,1]/M[1,2] + 1/p |
| | 654 | M[1,2]=M[1,2]*p; |
| | 655 | if (interval<0) |
| | 656 | { if (i<>g) |
| | 657 | { Mstring=string(M); |
| | 658 | for (j=1;j<=size(Mstring);j=j+1) |
| | 659 | { if (Mstring[j]=="e") |
| | 660 | { interval=0; |
| | 661 | break; |
| | 662 | } |
| | 663 | } |
| | 664 | } |
| | 665 | if (interval<>0) |
| | 666 | { s=matrix(syz(ideal(M))); // once gcd() is faster than syz() |
| | 667 | M[1,1]=-s[2,1]; // these three lines should be |
| | 668 | M[1,2]=s[1,1]; // replaced by the following three |
| | 669 | // p=gcd(M[1,1],M[1,2]); |
| | 670 | // M[1,1]=M[1,1]/p; |
| | 671 | // M[1,2]=M[1,2]/p; |
| | 672 | } |
| | 673 | else |
| | 674 | { interval=-1; |
| | 675 | } |
| | 676 | } |
| | 677 | else |
| | 678 | { if (interval<>0) // canceling common terms of denominator |
| | 679 | { if ((i/interval)*interval==i or i==g) // and enumerator |
| | 680 | { s=matrix(syz(ideal(M))); // once gcd() is faster than syz() |
| | 681 | M[1,1]=-s[2,1]; // these three lines should be |
| | 682 | M[1,2]=s[1,1]; // replaced by the following three |
| | 683 | // p=gcd(M[1,1],M[1,2]); |
| | 684 | // M[1,1]=M[1,1]/p; |
| | 685 | // M[1,2]=M[1,2]/p; |
| | 686 | } |
| | 687 | } |
| | 688 | } |
| | 689 | p=1; |
| | 690 | setring bre; |
| | 691 | if (v) |
| | 692 | { " Term "+string(i)+" of the Molien series has been computed."; |
| | 693 | } |
| | 694 | } |
| | 695 | if (v) |
| | 696 | { ""; |
| | 697 | } |
| | 698 | setring `newring`; |
| | 699 | if (interval==0) // canceling common terms of denominator |
| | 700 | { // and enumerator - |
| | 701 | s=matrix(syz(ideal(M))); // once gcd() is faster than syz() these |
| | 702 | M[1,1]=-s[2,1]; // three lines should be replaced by the |
| | 703 | M[1,2]=s[1,1]; // following three |
| | 704 | // p=gcd(M[1,1],M[1,2]); |
| | 705 | // M[1,1]=M[1,1]/p; |
| | 706 | // M[1,2]=M[1,2]/p; |
| | 707 | } |
| | 708 | map slead=`newring`,ideal(0); |
| | 709 | s=slead(M); // forcing the constant term of numerator |
| | 710 | M[1,1]=1/s[1,1]*M[1,1]; // and denominator to be 1 |
| | 711 | M[1,2]=1/s[1,2]*M[1,2]; |
| | 712 | kill slead; |
| | 713 | kill s; |
| | 714 | kill p; |
| | 715 | } |
| | 716 | else // if the group only contains an identity |
| | 717 | { ring `newring`=0,x,dp; // element, it is very easy to calculate |
| | 718 | matrix M[1][2]=1,(1-x)^n; // the Molien series |
| | 719 | } |
| | 720 | export `newring`; // we keep the ring where we computed the |
| | 721 | export M; // Molien series in such that we can |
| | 722 | setring br; // keep it |
| | 723 | if (v) |
| | 724 | { " We are done calculating the Molien series."; |
| | 725 | ""; |
| | 726 | } |
| | 727 | } |
| | 728 | else // i.e. char<>0 and mol_flag<>0, the user |
| | 729 | { // has specified that we are dealing with |
| | 730 | // a ring of large characteristic which |
| | 731 | // can be treated like a ring of |
| | 732 | // characteristic 0; we'll avoid the |
| | 733 | // Brauer lifts |
| | 734 | //----------------------- simulating characteristic 0 ------------------------ |
| | 735 | string chst=charstr(br); |
| | 736 | for (int i=1;i<=size(chst);i=i+1) |
| | 737 | { if (chst[i]==",") |
| | 738 | { break; |
| | 739 | } |
| | 740 | } |
| | 741 | //----------------- generating ring of characteristic 0 ---------------------- |
| | 742 | if (minpoly==0) |
| | 743 | { if (i>size(chst)) |
| | 744 | { execute "ring "+newring+"=0,("+varstr(br)+"),("+ordstr(br)+")"; |
| | 745 | } |
| | 746 | else |
| | 747 | { chst=chst[i..size(chst)]; |
| | 748 | execute "ring "+newring+"=(0"+chst+"),("+varstr(br)+"),("+ordstr(br)+")"; |
| | 749 | } |
| | 750 | } |
| | 751 | else |
| | 752 | { string minp=string(minpoly); |
| | 753 | minp=minp[2..size(minp)-1]; |
| | 754 | chst=chst[i..size(chst)]; |
| | 755 | execute "ring "+newring+"=(0"+chst+"),("+varstr(br)+"),("+ordstr(br)+")"; |
| | 756 | execute "minpoly="+minp; |
| | 757 | } |
| | 758 | matrix I=diag(1,n); |
| | 759 | poly v1=maxideal(1)[1]; // the Molien series will be in terms of |
| | 760 | // the first variable of the current |
| | 761 | // ring - |
| | 762 | matrix M[1][2]; // M will contain the Molien series - |
| | 763 | M[1,1]=0; // M[1,1] will be the numerator - |
| | 764 | M[1,2]=1; // M[1,2] will be the denominator - |
| | 765 | matrix s; // will help us canceling in the |
| | 766 | // fraction |
| | 767 | poly p; // will contain the denominator of the |
| | 768 | // new term of the Molien series |
| | 769 | int j; |
| | 770 | string links, rechts; |
| | 771 | //----------------- finding all terms of the Molien series ------------------- |
| | 772 | for (i=1;i<=g;i=i+1) |
| | 773 | { setring br; |
| | 774 | if (not(typeof(#[i])=="matrix")) |
| | 775 | { "ERROR: the parameters must be a list of matrices and maybe an <intvec>"; |
| | 776 | return(); |
| | 777 | } |
| | 778 | if ((n<>nrows(#[i])) or (n<>ncols(#[i]))) |
| | 779 | { "ERROR: matrices need to be square and of the same dimensions"; |
| | 780 | return(); |
| | 781 | } |
| | 782 | string stM(i)=string(#[i]); |
| | 783 | for (j=1;j<=size(stM(i));j=j+1) |
| | 784 | { if (stM(i)[j]==" |
| | 785 | ") |
| | 786 | { links=stM(i)[1..j-1]; |
| | 787 | rechts=stM(i)[j+1..size(stM(i))]; |
| | 788 | stM(i)=links+rechts; |
| | 789 | } |
| | 790 | } |
| | 791 | setring `newring`; |
| | 792 | execute "matrix G(i)["+string(n)+"]["+string(n)+"]="+stM(i); |
| | 793 | p=det(I-v1*G(i)); // denominator of new term - |
| | 794 | M[1,1]=M[1,1]*p+M[1,2]; // expanding M[1,1]/M[1,2] + 1/p |
| | 795 | M[1,2]=M[1,2]*p; |
| | 796 | if (interval<>0) // canceling common terms of denominator |
| | 797 | { if ((i/interval)*interval==i or i==g) // and enumerator |
| | 798 | { |
| | 799 | s=matrix(syz(ideal(M))); // once gcd() is faster than syz() these |
| | 800 | M[1,1]=-s[2,1]; // three lines should be replaced by the |
| | 801 | M[1,2]=s[1,1]; // following three |
| | 802 | // p=gcd(M[1,1],M[1,2]); |
| | 803 | // M[1,1]=M[1,1]/p; |
| | 804 | // M[1,2]=M[1,2]/p; |
| | 805 | } |
| | 806 | } |
| | 807 | if (v) |
| | 808 | { " Term "+string(i)+" of the Molien series has been computed."; |
| | 809 | } |
| | 810 | } |
| | 811 | if (interval==0) // canceling common terms of denominator |
| | 812 | { // and enumerator - |
| | 813 | s=matrix(syz(ideal(M))); // once gcd() is faster than syz() these |
| | 814 | M[1,1]=-s[2,1]; // three lines should be replaced by the |
| | 815 | M[1,2]=s[1,1]; // following three |
| | 816 | // p=gcd(M[1,1],M[1,2]); |
| | 817 | // M[1,1]=M[1,1]/p; |
| | 818 | // M[1,2]=M[1,2]/p; |
| | 819 | } |
| | 820 | map slead=`newring`,ideal(0); |
| | 821 | s=slead(M); |
| | 822 | M[1,1]=1/s[1,1]*M[1,1]; // numerator and denominator have to have |
| | 823 | M[1,2]=1/s[1,2]*M[1,2]; // a constant term of 1 |
| | 824 | if (v) |
| | 825 | { ""; |
| | 826 | " We are done calculating the Molien series."; |
| | 827 | ""; |
| | 828 | } |
| | 829 | kill G(1..g), s, slead, p, v1, I; |
| | 830 | export `newring`; // we keep the ring where we computed the |
| | 831 | export M; // the Molien series such that we can |
| | 832 | setring br; // keep it |
| | 833 | } |
| | 834 | } |
| | 835 | example |
| | 836 | { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; |
| | 837 | " note the case of prime characteristic"; |
| | 838 | echo=2; |
| | 839 | ring R=0,(x,y,z),dp; |
| | 840 | matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; |
| | 841 | list L=group_reynolds(A); |
| | 842 | matrix M=molien(L[2..size(L)]); |
| | 843 | print(M); |
| | 844 | ring S=3,(x,y,z),dp; |
| | 845 | string newring="alksdfjlaskdjf"; |
| | 846 | matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; |
| | 847 | list L=group_reynolds(A); |
| | 848 | molien(L[2..size(L)],newring); |
| | 849 | setring alksdfjlaskdjf; |
| | 850 | print(M); |
| | 851 | setring S; |
| | 852 | kill alksdfjlaskdjf; |
| | 853 | } |
| | 854 | |
| | 855 | proc reynolds_molien (list #) |
| | 856 | USAGE: reynolds_molien(G1,G2,...[,ringname,flags]); |
| | 857 | G1,G2,...: nxn <matrices> generating a finite matrix group, ringname: |
| | 858 | a <string> giving a name for a new ring of characteristic 0 for the |
| | 859 | Molien series in case of prime characteristic, flags: an optional |
| | 860 | <intvec> with three components: if the first element is not equal to 0 |
| | 861 | characteristic 0 is simulated, i.e. the Molien series is computed as |
| | 862 | if the base field were characteristic 0 (the user must choose a field |
| | 863 | of large prime characteristic, e.g. 32003) the second component should |
| | 864 | give the size of intervals between canceling common factors in the |
| | 865 | expansion of the Molien series, 0 (the default) means only once after |
| | 866 | generating all terms, in prime characteristic also a negative number |
| | 867 | can be given to indicate that common factors should always be canceled |
| | 868 | when the expansion is simple (the root of the extension field does not |
| | 869 | occur among the coefficients) |
| | 870 | ASSUME: n is the number of variables of the basering, G1,G2... are the group |
| | 871 | elements generated by group_reynolds(), g is the size of the group |
| | 872 | RETURN: a gxn <matrix> representing the Reynolds operator is the first return |
| | 873 | value and in case of characteristic 0 a 1x2 <matrix> giving enumerator |
| | 874 | and denominator of Molien series is the second one; in case of prime |
| | 875 | characteristic a ring with the name `ringname` of characteristic 0 is |
| | 876 | created where the same Molien series (named M) is stored |
| | 877 | DISPLAY: information if the third component of flags does not equal 0 |
| | 878 | EXAMPLE: example reynolds_molien; shows an example |
| | 879 | THEORY: The entire matrix group is generated by getting all left products of |
| | 880 | the generators with the new elements from the last run through the loop |
| | 881 | (or the generators themselves during the first run). All the ones that |
| | 882 | have been generated before are thrown out and the program terminates |
| | 883 | when there are no new elements found in one run. Additionally each time |
| | 884 | a new group element is found the corresponding ring mapping of which |
| | 885 | the Reynolds operator is made up is generated. They are stored in the |
| | 886 | rows of the first return value. In characteristic 0 the terms |
| | 887 | 1/det(1-xE) is computed whenever a new element E is found. In prime |
| | 888 | characteristic a Brauer lift is involved and the terms are only |
| | 889 | computed after the entire matrix group is generated (to avoid the |
| | 890 | modular case). The returned matrix gives enumerator and denominator of |
| | 891 | the expanded version where common factors have been canceled. |
| | 892 | { def br=basering; // the Molien series depends on the |
| | 893 | int ch=char(br); // characteristic of the coefficient |
| | 894 | // field |
| | 895 | int gen_num; |
| | 896 | //------------------- making sure the input is okay -------------------------- |
| | 897 | if (typeof(#[size(#)])=="intvec") |
| | 898 | { if (size(#[size(#)])==3) |
| | 899 | { int mol_flag=#[size(#)][1]; |
| | 900 | if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0))) |
| | 901 | { "ERROR: the second component of the <intvec> should be >=0"; |
| | 902 | return(); |
| | 903 | } |
| | 904 | int interval=#[size(#)][2]; |
| | 905 | int v=#[size(#)][3]; |
| | 906 | } |
| | 907 | else |
| | 908 | { "ERROR: <intvec> should have three components"; |
| | 909 | return(); |
| | 910 | } |
| | 911 | if (ch<>0) |
| | 912 | { if (typeof(#[size(#)-1])<>"string") |
| | 913 | { "ERROR: in characteristic p a <string> must be given for the name"; |
| | 914 | " of a new ring where the Molien series can be stored"; |
| | 915 | return(); |
| | 916 | } |
| | 917 | else |
| | 918 | { if (#[size(#)-1]=="") |
| | 919 | { "ERROR: <string> may not be empty"; |
| | 920 | return(); |
| | 921 | } |
| | 922 | string newring=#[size(#)-1]; |
| | 923 | gen_num=size(#)-2; |
| | 924 | } |
| | 925 | } |
| | 926 | else // then <string> ist not needed |
| | 927 | { gen_num=size(#)-1; |
| | 928 | } |
| | 929 | } |
| | 930 | else // last parameter is not <intvec> |
| | 931 | { int v=0; // no information is default |
| | 932 | int interval; |
| | 933 | int mol_flag=0; // computing of Molien series is default |
| | 934 | if (ch<>0) |
| | 935 | { if (typeof(#[size(#)])<>"string") |
| | 936 | { "ERROR: in characteristic p a <string> must be given for the name"; |
| | 937 | " of a new ring where the Molien series can be stored"; |
| | 938 | return(); |
| | 939 | } |
| | 940 | else |
| | 941 | { if (#[size(#)]=="") |
| | 942 | { "ERROR: <string> may not be empty"; |
| | 943 | return(); |
| | 944 | } |
| | 945 | string newring=#[size(#)]; |
| | 946 | gen_num=size(#)-1; |
| | 947 | } |
| | 948 | } |
| | 949 | else |
| | 950 | { gen_num=size(#); |
| | 951 | } |
| | 952 | } |
| | 953 | // ----------------- computing the terms with Brauer lift -------------------- |
| | 954 | if (ch<>0 && mol_flag==0) |
| | 955 | { list L=group_reynolds(#[1..gen_num],v); |
| | 956 | if (L[1]==0) |
| | 1151 | //------------------------ simulating characteristic 0 ----------------------- |
| | 1152 | else // if ch<>0 and mol_flag<>0 |
| | 1153 | { if (typeof(#[1])<>"matrix") |
| | 1154 | { "ERROR: the parameters must be a list of matrices and maybe an <intvec>"; |
| | 1155 | return(); |
| | 1156 | } |
| | 1157 | int n=nrows(#[1]); |
| | 1158 | if (n<>nvars(br)) |
| | 1159 | { "ERROR: the number of variables of the basering needs to be the same"; |
| | 1160 | " as the dimension of the matrices"; |
| | 1161 | return(); |
| | 1162 | } |
| | 1163 | if (n<>ncols(#[1])) |
| | 1164 | { "ERROR: matrices need to be square and of the same dimensions"; |
| | 1165 | return(); |
| | 1166 | } |
| | 1167 | matrix vars=matrix(maxideal(1)); // creating an nx1-matrix containing the |
| | 1168 | vars=transpose(vars); // variables of the ring - |
| | 1169 | matrix A(1)=#[1]*vars; // calculating the first ring mapping - |
| | 1170 | // A(1) will contain the Reynolds |
| | 1171 | // operator |
| | 1172 | string chst=charstr(br); |
| | 1173 | for (int i=1;i<=size(chst);i=i+1) |
| | 1174 | { if (chst[i]==",") |
| | 1175 | { break; |
| | 1176 | } |
| | 1177 | } |
| | 1178 | if (minpoly==0) |
| | 1179 | { if (i>size(chst)) |
| | 1180 | { execute "ring "+newring+"=0,("+varstr(br)+"),("+ordstr(br)+")"; |
| | 1181 | } |
| | 1182 | else |
| | 1183 | { chst=chst[i..size(chst)]; |
| | 1184 | execute "ring "+newring+"=(0"+chst+"),("+varstr(br)+"),("+ordstr(br)+")"; |
| | 1185 | } |
| | 1186 | } |
| | 1187 | else |
| | 1188 | { string minp=string(minpoly); |
| | 1189 | minp=minp[2..size(minp)-1]; |
| | 1190 | chst=chst[i..size(chst)]; |
| | 1191 | execute "ring "+newring+"=(0"+chst+"),("+varstr(br)+"),("+ordstr(br)+")"; |
| | 1192 | execute "minpoly="+minp; |
| | 1193 | } |
| | 1194 | poly v1=var(1); // the Molien series will be in terms of |
| | 1195 | // the first variable of the current |
| | 1196 | // ring |
| | 1197 | matrix I=diag(1,n); |
| | 1198 | int o; |
| | 1199 | setring br; |
| | 1200 | matrix G(1)=#[1]; |
| | 1201 | string links, rechts; |
| | 1202 | string stM(1)=string(#[1]); |
| | 1203 | for (o=1;o<=size(stM(1));o=o+1) |
| | 1204 | { if (stM(1)[o]==" |
| | 1205 | ") |
| | 1206 | { links=stM(1)[1..o-1]; |
| | 1207 | rechts=stM(1)[o+1..size(stM(1))]; |
| | 1208 | stM(1)=links+rechts; |
| | 1209 | } |
| | 1210 | } |
| | 1211 | setring `newring`; |
| | 1212 | execute "matrix G(1)["+string(n)+"]["+string(n)+"]="+stM(1); |
| | 1213 | matrix A(2)[1][2]; // A(2) will contain the Molien series - |
| | 1214 | A(2)[1,1]=1; // A(2)[1,1] will be the numerator |
| | 1215 | A(2)[1,2]=det(I-v1*(G(1))); // A(2)[1,2] will be the denominator - |
| | 1216 | matrix s; // will help us canceling in the |
| | 1217 | // fraction |
| | 1218 | poly p; // will contain the denominator of the |
| | 1219 | // new term of the Molien series |
| | 1220 | i=1; |
| | 1221 | for (int j=2;j<=gen_num;j=j+1) // this loop adds the parameters to the |
| | 1222 | { // group, leaving out doubles and |
| | 1223 | // checking whether the parameters are |
| | 1224 | // compatible with the task of the |
| | 1225 | // procedure |
| | 1226 | setring br; |
| | 1227 | if (not(typeof(#[j])=="matrix")) |
| | 1228 | { "ERROR: the parameters must be a list of matrices and maybe an <intvec>"; |
| | 1229 | return(); |
| | 1230 | } |
| | 1231 | if ((n!=nrows(#[j])) or (n!=ncols(#[j]))) |
| | 1232 | { "ERROR: matrices need to be square and of the same dimensions"; |
| | 1233 | return(); |
| | 1234 | } |
| | 1235 | if (unique(G(1..i),#[j])) |
| | 1236 | { i=i+1; |
| | 1237 | matrix G(i)=#[j]; |
| | 1238 | A(1)=concat(A(1),G(i)*vars); // adding ring homomorphisms to A(1) |
| | 1239 | string stM(i)=string(G(i)); |
| | 1240 | for (o=1;o<=size(stM(i));o=o+1) |
| | 1241 | { if (stM(i)[o]==" |
| | 1242 | ") |
| | 1243 | { links=stM(i)[1..o-1]; |
| | 1244 | rechts=stM(i)[o+1..size(stM(i))]; |
| | 1245 | stM(i)=links+rechts; |
| | 1246 | } |
| | 1247 | } |
| | 1248 | setring `newring`; |
| | 1249 | execute "matrix G(i)["+string(n)+"]["+string(n)+"]="+stM(i); |
| | 1250 | p=det(I-v1*G(i)); // denominator of new term - |
| | 1251 | A(2)[1,1]=A(2)[1,1]*p+A(2)[1,2]; // expanding A(2)[1,1]/A(2)[1,2] +1/p |
| | 1252 | A(2)[1,2]=A(2)[1,2]*p; |
| | 1253 | if (interval<>0) // canceling common terms of denominator |
| | 1254 | { if ((i/interval)*interval==i) // and enumerator |
| | 1255 | { |
| | 1256 | s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz() these |
| | 1257 | A(2)[1,1]=-s[2,1]; // three lines should be replaced by the |
| | 1258 | A(2)[1,2]=s[1,1]; // following three |
| | 1259 | // p=gcd(A(2)[1,1],A(2)[1,2]); |
| | 1260 | // A(2)[1,1]=A(2)[1,1]/p; |
| | 1261 | // A(2)[1,2]=A(2)[1,2]/p; |
| | 1262 | } |
| | 1263 | } |
| | 1264 | setring br; |
| | 1265 | } |
| | 1266 | } |
| | 1267 | int g=i; // G(1)..G(i) are generators without |
| | 1268 | // doubles - g generally is the number |
| | 1269 | // of elements in the group so far - |
| | 1270 | j=i; // j is the number of new elements that |
| | 1271 | // we use as factors |
| | 1272 | int k, m, l; |
| | 1273 | if (v) |
| | 1274 | { ""; |
| | 1275 | " Generating the entire matrix group. Whenever a new group element is found,"; |
| | 1276 | " the coressponding ring homomorphism of the Reynolds operator and the"; |
| | 1277 | " corresponding term of the Molien series is generated."; |
| | 1278 | ""; |
| | 1279 | } |
| | 1280 | // taking all elements in a ring of characteristic 0 and computing the terms |
| | 1281 | // of the Molien series there |
| | 1282 | while (1) |
| | 1283 | { l=0; // l is the number of products we get in |
| | 1284 | // one going |
| | 1285 | for (m=g-j+1;m<=g;m=m+1) |
| | 1286 | { for (k=1;k<=i;k=k+1) |
| | 1287 | { l=l+1; |
| | 1288 | matrix P(l)=G(k)*G(m); // possible new element |
| | 1289 | } |
| | 1290 | } |
| | 1291 | j=0; |
| | 1292 | for (k=1;k<=l;k=k+1) |
| | 1293 | { if (unique(G(1..g),P(k))) |
| | 1294 | { j=j+1; // a new factor for next run |
| | 1295 | g=g+1; |
| | 1296 | matrix G(g)=P(k); // a new group element - |
| | 1297 | A(1)=concat(A(1),G(g)*vars); // adding new mapping to A(1) |
| | 1298 | string stM(g)=string(G(g)); |
| | 1299 | for (o=1;o<=size(stM(g));o=o+1) |
| | 1300 | { if (stM(g)[o]==" |
| | 1301 | ") |
| | 1302 | { links=stM(g)[1..o-1]; |
| | 1303 | rechts=stM(g)[o+1..size(stM(g))]; |
| | 1304 | stM(g)=links+rechts; |
| | 1305 | } |
| | 1306 | } |
| | 1307 | setring `newring`; |
| | 1308 | execute "matrix G(g)["+string(n)+"]["+string(n)+"]="+stM(g); |
| | 1309 | p=det(I-v1*G(g)); // denominator of new term |
| | 1310 | A(2)[1,1]=A(2)[1,1]*p+A(2)[1,2]; |
| | 1311 | A(2)[1,2]=A(2)[1,2]*p; // expanding A(2)[1,1]/A(2)[1,2] + 1/p - |
| | 1312 | if (interval<>0) // canceling common terms of denominator |
| | 1313 | { if ((g/interval)*interval==g) // and enumerator |
| | 1314 | { |
| | 1315 | s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz() |
| | 1316 | A(2)[1,1]=-s[2,1]; // these three lines should be replaced |
| | 1317 | A(2)[1,2]=s[1,1]; // by the following three |
| | 1318 | // p=gcd(A(2)[1,1],A(2)[1,2]); |
| | 1319 | // A(2)[1,1]=A(2)[1,1]/p; |
| | 1320 | // A(2)[1,2]=A(2)[1,2]/p; |
| | 1321 | } |
| | 1322 | } |
| | 1323 | if (v) |
| | 1324 | { " Group element "+string(g)+" has been found."; |
| | 1325 | } |
| | 1326 | setring br; |
| | 1327 | } |
| | 1328 | kill P(k); |
| | 1329 | } |
| | 1330 | if (j==0) // when we didn't add any new elements |
| | 1331 | { break; // in one run through the while loop |
| | 1332 | } // we are done |
| | 1333 | } |
| | 1334 | if (v) |
| | 1335 | { if (g<=i) |
| | 1336 | { " There are only "+string(g)+" group elements."; |
| | 1337 | } |
| | 1338 | ""; |
| | 1339 | } |
| | 1340 | A(1)=transpose(A(1)); // when we evaluate the Reynolds operator |
| | 1341 | // later on, we actually want 1xn |
| | 1342 | // matrices |
| | 1343 | setring `newring`; |
| | 1344 | if (interval==0) // canceling common terms of denominator |
| | 1345 | { // and enumerator - |
| | 1346 | s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz() |
| | 1347 | A(2)[1,1]=-s[2,1]; // these three lines should be replaced |
| | 1348 | A(2)[1,2]=s[1,1]; // by the following three |
| | 1349 | // p=gcd(A(2)[1,1],A(2)[1,2]); |
| | 1350 | // A(2)[1,1]=A(2)[1,1]/p; |
| | 1351 | // A(2)[1,2]=A(2)[1,2]/p; |
| | 1352 | } |
| | 1353 | if (interval<>0) // canceling common terms of denominator |
| | 1354 | { if ((g/interval)*interval<>g) // and enumerator |
| | 1355 | { |
| | 1356 | s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz() |
| | 1357 | A(2)[1,1]=-s[2,1]; // these three lines should be replaced |
| | 1358 | A(2)[1,2]=s[1,1]; // by the following three |
| | 1359 | // p=gcd(A(2)[1,1],A(2)[1,2]); |
| | 1360 | // A(2)[1,1]=A(2)[1,1]/p; |
| | 1361 | // A(2)[1,2]=A(2)[1,2]/p; |
| | 1362 | } |
| | 1363 | } |
| | 1364 | map slead=`newring`,ideal(0); |
| | 1365 | s=slead(A(2)); |
| | 1366 | A(2)[1,1]=1/s[1,1]*A(2)[1,1]; // numerator and denominator have to have |
| | 1367 | A(2)[1,2]=1/s[1,2]*A(2)[1,2]; // a constant term of 1 |
| | 1368 | if (v) |
| | 1369 | { " Now we are done calculating Molien series and Reynolds operator."; |
| | 1370 | ""; |
| | 1371 | } |
| | 1372 | matrix M=A(2); |
| | 1373 | kill G(1..g), s, slead, p, v1, I, A(2); |
| | 1374 | export `newring`; // we keep the ring where we computed the |
| | 1375 | export M; // the Molien series such that we can |
| | 1376 | setring br; // keep it |
| | 1377 | return(A(1)); |
| | 1378 | } |
| 1077 | | proc group (list #) |
| 1078 | | { matrix G(1)=#[1]; // first group element |
| 1079 | | int i=1; |
| 1080 | | for (int j=2;j<=size(#);j=j+1) // throwing out doubles among the |
| 1081 | | { if (unique(G(1..i),#[j])) // generators |
| 1082 | | { i=i+1; |
| 1083 | | matrix G(i)=#[j]; |
| 1084 | | } |
| 1085 | | } |
| 1086 | | int g=i; // g: elements in the group so far, i: |
| 1087 | | j=i; // generators, j: new ones used as |
| 1088 | | int m, k, l; // as factors, l: counting possible new |
| 1089 | | // new elements |
| 1090 | | while (1) |
| 1091 | | { l=0; |
| 1092 | | for (m=g-j+1;m<=g;m=m+1) |
| 1093 | | { for (k=1;k<=i;k=k+1) |
| 1094 | | { l=l+1; |
| 1095 | | matrix P(l)=G(k)*G(m); // possible new element |
| 1096 | | } |
| 1097 | | } |
| 1098 | | j=0; |
| 1099 | | for (k=1;k<=l;k=k+1) // checking whether the P(k) are new |
| 1100 | | { if (unique(G(1..g),P(k))) |
| 1101 | | { j=j+1; |
| 1102 | | g=g+1; |
| 1103 | | matrix G(g)=P(k); // adding new elements - |
| 1104 | | } |
| 1105 | | kill P(k); |
| 1106 | | } |
| 1107 | | if (j==0) // when we didn't add any new elements |
| 1108 | | { break; // in one run through the while loop, we |
| 1109 | | } // are done |
| 1110 | | } |
| 1111 | | return(g); |
| | 1938 | proc p_search (int n,int d,ideal B,int cd,ideal P,ideal sP,int i,int dif,int dB,ideal CI) |
| | 1939 | { def br=basering; |
| | 1940 | degBound=0; |
| | 1941 | matrix vec(1)[1][cd]; // starting with 0-vector - |
| | 1942 | intmat new[1][cd]; // the coefficients for the next |
| | 1943 | // combination - |
| | 1944 | matrix pnew[1][cd]; // new needs to be mapped into br - |
| | 1945 | int counter=1; // counts the vectors |
| | 1946 | int j; |
| | 1947 | int p=char(br); |
| | 1948 | if (minpoly<>0) |
| | 1949 | { int ext_deg=pardeg(minpoly); // field has p^d elements |
| | 1950 | } |
| | 1951 | else |
| | 1952 | { int ext_deg=1; // field has p^d elements |
| | 1953 | } |
| | 1954 | poly test_poly; // the linear combination to test |
| | 1955 | int test_dim; |
| | 1956 | ring R=0,x,dp; // just to calculate next variable |
| | 1957 | // bound - |
| | 1958 | number bound=(number(p)^(ext_deg*cd)-1)/(number(p)^ext_deg-1)+1; // this is |
| | 1959 | // how many linearly independent vectors |
| | 1960 | // of size cd exist having entries in the |
| | 1961 | // base field of br |
| | 1962 | setring br; |
| | 1963 | intvec h; // Hilbert series |
| | 1964 | int k=i+1; |
| | 1965 | matrix tB=transpose(B); |
| | 1966 | ideal TEST; |
| | 1967 | while (k<=i+dif) |
| | 1968 | { CI=CI+ideal(var(k)^d); // homogeneous polynomial of the same |
| | 1969 | // degree as the one we're looking for is |
| | 1970 | // added |
| | 1971 | // h=hilb(std(CI),1); |
| | 1972 | dB=dB+d-1; // used as degBound |
| | 1973 | setring R; |
| | 1974 | while (number(counter)<>bound) // otherwise, we are done |
| | 1975 | { setring br; |
| | 1976 | new=next_vector(new); |
| | 1977 | for (j=1;j<=cd;j=j+1) |
| | 1978 | { pnew[1,j]=int_number_map(new[1,j]); // mapping an integer into br |
| | 1979 | } |
| | 1980 | if (unique(vec(1..counter),pnew)) // checking whether we tried pnew before |
| | 1981 | { counter=counter+1; |
| | 1982 | matrix vec(counter)=pnew; // keeping track of the ones we tried - |
| | 1983 | test_poly=(vec(counter)*tB)[1,1]; // linear combination - |
| | 1984 | // degBound=dB; |
| | 1985 | TEST=sP+ideal(test_poly); |
| | 1986 | attrib(TEST,"isSB",1); |
| | 1987 | test_dim=dim(TEST); |
| | 1988 | // degBound=0; |
| | 1989 | if (n-test_dim==k) // the dimension has been lowered by one |
| | 1990 | { sP=TEST; |
| | 1991 | setring R; |
| | 1992 | break; |
| | 1993 | } |
| | 1994 | // degBound=dB; |
| | 1995 | TEST=std(sP+ideal(test_poly)); // should soon to be replaced by next |
| | 1996 | // line |
| | 1997 | // TEST=std(sP,test_poly,h); // Hilbert driven std-calculation |
| | 1998 | test_dim=dim(TEST); |
| | 1999 | // degBound=0; |
| | 2000 | if (n-test_dim==k) // the dimension has been lowered by one |
| | 2001 | { sP=TEST; |
| | 2002 | setring R; |
| | 2003 | break; |
| | 2004 | } |
| | 2005 | } |
| | 2006 | setring R; |
| | 2007 | } |
| | 2008 | if (number(counter)<=bound) |
| | 2009 | { setring br; |
| | 2010 | P[k]=test_poly; // test_poly ist added to primary |
| | 2011 | } // invariants |
| | 2012 | else |
| | 2013 | { setring br; |
| | 2014 | CI=CI[1..size(CI)-1]; |
| | 2015 | return(P,sP,CI,dB-d+1); |
| | 2016 | } |
| | 2017 | k=k+1; |
| | 2018 | } |
| | 2019 | return(P,sP,CI,dB); |
| | 2020 | } |
| | 2021 | |
| | 2022 | proc primary_char0 (matrix REY,matrix M,list #) |
| | 2023 | USAGE: primary_char0(REY,M[,v]); |
| | 2024 | REY: a <matrix> representing the Reynolds operator, M: a 1x2 <matrix> |
| | 2025 | representing the Molien series, v: an optional <int> |
| | 2026 | ASSUME: REY is the first return value of group_reynolds or reynolds_molien and |
| | 2027 | M the one of molien or the second one of reynolds_molien |
| | 2028 | DISPLAY: information about the various stages of the programme if v does not |
| | 2029 | equal 0 |
| | 2030 | RETURN: primary invariants (type <matrix>) of the invariant ring |
| | 2031 | EXAMPLE: example primary_char0; shows an example |
| | 2032 | THEORY: Bases of homogeneous invariants are generated successively and those |
| | 2033 | are chosen as primary invariants that lower the dimension of the ideal |
| | 2034 | generated by the previously found invariants (see paper "Generating a |
| | 2035 | Noetherian Normalization of the Invariant Ring of a Finite Group" by |
| | 2036 | Decker, Heydtmann, Schreyer (1997) to appear in JSC). |
| | 2037 | { degBound=0; |
| | 2038 | if (char(basering)<>0) |
| | 2039 | { "ERROR: primary_char0 should only be used with rings of characteristic 0."; |
| | 2040 | return(); |
| | 2041 | } |
| | 2042 | //----------------- checking input and setting verbose mode ------------------ |
| | 2043 | if (size(#)>1) |
| | 2044 | { "ERROR: primary_char0 can only have three parameters."; |
| | 2045 | return(); |
| | 2046 | } |
| | 2047 | if (size(#)==1) |
| | 2048 | { if (typeof(#[1])<>"int") |
| | 2049 | { "ERROR: The third parameter should be of type <int>."; |
| | 2050 | return(); |
| | 2051 | } |
| | 2052 | else |
| | 2053 | { int v=#[1]; |
| | 2054 | } |
| | 2055 | } |
| | 2056 | else |
| | 2057 | { int v=0; |
| | 2058 | } |
| | 2059 | int n=nvars(basering); // n is the number of variables, as well |
| | 2060 | // as the size of the matrices, as well |
| | 2061 | // as the number of primary invariants, |
| | 2062 | // we should get |
| | 2063 | if (ncols(REY)<>n) |
| | 2064 | { "ERROR: First parameter ought to be the Reynolds operator." |
| | 2065 | return(); |
| | 2066 | } |
| | 2067 | if (ncols(M)<>2 or nrows(M)<>1) |
| | 2068 | { "ERROR: Second parameter ought to be the Molien series." |
| | 2069 | return(); |
| | 2070 | } |
| | 2071 | //---------------------------------------------------------------------------- |
| | 2072 | if (v && voice<>2) |
| | 2073 | { " We can start looking for primary invariants..."; |
| | 2074 | ""; |
| | 2075 | } |
| | 2076 | if (v && voice==2) |
| | 2077 | { ""; |
| | 2078 | } |
| | 2079 | //------------------------- initializing variables --------------------------- |
| | 2080 | int dB; |
| | 2081 | poly p(1..2); // p(1) will be used for single terms of |
| | 2082 | // the partial expansion, p(2) to store |
| | 2083 | p(1..2)=partial_molien(M,1); // the intermediate result - |
| | 2084 | poly v1=var(1); // we need v1 to split off coefficients |
| | 2085 | // in the partial expansion of M (which |
| | 2086 | // is in terms of the first variable) - |
| | 2087 | int j,d,cd,newdim,dif; // d: current degree, cd: dimension of |
| | 2088 | // space of invariants of degree d, |
| | 2089 | // newdim: dimension the ideal generated |
| | 2090 | // the primary invariants plus basis |
| | 2091 | // elements, dif=n-i-newdim, i.e. the |
| | 2092 | // number of new primary invairants that |
| | 2093 | // should be added in this degree - |
| | 2094 | ideal P,Pplus,sPplus,CI,B; // P: will contain primary invariants, |
| | 2095 | // Pplus: P+B, CI: a complete |
| | 2096 | // intersection with the same Hilbert |
| | 2097 | // function as P |
| | 2098 | ideal sP=std(P); |
| | 2099 | dB=1; // used as degree bound |
| | 2100 | int i=0; |
| | 2101 | //-------------- loop that searches for primary invariants ------------------ |
| | 2102 | while(1) // repeat until n primary invariants are |
| | 2103 | { // found - |
| | 2104 | p(1..2)=partial_molien(M,1,p(2)); // next term of the partial expansion - |
| | 2105 | d=deg(p(1)); // degree where we'll search - |
| | 2106 | cd=int(coef(p(1),v1)[2,1]); // dimension of the homogeneous space of |
| | 2107 | // inviarants of degree d |
| | 2108 | if (v) |
| | 2109 | { " Computing primary invariants in degree "+string(d)+":"; |
| | 2110 | } |
| | 2111 | B=invariant_basis_reynolds(REY,d,intvec(cd,6)); // basis of invariants of |
| | 2112 | // degree d |
| | 2113 | if (B[1]<>0) |
| | 2114 | { Pplus=P+B; |
| | 2115 | sPplus=std(Pplus); |
| | 2116 | newdim=dim(sPplus); |
| | 2117 | dif=n-i-newdim; |
| | 2118 | } |
| | 2119 | else |
| | 2120 | { dif=0; |
| | 2121 | } |
| | 2122 | if (dif<>0) // we have to find dif new primary |
| | 2123 | { // invariants |
| | 2124 | if (cd<>dif) |
| | 2125 | { P,sP,CI,dB=search(n,d,B,cd,P,sP,i,dif,dB,CI); // searching for dif invariants |
| | 2126 | } // i.e. we can take all of B |
| | 2127 | else |
| | 2128 | { for(j=i+1;j>i+dif;j=j+1) |
| | 2129 | { CI=CI+ideal(var(j)^d); |
| | 2130 | } |
| | 2131 | dB=dB+dif*(d-1); |
| | 2132 | P=Pplus; |
| | 2133 | sP=sPplus; |
| | 2134 | } |
| | 2135 | if (v) |
| | 2136 | { for (j=1;j<=dif;j=j+1) |
| | 2137 | { " We find: "+string(P[i+j]); |
| | 2138 | } |
| | 2139 | } |
| | 2140 | i=i+dif; |
| | 2141 | if (i==n) // found all primary invariants |
| | 2142 | { if (v) |
| | 2143 | { ""; |
| | 2144 | " We found all primary invariants."; |
| | 2145 | ""; |
| | 2146 | } |
| | 2147 | return(matrix(P)); |
| | 2148 | } |
| | 2149 | } // done with degree d |
| | 2150 | } |
| | 2151 | } |
| | 2152 | example |
| | 2153 | { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; |
| | 2154 | echo=2; |
| | 2155 | ring R=0,(x,y,z),dp; |
| | 2156 | matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; |
| | 2157 | matrix REY,M=reynolds_molien(A); |
| | 2158 | matrix P=primary_char0(REY,M); |
| | 2159 | print(P); |
| | 2160 | } |
| | 2161 | |
| | 2162 | proc primary_charp (matrix REY,string ring_name,list #) |
| | 2163 | USAGE: primary_charp(REY,ringname[,v]); |
| | 2164 | REY: a <matrix> representing the Reynolds operator, ringname: a |
| | 2165 | <string> giving the name of a ring where the Molien series is stored, |
| | 2166 | v: an optional <int> |
| | 2167 | ASSUME: REY is the first return value of group_reynolds or reynolds_molien and |
| | 2168 | ringname gives the name of a ring of characteristic 0 that has been |
| | 2169 | created by molien or reynolds_molien |
| | 2170 | DISPLAY: information about the various stages of the programme if v does not |
| | 2171 | equal 0 |
| | 2172 | RETURN: primary invariants (type <matrix>) of the invariant ring |
| | 2173 | EXAMPLE: example primary_charp; shows an example |
| | 2174 | THEORY: Bases of homogeneous invariants are generated successively and those |
| | 2175 | are chosen as primary invariants that lower the dimension of the ideal |
| | 2176 | generated by the previously found invariants (see paper "Generating a |
| | 2177 | Noetherian Normalization of the Invariant Ring of a Finite Group" by |
| | 2178 | Decker, Heydtmann, Schreyer (1997) to appear in JSC). |
| | 2179 | { degBound=0; |
| | 2180 | // ---------------- checking input and setting verbose mode ------------------- |
| | 2181 | if (char(basering)==0) |
| | 2182 | { "ERROR: primary_charp should only be used with rings of characteristic p>0."; |
| | 2183 | return(); |
| | 2184 | } |
| | 2185 | if (size(#)>1) |
| | 2186 | { "ERROR: primary_charp can only have three parameters."; |
| | 2187 | return(); |
| | 2188 | } |
| | 2189 | if (size(#)==1) |
| | 2190 | { if (typeof(#[1])<>"int") |
| | 2191 | { "ERROR: The third parameter should be of type <int>."; |
| | 2192 | return(); |
| | 2193 | } |
| | 2194 | else |
| | 2195 | { int v=#[1]; |
| | 2196 | } |
| | 2197 | } |
| | 2198 | else |
| | 2199 | { int v=0; |
| | 2200 | } |
| | 2201 | def br=basering; |
| | 2202 | int n=nvars(br); // n is the number of variables, as well |
| | 2203 | // as the size of the matrices, as well |
| | 2204 | // as the number of primary invariants, |
| | 2205 | // we should get |
| | 2206 | if (ncols(REY)<>n) |
| | 2207 | { "ERROR: First parameter ought to be the Reynolds operator." |
| | 2208 | return(); |
| | 2209 | } |
| | 2210 | if (typeof(`ring_name`)<>"ring") |
| | 2211 | { "ERROR: Second parameter ought to the name of a ring where the Molien"; |
| | 2212 | " is stored."; |
| | 2213 | return(); |
| | 2214 | } |
| | 2215 | //---------------------------------------------------------------------------- |
| | 2216 | if (v && voice<>2) |
| | 2217 | { " We can start looking for primary invariants..."; |
| | 2218 | ""; |
| | 2219 | } |
| | 2220 | if (v && voice==2) |
| | 2221 | { ""; |
| | 2222 | } |
| | 2223 | //----------------------- initializing variables ----------------------------- |
| | 2224 | int dB; |
| | 2225 | setring `ring_name`; // the Molien series is stores here - |
| | 2226 | poly p(1..2); // p(1) will be used for single terms of |
| | 2227 | // the partial expansion, p(2) to store |
| | 2228 | p(1..2)=partial_molien(M,1); // the intermediate result - |
| | 2229 | poly v1=var(1); // we need v1 to split off coefficients |
| | 2230 | // in the partial expansion of M (which |
| | 2231 | // is in terms of the first variable) |
| | 2232 | setring br; |
| | 2233 | int j,d,cd,newdim,dif; // d: current degree, cd: dimension of |
| | 2234 | // space of invariants of degree d, |
| | 2235 | // newdim: dimension the ideal generated |
| | 2236 | // the primary invariants plus basis |
| | 2237 | // elements, dif=n-i-newdim, i.e. the |
| | 2238 | // number of new primary invairants that |
| | 2239 | // should be added in this degree - |
| | 2240 | ideal P,Pplus,sPplus,CI,B; // P: will contain primary invariants, |
| | 2241 | // Pplus: P+B, CI: a complete |
| | 2242 | // intersection with the same Hilbert |
| | 2243 | // function as P |
| | 2244 | ideal sP=std(P); |
| | 2245 | dB=1; // used as degree bound |
| | 2246 | int i=0; |
| | 2247 | //---------------- loop that searches for primary invariants ----------------- |
| | 2248 | while(1) // repeat until n primary invariants are |
| | 2249 | { // found |
| | 2250 | setring `ring_name`; |
| | 2251 | p(1..2)=partial_molien(M,1,p(2)); // next term of the partial expansion - |
| | 2252 | d=deg(p(1)); // degree where we'll search - |
| | 2253 | cd=int(coef(p(1),v1)[2,1]); // dimension of the homogeneous space of |
| | 2254 | // inviarants of degree d |
| | 2255 | setring br; |
| | 2256 | if (v) |
| | 2257 | { " Computing primary invariants in degree "+string(d)+":"; |
| | 2258 | } |
| | 2259 | B=invariant_basis_reynolds(REY,d,intvec(cd,6)); // basis of invariants of |
| | 2260 | // degree d |
| | 2261 | if (B[1]<>0) |
| | 2262 | { Pplus=P+B; |
| | 2263 | sPplus=std(Pplus); |
| | 2264 | newdim=dim(sPplus); |
| | 2265 | dif=n-i-newdim; |
| | 2266 | } |
| | 2267 | else |
| | 2268 | { dif=0; |
| | 2269 | } |
| | 2270 | if (dif<>0) // we have to find dif new primary |
| | 2271 | { // invariants |
| | 2272 | if (cd<>dif) |
| | 2273 | { P,sP,CI,dB=p_search(n,d,B,cd,P,sP,i,dif,dB,CI); |
| | 2274 | } |
| | 2275 | else // i.e. we can take all of B |
| | 2276 | { for(j=i+1;j>i+dif;j=j+1) |
| | 2277 | { CI=CI+ideal(var(j)^d); |
| | 2278 | } |
| | 2279 | dB=dB+dif*(d-1); |
| | 2280 | P=Pplus; |
| | 2281 | sP=sPplus; |
| | 2282 | } |
| | 2283 | if (v) |
| | 2284 | { for (j=1;j<=size(P)-i;j=j+1) |
| | 2285 | { " We find: "+string(P[i+j]); |
| | 2286 | } |
| | 2287 | } |
| | 2288 | i=size(P); |
| | 2289 | if (i==n) // found all primary invariants |
| | 2290 | { if (v) |
| | 2291 | { ""; |
| | 2292 | " We found all primary invariants."; |
| | 2293 | ""; |
| | 2294 | } |
| | 2295 | return(matrix(P)); |
| | 2296 | } |
| | 2297 | } // done with degree d |
| | 2298 | } |
| | 2299 | } |
| | 2300 | example |
| | 2301 | { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed into"; |
| | 2302 | " characteristic 3)"; |
| | 2303 | echo=2; |
| | 2304 | ring R=3,(x,y,z),dp; |
| | 2305 | matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; |
| | 2306 | list L=group_reynolds(A); |
| | 2307 | string newring="alskdfj"; |
| | 2308 | molien(L[2..size(L)],newring); |
| | 2309 | matrix P=primary_charp(L[1],newring); |
| | 2310 | kill `newring`; |
| | 2311 | print(P); |
| | 2312 | } |
| | 2313 | |
| | 2314 | proc primary_char0_no_molien (matrix REY, list #) |
| | 2315 | USAGE: primary_char0_no_molien(REY[,v]); |
| | 2316 | REY: a <matrix> representing the Reynolds operator, v: an optional |
| | 2317 | <int> |
| | 2318 | ASSUME: REY is the first return value of group_reynolds or reynolds_molien |
| | 2319 | DISPLAY: information about the various stages of the programme if v does not |
| | 2320 | equal 0 |
| | 2321 | RETURN: primary invariants (type <matrix>) of the invariant ring and an |
| | 2322 | <intvec> listing some of the degrees where no non-trivial homogeneous |
| | 2323 | invariants are to be found |
| | 2324 | EXAMPLE: example primary_char0_no_molien; shows an example |
| | 2325 | THEORY: Bases of homogeneous invariants are generated successively and those |
| | 2326 | are chosen as primary invariants that lower the dimension of the ideal |
| | 2327 | generated by the previously found invariants (see paper "Generating a |
| | 2328 | Noetherian Normalization of the Invariant Ring of a Finite Group" by |
| | 2329 | Decker, Heydtmann, Schreyer (1997) to appear in JSC). |
| | 2330 | { degBound=0; |
| | 2331 | //-------------- checking input and setting verbose mode --------------------- |
| | 2332 | if (char(basering)<>0) |
| | 2333 | { "ERROR: primary_char0_no_molien should only be used with rings of"; |
| | 2334 | " characteristic 0."; |
| | 2335 | return(); |
| | 2336 | } |
| | 2337 | if (size(#)>1) |
| | 2338 | { "ERROR: primary_char0_no_molien can only have two parameters."; |
| | 2339 | return(); |
| | 2340 | } |
| | 2341 | if (size(#)==1) |
| | 2342 | { if (typeof(#[1])<>"int") |
| | 2343 | { "ERROR: The second parameter should be of type <int>."; |
| | 2344 | return(); |
| | 2345 | } |
| | 2346 | else |
| | 2347 | { int v=#[1]; |
| | 2348 | } |
| | 2349 | } |
| | 2350 | else |
| | 2351 | { int v=0; |
| | 2352 | } |
| | 2353 | int n=nvars(basering); // n is the number of variables, as well |
| | 2354 | // as the size of the matrices, as well |
| | 2355 | // as the number of primary invariants, |
| | 2356 | // we should get |
| | 2357 | if (ncols(REY)<>n) |
| | 2358 | { "ERROR: First parameter ought to be the Reynolds operator." |
| | 2359 | return(); |
| | 2360 | } |
| | 2361 | //---------------------------------------------------------------------------- |
| | 2362 | if (v && voice<>2) |
| | 2363 | { " We can start looking for primary invariants..."; |
| | 2364 | ""; |
| | 2365 | } |
| | 2366 | if (v && voice==2) |
| | 2367 | { ""; |
| | 2368 | } |
| | 2369 | //----------------------- initializing variables ----------------------------- |
| | 2370 | int dB; |
| | 2371 | int j,d,cd,newdim,dif; // d: current degree, cd: dimension of |
| | 2372 | // space of invariants of degree d, |
| | 2373 | // newdim: dimension the ideal generated |
| | 2374 | // the primary invariants plus basis |
| | 2375 | // elements, dif=n-i-newdim, i.e. the |
| | 2376 | // number of new primary invairants that |
| | 2377 | // should be added in this degree - |
| | 2378 | ideal P,Pplus,CI,B; // P: will contain primary invariants, |
| | 2379 | // Pplus: P+B, CI: a complete |
| | 2380 | // intersection with the same Hilbert |
| | 2381 | // function as P |
| | 2382 | ideal sP=std(P); |
| | 2383 | dB=1; // used as degree bound - |
| | 2384 | d=0; // initializing |
| | 2385 | int i=0; |
| | 2386 | intvec deg_vector; |
| | 2387 | //------------------ loop that searches for primary invariants --------------- |
| | 2388 | while(1) // repeat until n primary invariants are |
| | 2389 | { // found - |
| | 2390 | d=d+1; // degree where we'll search |
| | 2391 | if (v) |
| | 2392 | { " Computing primary invariants in degree "+string(d)+":"; |
| | 2393 | } |
| | 2394 | B=invariant_basis_reynolds(REY,d,intvec(-1,6)); // basis of invariants of |
| | 2395 | // degree d |
| | 2396 | if (B[1]<>0) |
| | 2397 | { Pplus=P+B; |
| | 2398 | newdim=dim(std(Pplus)); |
| | 2399 | dif=n-i-newdim; |
| | 2400 | } |
| | 2401 | else |
| | 2402 | { dif=0; |
| | 2403 | deg_vector=deg_vector,d; |
| | 2404 | } |
| | 2405 | if (dif<>0) // we have to find dif new primary |
| | 2406 | { // invariants |
| | 2407 | cd=size(B); |
| | 2408 | if (cd<>dif) |
| | 2409 | { P,sP,CI,dB=search(n,d,B,cd,P,sP,i,dif,dB,CI); |
| | 2410 | } |
| | 2411 | else // i.e. we can take all of B |
| | 2412 | { for(j=i+1;j<=i+dif;j=j+1) |
| | 2413 | { CI=CI+ideal(var(j)^d); |
| | 2414 | } |
| | 2415 | dB=dB+dif*(d-1); |
| | 2416 | P=Pplus; |
| | 2417 | sP=std(P); |
| | 2418 | } |
| | 2419 | if (v) |
| | 2420 | { for (j=1;j<=dif;j=j+1) |
| | 2421 | { " We find: "+string(P[i+j]); |
| | 2422 | } |
| | 2423 | } |
| | 2424 | i=i+dif; |
| | 2425 | if (i==n) // found all primary invariants |
| | 2426 | { if (v) |
| | 2427 | { ""; |
| | 2428 | " We found all primary invariants."; |
| | 2429 | ""; |
| | 2430 | } |
| | 2431 | if (deg_vector==0) |
| | 2432 | { return(matrix(P)); |
| | 2433 | } |
| | 2434 | else |
| | 2435 | { return(matrix(P),compress(deg_vector)); |
| | 2436 | } |
| | 2437 | } |
| | 2438 | } // done with degree d |
| | 2439 | else |
| | 2440 | { if (v) |
| | 2441 | { " None here..."; |
| | 2442 | } |
| | 2443 | } |
| | 2444 | } |
| | 2445 | } |
| | 2446 | example |
| | 2447 | { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; |
| | 2448 | echo=2; |
| | 2449 | ring R=0,(x,y,z),dp; |
| | 2450 | matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; |
| | 2451 | list L=group_reynolds(A); |
| | 2452 | list l=primary_char0_no_molien(L[1]); |
| | 2453 | print(l[1]); |
| | 2454 | } |
| | 2455 | |
| | 2456 | proc primary_charp_no_molien (matrix REY, list #) |
| | 2457 | USAGE: primary_charp_no_molien(REY[,v]); |
| | 2458 | REY: a <matrix> representing the Reynolds operator, v: an optional |
| | 2459 | <int> |
| | 2460 | ASSUME: REY is the first return value of group_reynolds or reynolds_molien |
| | 2461 | DISPLAY: information about the various stages of the programme if v does not |
| | 2462 | equal 0 |
| | 2463 | RETURN: primary invariants (type <matrix>) of the invariant ring and an |
| | 2464 | <intvec> listing some of the degrees where no non-trivial homogeneous |
| | 2465 | invariants are to be found |
| | 2466 | EXAMPLE: example primary_charp_no_molien; shows an example |
| | 2467 | THEORY: Bases of homogeneous invariants are generated successively and those |
| | 2468 | are chosen as primary invariants that lower the dimension of the ideal |
| | 2469 | generated by the previously found invariants (see paper "Generating a |
| | 2470 | Noetherian Normalization of the Invariant Ring of a Finite Group" by |
| | 2471 | Decker, Heydtmann, Schreyer (1997) to appear in JSC). |
| | 2472 | { degBound=0; |
| | 2473 | //----------------- checking input and setting verbose mode ------------------ |
| | 2474 | if (char(basering)==0) |
| | 2475 | { "ERROR: primary_charp_no_molien should only be used with rings of"; |
| | 2476 | " characteristic p>0."; |
| | 2477 | return(); |
| | 2478 | } |
| | 2479 | if (size(#)>1) |
| | 2480 | { "ERROR: primary_charp_no_molien can only have two parameters."; |
| | 2481 | return(); |
| | 2482 | } |
| | 2483 | if (size(#)==1) |
| | 2484 | { if (typeof(#[1])<>"int") |
| | 2485 | { "ERROR: The second parameter should be of type <int>."; |
| | 2486 | return(); |
| | 2487 | } |
| | 2488 | else |
| | 2489 | { int v=#[1]; |
| | 2490 | } |
| | 2491 | } |
| | 2492 | else |
| | 2493 | { int v=0; |
| | 2494 | } |
| | 2495 | int n=nvars(basering); // n is the number of variables, as well |
| | 2496 | // as the size of the matrices, as well |
| | 2497 | // as the number of primary invariants, |
| | 2498 | // we should get |
| | 2499 | if (ncols(REY)<>n) |
| | 2500 | { "ERROR: First parameter ought to be the Reynolds operator." |
| | 2501 | return(); |
| | 2502 | } |
| | 2503 | //---------------------------------------------------------------------------- |
| | 2504 | if (v && voice<>2) |
| | 2505 | { " We can start looking for primary invariants..."; |
| | 2506 | ""; |
| | 2507 | } |
| | 2508 | if (v && voice==2) |
| | 2509 | { ""; |
| | 2510 | } |
| | 2511 | //-------------------- initializing variables -------------------------------- |
| | 2512 | int dB; |
| | 2513 | int j,d,cd,newdim,dif; // d: current degree, cd: dimension of |
| | 2514 | // space of invariants of degree d, |
| | 2515 | // newdim: dimension the ideal generated |
| | 2516 | // the primary invariants plus basis |
| | 2517 | // elements, dif=n-i-newdim, i.e. the |
| | 2518 | // number of new primary invairants that |
| | 2519 | // should be added in this degree - |
| | 2520 | ideal P,Pplus,sPplus,CI,B; // P: will contain primary invariants, |
| | 2521 | // Pplus: P+B, CI: a complete |
| | 2522 | // intersection with the same Hilbert |
| | 2523 | // function as P |
| | 2524 | ideal sP=std(P); |
| | 2525 | dB=1; // used as degree bound - |
| | 2526 | d=0; // initializing |
| | 2527 | int i=0; |
| | 2528 | intvec deg_vector; |
| | 2529 | //------------------ loop that searches for primary invariants --------------- |
| | 2530 | while(1) // repeat until n primary invariants are |
| | 2531 | { // found - |
| | 2532 | d=d+1; // degree where we'll search |
| | 2533 | if (v) |
| | 2534 | { " Computing primary invariants in degree "+string(d)+":"; |
| | 2535 | } |
| | 2536 | B=invariant_basis_reynolds(REY,d,intvec(-1,6)); // basis of invariants of |
| | 2537 | // degree d |
| | 2538 | if (B[1]<>0) |
| | 2539 | { Pplus=P+B; |
| | 2540 | sPplus=std(Pplus); |
| | 2541 | newdim=dim(sPplus); |
| | 2542 | dif=n-i-newdim; |
| | 2543 | } |
| | 2544 | else |
| | 2545 | { dif=0; |
| | 2546 | deg_vector=deg_vector,d; |
| | 2547 | } |
| | 2548 | if (dif<>0) // we have to find dif new primary |
| | 2549 | { // invariants |
| | 2550 | cd=size(B); |
| | 2551 | if (cd<>dif) |
| | 2552 | { P,sP,CI,dB=p_search(n,d,B,cd,P,sP,i,dif,dB,CI); |
| | 2553 | } |
| | 2554 | else // i.e. we can take all of B |
| | 2555 | { for(j=i+1;j<=i+dif;j=j+1) |
| | 2556 | { CI=CI+ideal(var(j)^d); |
| | 2557 | } |
| | 2558 | dB=dB+dif*(d-1); |
| | 2559 | P=Pplus; |
| | 2560 | sP=sPplus; |
| | 2561 | } |
| | 2562 | if (v) |
| | 2563 | { for (j=1;j<=size(P)-i;j=j+1) |
| | 2564 | { " We find: "+string(P[i+j]); |
| | 2565 | } |
| | 2566 | } |
| | 2567 | i=size(P); |
| | 2568 | if (i==n) // found all primary invariants |
| | 2569 | { if (v) |
| | 2570 | { ""; |
| | 2571 | " We found all primary invariants."; |
| | 2572 | ""; |
| | 2573 | } |
| | 2574 | if (deg_vector==0) |
| | 2575 | { return(matrix(P)); |
| | 2576 | } |
| | 2577 | else |
| | 2578 | { return(matrix(P),compress(deg_vector)); |
| | 2579 | } |
| | 2580 | } |
| | 2581 | } // done with degree d |
| | 2582 | else |
| | 2583 | { if (v) |
| | 2584 | { " None here..."; |
| | 2585 | } |
| | 2586 | } |
| | 2587 | } |
| | 2588 | } |
| | 2589 | example |
| | 2590 | { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed into"; |
| | 2591 | " characteristic 3)"; |
| | 2592 | echo=2; |
| | 2593 | ring R=3,(x,y,z),dp; |
| | 2594 | matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; |
| | 2595 | list L=group_reynolds(A); |
| | 2596 | list l=primary_charp_no_molien(L[1]); |
| | 2597 | print(l[1]); |
| | 2598 | } |
| | 2599 | |
| | 2600 | proc primary_charp_without (list #) |
| | 2601 | USAGE: primary_charp_without(G1,G2,...[,v]); |
| | 2602 | G1,G2,...: <matrices> generating a finite matrix group, v: an optional |
| | 2603 | <int> |
| | 2604 | DISPLAY: information about the various stages of the programme if v does not |
| | 2605 | equal 0 |
| | 2606 | RETURN: primary invariants (type <matrix>) of the invariant ring |
| | 2607 | EXAMPLE: example primary_charp_without; shows an example |
| | 2608 | THEORY: Bases of homogeneous invariants are generated successively and those |
| | 2609 | are chosen as primary invariants that lower the dimension of the ideal |
| | 2610 | generated by the previously found invariants (see paper "Generating a |
| | 2611 | Noetherian Normalization of the Invariant Ring of a Finite Group" by |
| | 2612 | Decker, Heydtmann, Schreyer (1997) to appear in JSC). No Reynolds |
| | 2613 | operator or Molien series is used. |
| | 2614 | { degBound=0; |
| | 2615 | //--------------------- checking input and setting verbose mode -------------- |
| | 2616 | if (char(basering)==0) |
| | 2617 | { "ERROR: primary_charp_without should only be used with rings of"; |
| | 2618 | " characteristic 0."; |
| | 2619 | return(); |
| | 2620 | } |
| | 2621 | if (size(#)==0) |
| | 2622 | { "ERROR: There are no parameters."; |
| | 2623 | return(); |
| | 2624 | } |
| | 2625 | if (typeof(#[size(#)])=="int") |
| | 2626 | { int v=#[size(#)]; |
| | 2627 | int gen_num=size(#)-1; |
| | 2628 | if (gen_num==0) |
| | 2629 | { "ERROR: There are no generators of a finite matrix group given."; |
| | 2630 | return(); |
| | 2631 | } |
| | 2632 | } |
| | 2633 | else |
| | 2634 | { int v=0; |
| | 2635 | int gen_num=size(#); |
| | 2636 | } |
| | 2637 | int n=nvars(basering); // n is the number of variables, as well |
| | 2638 | // as the size of the matrices, as well |
| | 2639 | // as the number of primary invariants, |
| | 2640 | // we should get |
| | 2641 | for (int i=1;i<=gen_num;i=i+1) |
| | 2642 | { if (typeof(#[i])=="matrix") |
| | 2643 | { if (nrows(#[i])<>n or ncols(#[i])<>n) |
| | 2644 | { "ERROR: The number of variables of the base ring needs to be the same"; |
| | 2645 | " as the dimension of the square matrices"; |
| | 2646 | return(); |
| | 2647 | } |
| | 2648 | } |
| | 2649 | else |
| | 2650 | { "ERROR: The first parameters should be a list of matrices"; |
| | 2651 | return(); |
| | 2652 | } |
| | 2653 | } |
| | 2654 | //---------------------------------------------------------------------------- |
| | 2655 | if (v && voice==2) |
| | 2656 | { ""; |
| | 2657 | } |
| | 2658 | //---------------------------- initializing variables ------------------------ |
| | 2659 | int dB; |
| | 2660 | int j,d,cd,newdim,dif; // d: current degree, cd: dimension of |
| | 2661 | // space of invariants of degree d, |
| | 2662 | // newdim: dimension the ideal generated |
| | 2663 | // the primary invariants plus basis |
| | 2664 | // elements, dif=n-i-newdim, i.e. the |
| | 2665 | // number of new primary invairants that |
| | 2666 | // should be added in this degree - |
| | 2667 | ideal P,Pplus,sPplus,CI,B; // P: will contain primary invariants, |
| | 2668 | // Pplus: P+B, CI: a complete |
| | 2669 | // intersection with the same Hilbert |
| | 2670 | // function as P |
| | 2671 | ideal sP=std(P); |
| | 2672 | dB=1; // used as degree bound - |
| | 2673 | d=0; // initializing |
| | 2674 | i=0; |
| | 2675 | intvec deg_vector; |
| | 2676 | //-------------------- loop that searches for primary invariants ------------- |
| | 2677 | while(1) // repeat until n primary invariants are |
| | 2678 | { // found - |
| | 2679 | d=d+1; // degree where we'll search |
| | 2680 | if (v) |
| | 2681 | { " Computing primary invariants in degree "+string(d)+":"; |
| | 2682 | } |
| | 2683 | B=invariant_basis(d,#[1..gen_num]); // basis of invariants of degree d |
| | 2684 | if (B[1]<>0) |
| | 2685 | { Pplus=P+B; |
| | 2686 | sPplus=std(Pplus); |
| | 2687 | newdim=dim(sPplus); |
| | 2688 | dif=n-i-newdim; |
| | 2689 | } |
| | 2690 | else |
| | 2691 | { dif=0; |
| | 2692 | deg_vector=deg_vector,d; |
| | 2693 | } |
| | 2694 | if (dif<>0) // we have to find dif new primary |
| | 2695 | { // invariants |
| | 2696 | cd=size(B); |
| | 2697 | if (cd<>dif) |
| | 2698 | { P,sP,CI,dB=p_search(n,d,B,cd,P,sP,i,dif,dB,CI); |
| | 2699 | } |
| | 2700 | else // i.e. we can take all of B |
| | 2701 | { for(j=i+1;j<=i+dif;j=j+1) |
| | 2702 | { CI=CI+ideal(var(j)^d); |
| | 2703 | } |
| | 2704 | dB=dB+dif*(d-1); |
| | 2705 | P=Pplus; |
| | 2706 | sP=sPplus; |
| | 2707 | } |
| | 2708 | if (v) |
| | 2709 | { for (j=1;j<=size(P)-i;j=j+1) |
| | 2710 | { " We find: "+string(P[i+j]); |
| | 2711 | } |
| | 2712 | } |
| | 2713 | i=size(P); |
| | 2714 | if (i==n) // found all primary invariants |
| | 2715 | { if (v) |
| | 2716 | { ""; |
| | 2717 | " We found all primary invariants."; |
| | 2718 | ""; |
| | 2719 | } |
| | 2720 | return(matrix(P)); |
| | 2721 | } |
| | 2722 | } // done with degree d |
| | 2723 | else |
| | 2724 | { if (v) |
| | 2725 | { " None here..."; |
| | 2726 | } |
| | 2727 | } |
| | 2728 | } |
| | 2729 | } |
| | 2730 | example |
| | 2731 | { echo=2; |
| | 2732 | ring R=2,(x,y,z),dp; |
| | 2733 | matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; |
| | 2734 | matrix P=primary_charp_without(A); |
| | 2735 | print(P); |
| | 2736 | } |
| | 2737 | |
| | 2738 | proc primary_invariants (list #) |
| | 2739 | USAGE: primary_invariants(G1,G2,...[,flags]); |
| | 2740 | G1,G2,...: <matrices> generating a finite matrix group, flags: an |
| | 2741 | optional <intvec> with three entries, if the first one equals 0 (also |
| | 2742 | the default), the programme attempts to compute the Molien series and |
| | 2743 | Reynolds operator, if it equals 1, the programme is told that the |
| | 2744 | Molien series should not be computed, if it equals -1 characteristic 0 |
| | 2745 | is simulated, i.e. the Molien series is computed as if the base field |
| | 2746 | were characteristic 0 (the user must choose a field of large prime |
| | 2747 | characteristic, e.g. 32003) and if the first one is anything else, it |
| | 2748 | means that the characteristic of the base field divides the group |
| | 2749 | order, the second component should give the size of intervals between |
| | 2750 | canceling common factors in the expansion of the Molien series, 0 (the |
| | 2751 | default) means only once after generating all terms, in prime |
| | 2752 | characteristic also a negative number can be given to indicate that |
| | 2753 | common factors should always be canceled when the expansion is simple |
| | 2754 | (the root of the extension field does not occur among the coefficients) |
| | 2755 | DISPLAY: information about the various stages of the programme if the third |
| | 2756 | flag does not equal 0 |
| | 2757 | RETURN: primary invariants (type <matrix>) of the invariant ring and if |
| | 2758 | computable Reynolds operator (type <matrix>) and Molien series (type |
| | 2759 | <matrix>), if the first flag is 1 and we are in the non-modular case |
| | 2760 | then an <intvec> is returned giving some of the degrees where no |
| | 2761 | non-trivial homogeneous invariants can be found |
| | 2762 | EXAMPLE: example primary_invariants; shows an example |
| | 2763 | THEORY: Bases of homogeneous invariants are generated successively and those |
| | 2764 | are chosen as primary invariants that lower the dimension of the ideal |
| | 2765 | generated by the previously found invariants (see paper "Generating a |
| | 2766 | Noetherian Normalization of the Invariant Ring of a Finite Group" by |
| | 2767 | Decker, Heydtmann, Schreyer (1997) to appear in JSC). |
| | 2768 | { |
| | 2769 | // ----------------- checking input and setting flags ------------------------ |
| | 2770 | if (size(#)==0) |
| | 2771 | { "ERROR: There are no parameters."; |
| | 2772 | return(); |
| | 2773 | } |
| | 2774 | int ch=char(basering); // the algorithms depend very much on the |
| | 2775 | // characteristic of the ground field |
| | 2776 | int n=nvars(basering); // n is the number of variables, as well |
| | 2777 | // as the size of the matrices, as well |
| | 2778 | // as the number of primary invariants, |
| | 2779 | // we should get |
| | 2780 | int gen_num; |
| | 2781 | int mol_flag,v; |
| | 2782 | if (typeof(#[size(#)])=="intvec") |
| | 2783 | { if (size(#[size(#)])<>3) |
| | 2784 | { "ERROR: <intvec> should have three entries."; |
| | 2785 | return(); |
| | 2786 | } |
| | 2787 | gen_num=size(#)-1; |
| | 2788 | mol_flag=#[size(#)][1]; |
| | 2789 | if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag==-1))) |
| | 2790 | { "ERROR: the second component of <intvec> should be >=0"; |
| | 2791 | return(); |
| | 2792 | } |
| | 2793 | int interval=#[size(#)][2]; |
| | 2794 | v=#[size(#)][3]; |
| | 2795 | if (gen_num==0) |
| | 2796 | { "ERROR: There are no generators of a finite matrix group given."; |
| | 2797 | return(); |
| | 2798 | } |
| | 2799 | } |
| | 2800 | else |
| | 2801 | { gen_num=size(#); |
| | 2802 | mol_flag=0; |
| | 2803 | int interval=0; |
| | 2804 | v=0; |
| | 2805 | } |
| | 2806 | for (int i=1;i<=gen_num;i=i+1) |
| | 2807 | { if (typeof(#[i])=="matrix") |
| | 2808 | { if (nrows(#[i])<>n or ncols(#[i])<>n) |
| | 2809 | { "ERROR: The number of variables of the base ring needs to be the same"; |
| | 2810 | " as the dimension of the square matrices"; |
| | 2811 | return(); |
| | 2812 | } |
| | 2813 | } |
| | 2814 | else |
| | 2815 | { "ERROR: The first parameters should be a list of matrices"; |
| | 2816 | return(); |
| | 2817 | } |
| | 2818 | } |
| | 2819 | //---------------------------------------------------------------------------- |
| | 2820 | if (mol_flag==0) |
| | 2821 | { if (ch==0) |
| | 2822 | { matrix REY,M=reynolds_molien(#[1..gen_num],intvec(mol_flag,interval,v)); |
| | 2823 | // one will contain Reynolds operator and |
| | 2824 | // the other enumerator and denominator |
| | 2825 | // of Molien series |
| | 2826 | matrix P=primary_char0(REY,M,v); |
| | 2827 | return(P,REY,M); |
| | 2828 | } |
| | 2829 | else |
| | 2830 | { list L=group_reynolds(#[1..gen_num],v); |
| | 2831 | if (L[1]<>0) // testing whether we are in the modular |
| | 2832 | { string newring="aksldfalkdsflkj"; // case |
| | 2833 | if (minpoly==0) |
| | 2834 | { if (v) |
| | 2835 | { " We are dealing with the non-modular case."; |
| | 2836 | } |
| | 2837 | if (typeof(L[2])=="int") |
| | 2838 | { molien(L[3..size(L)],newring,L[2],intvec(mol_flag,interval,v)); |
| | 2839 | } |
| | 2840 | else |
| | 2841 | { molien(L[2..size(L)],newring,intvec(mol_flag,interval,v)); |
| | 2842 | } |
| | 2843 | matrix P=primary_charp(L[1],newring,v); |
| | 2844 | return(P,L[1],newring); |
| | 2845 | } |
| | 2846 | else |
| | 2847 | { if (v) |
| | 2848 | { " Since it is impossible for this programme to calculate the Molien series for"; |
| | 2849 | " invariant rings over extension fields of prime characteristic, we have to"; |
| | 2850 | " continue without it."; |
| | 2851 | ""; |
| | 2852 | |
| | 2853 | } |
| | 2854 | list l=primary_charp_no_molien(L[1],v); |
| | 2855 | if (size(l)==2) |
| | 2856 | { return(l[1],L[1],l[2]); |
| | 2857 | } |
| | 2858 | else |
| | 2859 | { return(l[1],L[1]); |
| | 2860 | } |
| | 2861 | } |
| | 2862 | } |
| | 2863 | else // the modular case |
| | 2864 | { if (v) |
| | 2865 | { " There is also no Molien series, we can make use of..."; |
| | 2866 | ""; |
| | 2867 | " We can start looking for primary invariants..."; |
| | 2868 | ""; |
| | 2869 | } |
| | 2870 | return(primary_charp_without(#[1..gen_num],v)); |
| | 2871 | } |
| | 2872 | } |
| | 2873 | } |
| | 2874 | if (mol_flag==1) // the user wants no calculation of the |
| | 2875 | { list L=group_reynolds(#[1..gen_num],v); // Molien series |
| | 2876 | if (ch==0) |
| | 2877 | { list l=primary_char0_no_molien(L[1],v); |
| | 2878 | if (size(l)==2) |
| | 2879 | { return(l[1],L[1],l[2]); |
| | 2880 | } |
| | 2881 | else |
| | 2882 | { return(l[1],L[1]); |
| | 2883 | } |
| | 2884 | } |
| | 2885 | else |
| | 2886 | { if (L[1]<>0) // testing whether we are in the modular |
| | 2887 | { list l=primary_charp_no_molien(L[1],v); // case |
| | 2888 | if (size(l)==2) |
| | 2889 | { return(l[1],L[1],l[2]); |
| | 2890 | } |
| | 2891 | else |
| | 2892 | { return(l[1],L[1]); |
| | 2893 | } |
| | 2894 | } |
| | 2895 | else // the modular case |
| | 2896 | { if (v) |
| | 2897 | { " We can start looking for primary invariants..."; |
| | 2898 | ""; |
| | 2899 | } |
| | 2900 | return(primary_charp_without(#[1..gen_num],v)); |
| | 2901 | } |
| | 2902 | } |
| | 2903 | } |
| | 2904 | if (mol_flag==-1) |
| | 2905 | { if (ch==0) |
| | 2906 | { "ERROR: Characteristic 0 can only be simulated in characteristic p>>0."; |
| | 2907 | return(); |
| | 2908 | } |
| | 2909 | list L=group_reynolds(#[1..gen_num],v); |
| | 2910 | string newring="aksldfalkdsflkj"; |
| | 2911 | molien(L[2..size(L)],newring,intvec(1,interval,v)); |
| | 2912 | matrix P=primary_charp(L[1],newring,v); |
| | 2913 | return(P,L[1],newring); |
| | 2914 | } |
| | 2915 | else // the user specified that the |
| | 2916 | { if (ch==0) // characteristic divides the group order |
| | 2917 | { "ERROR: The characteristic cannot divide the group order when it is 0."; |
| | 2918 | return(); |
| | 2919 | } |
| | 2920 | if (v) |
| | 2921 | { ""; |
| | 2922 | } |
| | 2923 | return(primary_charp_without(#[1..gen_num],v)); |
| | 2924 | } |
| | 2925 | } |
| | 2926 | example |
| | 2927 | { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; |
| | 2928 | echo=2; |
| | 2929 | ring R=0,(x,y,z),dp; |
| | 2930 | matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; |
| | 2931 | list L=primary_invariants(A); |
| | 2932 | print(L[1]); |
| 1134 | | proc sec_minus_mol (ideal P, ideal sP, int g, int v, list#) |
| | 2940 | proc search_random (int n,int d,ideal B,int cd,ideal P,int i,int dif,int dB,ideal CI,int max) |
| | 2941 | { string answer; |
| | 2942 | degBound=0; |
| | 2943 | int j,k,test_dim,flag; |
| | 2944 | matrix test_matrix[1][dif]; // the linear combination to test |
| | 2945 | intvec h; // Hilbert series |
| | 2946 | for (j=i+1;j<=i+dif;j=j+1) |
| | 2947 | { CI=CI+ideal(var(j)^d); // homogeneous polynomial of the same |
| | 2948 | // degree as the one we're looking for |
| | 2949 | // is added |
| | 2950 | } |
| | 2951 | ideal TEST; |
| | 2952 | // h=hilb(std(CI),1); |
| | 2953 | dB=dB+dif*(d-1); // used as degBound |
| | 2954 | while (1) |
| | 2955 | { test_matrix=matrix(B)*random(max,cd,dif); |
| | 2956 | // degBound=dB; |
| | 2957 | TEST=P+ideal(test_matrix); |
| | 2958 | attrib(TEST,"isSB",1); |
| | 2959 | test_dim=dim(TEST); |
| | 2960 | // degBound=0; |
| | 2961 | if (n-test_dim==i+dif) |
| | 2962 | { break; |
| | 2963 | } |
| | 2964 | // degBound=dB; |
| | 2965 | test_dim=dim(std(TEST)); |
| | 2966 | // test_dim=dim(std(TEST,h)); // Hilbert driven std-calculation |
| | 2967 | // degBound=0; |
| | 2968 | if (n-test_dim==i+dif) |
| | 2969 | { break; |
| | 2970 | } |
| | 2971 | else |
| | 2972 | { "HELP: The "+string(dif)+" random combination(s) of the "+string(cd)+" basis elements with"; |
| | 2973 | " coefficients in the range from -"+string(max)+" to "+string(max)+" did not lower the"; |
| | 2974 | " dimension by "+string(dif)+". You can abort, try again or give a new range:"; |
| | 2975 | answer=""; |
| | 2976 | while (answer<>"n |
| | 2977 | " && answer<>"y |
| | 2978 | ") |
| | 2979 | { " Do you want to abort (y/n)?"; |
| | 2980 | answer=read(""); |
| | 2981 | } |
| | 2982 | if (answer=="y |
| | 2983 | ") |
| | 2984 | { flag=1; |
| | 2985 | break; |
| | 2986 | } |
| | 2987 | answer=""; |
| | 2988 | while (answer<>"n |
| | 2989 | " && answer<>"y |
| | 2990 | ") |
| | 2991 | { " Do you want to try again (y/n)?"; |
| | 2992 | answer=read(""); |
| | 2993 | } |
| | 2994 | if (answer=="n |
| | 2995 | ") |
| | 2996 | { flag=1; |
| | 2997 | while (flag) |
| | 2998 | { " Give a new <int> > "+string(max)+" that bounds the range of coefficients:"; |
| | 2999 | answer=read(""); |
| | 3000 | for (j=1;j<=size(answer)-1;j=j+1) |
| | 3001 | { for (k=0;k<=9;k=k+1) |
| | 3002 | { if (answer[j]==string(k)) |
| | 3003 | { break; |
| | 3004 | } |
| | 3005 | } |
| | 3006 | if (k>9) |
| | 3007 | { flag=1; |
| | 3008 | break; |
| | 3009 | } |
| | 3010 | flag=0; |
| | 3011 | } |
| | 3012 | if (not(flag)) |
| | 3013 | { execute "test_dim="+string(answer[1..size(answer)]); |
| | 3014 | if (test_dim<=max) |
| | 3015 | { flag=1; |
| | 3016 | } |
| | 3017 | else |
| | 3018 | { max=test_dim; |
| | 3019 | } |
| | 3020 | } |
| | 3021 | } |
| | 3022 | } |
| | 3023 | } |
| | 3024 | } |
| | 3025 | if (not(flag)) |
| | 3026 | { P[(i+1)..(i+dif)]=test_matrix[1,1..dif]; |
| | 3027 | } |
| | 3028 | return(P,CI,dB); |
| | 3029 | } |
| | 3030 | |
| | 3031 | //////////////////////////////////////////////////////////////////////////////// |
| | 3032 | // This procedure finds at most dif primary invariants in degree d. It returns |
| | 3033 | // all primary invariants found so far. The coefficients lie in the field of |
| | 3034 | // characteristic p>0. |
| | 3035 | //////////////////////////////////////////////////////////////////////////////// |
| | 3036 | proc p_search_random (int n,int d,ideal B,int cd,ideal P,int i,int dif,int dB,ideal CI,int max) |
| | 3037 | { string answer; |
| | 3038 | degBound=0; |
| | 3039 | int j,k,test_dim,flag; |
| | 3040 | matrix test_matrix[1][dif]; // the linear combination to test |
| | 3041 | intvec h; // Hilbert series |
| | 3042 | ideal TEST; |
| | 3043 | while (dif>0) |
| | 3044 | { for (j=i+1;j<=i+dif;j=j+1) |
| | 3045 | { CI=CI+ideal(var(j)^d); // homogeneous polynomial of the same |
| | 3046 | // degree as the one we're looking for |
| | 3047 | // is added |
| | 3048 | } |
| | 3049 | // h=hilb(std(CI),1); |
| | 3050 | dB=dB+dif*(d-1); // used as degBound |
| | 3051 | test_matrix=matrix(B)*random(max,cd,dif); |
| | 3052 | // degBound=dB; |
| | 3053 | TEST=P+ideal(test_matrix); |
| | 3054 | attrib(TEST,"isSB",1); |
| | 3055 | test_dim=dim(TEST); |
| | 3056 | // degBound=0; |
| | 3057 | if (n-test_dim==i+dif) |
| | 3058 | { break; |
| | 3059 | } |
| | 3060 | // degBound=dB; |
| | 3061 | test_dim=dim(std(TEST)); |
| | 3062 | // test_dim=dim(std(TEST,h)); // Hilbert driven std-calculation |
| | 3063 | // degBound=0; |
| | 3064 | if (n-test_dim==i+dif) |
| | 3065 | { break; |
| | 3066 | } |
| | 3067 | else |
| | 3068 | { "HELP: The "+string(dif)+" random combination(s) of the "+string(cd)+" basis elements with"; |
| | 3069 | " coefficients in the range from -"+string(max)+" to "+string(max)+" did not lower the"; |
| | 3070 | " dimension by "+string(dif)+". You can abort, try again, lower the number of"; |
| | 3071 | " combinations searched for by 1 or give a larger coefficient range:"; |
| | 3072 | answer=""; |
| | 3073 | while (answer<>"n |
| | 3074 | " && answer<>"y |
| | 3075 | ") |
| | 3076 | { " Do you want to abort (y/n)?"; |
| | 3077 | answer=read(""); |
| | 3078 | } |
| | 3079 | if (answer=="y |
| | 3080 | ") |
| | 3081 | { flag=1; |
| | 3082 | break; |
| | 3083 | } |
| | 3084 | answer=""; |
| | 3085 | while (answer<>"n |
| | 3086 | " && answer<>"y |
| | 3087 | ") |
| | 3088 | { " Do you want to try again (y/n)?"; |
| | 3089 | answer=read(""); |
| | 3090 | } |
| | 3091 | if (answer=="n |
| | 3092 | ") |
| | 3093 | { answer=""; |
| | 3094 | while (answer<>"n |
| | 3095 | " && answer<>"y |
| | 3096 | ") |
| | 3097 | { " Do you want to lower the number of combinations by 1 (y/n)?"; |
| | 3098 | answer=read(""); |
| | 3099 | } |
| | 3100 | if (answer=="y |
| | 3101 | ") |
| | 3102 | { dif=dif-1; |
| | 3103 | } |
| | 3104 | else |
| | 3105 | { flag=1; |
| | 3106 | while (flag) |
| | 3107 | { " Give a new <int> > "+string(max)+" that bounds the range of coefficients:"; |
| | 3108 | answer=read(""); |
| | 3109 | for (j=1;j<=size(answer)-1;j=j+1) |
| | 3110 | { for (k=0;k<=9;k=k+1) |
| | 3111 | { if (answer[j]==string(k)) |
| | 3112 | { break; |
| | 3113 | } |
| | 3114 | } |
| | 3115 | if (k>9) |
| | 3116 | { flag=1; |
| | 3117 | break; |
| | 3118 | } |
| | 3119 | flag=0; |
| | 3120 | } |
| | 3121 | if (not(flag)) |
| | 3122 | { execute "test_dim="+string(answer[1..size(answer)]); |
| | 3123 | if (test_dim<=max) |
| | 3124 | { flag=1; |
| | 3125 | } |
| | 3126 | else |
| | 3127 | { max=test_dim; |
| | 3128 | } |
| | 3129 | } |
| | 3130 | } |
| | 3131 | } |
| | 3132 | } |
| | 3133 | } |
| | 3134 | CI=CI[1..i]; |
| | 3135 | dB=dB-dif*(d-1); |
| | 3136 | } |
| | 3137 | if (dif && not(flag)) |
| | 3138 | { P[(i+1)..(i+dif)]=test_matrix[1,1..dif]; |
| | 3139 | } |
| | 3140 | if (dif && flag) |
| | 3141 | { P[n+1]=0; |
| | 3142 | } |
| | 3143 | return(P,CI,dB); |
| | 3144 | } |
| | 3145 | |
| | 3146 | proc primary_char0_random (matrix REY,matrix M,int max,list #) |
| | 3147 | USAGE: primary_char0_random(REY,M,r[,v]); |
| | 3148 | REY: a <matrix> representing the Reynolds operator, M: a 1x2 <matrix> |
| | 3149 | representing the Molien series, r: an <int> where -|r| to |r| is the |
| | 3150 | range of coefficients of the random combinations of bases elements, |
| | 3151 | v: an optional <int> |
| | 3152 | ASSUME: REY is the first return value of group_reynolds or reynolds_molien and |
| | 3153 | M the one of molien or the second one of reynolds_molien |
| | 3154 | DISPLAY: information about the various stages of the programme if v does not |
| | 3155 | equal 0 |
| | 3156 | RETURN: primary invariants (type <matrix>) of the invariant ring |
| | 3157 | EXAMPLE: example primary_char0_random; shows an example |
| | 3158 | THEORY: Bases of homogeneous invariants are generated successively and random |
| | 3159 | linear combinations are chosen as primary invariants that lower the |
| | 3160 | dimension of the ideal generated by the previously found invariants |
| | 3161 | (see paper "Generating a Noetherian Normalization of the Invariant Ring |
| | 3162 | of a Finite Group" by Decker, Heydtmann, Schreyer (1997) to appear in |
| | 3163 | JSC). |
| | 3164 | { degBound=0; |
| | 3165 | if (char(basering)<>0) |
| | 3166 | { "ERROR: primary_char0_random should only be used with rings of"; |
| | 3167 | " characteristic 0."; |
| | 3168 | return(); |
| | 3169 | } |
| | 3170 | //----------------- checking input and setting verbose mode ------------------ |
| | 3171 | if (size(#)>1) |
| | 3172 | { "ERROR: primary_char0_random can only have four parameters."; |
| | 3173 | return(); |
| | 3174 | } |
| | 3175 | if (size(#)==1) |
| | 3176 | { if (typeof(#[1])<>"int") |
| | 3177 | { "ERROR: The fourth parameter should be of type <int>."; |
| | 3178 | return(); |
| | 3179 | } |
| | 3180 | else |
| | 3181 | { int v=#[1]; |
| | 3182 | } |
| | 3183 | } |
| | 3184 | else |
| | 3185 | { int v=0; |
| | 3186 | } |
| | 3187 | int n=nvars(basering); // n is the number of variables, as well |
| | 3188 | // as the size of the matrices, as well |
| | 3189 | // as the number of primary invariants, |
| | 3190 | // we should get |
| | 3191 | if (ncols(REY)<>n) |
| | 3192 | { "ERROR: First parameter ought to be the Reynolds operator." |
| | 3193 | return(); |
| | 3194 | } |
| | 3195 | if (ncols(M)<>2 or nrows(M)<>1) |
| | 3196 | { "ERROR: Second parameter ought to be the Molien series." |
| | 3197 | return(); |
| | 3198 | } |
| | 3199 | //---------------------------------------------------------------------------- |
| | 3200 | if (v && voice<>2) |
| | 3201 | { " We can start looking for primary invariants..."; |
| | 3202 | ""; |
| | 3203 | } |
| | 3204 | if (v && voice==2) |
| | 3205 | { ""; |
| | 3206 | } |
| | 3207 | //------------------------- initializing variables --------------------------- |
| | 3208 | int dB; |
| | 3209 | poly p(1..2); // p(1) will be used for single terms of |
| | 3210 | // the partial expansion, p(2) to store |
| | 3211 | p(1..2)=partial_molien(M,1); // the intermediate result - |
| | 3212 | poly v1=var(1); // we need v1 to split off coefficients |
| | 3213 | // in the partial expansion of M (which |
| | 3214 | // is in terms of the first variable) - |
| | 3215 | int j,d,cd,newdim,dif; // d: current degree, cd: dimension of |
| | 3216 | // space of invariants of degree d, |
| | 3217 | // newdim: dimension the ideal generated |
| | 3218 | // the primary invariants plus basis |
| | 3219 | // elements, dif=n-i-newdim, i.e. the |
| | 3220 | // number of new primary invairants that |
| | 3221 | // should be added in this degree - |
| | 3222 | ideal P,Pplus,CI,B; // P: will contain primary invariants, |
| | 3223 | // Pplus: P+B,CI: a complete |
| | 3224 | // intersection with the same Hilbert |
| | 3225 | // function as P - |
| | 3226 | dB=1; // used as degree bound |
| | 3227 | int i=0; |
| | 3228 | //-------------- loop that searches for primary invariants ------------------ |
| | 3229 | while(1) // repeat until n primary invariants are |
| | 3230 | { // found - |
| | 3231 | p(1..2)=partial_molien(M,1,p(2)); // next term of the partial expansion - |
| | 3232 | d=deg(p(1)); // degree where we'll search - |
| | 3233 | cd=int(coef(p(1),v1)[2,1]); // dimension of the homogeneous space of |
| | 3234 | // inviarants of degree d |
| | 3235 | if (v) |
| | 3236 | { " Computing primary invariants in degree "+string(d)+":"; |
| | 3237 | } |
| | 3238 | B=invariant_basis_reynolds(REY,d,intvec(cd,6)); // basis of invariants of |
| | 3239 | // degree d |
| | 3240 | if (B[1]<>0) |
| | 3241 | { Pplus=P+B; |
| | 3242 | newdim=dim(std(Pplus)); |
| | 3243 | dif=n-i-newdim; |
| | 3244 | } |
| | 3245 | else |
| | 3246 | { dif=0; |
| | 3247 | } |
| | 3248 | if (dif<>0) // we have to find dif new primary |
| | 3249 | { // invariants |
| | 3250 | if (cd<>dif) |
| | 3251 | { P,CI,dB=search_random(n,d,B,cd,P,i,dif,dB,CI,max); // searching for |
| | 3252 | } // dif invariants - |
| | 3253 | else // i.e. we can take all of B |
| | 3254 | { for(j=i+1;j>i+dif;j=j+1) |
| | 3255 | { CI=CI+ideal(var(j)^d); |
| | 3256 | } |
| | 3257 | dB=dB+dif*(d-1); |
| | 3258 | P=Pplus; |
| | 3259 | } |
| | 3260 | if (ncols(P)==i) |
| | 3261 | { "WARNING: The return value is not a set of primary invariants, but"; |
| | 3262 | " polynomials qualifying as the first "+string(i)+" primary invariants."; |
| | 3263 | return(matrix(P)); |
| | 3264 | } |
| | 3265 | if (v) |
| | 3266 | { for (j=1;j<=dif;j=j+1) |
| | 3267 | { " We find: "+string(P[i+j]); |
| | 3268 | } |
| | 3269 | } |
| | 3270 | i=i+dif; |
| | 3271 | if (i==n) // found all primary invariants |
| | 3272 | { if (v) |
| | 3273 | { ""; |
| | 3274 | " We found all primary invariants."; |
| | 3275 | ""; |
| | 3276 | } |
| | 3277 | return(matrix(P)); |
| | 3278 | } |
| | 3279 | } // done with degree d |
| | 3280 | } |
| | 3281 | } |
| | 3282 | example |
| | 3283 | { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; |
| | 3284 | echo=2; |
| | 3285 | ring R=0,(x,y,z),dp; |
| | 3286 | matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; |
| | 3287 | matrix REY,M=reynolds_molien(A); |
| | 3288 | matrix P=primary_char0_random(REY,M,1); |
| | 3289 | print(P); |
| | 3290 | } |
| | 3291 | |
| | 3292 | proc primary_charp_random (matrix REY,string ring_name,int max,list #) |
| | 3293 | USAGE: primary_charp_random(REY,ringname,r[,v]); |
| | 3294 | REY: a <matrix> representing the Reynolds operator, ringname: a |
| | 3295 | <string> giving the name of a ring where the Molien series is stored, |
| | 3296 | r: an <int> where -|r| to |r| is the range of coefficients of the |
| | 3297 | random combinations of bases elements, v: an optional <int> |
| | 3298 | ASSUME: REY is the first return value of group_reynolds or reynolds_molien and |
| | 3299 | ringname gives the name of a ring of characteristic 0 that has been |
| | 3300 | created by molien or reynolds_molien |
| | 3301 | DISPLAY: information about the various stages of the programme if v does not |
| | 3302 | equal 0 |
| | 3303 | RETURN: primary invariants (type <matrix>) of the invariant ring |
| | 3304 | EXAMPLE: example primary_charp_random; shows an example |
| | 3305 | THEORY: Bases of homogeneous invariants are generated successively and random |
| | 3306 | linear combinations are chosen as primary invariants that lower the |
| | 3307 | dimension of the ideal generated by the previously found invariants |
| | 3308 | (see paper "Generating a Noetherian Normalization of the Invariant Ring |
| | 3309 | of a Finite Group" by Decker, Heydtmann, Schreyer (1997) to appear in |
| | 3310 | JSC). |
| | 3311 | { degBound=0; |
| | 3312 | // ---------------- checking input and setting verbose mode ------------------ |
| | 3313 | if (char(basering)==0) |
| | 3314 | { "ERROR: primary_charp_random should only be used with rings of"; |
| | 3315 | " characteristic p>0."; |
| | 3316 | return(); |
| | 3317 | } |
| | 3318 | if (size(#)>1) |
| | 3319 | { "ERROR: primary_charp_random can only have four parameters."; |
| | 3320 | return(); |
| | 3321 | } |
| | 3322 | if (size(#)==1) |
| | 3323 | { if (typeof(#[1])<>"int") |
| | 3324 | { "ERROR: The fourth parameter should be of type <int>."; |
| | 3325 | return(); |
| | 3326 | } |
| | 3327 | else |
| | 3328 | { int v=#[1]; |
| | 3329 | } |
| | 3330 | } |
| | 3331 | else |
| | 3332 | { int v=0; |
| | 3333 | } |
| | 3334 | def br=basering; |
| | 3335 | int n=nvars(br); // n is the number of variables, as well |
| | 3336 | // as the size of the matrices, as well |
| | 3337 | // as the number of primary invariants, |
| | 3338 | // we should get |
| | 3339 | if (ncols(REY)<>n) |
| | 3340 | { "ERROR: First parameter ought to be the Reynolds operator." |
| | 3341 | return(); |
| | 3342 | } |
| | 3343 | if (typeof(`ring_name`)<>"ring") |
| | 3344 | { "ERROR: Second parameter ought to the name of a ring where the Molien"; |
| | 3345 | " is stored."; |
| | 3346 | return(); |
| | 3347 | } |
| | 3348 | //---------------------------------------------------------------------------- |
| | 3349 | if (v && voice<>2) |
| | 3350 | { " We can start looking for primary invariants..."; |
| | 3351 | ""; |
| | 3352 | } |
| | 3353 | if (v && voice==2) |
| | 3354 | { ""; |
| | 3355 | } |
| | 3356 | //----------------------- initializing variables ----------------------------- |
| | 3357 | int dB; |
| | 3358 | setring `ring_name`; // the Molien series is stores here - |
| | 3359 | poly p(1..2); // p(1) will be used for single terms of |
| | 3360 | // the partial expansion, p(2) to store |
| | 3361 | p(1..2)=partial_molien(M,1); // the intermediate result - |
| | 3362 | poly v1=var(1); // we need v1 to split off coefficients |
| | 3363 | // in the partial expansion of M (which |
| | 3364 | // is in terms of the first variable) |
| | 3365 | setring br; |
| | 3366 | int j,d,cd,newdim,dif; // d: current degree, cd: dimension of |
| | 3367 | // space of invariants of degree d, |
| | 3368 | // newdim: dimension the ideal generated |
| | 3369 | // the primary invariants plus basis |
| | 3370 | // elements, dif=n-i-newdim, i.e. the |
| | 3371 | // number of new primary invairants that |
| | 3372 | // should be added in this degree - |
| | 3373 | ideal P,Pplus,CI,B; // P: will contain primary invariants, |
| | 3374 | // Pplus: P+B, CI: a complete |
| | 3375 | // intersection with the same Hilbert |
| | 3376 | // function as P - |
| | 3377 | dB=1; // used as degree bound |
| | 3378 | int i=0; |
| | 3379 | //---------------- loop that searches for primary invariants ----------------- |
| | 3380 | while(1) // repeat until n primary invariants are |
| | 3381 | { // found |
| | 3382 | setring `ring_name`; |
| | 3383 | p(1..2)=partial_molien(M,1,p(2)); // next term of the partial expansion - |
| | 3384 | d=deg(p(1)); // degree where we'll search - |
| | 3385 | cd=int(coef(p(1),v1)[2,1]); // dimension of the homogeneous space of |
| | 3386 | // inviarants of degree d |
| | 3387 | setring br; |
| | 3388 | if (v) |
| | 3389 | { " Computing primary invariants in degree "+string(d)+":"; |
| | 3390 | } |
| | 3391 | B=invariant_basis_reynolds(REY,d,intvec(cd,6)); // basis of invariants of |
| | 3392 | // degree d |
| | 3393 | if (B[1]<>0) |
| | 3394 | { Pplus=P+B; |
| | 3395 | newdim=dim(std(Pplus)); |
| | 3396 | dif=n-i-newdim; |
| | 3397 | } |
| | 3398 | else |
| | 3399 | { dif=0; |
| | 3400 | } |
| | 3401 | if (dif<>0) // we have to find dif new primary |
| | 3402 | { // invariants |
| | 3403 | if (cd<>dif) |
| | 3404 | { P,CI,dB=p_search_random(n,d,B,cd,P,i,dif,dB,CI,max); |
| | 3405 | } |
| | 3406 | else // i.e. we can take all of B |
| | 3407 | { for(j=i+1;j>i+dif;j=j+1) |
| | 3408 | { CI=CI+ideal(var(j)^d); |
| | 3409 | } |
| | 3410 | dB=dB+dif*(d-1); |
| | 3411 | P=Pplus; |
| | 3412 | } |
| | 3413 | if (ncols(P)==n+1) |
| | 3414 | { "WARNING: The first return value is not a set of primary invariants,"; |
| | 3415 | " but polynomials qualifying as the first "+string(i)+" primary invariants."; |
| | 3416 | return(matrix(P)); |
| | 3417 | } |
| | 3418 | if (v) |
| | 3419 | { for (j=1;j<=size(P)-i;j=j+1) |
| | 3420 | { " We find: "+string(P[i+j]); |
| | 3421 | } |
| | 3422 | } |
| | 3423 | i=size(P); |
| | 3424 | if (i==n) // found all primary invariants |
| | 3425 | { if (v) |
| | 3426 | { ""; |
| | 3427 | " We found all primary invariants."; |
| | 3428 | ""; |
| | 3429 | } |
| | 3430 | return(matrix(P)); |
| | 3431 | } |
| | 3432 | } // done with degree d |
| | 3433 | } |
| | 3434 | } |
| | 3435 | example |
| | 3436 | { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed into"; |
| | 3437 | " characteristic 3)"; |
| | 3438 | echo=2; |
| | 3439 | ring R=3,(x,y,z),dp; |
| | 3440 | matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; |
| | 3441 | list L=group_reynolds(A); |
| | 3442 | string newring="alskdfj"; |
| | 3443 | molien(L[2..size(L)],newring); |
| | 3444 | matrix P=primary_charp_random(L[1],newring,1); |
| | 3445 | kill `newring`; |
| | 3446 | print(P); |
| | 3447 | } |
| | 3448 | |
| | 3449 | proc primary_char0_no_molien_random (matrix REY, int max, list #) |
| | 3450 | USAGE: primary_char0_no_molien_random(REY,r[,v]); |
| | 3451 | REY: a <matrix> representing the Reynolds operator, r: an <int> where |
| | 3452 | -|r| to |r| is the range of coefficients of the random combinations of |
| | 3453 | bases elements, v: an optional <int> |
| | 3454 | ASSUME: REY is the first return value of group_reynolds or reynolds_molien |
| | 3455 | DISPLAY: information about the various stages of the programme if v does not |
| | 3456 | equal 0 |
| | 3457 | RETURN: primary invariants (type <matrix>) of the invariant ring and an |
| | 3458 | <intvec> listing some of the degrees where no non-trivial homogeneous |
| | 3459 | invariants are to be found |
| | 3460 | EXAMPLE: example primary_char0_no_molien_random; shows an example |
| | 3461 | THEORY: Bases of homogeneous invariants are generated successively and random |
| | 3462 | linear combinations are chosen as primary invariants that lower the |
| | 3463 | dimension of the ideal generated by the previously found invariants |
| | 3464 | (see paper "Generating a Noetherian Normalization of the Invariant Ring |
| | 3465 | of a Finite Group" by Decker, Heydtmann, Schreyer (1997) to appear in |
| | 3466 | JSC). |
| | 3467 | { degBound=0; |
| | 3468 | //-------------- checking input and setting verbose mode --------------------- |
| | 3469 | if (char(basering)<>0) |
| | 3470 | { "ERROR: primary_char0_no_molien_random should only be used with rings of"; |
| | 3471 | " characteristic 0."; |
| | 3472 | return(); |
| | 3473 | } |
| | 3474 | if (size(#)>1) |
| | 3475 | { "ERROR: primary_char0_no_molien_random can only have three parameters."; |
| | 3476 | return(); |
| | 3477 | } |
| | 3478 | if (size(#)==1) |
| | 3479 | { if (typeof(#[1])<>"int") |
| | 3480 | { "ERROR: The third parameter should be of type <int>."; |
| | 3481 | return(); |
| | 3482 | } |
| | 3483 | else |
| | 3484 | { int v=#[1]; |
| | 3485 | } |
| | 3486 | } |
| | 3487 | else |
| | 3488 | { int v=0; |
| | 3489 | } |
| | 3490 | int n=nvars(basering); // n is the number of variables, as well |
| | 3491 | // as the size of the matrices, as well |
| | 3492 | // as the number of primary invariants, |
| | 3493 | // we should get |
| | 3494 | if (ncols(REY)<>n) |
| | 3495 | { "ERROR: First parameter ought to be the Reynolds operator." |
| | 3496 | return(); |
| | 3497 | } |
| | 3498 | //---------------------------------------------------------------------------- |
| | 3499 | if (v && voice<>2) |
| | 3500 | { " We can start looking for primary invariants..."; |
| | 3501 | ""; |
| | 3502 | } |
| | 3503 | if (v && voice==2) |
| | 3504 | { ""; |
| | 3505 | } |
| | 3506 | //----------------------- initializing variables ----------------------------- |
| | 3507 | int dB; |
| | 3508 | int j,d,cd,newdim,dif; // d: current degree, cd: dimension of |
| | 3509 | // space of invariants of degree d, |
| | 3510 | // newdim: dimension the ideal generated |
| | 3511 | // the primary invariants plus basis |
| | 3512 | // elements, dif=n-i-newdim, i.e. the |
| | 3513 | // number of new primary invairants that |
| | 3514 | // should be added in this degree - |
| | 3515 | ideal P,Pplus,CI,B; // P: will contain primary invariants, |
| | 3516 | // Pplus: P+B, CI: a complete |
| | 3517 | // intersection with the same Hilbert |
| | 3518 | // function as P - |
| | 3519 | dB=1; // used as degree bound - |
| | 3520 | d=0; // initializing |
| | 3521 | int i=0; |
| | 3522 | intvec deg_vector; |
| | 3523 | //------------------ loop that searches for primary invariants --------------- |
| | 3524 | while(1) // repeat until n primary invariants are |
| | 3525 | { // found - |
| | 3526 | d=d+1; // degree where we'll search |
| | 3527 | if (v) |
| | 3528 | { " Computing primary invariants in degree "+string(d)+":"; |
| | 3529 | } |
| | 3530 | B=invariant_basis_reynolds(REY,d,intvec(-1,6)); // basis of invariants of |
| | 3531 | // degree d |
| | 3532 | if (B[1]<>0) |
| | 3533 | { Pplus=P+B; |
| | 3534 | newdim=dim(std(Pplus)); |
| | 3535 | dif=n-i-newdim; |
| | 3536 | } |
| | 3537 | else |
| | 3538 | { dif=0; |
| | 3539 | deg_vector=deg_vector,d; |
| | 3540 | } |
| | 3541 | if (dif<>0) // we have to find dif new primary |
| | 3542 | { // invariants |
| | 3543 | cd=size(B); |
| | 3544 | if (cd<>dif) |
| | 3545 | { P,CI,dB=search_random(n,d,B,cd,P,i,dif,dB,CI,max); |
| | 3546 | } |
| | 3547 | else // i.e. we can take all of B |
| | 3548 | { for(j=i+1;j<=i+dif;j=j+1) |
| | 3549 | { CI=CI+ideal(var(j)^d); |
| | 3550 | } |
| | 3551 | dB=dB+dif*(d-1); |
| | 3552 | P=Pplus; |
| | 3553 | } |
| | 3554 | if (ncols(P)==i) |
| | 3555 | { "WARNING: The first return value is not a set of primary invariants,"; |
| | 3556 | " but polynomials qualifying as the first "+string(i)+" primary invariants."; |
| | 3557 | return(matrix(P)); |
| | 3558 | } |
| | 3559 | if (v) |
| | 3560 | { for (j=1;j<=dif;j=j+1) |
| | 3561 | { " We find: "+string(P[i+j]); |
| | 3562 | } |
| | 3563 | } |
| | 3564 | i=i+dif; |
| | 3565 | if (i==n) // found all primary invariants |
| | 3566 | { if (v) |
| | 3567 | { ""; |
| | 3568 | " We found all primary invariants."; |
| | 3569 | ""; |
| | 3570 | } |
| | 3571 | if (deg_vector==0) |
| | 3572 | { return(matrix(P)); |
| | 3573 | } |
| | 3574 | else |
| | 3575 | { return(matrix(P),compress(deg_vector)); |
| | 3576 | } |
| | 3577 | } |
| | 3578 | } // done with degree d |
| | 3579 | else |
| | 3580 | { if (v) |
| | 3581 | { " None here..."; |
| | 3582 | } |
| | 3583 | } |
| | 3584 | } |
| | 3585 | } |
| | 3586 | example |
| | 3587 | { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; |
| | 3588 | echo=2; |
| | 3589 | ring R=0,(x,y,z),dp; |
| | 3590 | matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; |
| | 3591 | list L=group_reynolds(A); |
| | 3592 | list l=primary_char0_no_molien_random(L[1],1); |
| | 3593 | print(l[1]); |
| | 3594 | } |
| | 3595 | |
| | 3596 | proc primary_charp_no_molien_random (matrix REY, int max, list #) |
| | 3597 | USAGE: primary_charp_no_molien_random(REY,r[,v]); |
| | 3598 | REY: a <matrix> representing the Reynolds operator, r: an <int> where |
| | 3599 | -|r| to |r| is the range of coefficients of the random combinations of |
| | 3600 | bases elements, v: an optional <int> |
| | 3601 | ASSUME: REY is the first return value of group_reynolds or reynolds_molien |
| | 3602 | DISPLAY: information about the various stages of the programme if v does not |
| | 3603 | equal 0 |
| | 3604 | RETURN: primary invariants (type <matrix>) of the invariant ring and an |
| | 3605 | <intvec> listing some of the degrees where no non-trivial homogeneous |
| | 3606 | invariants are to be found |
| | 3607 | EXAMPLE: example primary_charp_no_molien_random; shows an example |
| | 3608 | THEORY: Bases of homogeneous invariants are generated successively and random |
| | 3609 | linear combinations are chosen as primary invariants that lower the |
| | 3610 | dimension of the ideal generated by the previously found invariants |
| | 3611 | (see paper "Generating a Noetherian Normalization of the Invariant Ring |
| | 3612 | of a Finite Group" by Decker, Heydtmann, Schreyer (1997) to appear in |
| | 3613 | JSC). |
| | 3614 | { degBound=0; |
| | 3615 | //----------------- checking input and setting verbose mode ------------------ |
| | 3616 | if (char(basering)==0) |
| | 3617 | { "ERROR: primary_charp_no_molien_random should only be used with rings of"; |
| | 3618 | " characteristic p>0."; |
| | 3619 | return(); |
| | 3620 | } |
| | 3621 | if (size(#)>1) |
| | 3622 | { "ERROR: primary_charp_no_molien_random can only have three parameters."; |
| | 3623 | return(); |
| | 3624 | } |
| | 3625 | if (size(#)==1) |
| | 3626 | { if (typeof(#[1])<>"int") |
| | 3627 | { "ERROR: The third parameter should be of type <int>."; |
| | 3628 | return(); |
| | 3629 | } |
| | 3630 | else |
| | 3631 | { int v=#[1]; |
| | 3632 | } |
| | 3633 | } |
| | 3634 | else |
| | 3635 | { int v=0; |
| | 3636 | } |
| | 3637 | int n=nvars(basering); // n is the number of variables, as well |
| | 3638 | // as the size of the matrices, as well |
| | 3639 | // as the number of primary invariants, |
| | 3640 | // we should get |
| | 3641 | if (ncols(REY)<>n) |
| | 3642 | { "ERROR: First parameter ought to be the Reynolds operator." |
| | 3643 | return(); |
| | 3644 | } |
| | 3645 | //---------------------------------------------------------------------------- |
| | 3646 | if (v && voice<>2) |
| | 3647 | { " We can start looking for primary invariants..."; |
| | 3648 | ""; |
| | 3649 | } |
| | 3650 | if (v && voice==2) |
| | 3651 | { ""; |
| | 3652 | } |
| | 3653 | //-------------------- initializing variables -------------------------------- |
| | 3654 | int dB; |
| | 3655 | int j,d,cd,newdim,dif; // d: current degree, cd: dimension of |
| | 3656 | // space of invariants of degree d, |
| | 3657 | // newdim: dimension the ideal generated |
| | 3658 | // the primary invariants plus basis |
| | 3659 | // elements, dif=n-i-newdim, i.e. the |
| | 3660 | // number of new primary invairants that |
| | 3661 | // should be added in this degree - |
| | 3662 | ideal P,Pplus,CI,B; // P: will contain primary invariants, |
| | 3663 | // Pplus: P+B, CI: a complete |
| | 3664 | // intersection with the same Hilbert |
| | 3665 | // function as P - |
| | 3666 | dB=1; // used as degree bound - |
| | 3667 | d=0; // initializing |
| | 3668 | int i=0; |
| | 3669 | intvec deg_vector; |
| | 3670 | //------------------ loop that searches for primary invariants --------------- |
| | 3671 | while(1) // repeat until n primary invariants are |
| | 3672 | { // found - |
| | 3673 | d=d+1; // degree where we'll search |
| | 3674 | if (v) |
| | 3675 | { " Computing primary invariants in degree "+string(d)+":"; |
| | 3676 | } |
| | 3677 | B=invariant_basis_reynolds(REY,d,intvec(-1,6)); // basis of invariants of |
| | 3678 | // degree d |
| | 3679 | if (B[1]<>0) |
| | 3680 | { Pplus=P+B; |
| | 3681 | newdim=dim(std(Pplus)); |
| | 3682 | dif=n-i-newdim; |
| | 3683 | } |
| | 3684 | else |
| | 3685 | { dif=0; |
| | 3686 | deg_vector=deg_vector,d; |
| | 3687 | } |
| | 3688 | if (dif<>0) // we have to find dif new primary |
| | 3689 | { // invariants |
| | 3690 | cd=size(B); |
| | 3691 | if (cd<>dif) |
| | 3692 | { P,CI,dB=p_search_random(n,d,B,cd,P,i,dif,dB,CI,max); |
| | 3693 | } |
| | 3694 | else // i.e. we can take all of B |
| | 3695 | { for(j=i+1;j<=i+dif;j=j+1) |
| | 3696 | { CI=CI+ideal(var(j)^d); |
| | 3697 | } |
| | 3698 | dB=dB+dif*(d-1); |
| | 3699 | P=Pplus; |
| | 3700 | } |
| | 3701 | if (ncols(P)==n+1) |
| | 3702 | { "WARNING: The first return value is not a set of primary invariants,"; |
| | 3703 | " but polynomials qualifying as the first "+string(i)+" primary invariants."; |
| | 3704 | return(matrix(P)); |
| | 3705 | } |
| | 3706 | if (v) |
| | 3707 | { for (j=1;j<=size(P)-i;j=j+1) |
| | 3708 | { " We find: "+string(P[i+j]); |
| | 3709 | } |
| | 3710 | } |
| | 3711 | i=size(P); |
| | 3712 | if (i==n) // found all primary invariants |
| | 3713 | { if (v) |
| | 3714 | { ""; |
| | 3715 | " We found all primary invariants."; |
| | 3716 | ""; |
| | 3717 | } |
| | 3718 | if (deg_vector==0) |
| | 3719 | { return(matrix(P)); |
| | 3720 | } |
| | 3721 | else |
| | 3722 | { return(matrix(P),compress(deg_vector)); |
| | 3723 | } |
| | 3724 | } |
| | 3725 | } // done with degree d |
| | 3726 | else |
| | 3727 | { if (v) |
| | 3728 | { " None here..."; |
| | 3729 | } |
| | 3730 | } |
| | 3731 | } |
| | 3732 | } |
| | 3733 | example |
| | 3734 | { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed into"; |
| | 3735 | " characteristic 3)"; |
| | 3736 | echo=2; |
| | 3737 | ring R=3,(x,y,z),dp; |
| | 3738 | matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; |
| | 3739 | list L=group_reynolds(A); |
| | 3740 | list l=primary_charp_no_molien_random(L[1],1); |
| | 3741 | print(l[1]); |
| | 3742 | } |
| | 3743 | |
| | 3744 | proc primary_charp_without_random (list #) |
| | 3745 | USAGE: primary_charp_without_random(G1,G2,...,r[,v]); |
| | 3746 | G1,G2,...: <matrices> generating a finite matrix group, r: an <int> |
| | 3747 | where -|r| to |r| is the range of coefficients of the random |
| | 3748 | combinations of bases elements, v: an optional <int> |
| | 3749 | DISPLAY: information about the various stages of the programme if v does not |
| | 3750 | equal 0 |
| | 3751 | RETURN: primary invariants (type <matrix>) of the invariant ring |
| | 3752 | EXAMPLE: example primary_charp_without_random; shows an example |
| | 3753 | THEORY: Bases of homogeneous invariants are generated successively and random |
| | 3754 | linear combinations are chosen as primary invariants that lower the |
| | 3755 | dimension of the ideal generated by the previously found invariants |
| | 3756 | (see paper "Generating a Noetherian Normalization of the Invariant Ring |
| | 3757 | of a Finite Group" by Decker, Heydtmann, Schreyer (1997) to appear in |
| | 3758 | JSC). No Reynolds operator or Molien series is used. |
| | 3759 | { degBound=0; |
| | 3760 | //--------------------- checking input and setting verbose mode -------------- |
| | 3761 | if (char(basering)==0) |
| | 3762 | { "ERROR: primary_charp_without_random should only be used with rings of"; |
| | 3763 | " characteristic 0."; |
| | 3764 | return(); |
| | 3765 | } |
| | 3766 | if (size(#)<2) |
| | 3767 | { "ERROR: There are too few parameters."; |
| | 3768 | return(); |
| | 3769 | } |
| | 3770 | if (typeof(#[size(#)])=="int" && typeof(#[size(#)-1])=="int") |
| | 3771 | { int v=#[size(#)]; |
| | 3772 | int max=#[size(#)-1]; |
| | 3773 | int gen_num=size(#)-2; |
| | 3774 | if (gen_num==0) |
| | 3775 | { "ERROR: There are no generators of a finite matrix group given."; |
| | 3776 | return(); |
| | 3777 | } |
| | 3778 | } |
| | 3779 | else |
| | 3780 | { if (typeof(#[size(#)])=="int") |
| | 3781 | { int max=#[size(#)]; |
| | 3782 | int v=0; |
| | 3783 | int gen_num=size(#)-1; |
| | 3784 | } |
| | 3785 | else |
| | 3786 | { "ERROR: The last parameter should be an <int>."; |
| | 3787 | return(); |
| | 3788 | } |
| | 3789 | } |
| | 3790 | int n=nvars(basering); // n is the number of variables, as well |
| | 3791 | // as the size of the matrices, as well |
| | 3792 | // as the number of primary invariants, |
| | 3793 | // we should get |
| | 3794 | for (int i=1;i<=gen_num;i=i+1) |
| | 3795 | { if (typeof(#[i])=="matrix") |
| | 3796 | { if (nrows(#[i])<>n or ncols(#[i])<>n) |
| | 3797 | { "ERROR: The number of variables of the base ring needs to be the same"; |
| | 3798 | " as the dimension of the square matrices"; |
| | 3799 | return(); |
| | 3800 | } |
| | 3801 | } |
| | 3802 | else |
| | 3803 | { "ERROR: The first parameters should be a list of matrices"; |
| | 3804 | return(); |
| | 3805 | } |
| | 3806 | } |
| | 3807 | //---------------------------------------------------------------------------- |
| | 3808 | if (v && voice==2) |
| | 3809 | { ""; |
| | 3810 | } |
| | 3811 | //---------------------------- initializing variables ------------------------ |
| | 3812 | int dB; |
| | 3813 | int j,d,cd,newdim,dif; // d: current degree, cd: dimension of |
| | 3814 | // space of invariants of degree d, |
| | 3815 | // newdim: dimension the ideal generated |
| | 3816 | // the primary invariants plus basis |
| | 3817 | // elements, dif=n-i-newdim, i.e. the |
| | 3818 | // number of new primary invairants that |
| | 3819 | // should be added in this degree - |
| | 3820 | ideal P,Pplus,CI,B; // P: will contain primary invariants, |
| | 3821 | // Pplus: P+B, CI: a complete |
| | 3822 | // intersection with the same Hilbert |
| | 3823 | // function as P - |
| | 3824 | dB=1; // used as degree bound - |
| | 3825 | d=0; // initializing |
| | 3826 | i=0; |
| | 3827 | intvec deg_vector; |
| | 3828 | //-------------------- loop that searches for primary invariants ------------- |
| | 3829 | while(1) // repeat until n primary invariants are |
| | 3830 | { // found - |
| | 3831 | d=d+1; // degree where we'll search |
| | 3832 | if (v) |
| | 3833 | { " Computing primary invariants in degree "+string(d)+":"; |
| | 3834 | } |
| | 3835 | B=invariant_basis(d,#[1..gen_num]); // basis of invariants of degree d |
| | 3836 | if (B[1]<>0) |
| | 3837 | { Pplus=P+B; |
| | 3838 | newdim=dim(std(Pplus)); |
| | 3839 | dif=n-i-newdim; |
| | 3840 | } |
| | 3841 | else |
| | 3842 | { dif=0; |
| | 3843 | deg_vector=deg_vector,d; |
| | 3844 | } |
| | 3845 | if (dif<>0) // we have to find dif new primary |
| | 3846 | { // invariants |
| | 3847 | cd=size(B); |
| | 3848 | if (cd<>dif) |
| | 3849 | { P,CI,dB=p_search_random(n,d,B,cd,P,i,dif,dB,CI,max); |
| | 3850 | } |
| | 3851 | else // i.e. we can take all of B |
| | 3852 | { for(j=i+1;j<=i+dif;j=j+1) |
| | 3853 | { CI=CI+ideal(var(j)^d); |
| | 3854 | } |
| | 3855 | dB=dB+dif*(d-1); |
| | 3856 | P=Pplus; |
| | 3857 | } |
| | 3858 | if (ncols(P)==n+1) |
| | 3859 | { "WARNING: The first return value is not a set of primary invariants,"; |
| | 3860 | " but polynomials qualifying as the first "+string(i)+" primary invariants."; |
| | 3861 | return(matrix(P)); |
| | 3862 | } |
| | 3863 | if (v) |
| | 3864 | { for (j=1;j<=size(P)-i;j=j+1) |
| | 3865 | { " We find: "+string(P[i+j]); |
| | 3866 | } |
| | 3867 | } |
| | 3868 | i=size(P); |
| | 3869 | if (i==n) // found all primary invariants |
| | 3870 | { if (v) |
| | 3871 | { ""; |
| | 3872 | " We found all primary invariants."; |
| | 3873 | ""; |
| | 3874 | } |
| | 3875 | return(matrix(P)); |
| | 3876 | } |
| | 3877 | } // done with degree d |
| | 3878 | else |
| | 3879 | { if (v) |
| | 3880 | { " None here..."; |
| | 3881 | } |
| | 3882 | } |
| | 3883 | } |
| | 3884 | } |
| | 3885 | example |
| | 3886 | { echo=2; |
| | 3887 | ring R=2,(x,y,z),dp; |
| | 3888 | matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; |
| | 3889 | matrix P=primary_charp_without_random(A,1); |
| | 3890 | print(P); |
| | 3891 | } |
| | 3892 | |
| | 3893 | proc primary_invariants_random (list #) |
| | 3894 | USAGE: primary_invariants_random(G1,G2,...,r[,flags]); |
| | 3895 | G1,G2,...: <matrices> generating a finite matrix group, r: an <int> |
| | 3896 | where -|r| to |r| is the range of coefficients of the random |
| | 3897 | combinations of bases elements, flags: an optional <intvec> with three |
| | 3898 | entries, if the first one equals 0 (also the default), the programme |
| | 3899 | attempts to compute the Molien series and Reynolds operator, if it |
| | 3900 | equals 1, the programme is told that the Molien series should not be |
| | 3901 | computed, if it equals -1 characteristic 0 is simulated, i.e. the |
| | 3902 | Molien series is computed as if the base field were characteristic 0 |
| | 3903 | (the user must choose a field of large prime characteristic, e.g. |
| | 3904 | 32003) and if the first one is anything else, it means that the |
| | 3905 | characteristic of the base field divides the group order, the second |
| | 3906 | component should give the size of intervals between canceling common |
| | 3907 | factors in the expansion of the Molien series, 0 (the default) means |
| | 3908 | only once after generating all terms, in prime characteristic also a |
| | 3909 | negative number can be given to indicate that common factors should |
| | 3910 | always be canceled when the expansion is simple (the root of the |
| | 3911 | extension field does not occur among the coefficients) |
| | 3912 | DISPLAY: information about the various stages of the programme if the third |
| | 3913 | flag does not equal 0 |
| | 3914 | RETURN: primary invariants (type <matrix>) of the invariant ring and if |
| | 3915 | computable Reynolds operator (type <matrix>) and Molien series (type |
| | 3916 | <matrix>), if the first flag is 1 and we are in the non-modular case |
| | 3917 | then an <intvec> is returned giving some of the degrees where no |
| | 3918 | non-trivial homogeneous invariants can be found |
| | 3919 | EXAMPLE: example primary_invariants_random; shows an example |
| | 3920 | THEORY: Bases of homogeneous invariants are generated successively and random |
| | 3921 | linear combinations are chosen as primary invariants that lower the |
| | 3922 | dimension of the ideal generated by the previously found invariants |
| | 3923 | (see paper "Generating a Noetherian Normalization of the Invariant Ring |
| | 3924 | of a Finite Group" by Decker, Heydtmann, Schreyer (1997) to appear in |
| | 3925 | JSC). |
| | 3926 | { |
| | 3927 | // ----------------- checking input and setting flags ------------------------ |
| | 3928 | if (size(#)<2) |
| | 3929 | { "ERROR: There are too few parameters."; |
| | 3930 | return(); |
| | 3931 | } |
| | 3932 | int ch=char(basering); // the algorithms depend very much on the |
| | 3933 | // characteristic of the ground field |
| | 3934 | int n=nvars(basering); // n is the number of variables, as well |
| | 3935 | // as the size of the matrices, as well |
| | 3936 | // as the number of primary invariants, |
| | 3937 | // we should get |
| | 3938 | int gen_num; |
| | 3939 | int mol_flag,v; |
| | 3940 | if (typeof(#[size(#)])=="intvec" && typeof(#[size(#)-1])=="int") |
| | 3941 | { if (size(#[size(#)])<>3) |
| | 3942 | { "ERROR: <intvec> should have three entries."; |
| | 3943 | return(); |
| | 3944 | } |
| | 3945 | gen_num=size(#)-2; |
| | 3946 | mol_flag=#[size(#)][1]; |
| | 3947 | if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0))) |
| | 3948 | { "ERROR: the second component of <intvec> should be >=0"; |
| | 3949 | return(); |
| | 3950 | } |
| | 3951 | int interval=#[size(#)][2]; |
| | 3952 | v=#[size(#)][3]; |
| | 3953 | int max=#[size(#)-1]; |
| | 3954 | if (gen_num==0) |
| | 3955 | { "ERROR: There are no generators of a finite matrix group given."; |
| | 3956 | return(); |
| | 3957 | } |
| | 3958 | } |
| | 3959 | else |
| | 3960 | { if (typeof(#[size(#)])=="int") |
| | 3961 | { gen_num=size(#)-1; |
| | 3962 | mol_flag=0; |
| | 3963 | int interval=0; |
| | 3964 | v=0; |
| | 3965 | int max=#[size(#)]; |
| | 3966 | } |
| | 3967 | else |
| | 3968 | { "ERROR: If the two last parameters are not <int> and <intvec>, the last"; |
| | 3969 | " parameter should be an <int>."; |
| | 3970 | return(); |
| | 3971 | } |
| | 3972 | } |
| | 3973 | for (int i=1;i<=gen_num;i=i+1) |
| | 3974 | { if (typeof(#[i])=="matrix") |
| | 3975 | { if (nrows(#[i])<>n or ncols(#[i])<>n) |
| | 3976 | { "ERROR: The number of variables of the base ring needs to be the same"; |
| | 3977 | " as the dimension of the square matrices"; |
| | 3978 | return(); |
| | 3979 | } |
| | 3980 | } |
| | 3981 | else |
| | 3982 | { "ERROR: The first parameters should be a list of matrices"; |
| | 3983 | return(); |
| | 3984 | } |
| | 3985 | } |
| | 3986 | //---------------------------------------------------------------------------- |
| | 3987 | if (mol_flag==0) |
| | 3988 | { if (ch==0) |
| | 3989 | { matrix REY,M=reynolds_molien(#[1..gen_num],intvec(0,interval,v)); |
| | 3990 | // one will contain Reynolds operator and |
| | 3991 | // the other enumerator and denominator |
| | 3992 | // of Molien series |
| | 3993 | matrix P=primary_char0_random(REY,M,max,v); |
| | 3994 | return(P,REY,M); |
| | 3995 | } |
| | 3996 | else |
| | 3997 | { list L=group_reynolds(#[1..gen_num],v); |
| | 3998 | if (L[1]<>0) // testing whether we are in the modular |
| | 3999 | { string newring="aksldfalkdsflkj"; // case |
| | 4000 | if (minpoly==0) |
| | 4001 | { if (v) |
| | 4002 | { " We are dealing with the non-modular case."; |
| | 4003 | } |
| | 4004 | molien(L[2..size(L)],newring,intvec(0,interval,v)); |
| | 4005 | matrix P=primary_charp_random(L[1],newring,max,v); |
| | 4006 | return(P,L[1],newring); |
| | 4007 | } |
| | 4008 | else |
| | 4009 | { if (v) |
| | 4010 | { " Since it is impossible for this programme to calculate the Molien series for"; |
| | 4011 | " invariant rings over extension fields of prime characteristic, we have to"; |
| | 4012 | " continue without it."; |
| | 4013 | ""; |
| | 4014 | |
| | 4015 | } |
| | 4016 | list l=primary_charp_no_molien_random(L[1],max,v); |
| | 4017 | if (size(l)==2) |
| | 4018 | { return(l[1],L[1],l[2]); |
| | 4019 | } |
| | 4020 | else |
| | 4021 | { return(l[1],L[1]); |
| | 4022 | } |
| | 4023 | } |
| | 4024 | } |
| | 4025 | else // the modular case |
| | 4026 | { if (v) |
| | 4027 | { " There is also no Molien series, we can make use of..."; |
| | 4028 | ""; |
| | 4029 | " We can start looking for primary invariants..."; |
| | 4030 | ""; |
| | 4031 | } |
| | 4032 | return(primary_charp_without_random(#[1..gen_num],max,v)); |
| | 4033 | } |
| | 4034 | } |
| | 4035 | } |
| | 4036 | if (mol_flag==1) // the user wants no calculation of the |
| | 4037 | { list L=group_reynolds(#[1..gen_num],v); // Molien series |
| | 4038 | if (ch==0) |
| | 4039 | { list l=primary_char0_no_molien_random(L[1],max,v); |
| | 4040 | if (size(l)==2) |
| | 4041 | { return(l[1],L[1],l[2]); |
| | 4042 | } |
| | 4043 | else |
| | 4044 | { return(l[1],L[1]); |
| | 4045 | } |
| | 4046 | } |
| | 4047 | else |
| | 4048 | { if (L[1]<>0) // testing whether we are in the modular |
| | 4049 | { list l=primary_charp_no_molien_random(L[1],max,v); // case |
| | 4050 | if (size(l)==2) |
| | 4051 | { return(l[1],L[1],l[2]); |
| | 4052 | } |
| | 4053 | else |
| | 4054 | { return(l[1],L[1]); |
| | 4055 | } |
| | 4056 | } |
| | 4057 | else // the modular case |
| | 4058 | { if (v) |
| | 4059 | { " We can start looking for primary invariants..."; |
| | 4060 | ""; |
| | 4061 | } |
| | 4062 | return(primary_charp_without_random(#[1..gen_num],max,v)); |
| | 4063 | } |
| | 4064 | } |
| | 4065 | } |
| | 4066 | if (mol_flag==-1) |
| | 4067 | { if (ch==0) |
| | 4068 | { "ERROR: Characteristic 0 can only be simulated in characteristic p>>0."; |
| | 4069 | return(); |
| | 4070 | } |
| | 4071 | list L=group_reynolds(#[1..gen_num],v); |
| | 4072 | string newring="aksldfalkdsflkj"; |
| | 4073 | molien(L[2..size(L)],newring,intvec(1,interval,v)); |
| | 4074 | matrix P=primary_charp_random(L[1],newring,max,v); |
| | 4075 | return(P,L[1],newring); |
| | 4076 | } |
| | 4077 | else // the user specified that the |
| | 4078 | { if (ch==0) // characteristic divides the group order |
| | 4079 | { "ERROR: The characteristic cannot divide the group order when it is 0."; |
| | 4080 | return(); |
| | 4081 | } |
| | 4082 | if (v) |
| | 4083 | { ""; |
| | 4084 | } |
| | 4085 | return(primary_charp_without_random(#[1..gen_num],max,v)); |
| | 4086 | } |
| | 4087 | } |
| | 4088 | example |
| | 4089 | { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; |
| | 4090 | echo=2; |
| | 4091 | ring R=0,(x,y,z),dp; |
| | 4092 | matrix A[3][3]=0,1,0,-1,0,0,0,0,-1; |
| | 4093 | list L=primary_invariants_random(A,1); |
| | 4094 | print(L[1]); |
| | 4095 | } |
| | 4096 | |
| | 4097 | proc concat_intmat(intmat A,intmat B) |
| | 4098 | { int n=nrows(A); |
| | 4099 | int m1=ncols(A); |
| | 4100 | int m2=ncols(B); |
| | 4101 | intmat C[n][m1+m2]; |
| | 4102 | C[1..n,1..m1]=A[1..n,1..m1]; |
| | 4103 | C[1..n,m1+1..m1+m2]=B[1..n,1..m2]; |
| | 4104 | return(C); |
| | 4105 | } |
| | 4106 | |
| | 4107 | proc power_products(intvec deg_vec,int d) |
| | 4108 | USAGE: power_products(dv,d); |
| | 4109 | dv: an <intvec> giving the degrees of homogeneous polynomials, d: the |
| | 4110 | degree of the desired power products |
| | 4111 | RETURN: a size(dv)*m <intmat> where each column ought to be interpreted as |
| | 4112 | containing the exponents of the corresponding polynomials. The product |
| | 4113 | of the powers is then homogeneous of degree d. |
| | 4114 | EXAMPLE: example power_products; gives an example |
| | 4115 | { if (d<=0) |
| | 4116 | { "ERROR: The <int> may not be <= 0"; |
| | 4117 | return(); |
| | 4118 | } |
| | 4119 | int d_neu,j,nc; |
| | 4120 | int s=size(deg_vec); |
| | 4121 | intmat PP[s][1]; |
| | 4122 | intmat TEST[s][1]; |
| | 4123 | for (int i=1;i<=s;i=i+1) |
| | 4124 | { if (i<0) |
| | 4125 | { "ERROR: The entries of <intvec> may not be <= 0"; |
| | 4126 | return(); |
| | 4127 | } |
| | 4128 | d_neu=d-deg_vec[i]; |
| | 4129 | if (d_neu>0) |
| | 4130 | { intmat PPd_neu=power_products(intvec(deg_vec[i..s]),d_neu); |
| | 4131 | if (size(ideal(PPd_neu))<>0) |
| | 4132 | { nc=ncols(PPd_neu); |
| | 4133 | intmat PPd_neu_gross[s][nc]; |
| | 4134 | PPd_neu_gross[i..s,1..nc]=PPd_neu[1..s-i+1,1..nc]; |
| | 4135 | for (j=1;j<=nc;j=j+1) |
| | 4136 | { PPd_neu_gross[i,j]=PPd_neu_gross[i,j]+1; |
| | 4137 | } |
| | 4138 | PP=concat_intmat(PP,PPd_neu_gross); |
| | 4139 | kill PPd_neu_gross; |
| | 4140 | } |
| | 4141 | kill PPd_neu; |
| | 4142 | } |
| | 4143 | if (d_neu==0) |
| | 4144 | { intmat PPd_neu[s][1]; |
| | 4145 | PPd_neu[i,1]=1; |
| | 4146 | PP=concat_intmat(PP,PPd_neu); |
| | 4147 | kill PPd_neu; |
| | 4148 | } |
| | 4149 | } |
| | 4150 | if (matrix(PP)<>matrix(TEST)) |
| | 4151 | { PP=compress(PP); |
| | 4152 | } |
| | 4153 | return(PP); |
| | 4154 | } |
| | 4155 | example |
| | 4156 | { echo=2; |
| | 4157 | intvec dv=5,5,5,10,10; |
| | 4158 | print(power_products(dv,10)); |
| | 4159 | print(power_products(dv,7)); |
| | 4160 | } |
| | 4161 | |
| | 4162 | proc secondary_char0 (matrix P, matrix REY, matrix M, list #) |
| | 4163 | USAGE: secondary_char0(P,REY,M[,v]); |
| | 4164 | P: a 1xn <matrix> with primary invariants, REY: a gxn <matrix> |
| | 4165 | representing the Reynolds operator, M: a 1x2 <matrix> giving enumerator |
| | 4166 | and denominator of the Molien series, v: an optional <int> |
| | 4167 | ASSUME: n is the number of variables of the basering, g the size of the group, |
| | 4168 | REY is the 1st return value of group_reynolds(), reynolds_molien() or |
| | 4169 | the second one of primary_invariants(), M the return value of molien() |
| | 4170 | or the second one of reynolds_molien() or the third one of |
| | 4171 | primary_invariants() |
| | 4172 | RETURN: secondary invariants of the invariant ring (type <matrix>) and |
| | 4173 | irreducible secondary invariants (type <matrix>) |
| | 4174 | DISPLAY: information if v does not equal 0 |
| | 4175 | EXAMPLE: example secondary_char0; shows an example |
| | 4176 | THEORY: The secondary invariants are calculated by finding a basis (in terms of |
| | 4177 | monomials) of the basering modulo the primary invariants, mapping those |
| | 4178 | to invariants with the Reynolds operator and using these images or |
| | 4179 | their power products such that they are linearly independent modulo the |
| | 4180 | primary invariants (see paper "Some Algorithms in Invariant Theory of |
| | 4181 | Finite Groups" by Kemper and Steel (1997)). |
