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)). |