source: git/Singular/LIB/finvar.lib @ 56f1845

spielwiese
Last change on this file since 56f1845 was 56f1845, checked in by Hans Schönemann <hannes@…>, 26 years ago
* agnes: new version of finvar.lib - complete rewrite git-svn-id: file:///usr/local/Singular/svn/trunk@966 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 214.2 KB
Line 
1// $Header: /usr/local/Singular/cvsroot/Singular/LIB/finvar.lib,v 1.1.4.1 1997/0
27/03 15:32:03 Singular Exp $
3////////////////////////////////////////////////////////////////////////////////
4
5LIBRARY:  finvar.lib         LIBRARY TO CALCULATE INVARIANT RINGS & MORE
6                             (c) Agnes Eileen Heydtmann,
7                             send bugs and comments to agnes@math.uni-sb.de
8
9 cyclotomic                                cyclotomic polynomial
10 group_reynolds                            finite group and Reynolds operator
11 molien                                    Molien series
12 reynolds_molien                           Reynolds operator and Molien series
13 partial_molien                            partial expansion of Molien series
14 evaluate_reynolds                         image under the Reynolds operator
15 invariant_basis                           basis of homogeneous invariants
16 invariant_basis_reynolds                  basis of homogeneous invariants
17 primary_char0                             primary invariants
18 primary_charp                             primary invariants
19 primary_char0_no_molien                   primary invariants
20 primary_charp_no_molien                   primary invariants
21 primary_charp_without                     primary invariants
22 primary_invariants                        primary invariants
23 primary_char0_random                      primary invariants
24 primary_charp_random                      primary invariants
25 primary_char0_no_molien_random            primary invariants
26 primary_charp_no_molien_random            primary invariants
27 primary_charp_without_random              primary invariants
28 primary_invariants_random                 primary invariants
29 power_products                            exponents for power products
30 secondary_char0                           secondary invariants
31 secondary_charp                           secondary invariants
32 secondary_no_molien                       secondary invariants
33 secondary_with_irreducible_ones_no_molien secondary invariants
34 secondary_not_cohen_macaulay              secondary invariants
35 invariant_ring                            primary and secondary invariants
36 invariant_ring_random                     primary and secondary invariants
37 algebra_containment                       answers query of algebra containment
38 module_containment                        answers query of module containment
39 orbit_variety                             ideal of the orbit variety
40 relative_orbit_variety                    ideal of a relative orbit variety
41 image_of_variety                          ideal of the image of a variety
42
43////////////////////////////////////////////////////////////////////////////////
44LIB "matrix.lib";
45LIB "elim.lib";
46LIB "general.lib";
47////////////////////////////////////////////////////////////////////////////////
48
49////////////////////////////////////////////////////////////////////////////////
50// Checks whether the last parameter, being a matrix, is among the previous
51// parameters, also being matrices
52////////////////////////////////////////////////////////////////////////////////
53proc unique (list #)
54{ for (int i=1;i<size(#);i=i+1)
55  { if (#[i]==#[size(#)])
56    { return(0);
57    }
58  }
59  return(1);
60}
61
62proc cyclotomic (int i)
63USAGE:   cyclotomic(i);
64         i: an <int> > 0
65RETURNS: the i-th cyclotomic polynomial (type <poly>) as one in the first ring
66         variable
67EXAMPLE: example cyclotomic; shows an example
68THEORY:  x^i-1 is divided by the j-th cyclotomic polynomial where j takes on the
69         value of proper divisors of i
70{ if (i<=0)
71  { "ERROR:   the input should be > 0.";
72    return();
73  }
74  poly v1=var(1);
75  if (i==1)
76  { return(v1-1);                      // 1-st cyclotomic polynomial
77  }
78  poly min=v1^i-1;
79  matrix s[1][2];
80  min=min/(v1-1);                      // dividing by the 1-st cyclotomic
81                                       // polynomial
82  int j=2;
83  int n;
84  poly c;
85  int flag=1;
86  while(2*j<=i)                        // there are no proper divisors of i
87  { if ((i%j)==0)                      // greater than i/2
88    { if (flag==1)
89      { n=j;                           // n stores the first proper divisor of
90      }                                // i > 1
91      flag=0;
92      c=cyclotomic(j);                 // recursive computation
93      s=min,c;
94      s=matrix(syz(ideal(s)));         // dividing
95      min=s[2,1];
96    }
97    if (n*j==i)                        // the earliest possible point to break
98    { break;
99    }
100    j=j+1;
101  }
102  min=min/leadcoef(min);               // making sure that the leading
103  return(min);                         // coefficient is 1
104}
105example
106{ echo=2;
107          ring R=0,(x,y,z),dp;
108          print(cyclotomic(25));
109}
110
111proc group_reynolds (list #)
112USAGE:   group_reynolds(G1,G2,...[,v]);
113         G1,G2,...: nxn <matrices> generating a finite matrix group, v: an
114         optional <int>
115ASSUME:  n is the number of variables of the basering, g the number of group
116         elements
117RETURN:  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
123DISPLAY: information if v does not equal 0
124EXAMPLE: example group_reynolds; shows an example
125THEORY:  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";
141      return();
142    }
143    int v=#[size(#)];
144    gen_num=size(#)-1;
145  }
146  else                                 // last parameter is not <int>
147  { int v=0;                           // no information is default
148    gen_num=size(#);
149  }
150  if (typeof(#[1])<>"matrix")
151  { "ERROR:   the parameters must be a list of matrices and maybe an <int>";
152    return();
153  }
154  int n=nrows(#[1]);
155  if (n<>nvars(basering))
156  { "ERROR:   the number of variables of the basering needs to be the same";
157    "         as the dimension of the matrices";
158    return();
159  }
160  if (n<>ncols(#[1]))
161  { "ERROR:   matrices need to be square and of the same dimensions";
162    return();
163  }
164  matrix vars=matrix(maxideal(1));     // creating an nx1-matrix containing the
165  vars=transpose(vars);                // variables of the ring -
166  matrix REY=#[1]*vars;                // calculating the first ring mapping -
167                                       // REY will contain the Reynolds
168                                       // operator -
169  matrix G(1)=#[1];                    // G(k) are elements of the group -
170  if (ch<>0 && minpoly==0 && gen_num<>1) // finding out of which order the group
171  { matrix I=diag(1,n);                // element is
172    matrix TEST=G(1);
173    int o1=1;
174    int o2;
175    while (TEST<>I)
176    { TEST=TEST*G(1);
177      o1=o1+1;
178    }
179  }
180  int i=1;
181 // -------------- doubles among the generators should be avoided -------------
182  for (int j=2;j<=gen_num;j=j+1)       // this loop adds the parameters to the
183  {                                    // group, leaving out doubles and
184                                       // checking whether the parameters are
185                                       // compatible with the task of the
186                                       // procedure
187    if (not(typeof(#[j])=="matrix"))
188    { "ERROR:   the parameters must be a list of matrices and maybe an <int>";
189      return();
190    }
191    if ((n!=nrows(#[j])) or (n!=ncols(#[j])))
192    { "ERROR:   matrices need to be square and of the same dimensions";
193       return();
194    }
195    if (unique(G(1..i),#[j]))
196    { i=i+1;
197      matrix G(i)=#[j];
198      if (ch<>0 && minpoly==0)         // finding out of which order the group
199      { TEST=G(i);                     // element is
200        o2=1;
201        while (TEST<>I)
202        { TEST=TEST*G(i);
203          o2=o2+1;
204        }
205        o1=o1*o2/gcd(o1,o2);           // lowest common multiple of the element
206      }                                // orders -
207      REY=concat(REY,#[j]*vars);       // adding ring homomorphisms to REY
208    }
209  }
210  int g=i;                             // G(1)..G(i) are generators without
211                                       // doubles - g generally is the number
212                                       // of elements in the group so far -
213  j=i;                                 // j is the number of new elements that
214                                       // we use as factors
215  int k, m, l;
216  if (v)
217  { "";
218    "  Generating the entire matrix group and the Reynolds operator...";
219    "";
220  }
221 // -------------- main loop that finds all the group elements ----------------
222  while (1)
223  { l=0;                               // l is the number of products we get in
224                                       // one going
225    for (m=g-j+1;m<=g;m=m+1)
226    { for (k=1;k<=i;k=k+1)
227      { l=l+1;
228        matrix P(l)=G(k)*G(m);         // possible new element
229      }
230    }
231    j=0;
232    for (k=1;k<=l;k=k+1)
233    { if (unique(G(1..g),P(k)))
234      { j=j+1;                         // a new factor for next run
235        g=g+1;
236        matrix G(g)=P(k);              // a new group element -
237        if (ch<>0 && minpoly==0 && i<>1) // finding out of which order the group
238        { TEST=G(g);                   // element is
239          o2=1;
240          while (TEST<>I)
241          { TEST=TEST*G(g);
242            o2=o2+1;
243          }
244          o1=o1*o2/gcd(o1,o2);         // lowest common multiple of the element
245        }                              // orders -
246        REY=concat(REY,P(k)*vars);     // adding new mapping to REY
247        if (v)
248        { "  Group element "+string(g)+" has been found.";
249        }
250      }
251      kill P(k);
252    }
253    if (j==0)                          // when we didn't add any new elements
254    { break;                           // in one run through the while loop
255    }                                  // we are done
256  }
257  if (v)
258  { if (g<=i)
259    { "  There are only "+string(g)+" group elements.";
260    }
261    "";
262  }
263  REY=transpose(REY);                  // when we evaluate the Reynolds operator
264                                       // later on, we actually want 1xn
265                                       // matrices
266  if (ch<>0)
267  { if ((g%ch)==0)
268    { if (voice==2)
269      { "WARNING: The characteristic of the coefficient field divides the group order.";
270        "         Proceed without the Reynolds operator!";
271      }
272      else
273      { if (v)
274        { "  The characteristic of the base field divides the group order.";
275          "  We have to continue without Reynolds operator...";
276          "";
277        }
278      }
279      kill REY;
280      matrix REY[1][1]=0;
281      return(REY,G(1..g));
282    }
283    if (minpoly==0)
284    { if (i>1)
285      { return(REY,o1,G(1..g));
286      }
287      return(REY,G(1..g));
288    }
289  }
290  if (v)
291  { "  Done generating the group and the Reynolds operator.";
292    "";
293  }
294  return(REY,G(1..g));
295}
296example
297{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
298  echo=2;
299         ring R=0,(x,y,z),dp;
300         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
301         list L=group_reynolds(A);
302         print(L[1]);
303         print(L[2..size(L)]);
304}
305
306////////////////////////////////////////////////////////////////////////////////
307// Returns i such that root^i==n, i.e. it heavily relies on the right input.
308////////////////////////////////////////////////////////////////////////////////
309proc exponent(number n, number root)
310{ int i=0;
311   while((n/root^i)<>1)
312   { i=i+1;
313   }
314   return(i);
315}
316
317proc molien (list #)
318USAGE:   molien(G1,G2,...[,ringname,lcm,flags]);
319         G1,G2,...: nxn <matrices> generating a finite matrix group, ringname:
320         a <string> giving a name for a new ring of characteristic 0 for the
321         Molien series in case of prime characteristic, lcm: an <int> giving the
322         lowest common multiple of the elements' orders in case of prime
323         characteristic, minpoly==0 and a non-cyclic group, flags: an optional
324         <intvec> with three components: if the first element is not equal to 0
325         characteristic 0 is simulated, i.e. the Molien series is computed as
326         if the base field were characteristic 0 (the user must choose a field
327         of large prime characteristic, e.g. 32003), the second component should
328         give the size of intervals between canceling common factors in the
329         expansion of the Molien series, 0 (the default) means only once after
330         generating all terms, in prime characteristic also a negative number
331         can be given to indicate that common factors should always be canceled
332         when the expansion is simple (the root of the extension field does not
333         occur among the coefficients)
334ASSUME:  n is the number of variables of the basering, G1,G2... are the group
335         elements generated by group_reynolds(), lcm is the second return value
336         of group_reynolds()
337RETURN:  in case of characteristic 0 a 1x2 <matrix> giving enumerator and
338         denominator of Molien series; in case of prime characteristic a ring
339         with the name `ringname` of characteristic 0 is created where the same
340         Molien series (named M) is stored
341DISPLAY: information if the third component of flags does not equal 0
342EXAMPLE: example molien; shows an example
343THEORY:  In characteristic 0 the terms 1/det(1-xE) for all group elements of the
344         Molien series are computed in a straight forward way. In prime
345         characteristic a Brauer lift is involved. The returned matrix gives
346         enumerator and denominator of the expanded version where common factors
347         have been canceled.
348{ def br=basering;                     // the Molien series depends on the
349  int ch=char(br);                     // characteristic of the coefficient
350                                       // field -
351  int g;                               // size of the group
352 //---------------------- making sure the input is okay -----------------------
353  if (typeof(#[size(#)])=="intvec")
354  { if (size(#[size(#)])==3)
355    { int mol_flag=#[size(#)][1];
356      if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0)))
357      { "ERROR:   the second component of <intvec> should be >=0"
358        return();
359      }
360      int interval=#[size(#)][2];
361      int v=#[size(#)][3];
362    }
363    else
364    { "ERROR:   <intvec> should have three components";
365      return();
366    }
367    if (ch<>0)
368    { if (typeof(#[size(#)-1])=="int")
369      { int r=#[size(#)-1];
370        if (typeof(#[size(#)-2])<>"string")
371        { "ERROR:   in characteristic p>0 a <string> must be given for the name of a new";
372          "         ring where the Molien series can be stored";
373          return();
374        }
375        else
376        { if (#[size(#)-2]=="")
377          { "ERROR:   <string> may not be empty";
378            return();
379          }
380          string newring=#[size(#)-2];
381          g=size(#)-3;
382        }
383      }
384      else
385      { if (typeof(#[size(#)-1])<>"string")
386        { "ERROR:   in characteristic p>0 a <string> must be given for the name of a new";
387          "         ring where the Molien series can be stored";
388          return();
389        }
390        else
391        { if (#[size(#)-1]=="")
392          { "ERROR:   <string> may not be empty";
393            return();
394          }
395          string newring=#[size(#)-1];
396          g=size(#)-2;
397          int r=g;
398        }
399      }
400    }
401    else                               // then <string> ist not needed
402    { g=size(#)-1;
403    }
404  }
405  else                                 // last parameter is not <intvec>
406  { int v=0;                           // no information is default
407    int mol_flag=0;                    // computing of Molien series is default
408    int interval=0;
409    if (ch<>0)
410    { if (typeof(#[size(#)])=="int")
411      { int r=#[size(#)];
412        if (typeof(#[size(#)-1])<>"string")
413        { "ERROR:   in characteristic p>0 a <string> must be given for the name of a new";
414          "         ring where the Molien series can be stored";
415            return();
416        }
417        else
418        { if (#[size(#)-1]=="")
419          { "ERROR:   <string> may not be empty";
420            return();
421          }
422          string newring=#[size(#)-1];
423          g=size(#)-2;
424        }
425      }
426      else
427      { if (typeof(#[size(#)])<>"string")
428        { "ERROR:   in characteristic p>0 a <string> must be given for the name of a new";
429          "         ring where the Molien series can be stored";
430          return();
431        }
432        else
433        { if (#[size(#)]=="")
434          { "ERROR:   <string> may not be empty";
435            return();
436          }
437          string newring=#[size(#)];
438          g=size(#)-1;
439          int r=g;
440        }
441      }
442    }
443    else
444    { g=size(#);
445    }
446  }
447  if (ch<>0)
448  { if ((g/r)*r<>g)
449   { "ERROR:   <int> should divide the group order."
450      return();
451    }
452  }
453  if (ch<>0)
454  { if ((g%ch)==0)
455    { if (voice==2)
456      { "WARNING: The characteristic of the coefficient field divides the group";
457        "         order. Proceed without the Molien series!";
458      }
459      else
460      { if (v)
461        { "  The characteristic of the base field divides the group order.";
462          "  We have to continue without Molien series...";
463          "";
464        }
465      }
466    }
467    if (minpoly<>0 && mol_flag==0)
468    { if (voice==2)
469      { "WARNING: It is impossible for this program to calculate the Molien series";
470        "         for finite groups over extension fields of prime characteristic.";
471      }
472      else
473      { if (v)
474        { "  Since it is impossible for this program to calculate the Molien series for";
475          "  invariant rings over extension fields of prime characteristic, we have to";
476          "  continue without it.";
477          "";
478        }
479      }
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}
835example
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
855proc reynolds_molien (list #)
856USAGE:   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)
870ASSUME:  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
872RETURN:  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
877DISPLAY: information if the third component of flags does not equal 0
878EXAMPLE: example reynolds_molien; shows an example
879THEORY:  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)
957    { if (voice==2)
958      { "WARNING: The characteristic of the coefficient field divides the group order.";
959        "         Proceed without the Reynolds operator or the Molien series!";
960        return();
961      }
962      if (v)
963      { "  The characteristic of the base field divides the group order.";
964        "  We have to continue without Reynolds operator or the Molien series...";
965        return();
966      }
967    }
968    if (minpoly<>0)
969    { if (voice==2)
970      { "WARNING: It is impossible for this program to calculate the Molien series";
971        "         for finite groups over extension fields of prime characteristic.";
972        return(L[1]);
973      }
974      else
975      { if (v)
976        { "  Since it is impossible for this program to calculate the Molien series for";
977          "  invariant rings over extension fields of prime characteristic, we have to";
978          "  continue without it.";
979          return(L[1]);
980        }
981      }
982    }
983    if (typeof(L[2])=="int")
984    { molien(L[3..size(L)],newring,L[2],intvec(mol_flag,interval,v));
985    }
986    else
987    { molien(L[2..size(L)],newring,intvec(mol_flag,interval,v));
988    }
989    return(L[1]);
990  }
991 //----------- computing Molien series in the straight forward way ------------
992  if (ch==0)
993  { if (typeof(#[1])<>"matrix")
994    { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
995      return();
996    }
997    int n=nrows(#[1]);
998    if (n<>nvars(br))
999    { "ERROR:   the number of variables of the basering needs to be the same";
1000      "         as the dimension of the matrices";
1001      return();
1002    }
1003    if (n<>ncols(#[1]))
1004    { "ERROR:   matrices need to be square and of the same dimensions";
1005      return();
1006    }
1007    matrix vars=matrix(maxideal(1));   // creating an nx1-matrix containing the
1008    vars=transpose(vars);              // variables of the ring -
1009    matrix A(1)=#[1]*vars;             // calculating the first ring mapping -
1010                                       // A(1) will contain the Reynolds
1011                                       // operator -
1012    poly v1=vars[1,1];                 // the Molien series will be in terms of
1013                                       // the first variable of the current
1014                                       // ring
1015    matrix I=diag(1,n);
1016    matrix A(2)[1][2];                 // A(2) will contain the Molien series -
1017    A(2)[1,1]=1;                       // A(2)[1,1] will be the numerator
1018    matrix G(1)=#[1];                  // G(k) are elements of the group -
1019    A(2)[1,2]=det(I-v1*(G(1)));        // A(2)[1,2] will be the denominator -
1020    matrix s;                          // will help us canceling in the
1021                                       // fraction
1022    poly p;                            // will contain the denominator of the
1023                                       // new term of the Molien series
1024    int i=1;
1025    for (int j=2;j<=gen_num;j=j+1)     // this loop adds the parameters to the
1026    {                                  // group, leaving out doubles and
1027                                       // checking whether the parameters are
1028                                       // compatible with the task of the
1029                                       // procedure
1030      if (not(typeof(#[j])=="matrix"))
1031      { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
1032        return();
1033      }
1034      if ((n!=nrows(#[j])) or (n!=ncols(#[j])))
1035      { "ERROR:   matrices need to be square and of the same dimensions";
1036         return();
1037      }
1038      if (unique(G(1..i),#[j]))
1039      { i=i+1;
1040        matrix G(i)=#[j];
1041        A(1)=concat(A(1),#[j]*vars);   // adding ring homomorphisms to A(1) -
1042        p=det(I-v1*#[j]);              // denominator of new term -
1043        A(2)[1,1]=A(2)[1,1]*p+A(2)[1,2]; // expanding A(2)[1,1]/A(2)[1,2] +1/p
1044        A(2)[1,2]=A(2)[1,2]*p;
1045        if (interval<>0)               // canceling common terms of denominator
1046        { if ((i/interval)*interval==i) // and enumerator
1047          {
1048            s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz() these
1049            A(2)[1,1]=-s[2,1];         // three lines should be replaced by the
1050            A(2)[1,2]=s[1,1];          // following three
1051            // p=gcd(A(2)[1,1],A(2)[1,2]);
1052            // A(2)[1,1]=A(2)[1,1]/p;
1053            // A(2)[1,2]=A(2)[1,2]/p;
1054          }
1055        }
1056      }
1057    }
1058    int g=i;                           // G(1)..G(i) are generators without
1059                                       // doubles - g generally is the number
1060                                       // of elements in the group so far -
1061    j=i;                               // j is the number of new elements that
1062                                       // we use as factors
1063    int k, m, l;
1064    if (v)
1065    { "";
1066      "  Generating the entire matrix group. Whenever a new group element is found,";
1067      "  the coressponding ring homomorphism of the Reynolds operator and the";
1068      "  corresponding term of the Molien series is generated.";
1069      "";
1070    }
1071 //----------- computing 1/det(I-xE) whenever a new element E is found --------
1072    while (1)
1073    { l=0;                             // l is the number of products we get in
1074                                       // one going
1075      for (m=g-j+1;m<=g;m=m+1)
1076      { for (k=1;k<=i;k=k+1)
1077        { l=l+1;
1078          matrix P(l)=G(k)*G(m);       // possible new element
1079        }
1080      }
1081      j=0;
1082      for (k=1;k<=l;k=k+1)
1083      { if (unique(G(1..g),P(k)))
1084        { j=j+1;                       // a new factor for next run
1085          g=g+1;
1086          matrix G(g)=P(k);            // a new group element -
1087          A(1)=concat(A(1),P(k)*vars); // adding new mapping to A(1)
1088          p=det(I-v1*P(k));            // denominator of new term
1089          A(2)[1,1]=A(2)[1,1]*p+A(2)[1,2];
1090          A(2)[1,2]=A(2)[1,2]*p;       // expanding A(2)[1,1]/A(2)[1,2] + 1/p -
1091          if (interval<>0)             // canceling common terms of denominator
1092          { if ((g/interval)*interval==g) // and enumerator
1093            {
1094              s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz()
1095              A(2)[1,1]=-s[2,1];       // these three lines should be replaced
1096              A(2)[1,2]=s[1,1];        // by the following three
1097              // p=gcd(A(2)[1,1],A(2)[1,2]);
1098              // A(2)[1,1]=A(2)[1,1]/p;
1099              // A(2)[1,2]=A(2)[1,2]/p;
1100            }
1101          }
1102          if (v)
1103          { "  Group element "+string(g)+" has been found.";
1104          }
1105        }
1106        kill P(k);
1107      }
1108      if (j==0)                        // when we didn't add any new elements
1109      { break;                         // in one run through the while loop
1110      }                                // we are done
1111    }
1112    if (v)
1113    { if (g<=i)
1114      { "  There are only "+string(g)+" group elements.";
1115      }
1116      "";
1117    }
1118    A(1)=transpose(A(1));              // when we evaluate the Reynolds operator
1119                                       // later on, we actually want 1xn
1120                                       // matrices
1121    if (interval==0)                   // canceling common terms of denominator
1122    {                                  // and enumerator -
1123      s=matrix(syz(ideal(A(2))));      // once gcd() is faster than syz()
1124      A(2)[1,1]=-s[2,1];               // these three lines should be replaced
1125      A(2)[1,2]=s[1,1];                // by the following three
1126      // p=gcd(A(2)[1,1],A(2)[1,2]);
1127      // A(2)[1,1]=A(2)[1,1]/p;
1128      // A(2)[1,2]=A(2)[1,2]/p;
1129    }
1130    if (interval<>0)                   // canceling common terms of denominator
1131    { if ((g/interval)*interval<>g)    // and enumerator
1132      {
1133        s=matrix(syz(ideal(A(2))));    // once gcd() is faster than syz()
1134        A(2)[1,1]=-s[2,1];             // these three lines should be replaced
1135        A(2)[1,2]=s[1,1];              // by the following three
1136        // p=gcd(A(2)[1,1],A(2)[1,2]);
1137        // A(2)[1,1]=A(2)[1,1]/p;
1138        // A(2)[1,2]=A(2)[1,2]/p;
1139      }
1140    }
1141    map slead=br,ideal(0);
1142    s=slead(A(2));
1143    A(2)[1,1]=1/s[1,1]*A(2)[1,1];      // numerator and denominator have to have
1144    A(2)[1,2]=1/s[1,2]*A(2)[1,2];      // a constant term of 1
1145    if (v)
1146    { "  Now we are done calculating Molien series and Reynolds operator.";
1147      "";
1148    }
1149    return(A(1..2));
1150  }
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  }
1379}
1380example
1381{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
1382  "         note the case of prime characteristic";
1383  echo=2;
1384         ring R=0,(x,y,z),dp;
1385         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
1386         matrix REY,M=reynolds_molien(A);
1387         print(REY);
1388         print(M);
1389         ring S=3,(x,y,z),dp;
1390         string newring="Qadjoint";
1391         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
1392         matrix REY=reynolds_molien(A,newring);
1393         print(REY);
1394         setring Qadjoint;
1395         print(M);
1396         setring S;
1397         kill Qadjoint;
1398}
1399
1400proc partial_molien (matrix M, int n, list #)
1401USAGE:   partial_molien(M,n[,p]);
1402         M: a 1x2 <matrix>, n: an <int> indicating  number of terms in the
1403         expansion, p: an optional <poly>
1404ASSUME:  M is the return value of molien or the second return value of
1405         reynolds_molien, p ought to be the second return value of a previous
1406         run of partial_molien and avoids recalculating known terms
1407RETURN:  n terms (type <poly>) of the partial expansion of the Molien series
1408         (first n if there is no third parameter given, otherwise the next n
1409         terms depending on a previous calculation) and an intermediate result
1410         (type <poly>) of the calculation to be used as third parameter in a next
1411         run of partial_molien
1412THEORY:  The following calculation is implemented:
1413         (1+a1x+a2x^2+...+anx^n)/(1+b1x+b2x^2+...+bmx^m)=(1+(a1-b1)x+...
1414         (1+b1x+b2x^2+...+bmx^m)
1415         -----------------------
1416            (a1-b1)x+(a2-b2)x^2+...
1417            (a1-b1)x+b1(a1-b1)x^2+...
1418EXAMPLE: example partial_molien; shows an example
1419{ poly A(2);                           // A(2) will contain the return value of
1420                                       // the intermediate result
1421  if (char(basering)<>0)
1422  { "ERROR:   you have to change to a basering of characteristic 0, one in";
1423    "         which the Molien series is defined";
1424  }
1425  if (ncols(M)==2 && nrows(M)==1 && n>0 && size(#)<2)
1426  { def br=basering;                   // keeping track of the old ring
1427    map slead=br,ideal(0);
1428    matrix s=slead(M);
1429    if (s[1,1]<>1 || s[1,2]<>1)
1430    { "ERROR:   the constant terms of enumerator and denominator are not 1";
1431      return();
1432    }
1433
1434    if (size(#)==0)
1435    { A(2)=M[1,1];                     // if a third parameter is not given, the
1436                                       // intermediate result from the last run
1437                                       // corresponds to the numerator - we need
1438    }                                  // its smallest term
1439    else
1440    { if (typeof(#[1])=="poly")
1441      { A(2)=#[1];                     // if a third term is given we 'start'
1442      }                                // with its smallest term
1443      else
1444      { "ERROR:   <poly> as third parameter expected";
1445        return();
1446      }
1447    }
1448    poly A(1)=M[1,2];                  // denominator of Molien series (for now)
1449    string mp=string(minpoly);
1450    execute "ring R=("+charstr(br)+"),("+varstr(br)+"),ds;";
1451    execute "minpoly=number("+mp+");";
1452    poly A(1)=0;                       // A(1) will contain the sum of n terms -
1453    poly min;                          // min will be our smallest term -
1454    poly A(2)=fetch(br,A(2));          // fetching A(2) and M[1,2] into R
1455    poly den=fetch(br,A(1));
1456    for (int i=1; i<=n; i=i+1)         // getting n terms and adding them up
1457    { min=lead(A(2));
1458      A(1)=A(1)+min;
1459      A(2)=A(2)-min*den;
1460    }
1461    setring br;                        // moving A(1) and A(2) back in the
1462    A(1)=fetch(R,A(1));                // actual ring for output
1463    A(2)=fetch(R,A(2));
1464    return(A(1..2));
1465  }
1466  else
1467  { "ERROR:   the first parameter has to be a 1x2-matrix, i.e. the matrix";
1468    "         returned by the procedure 'reynolds_molien', the second one";
1469    "         should be > 0 and there should be no more than 3 parameters;"
1470    return();
1471  }
1472}
1473example
1474{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
1475  echo=2;
1476         ring R=0,(x,y,z),dp;
1477         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
1478         matrix REY,M=reynolds_molien(A);
1479         poly p(1..2);
1480         p(1..2)=partial_molien(M,5);
1481         p(1);
1482         p(1..2)=partial_molien(M,5,p(2));
1483         p(1);
1484}
1485
1486proc evaluate_reynolds (matrix REY, ideal I)
1487USAGE:   evaluate_reynolds(REY,I);
1488         REY: a <matrix> representing the Reynolds operator, I: an arbitrary
1489         <ideal>
1490ASSUME:  REY is the first return value of group_reynolds() or reynolds_molien()
1491RETURNS: image of the polynomials defining I under the Reynolds operator
1492         (type <ideal>)
1493NOTE:    the characteristic of the coefficient field of the polynomial ring
1494         should not divide the order of the finite matrix group
1495EXAMPLE: example evaluate_reynolds; shows an example
1496THEORY:  REY has been constructed in such a way that each row serves as a ring
1497         mapping of which the Reynolds operator is made up.
1498{ def br=basering;
1499  int n=nvars(br);
1500  if (ncols(REY)==n)
1501  { int m=nrows(REY);                  // we need m to 'cut' the ring
1502                                       // homomorphisms 'out' of REY and to
1503                                       // divide by the group order in the end
1504    int num_poly=ncols(I);
1505    matrix MI=matrix(I);
1506    matrix MiI[1][num_poly];
1507    map pREY;
1508    matrix rowREY[1][n];
1509    for (int i=1;i<=m;i=i+1)
1510    { rowREY=REY[i,1..n];
1511      pREY=br,ideal(rowREY);           // f is now the i-th ring homomorphism
1512      MiI=pREY(MI)+MiI;
1513    }
1514    MiI=(1/number(m))*MiI;
1515    return(ideal(MiI));
1516  }
1517  else
1518  { "ERROR:   the number of columns in the <matrix> should be the same as the";
1519    "         number of variables in the basering; in fact it should be first";
1520    "         return value of group_reynolds() or reynolds_molien().";
1521    return();
1522  }
1523}
1524example
1525{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
1526  echo=2;
1527         ring R=0,(x,y,z),dp;
1528         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
1529         list L=group_reynolds(A);
1530         ideal I=x2,y2,z2;
1531         print(evaluate_reynolds(L[1],I));
1532}
1533
1534proc invariant_basis (int g,list #)
1535USAGE:   invariant_basis(g,G1,G2,...);
1536         g: an <int> indicating of which degree (>0) the homogeneous basis
1537         shoud be, G1,G2,...: <matrices> generating a finite matrix group
1538RETURNS: the basis (type <ideal>) of the space of invariants of degree g
1539EXAMPLE: example invariant_basis; shows an example
1540THEORY:  A general polynomial of degree g is generated and the generators of the
1541         matrix group applied. The difference ought to be 0 and this way a
1542         system of linear equations is created. It is solved by computing
1543         syzygies.
1544{ if (g<=0)
1545  { "ERROR:   the first parameter should be > 0";
1546    return();
1547  }
1548  def br=basering;
1549  ideal mon=sort(maxideal(g))[1];      // needed for constructing a general
1550  int m=ncols(mon);                    // homogeneous polynomial of degree g
1551  mon=sort(mon,intvec(m..1))[1];
1552  int a=size(#);
1553  int i;
1554  int n=nvars(br);
1555 //---------------------- checking that the input is ok -----------------------
1556  for (i=1;i<=a;i=i+1)
1557  { if (typeof(#[i])=="matrix")
1558    { if (nrows(#[i])==n && ncols(#[i])==n)
1559      { matrix G(i)=#[i];
1560      }
1561      else
1562      { "ERROR:   the number of variables of the base ring needs to be the same";
1563        "         as the dimension of the square matrices";
1564        return();
1565      }
1566    }
1567    else
1568    { "ERROR:   the last parameters should be a list of matrices";
1569      return();
1570    }
1571  }
1572 //----------------------------------------------------------------------------
1573  execute "ring T=("+charstr(br)+"),("+varstr(br)+",p(1..m)),lp;";
1574  // p(1..m) are the general coefficients of the general polynomial of degree g
1575  execute "ideal vars="+varstr(br)+";";
1576  map f;
1577  ideal mon=imap(br,mon);
1578  poly P=0;
1579  for (i=m;i>=1;i=i-1)
1580  { P=P+p(i)*mon[i];                   // P is the general polynomial
1581  }
1582  ideal I;                             // will help substituting variables in P
1583                                       // by linear combinations of variables -
1584  poly Pnew,temp;                      // Pnew is P with substitutions -
1585  matrix S[m*a][m];                    // will contain system of linear
1586                                       // equations
1587  int j,k;
1588 //------------------- building the system of linear equations ----------------
1589  for (i=1;i<=a;i=i+1)
1590  { I=ideal(matrix(vars)*transpose(imap(br,G(i))));
1591    I=I,p(1..m);
1592    f=T,I;
1593    Pnew=f(P);
1594    for (j=1;j<=m;j=j+1)
1595    { temp=P/mon[j]-Pnew/mon[j];
1596      for (k=1;k<=m;k=k+1)
1597      { S[m*(i-1)+j,k]=temp/p(k);
1598      }
1599    }
1600  }
1601 //----------------------------------------------------------------------------
1602  setring br;
1603  map f=T,ideal(0);
1604  matrix S=f(S);
1605  matrix s=matrix(syz(S));             // s contains a basis of the space of
1606                                       // solutions -
1607  ideal I=ideal(matrix(mon)*s);        // I contains a basis of homogeneous
1608  if (I[1]<>0)                         // invariants of degree d
1609  { for (i=1;i<=ncols(I);i=i+1)
1610    { I[i]=I[i]/leadcoef(I[i]);        // setting leading coefficients to 1
1611    }
1612  }
1613  return(I);
1614}
1615example
1616{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
1617  echo=2;
1618           ring R=0,(x,y,z),dp;
1619           matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
1620           print(invariant_basis(2,A));
1621}
1622
1623proc invariant_basis_reynolds (matrix REY,int d,list #)
1624USAGE:   invariant_basis_reynolds(REY,d[,flags]);
1625         REY: a <matrix> representing the Reynolds operator, d: an <int>
1626         indicating of which degree (>0) the homogeneous basis shoud be, flags:
1627         an optional <intvec> with two entries: its first component gives the
1628         dimension of the space (default <0 meaning unknown) and its second
1629         component is used as the number of polynomials that should be mapped
1630         to invariants during one call of evaluate_reynolds if the dimension of
1631         the space is unknown or the number such that number x dimension
1632         polynomials are mapped to invariants during one call of
1633         evaluate_reynolds
1634ASSUME:  REY is the first return value of group_reynolds() or reynolds_molien()
1635         and flags[1] given by partial_molien
1636RETURN:  the basis (type <ideal>) of the space of invariants of degree d
1637EXAMPLE: example invariant_basis_reynolds; shows an example
1638THEORY:  Monomials of degree d are mapped to invariants with the Reynolds
1639         operator. A linearly independent set is generated with the help of
1640         minbase.
1641{
1642 //---------------------- checking that the input is ok -----------------------
1643  if (d<=0)
1644  { "  ERROR:   the second parameter should be > 0";
1645     return();
1646  }
1647  if (size(#)>1)
1648  { "  ERROR:   there should be at most three parameters";
1649    return();
1650  }
1651  if (size(#)==1)
1652  { if (typeof(#[1])<>"intvec")
1653    { "  ERROR: the third parameter should be of type <intvec>";
1654      return();
1655    }
1656    if (size(#[1])<>2)
1657    { "  ERROR: there should be two components in <intvec>";
1658      return();
1659    }
1660    else
1661    { int cd=#[1][1];
1662      int step_fac=#[1][2];
1663    }
1664    if (step_fac<=0)
1665    { "  ERROR: the second component of <intvec> should be > 0";
1666      return();
1667    }
1668    if (cd==0)
1669    { return(ideal(0));
1670    }
1671  }
1672  else
1673  { int step_fac=1;
1674    int cd=-1;
1675  }
1676  if (ncols(REY)<>nvars(basering))
1677  { "ERROR:   the number of columns in the <matrix> should be the same as the";
1678    "         number of variables in the basering; in fact it should be first";
1679    "         return value of group_reynolds() or reynolds_molien().";
1680    return();
1681  }
1682 //----------------------------------------------------------------------------
1683  ideal mon=sort(maxideal(d))[1];
1684  degBound=d;
1685  int j=ncols(mon);
1686  mon=sort(mon,intvec(j..1))[1];
1687  ideal B;                             // will contain the basis
1688  if (cd<0)
1689  { if (step_fac>j)                    // all of mon will be mapped to
1690    { B=evaluate_reynolds(REY,mon);    // invariants at once
1691      B=minbase(B);
1692      degBound=0;
1693      return(B);
1694    }
1695  }
1696  else
1697  { if (step_fac*cd>j)                 // all of mon will be mapped to
1698    { B=evaluate_reynolds(REY,mon);    // invariants at once
1699      B=minbase(B);
1700      degBound=0;
1701      return(B);
1702    }
1703  }
1704  int i,k;
1705  int upper_bound=0;
1706  int lower_bound=0;
1707  ideal part_mon;                      // a part of mon of size step_fac*cd
1708  while (1)
1709  { lower_bound=upper_bound+1;
1710    if (cd<0)
1711    { upper_bound=upper_bound+step_fac;
1712    }
1713    else
1714    { upper_bound=upper_bound+step_fac*cd;
1715    }
1716    if (upper_bound>j)
1717    { upper_bound=j;
1718    }
1719    part_mon=mon[lower_bound..upper_bound];
1720    B=minbase(B+evaluate_reynolds(REY,part_mon));
1721    if (ncols(B)==cd or upper_bound==j)
1722    { degBound=0;
1723      return(B);
1724    }
1725  }
1726}
1727example
1728{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
1729  echo=2;
1730           ring R=0,(x,y,z),dp;
1731           matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
1732           intvec flags=0,1,0;
1733           matrix REY,M=reynolds_molien(A,flags);
1734           flags=8,6;
1735           print(invariant_basis_reynolds(REY,6,flags));
1736}
1737
1738////////////////////////////////////////////////////////////////////////////////
1739// This procedure generates linearly independent invariant polynomials of degree
1740// d that do not reduce to 0 modulo the primary invariants. It does this by
1741// applying the Reynolds operator to the monomials returned by kbase(sP,d). The
1742// result is used when computing secondary invariants.
1743////////////////////////////////////////////////////////////////////////////////
1744proc sort_of_invariant_basis (ideal sP,matrix REY,int d,int step_fac)
1745{ ideal mon=kbase(sP,d);
1746  degBound=d;
1747  int j=ncols(mon);
1748  int i;
1749  mon=sort(mon,intvec(j..1))[1];
1750  ideal B;                             // will contain the "sort of basis"
1751  if (step_fac>j)
1752  { B=compress(evaluate_reynolds(REY,mon));
1753    for (i=1;i<=ncols(B);i=i+1)        // those are taken our that are o mod sP
1754    { if (reduce(B[i],sP)==0)
1755      { B[i]=0;
1756      }
1757    }
1758    B=minbase(B);                      // here are the linearly independent ones
1759    degBound=0;
1760    return(B);
1761  }
1762  int upper_bound=0;
1763  int lower_bound=0;
1764  ideal part_mon;                      // parts of mon
1765  while (1)
1766  { lower_bound=upper_bound+1;
1767    upper_bound=upper_bound+step_fac;
1768    if (upper_bound>j)
1769    { upper_bound=j;
1770    }
1771    part_mon=mon[lower_bound..upper_bound];
1772    part_mon=compress(evaluate_reynolds(REY,part_mon));
1773    for (i=1;i<=ncols(part_mon);i=i+1)
1774    { if (reduce(part_mon[i],sP)==0)
1775      { part_mon[i]=0;
1776      }
1777    }
1778    B=minbase(B+part_mon);             // here are the linearly independent ones
1779    if (upper_bound==j)
1780    { degBound=0;
1781      return(B);
1782    }
1783  }
1784}
1785
1786////////////////////////////////////////////////////////////////////////////////
1787// Procedure returning the succeeding vector after vec. It is used to list
1788// all the vectors of Z^n with first nonzero entry 1. They are listed by
1789// increasing sum of the absolute value of their entries.
1790////////////////////////////////////////////////////////////////////////////////
1791proc next_vector(intmat vec)
1792{ int n=ncols(vec);                    // p: >0, n: <0, p0: >=0, n0: <=0
1793  for (int i=1;i<=n;i=i+1)             // finding out which is the first
1794  { if (vec[1,i]<>0)                   // component <>0
1795    { break;
1796    }
1797  }
1798  intmat new[1][n];
1799  if (i>n)                             // 0,...,0 --> 1,0....,0
1800  { new[1,1]=1;
1801    return(new);
1802  }
1803  if (i==n)                            // 0,...,1 --> 1,1,0,...,0
1804  { new[1,1..2]=1,1;
1805    return(new);
1806  }
1807  if (i==n-1)
1808  { if (vec[1,n]==0)                   // 0,...,0,1,0 --> 0,...,0,1
1809    { new[1,n]=1;
1810      return(new);
1811    }
1812    if (vec[1,n]>0)                    // 0,..,0,1,p --> 0,...,0,1,-p
1813    { new[1,1..n]=vec[1,1..n-1],-vec[1,n];
1814      return(new);
1815    }
1816    new[1,1..2]=1,1-vec[1,n];          // 0,..,0,1,n --> 1,1-n,0,..,0
1817    return(new);
1818  }
1819  if (i>1)
1820  { intmat temp[1][n-i+1]=vec[1,i..n]; // 0,...,0,1,*,...,* --> 1,*,...,*
1821    temp=next_vector(temp);
1822    new[1,i..n]=temp[1,1..n-i+1];
1823    return(new);
1824  }                                    // case left: 1,*,...,*
1825  for (i=2;i<=n;i=i+1)
1826  { if (vec[1,i]>0)                    // make first positive negative and
1827    { vec[1,i]=-vec[1,i];              // return
1828      return(vec);
1829    }
1830    else
1831    { vec[1,i]=-vec[1,i];              // make all negatives before positives
1832    }                                  // positive
1833  }
1834  for (i=2;i<=n-1;i=i+1)               // case: 1,p,...,p after 1,n,...,n
1835  { if (vec[1,i]>0)
1836    { vec[1,2]=vec[1,i]-1;             // shuffleing things around...
1837      if (i>2)                         // same sum of absolute values of entries
1838      { vec[1,i]=0;
1839      }
1840      vec[1,i+1]=vec[1,i+1]+1;
1841      return(vec);
1842    }
1843  }                                    // case left: 1,0,...,0 --> 1,1,0,...,0
1844  new[1,2..3]=1,vec[1,n];              // and: 1,0,...,0,1 --> 0,1,1,0,...,0
1845  return(new);
1846}
1847
1848////////////////////////////////////////////////////////////////////////////////
1849// Maps integers to elements of the base field. It is only called if the base
1850// field is of prime characteristic. If the base field has q elements (depending
1851// on minpoly) 1..q is mapped to those q elements.
1852////////////////////////////////////////////////////////////////////////////////
1853proc int_number_map (int i)
1854{ int p=char(basering);
1855  if (minpoly==0)                      // if no minpoly is given, we have p
1856  { i=i%p;                             // elements in the field
1857    return(number(i));
1858  }
1859  int d=pardeg(minpoly);
1860  if (i<0)
1861  { int bool=1;
1862    i=(-1)*i;
1863  }
1864  i=i%p^d;                             // base field has p^d elements -
1865  number a=par(1);                     // a is the root of the minpoly - we have
1866  number out=0;                        // to construct a linear combination of
1867  int j=1;                             // a^k
1868  int k;
1869  while (1)
1870  { if (i<p^j)                         // finding an upper bound on i
1871    { for (k=0;k<j-1;k=k+1)
1872      { out=out+((i/p^k)%p)*a^k;       // finding how often p^k is contained in
1873      }                                // i
1874      out=out+(i/p^(j-1))*a^(j-1);
1875      if (defined(bool)==voice)
1876      { return((-1)*out);
1877      }
1878      return(out);
1879    }
1880    j=j+1;
1881  }
1882}
1883
1884////////////////////////////////////////////////////////////////////////////////
1885// This procedure finds dif primary invariants in degree d. It returns all
1886// primary invariants found so far. The coefficients lie in a field of
1887// characteristic 0.
1888////////////////////////////////////////////////////////////////////////////////
1889proc search (int n,int d,ideal B,int cd,ideal P,ideal sP,int i,int dif,int dB,ideal CI)
1890{ intmat vec[1][cd];                   // the coefficients for the next
1891                                       // combination -
1892  degBound=0;
1893  poly test_poly;                      // the linear combination to test
1894  int test_dim;
1895  intvec h;                            // Hilbert series
1896  int j=i+1;
1897  matrix tB=transpose(B);
1898  ideal TEST;
1899  while(j<=i+dif)
1900  { CI=CI+ideal(var(j)^d);             // homogeneous polynomial of the same
1901                                       // degree as the one we're looking for is
1902                                       // added
1903    // h=hilb(std(CI),1);
1904    dB=dB+d-1;                         // used as degBound
1905    while(1)
1906    { vec=next_vector(vec);            // next vector
1907      test_poly=(vec*tB)[1,1];
1908      // degBound=dB;
1909      TEST=sP+ideal(test_poly);
1910      attrib(TEST,"isSB",1);
1911      test_dim=dim(TEST);
1912      // degBound=0;
1913      if (n-test_dim==j)               // the dimension has been lowered by one
1914      { sP=TEST;
1915        break;
1916      }
1917      // degBound=dB;
1918      TEST=std(sP+ideal(test_poly));   // should soon be replaced by next line
1919      // TEST=std(sP,test_poly,h);        // Hilbert driven std-calculation
1920      test_dim=dim(TEST);
1921      // degBound=0;
1922      if (n-test_dim==j)               // the dimension has been lowered by one
1923      { sP=TEST;
1924        break;
1925      }
1926    }
1927    P[j]=test_poly;                    // test_poly ist added to primary
1928    j=j+1;                             // invariants
1929  }
1930  return(P,sP,CI,dB);
1931}
1932
1933////////////////////////////////////////////////////////////////////////////////
1934// This procedure finds at most dif primary invariants in degree d. It returns
1935// all primary invariants found so far. The coefficients lie in the field of
1936// characteristic p>0.
1937////////////////////////////////////////////////////////////////////////////////
1938proc 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
2022proc primary_char0 (matrix REY,matrix M,list #)
2023USAGE:   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>
2026ASSUME:  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
2028DISPLAY: information about the various stages of the programme if v does not
2029         equal 0
2030RETURN:  primary invariants (type <matrix>) of the invariant ring
2031EXAMPLE: example primary_char0; shows an example
2032THEORY:  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}
2152example
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
2162proc primary_charp (matrix REY,string ring_name,list #)
2163USAGE:   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>
2167ASSUME:  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
2170DISPLAY: information about the various stages of the programme if v does not
2171         equal 0
2172RETURN:  primary invariants (type <matrix>) of the invariant ring
2173EXAMPLE: example primary_charp; shows an example
2174THEORY:  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}
2300example
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
2314proc primary_char0_no_molien (matrix REY, list #)
2315USAGE:   primary_char0_no_molien(REY[,v]);
2316         REY: a <matrix> representing the Reynolds operator, v: an optional
2317         <int>
2318ASSUME:  REY is the first return value of group_reynolds or reynolds_molien
2319DISPLAY: information about the various stages of the programme if v does not
2320         equal 0
2321RETURN:  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
2324EXAMPLE: example primary_char0_no_molien; shows an example
2325THEORY:  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}
2446example
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
2456proc primary_charp_no_molien (matrix REY, list #)
2457USAGE:   primary_charp_no_molien(REY[,v]);
2458         REY: a <matrix> representing the Reynolds operator, v: an optional
2459         <int>
2460ASSUME:  REY is the first return value of group_reynolds or reynolds_molien
2461DISPLAY: information about the various stages of the programme if v does not
2462         equal 0
2463RETURN:  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
2466EXAMPLE: example primary_charp_no_molien; shows an example
2467THEORY:  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}
2589example
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
2600proc primary_charp_without (list #)
2601USAGE:   primary_charp_without(G1,G2,...[,v]);
2602         G1,G2,...: <matrices> generating a finite matrix group, v: an optional
2603         <int>
2604DISPLAY: information about the various stages of the programme if v does not
2605         equal 0
2606RETURN:  primary invariants (type <matrix>) of the invariant ring
2607EXAMPLE: example primary_charp_without; shows an example
2608THEORY:  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}
2730example
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
2738proc primary_invariants (list #)
2739USAGE:   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)
2755DISPLAY: information about the various stages of the programme if the third
2756         flag does not equal 0
2757RETURN:  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
2762EXAMPLE: example primary_invariants; shows an example
2763THEORY:  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}
2926example
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]);
2933}
2934
2935////////////////////////////////////////////////////////////////////////////////
2936// This procedure finds dif primary invariants in degree d. It returns all
2937// primary invariants found so far. The coefficients lie in a field of
2938// characteristic 0.
2939////////////////////////////////////////////////////////////////////////////////
2940proc 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////////////////////////////////////////////////////////////////////////////////
3036proc 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
3146proc primary_char0_random (matrix REY,matrix M,int max,list #)
3147USAGE:   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>
3152ASSUME:  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
3154DISPLAY: information about the various stages of the programme if v does not
3155         equal 0
3156RETURN:  primary invariants (type <matrix>) of the invariant ring
3157EXAMPLE: example primary_char0_random; shows an example
3158THEORY:  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}
3282example
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
3292proc primary_charp_random (matrix REY,string ring_name,int max,list #)
3293USAGE:   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>
3298ASSUME:  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
3301DISPLAY: information about the various stages of the programme if v does not
3302         equal 0
3303RETURN:  primary invariants (type <matrix>) of the invariant ring
3304EXAMPLE: example primary_charp_random; shows an example
3305THEORY:  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}
3435example
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
3449proc primary_char0_no_molien_random (matrix REY, int max, list #)
3450USAGE:   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>
3454ASSUME:  REY is the first return value of group_reynolds or reynolds_molien
3455DISPLAY: information about the various stages of the programme if v does not
3456         equal 0
3457RETURN:  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
3460EXAMPLE: example primary_char0_no_molien_random; shows an example
3461THEORY:  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}
3586example
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
3596proc primary_charp_no_molien_random (matrix REY, int max, list #)
3597USAGE:   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>
3601ASSUME:  REY is the first return value of group_reynolds or reynolds_molien
3602DISPLAY: information about the various stages of the programme if v does not
3603         equal 0
3604RETURN:  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
3607EXAMPLE: example primary_charp_no_molien_random; shows an example
3608THEORY:  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}
3733example
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
3744proc primary_charp_without_random (list #)
3745USAGE:   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>
3749DISPLAY: information about the various stages of the programme if v does not
3750         equal 0
3751RETURN:  primary invariants (type <matrix>) of the invariant ring
3752EXAMPLE: example primary_charp_without_random; shows an example
3753THEORY:  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}
3885example
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
3893proc primary_invariants_random (list #)
3894USAGE:   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)
3912DISPLAY: information about the various stages of the programme if the third
3913         flag does not equal 0
3914RETURN:  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
3919EXAMPLE: example primary_invariants_random; shows an example
3920THEORY:  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}
4088example
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
4097proc 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
4107proc power_products(intvec deg_vec,int d)
4108USAGE:   power_products(dv,d);
4109         dv: an <intvec> giving the degrees of homogeneous polynomials, d: the
4110         degree of the desired power products
4111RETURN:  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.
4114EXAMPLE: 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}
4155example
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
4162proc secondary_char0 (matrix P, matrix REY, matrix M, list #)
4163USAGE:   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>
4167ASSUME:  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()
4172RETURN:  secondary invariants of the invariant ring (type <matrix>) and
4173         irreducible secondary invariants (type <matrix>)
4174DISPLAY: information if v does not equal 0
4175EXAMPLE: example secondary_char0; shows an example
4176THEORY:  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)).
4182{ def br=basering;
4183  degBound=0;
4184 //----------------- checking input and setting verbose mode ------------------
4185  if (char(br)<>0)
4186  { "ERROR:   secondary_char0 should only be used with rings of characteristic 0.";
4187    return();
4188  }
4189  int i;
4190  if (size(#)>0)
4191  { if (typeof(#[size(#)])=="int")
4192    { int v=#[size(#)];
4193    }
4194    else
4195    { int v=0;
4196    }
4197  }
4198  else
4199  { int v=0;
4200  }
4201  int n=nvars(br);                     // n is the number of variables, as well
4202                                       // as the size of the matrices, as well
4203                                       // as the number of primary invariants,
4204                                       // we should get
4205  if (ncols(P)<>n)
4206  { "ERROR:   The first parameter ought to be the matrix of the primary";
4207    "         invariants."
4208    return();
4209  }
4210  if (ncols(REY)<>n)
4211  { "ERROR:   The second parameter ought to be the Reynolds operator."
4212    return();
4213  }
4214  if (ncols(M)<>2 or nrows(M)<>1)
4215  { "ERROR:   The third parameter ought to be the Molien series."
4216    return();
4217  }
4218  if (v && voice==2)
4219  { "";
4220  }
4221  int j, m, counter;
4222 //- finding the polynomial giving number and degrees of secondary invariants -
4223  poly p=1;
4224  for (j=1;j<=n;j=j+1)                 // calculating the denominator of the
4225  { p=p*(1-var(1)^deg(P[j]));          // Hilbert series of the ring generated
4226  }                                    // by the primary invariants -
4227  matrix s[1][2]=M[1,1]*p,M[1,2];      // s is used for canceling
4228  s=matrix(syz(ideal(s)));
4229  p=s[2,1];                            // the polynomial telling us where to
4230                                       // search for secondary invariants
4231  map slead=br,ideal(0);
4232  p=1/slead(p)*p;                      // smallest term of p needs to be 1
4233  if (v)
4234  { "  Polynomial telling us where to look for secondary invariants:";
4235    "   "+string(p);
4236    "";
4237  }
4238  matrix dimmat=coeffs(p,var(1));      // dimmat will contain the number of
4239                                       // secondary invariants, we need to find
4240                                       // of a certain degree -
4241  m=nrows(dimmat);                     // m-1 is the highest degree
4242  if (v)
4243  { "  In degree 0 we have: 1";
4244    "";
4245  }
4246 //-------------------------- initializing variables --------------------------
4247  intmat PP;
4248  poly pp;
4249  int k;
4250  intvec deg_vec;
4251  ideal sP=std(ideal(P));
4252  ideal TEST,B,IS;
4253  ideal S=1;                           // 1 is the first secondary invariant -
4254 //--------------------- generating secondary invariants ----------------------
4255  for (i=2;i<=m;i=i+1)                 // going through dimmat -
4256  { if (int(dimmat[i,1])<>0)           // when it is == 0 we need to find 0
4257    {                                  // elements in the current degree (i-1)
4258      if (v)
4259      { "  Searching in degree "+string(i-1)+", we need to find "+string(int(dimmat[i,1]))+" invariant(s)...";
4260      }
4261      TEST=sP;
4262      counter=0;                       // we'll count up to degvec[i]
4263      if (IS[1]<>0)
4264      { PP=power_products(deg_vec,i-1); // finding power products of irreducible
4265      }                                // secondary invariants
4266      if (size(ideal(PP))<>0)
4267      { for (j=1;j<=ncols(PP);j=j+1)   // going through all the power products
4268        { pp=1;
4269          for (k=1;k<=nrows(PP);k=k+1)
4270          { pp=pp*IS[1,k]^PP[k,j];
4271          }
4272          if (reduce(pp,TEST)<>0)
4273          { S=S,pp;
4274            counter=counter+1;
4275            if (v)
4276            { "  We find: "+string(pp);
4277            }
4278            if (int(dimmat[i,1])<>counter)
4279            { TEST=std(TEST+ideal(NF(pp,TEST))); // should be replaced by next
4280                                                 // line soon
4281            // TEST=std(TEST,NF(pp,TEST));
4282            }
4283            else
4284            { break;
4285            }
4286          }
4287        }
4288      }
4289      if (int(dimmat[i,1])<>counter)
4290      { B=sort_of_invariant_basis(sP,REY,i-1,int(dimmat[i,1])*6); // B contains
4291                                       // images of kbase(sP,i-1) under the
4292                                       // Reynolds operator that are linearly
4293                                       // independent and that don't reduce to
4294                                       // 0 modulo sP -
4295        if (counter==0 && ncols(B)==int(dimmat[i,1])) // then we can take all of
4296        { S=S,B;                       // B
4297          IS=IS+B;
4298          if (deg_vec[1]==0)
4299          { deg_vec=i-1;
4300            if (v)
4301            { "  We find: "+string(B[1]);
4302            }
4303            for (j=2;j<=int(dimmat[i,1]);j=j+1)
4304            { deg_vec=deg_vec,i-1;
4305              if (v)
4306              { "  We find: "+string(B[j]);
4307              }
4308            }
4309          }
4310          else
4311          { for (j=1;j<=int(dimmat[i,1]);j=j+1)
4312            { deg_vec=deg_vec,i-1;
4313              if (v)
4314              { "  We find: "+string(B[j]);
4315              }
4316            }
4317          }
4318        }
4319        else
4320        { j=0;                         // j goes through all of B -
4321          while (int(dimmat[i,1])<>counter) // need to find dimmat[i,1]
4322          {                            // invariants that are linearly
4323                                       // independent modulo TEST
4324            j=j+1;
4325            if (reduce(B[j],TEST)<>0)  // B[j] should be added
4326            { S=S,B[j];
4327              IS=IS+ideal(B[j]);
4328              if (deg_vec[1]==0)
4329              { deg_vec[1]=i-1;
4330              }
4331              else
4332              { deg_vec=deg_vec,i-1;
4333              }
4334              counter=counter+1;
4335              if (v)
4336              { "  We find: "+string(B[j]);
4337              }
4338              if (int(dimmat[i,1])<>counter)
4339              { TEST=std(TEST+ideal(NF(B[j],TEST))); // should be replaced by
4340                                                     // next line
4341              // TEST=std(TEST,NF(B[j],TEST));
4342              }
4343            }
4344          }
4345        }
4346      }
4347      if (v)
4348      { "";
4349      }
4350    }
4351  }
4352  if (v)
4353  { "  We're done!";
4354    "";
4355  }
4356  return(matrix(S),matrix(IS));
4357}
4358example
4359{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
4360  echo=2;
4361         ring R=0,(x,y,z),dp;
4362         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
4363         list L=primary_invariants(A);
4364         matrix S,IS=secondary_char0(L[1..3]);
4365         print(S);
4366         print(IS);
4367}
4368
4369proc secondary_charp (matrix P, matrix REY, string ring_name, list #)
4370USAGE:   secondary_charp(P,REY,ringname[,v]);
4371         P: a 1xn <matrix> with primary invariants, REY: a gxn <matrix>
4372         representing the Reynolds operator, ringname: a <string> giving the
4373         name of a ring of characteristic 0 where the Molien series is stored,
4374         v: an optional <int>
4375ASSUME:  n is the number of variables of the basering, g the size of the group,
4376         REY is the 1st return value of group_reynolds(), reynolds_molien() or
4377         the second one of primary_invariants(), `ringname` is ring of
4378         characteristic 0 that has been created by molien() or reynolds_molien()
4379         or primary_invariants()
4380RETURN:  secondary invariants of the invariant ring (type <matrix>) and
4381         irreducible secondary invariants (type <matrix>)
4382DISPLAY: information if v does not equal 0
4383EXAMPLE: example secondary_charp; shows an example
4384THEORY:  The secondary invariants are calculated by finding a basis (in terms of
4385         monomials) of the basering modulo the primary invariants, mapping those
4386         to invariants with the Reynolds operator and using these images or
4387         their power products such that they are linearly independent modulo the
4388         primary invariants (see paper "Some Algorithms in Invariant Theory of
4389         Finite Groups" by Kemper and Steel (1997)).
4390{ def br=basering;
4391  degBound=0;
4392 //---------------- checking input and setting verbose mode -------------------
4393  if (char(br)==0)
4394  { "ERROR:   secondary_charp should only be used with rings of characteristic p>0.";
4395    return();
4396  }
4397  int i;
4398  if (size(#)>0)
4399  { if (typeof(#[size(#)])=="int")
4400    { int v=#[size(#)];
4401    }
4402    else
4403    { int v=0;
4404    }
4405  }
4406  else
4407  { int v=0;
4408  }
4409  int n=nvars(br);                     // n is the number of variables, as well
4410                                       // as the size of the matrices, as well
4411                                       // as the number of primary invariants,
4412                                       // we should get
4413  if (ncols(P)<>n)
4414  { "ERROR:   The first parameter ought to be the matrix of the primary";
4415    "         invariants."
4416    return();
4417  }
4418  if (ncols(REY)<>n)
4419  { "ERROR:   The second parameter ought to be the Reynolds operator."
4420    return();
4421  }
4422  if (typeof(`ring_name`)<>"ring")
4423  { "ERROR:   The <string> should give the name of the ring where the Molien."
4424    "         series is stored.";
4425    return();
4426  }
4427  if (v && voice==2)
4428  { "";
4429  }
4430  int j, m, counter, d;
4431  intvec deg_dim_vec;
4432 //- finding the polynomial giving number and degrees of secondary invariants -
4433  for (j=1;j<=n;j=j+1)
4434  { deg_dim_vec[j]=deg(P[j]);
4435  }
4436  setring `ring_name`;
4437  poly p=1;
4438  for (j=1;j<=n;j=j+1)                 // calculating the denominator of the
4439  { p=p*(1-var(1)^deg_dim_vec[j]);     // Hilbert series of the ring generated
4440  }                                    // by the primary invariants -
4441  matrix s[1][2]=M[1,1]*p,M[1,2];      // s is used for canceling
4442  s=matrix(syz(ideal(s)));
4443  p=s[2,1];                            // the polynomial telling us where to
4444                                       // search for secondary invariants
4445  map slead=basering,ideal(0);
4446  p=1/slead(p)*p;                      // smallest term of p needs to be 1
4447  if (v)
4448  { "  Polynomial telling us where to look for secondary invariants:";
4449    "   "+string(p);
4450    "";
4451  }
4452  matrix dimmat=coeffs(p,var(1));      // dimmat will contain the number of
4453                                       // secondary invariants, we need to find
4454                                       // of a certain degree -
4455  m=nrows(dimmat);                     // m-1 is the highest degree
4456  deg_dim_vec=1;
4457  for (j=2;j<=m;j=j+1)
4458  { deg_dim_vec=deg_dim_vec,int(dimmat[j,1]);
4459  }
4460  if (v)
4461  { "  In degree 0 we have: 1";
4462    "";
4463  }
4464 //------------------------ initializing variables ----------------------------
4465  setring br;
4466  intmat PP;
4467  poly pp;
4468  int k;
4469  intvec deg_vec;
4470  ideal sP=std(ideal(P));
4471  ideal TEST,B,IS;
4472  ideal S=1;                           // 1 is the first secondary invariant
4473 //------------------- generating secondary invariants ------------------------
4474  for (i=2;i<=m;i=i+1)                 // going through deg_dim_vec -
4475  { if (deg_dim_vec[i]<>0)             // when it is == 0 we need to find 0
4476    {                                  // elements in the current degree (i-1)
4477      if (v)
4478      { "  Searching in degree "+string(i-1)+", we need to find "+string(deg_dim_vec[i])+" invariant(s)...";
4479      }
4480      TEST=sP;
4481      counter=0;                       // we'll count up to degvec[i]
4482      if (IS[1]<>0)
4483      { PP=power_products(deg_vec,i-1); // generating power products of
4484      }                                // irreducible secondary invariants
4485      if (size(ideal(PP))<>0)
4486      { for (j=1;j<=ncols(PP);j=j+1)   // going through all of those
4487        { pp=1;
4488          for (k=1;k<=nrows(PP);k=k+1)
4489          { pp=pp*IS[1,k]^PP[k,j];
4490          }
4491          if (reduce(pp,TEST)<>0)
4492          { S=S,pp;
4493            counter=counter+1;
4494            if (v)
4495            { "  We find: "+string(pp);
4496            }
4497            if (deg_dim_vec[i]<>counter)
4498            { TEST=std(TEST+ideal(NF(pp,TEST))); // should be soon replaced by
4499                                                 // next line
4500            // TEST=std(TEST,NF(pp,TEST));
4501            }
4502            else
4503            { break;
4504            }
4505          }
4506        }
4507      }
4508      if (deg_dim_vec[i]<>counter)
4509      { B=sort_of_invariant_basis(sP,REY,i-1,deg_dim_vec[i]*6); // B contains
4510                                       // images of kbase(sP,i-1) under the
4511                                       // Reynolds operator that are linearly
4512                                       // independent and that don't reduce to
4513                                       // 0 modulo sP -
4514        if (counter==0 && ncols(B)==deg_dim_vec[i]) // then we can add all of B
4515        { S=S,B;
4516          IS=IS+B;
4517          if (deg_vec[1]==0)
4518          { deg_vec=i-1;
4519            if (v)
4520            { "  We find: "+string(B[1]);
4521            }
4522            for (j=2;j<=deg_dim_vec[i];j=j+1)
4523            { deg_vec=deg_vec,i-1;
4524              if (v)
4525              { "  We find: "+string(B[j]);
4526              }
4527            }
4528          }
4529          else
4530          { for (j=1;j<=deg_dim_vec[i];j=j+1)
4531            { deg_vec=deg_vec,i-1;
4532              if (v)
4533              { "  We find: "+string(B[j]);
4534              }
4535            }
4536          }
4537        }
4538        else
4539        { j=0;                         // j goes through all of B -
4540          while (deg_dim_vec[i]<>counter) // need to find deg_dim_vec[i]
4541          {                            // invariants that are linearly
4542                                       // independent modulo TEST
4543            j=j+1;
4544            if (reduce(B[j],TEST)<>0)   // B[j] should be added
4545            { S=S,B[j];
4546              IS=IS+ideal(B[j]);
4547              if (deg_vec[1]==0)
4548              { deg_vec[1]=i-1;
4549              }
4550              else
4551              { deg_vec=deg_vec,i-1;
4552              }
4553              counter=counter+1;
4554              if (v)
4555              { "  We find: "+string(B[j]);
4556              }
4557              if (deg_dim_vec[i]<>counter)
4558              { TEST=std(TEST+ideal(NF(B[j],TEST))); // should be soon replaced
4559                                                     // by next line
4560              // TEST=std(TEST,NF(B[j],TEST));
4561              }
4562            }
4563          }
4564        }
4565      }
4566      if (v)
4567      { "";
4568      }
4569    }
4570  }
4571  if (v)
4572  { "  We're done!";
4573    "";
4574  }
4575  if (ring_name=="aksldfalkdsflkj")
4576  { kill `ring_name`;
4577  }
4578  return(matrix(S),matrix(IS));
4579}
4580example
4581{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed into";
4582  "         characteristic 3)";
4583  echo=2;
4584         ring R=3,(x,y,z),dp;
4585         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
4586         list L=primary_invariants(A);
4587         matrix S,IS=secondary_charp(L[1..size(L)]);
4588         print(S);
4589         print(IS);
4590}
4591
4592proc secondary_no_molien (matrix P, matrix REY, list #)
4593USAGE:   secondary_no_molien(P,REY[,deg_vec,v]);
4594         P: a 1xn <matrix> with primary invariants, REY: a gxn <matrix>
4595         representing the Reynolds operator, deg_vec: an optional <intvec>
4596         listing some degrees where no non-trivial homogeneous invariants can be
4597         found, v: an optional <int>
4598ASSUME:  n is the number of variables of the basering, g the size of the group,
4599         REY is the 1st return value of group_reynolds(), reynolds_molien() or
4600         the second one of primary_invariants(), deg_vec is the second return
4601         value of primary_char0_no_molien(), primary_charp_no_molien(),
4602         primary_char0_no_molien_random() or primary_charp_no_molien_random()
4603RETURN:  secondary invariants of the invariant ring (type <matrix>)
4604DISPLAY: information if v does not equal 0
4605EXAMPLE: example secondary_no_molien; shows an example
4606THEORY:  The secondary invariants are calculated by finding a basis (in terms of
4607         monomials) of the basering modulo the primary invariants, mapping those
4608         to invariants with the Reynolds operator and using these images as
4609         candidates for secondary invariants.
4610{ int i;
4611  degBound=0;
4612 //------------------ checking input and setting verbose ----------------------
4613  if (size(#)==1 or size(#)==2)
4614  { if (typeof(#[size(#)])=="int")
4615    { if (size(#)==2)
4616      { if (typeof(#[size(#)-1])=="intvec")
4617        { intvec deg_vec=#[size(#)-1];
4618        }
4619        else
4620        { "ERROR:   the third parameter should be an <intvec>";
4621          return();
4622        }
4623      }
4624      int v=#[size(#)];
4625    }
4626    else
4627    { if (size(#)==1)
4628      { if (typeof(#[size(#)])=="intvec")
4629        { intvec deg_vec=#[size(#)];
4630          int v=0;
4631        }
4632        else
4633        { "ERROR:   the third parameter should be an <intvec>";
4634          return();
4635        }
4636      }
4637      else
4638      { "ERROR:   wrong list of parameters";
4639        return();
4640      }
4641    }
4642  }
4643  else
4644  { if (size(#)>2)
4645    { "ERROR:   there are too many parameters";
4646      return();
4647    }
4648    int v=0;
4649  }
4650  int n=nvars(basering);               // n is the number of variables, as well
4651                                       // as the size of the matrices, as well
4652                                       // as the number of primary invariants,
4653                                       // we should get
4654  if (ncols(P)<>n)
4655  { "ERROR:   The first parameter ought to be the matrix of the primary";
4656    "         invariants."
4657    return();
4658  }
4659  if (ncols(REY)<>n)
4660  { "ERROR:   The second parameter ought to be the Reynolds operator."
4661    return();
4662  }
4663  if (v && voice==2)
4664  { "";
4665  }
4666  int j, m, d;
4667  int max=1;
4668  for (j=1;j<=n;j=j+1)
4669  { max=max*deg(P[j]);
4670  }
4671  max=max/nrows(REY);
4672  if (v)
4673  { "  We need to find "+string(max)+" secondary invariants.";
4674    "";
4675    "  In degree 0 we have: 1";
4676    "";
4677  }
4678 //------------------------- initializing variables ---------------------------
4679  ideal sP=std(ideal(P));
4680  ideal B, TEST;
4681  ideal S=1;                           // 1 is the first secondary invariant
4682  int counter=1;
4683  i=0;
4684  if (defined(deg_vec)<>voice)
4685  { intvec deg_vec;
4686  }
4687  int k=1;
4688 //--------------------- generating secondary invariants ----------------------
4689  while (counter<>max)
4690  { i=i+1;
4691    if (deg_vec[k]<>i)
4692    { if (v)
4693      { "  Searching in degree "+string(i)+"...";
4694      }
4695      B=sort_of_invariant_basis(sP,REY,i,max); // B contains images of
4696                                       // kbase(sP,i) under the Reynolds
4697                                       // operator that are linearly independent
4698                                       // and that don't reduce to 0 modulo sP
4699      TEST=sP;
4700      for (j=1;j<=ncols(B);j=j+1)
4701      {                                // that are linearly independent modulo
4702                                       // TEST
4703        if (reduce(B[j],TEST)<>0)      // B[j] should be added
4704        { S=S,B[j];
4705          counter=counter+1;
4706          if (v)
4707          { "  We find: "+string(B[j]);
4708          }
4709          if (counter==max)
4710          { break;
4711          }
4712          else
4713          { if (j<>ncols(B))
4714            { TEST=std(TEST+ideal(NF(B[j],TEST))); // should soon be replaced by
4715                                                   // next line
4716            // TEST=std(TEST,NF(B[j],TEST));
4717            }
4718          }
4719        }
4720      }
4721    }
4722    else
4723    { if (size(deg_vec)==k)
4724      { k=1;
4725      }
4726      else
4727      { k=k+1;
4728      }
4729    }
4730  }
4731  if (v)
4732  { "";
4733  }
4734  if (v)
4735  { "  We're done!";
4736    "";
4737  }
4738  return(matrix(S));
4739}
4740example
4741{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
4742  echo=2;
4743         ring R=0,(x,y,z),dp;
4744         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
4745         list L=primary_invariants(A,intvec(1,1,0));
4746         matrix S=secondary_no_molien(L[1..3]);
4747         print(S);
4748}
4749
4750proc secondary_with_irreducible_ones_no_molien (matrix P, matrix REY, list #)
4751USAGE:   secondary_with_irreducible_ones_no_molien(P,REY[,v]);
4752         P: a 1xn <matrix> with primary invariants, REY: a gxn <matrix>
4753         representing the Reynolds operator, v: an optional <int>
4754ASSUME:  n is the number of variables of the basering, g the size of the group,
4755         REY is the 1st return value of group_reynolds(), reynolds_molien() or
4756         the second one of primary_invariants()
4757RETURN:  secondary invariants of the invariant ring (type <matrix>) and
4758         irreducible secondary invariants (type <matrix>)
4759DISPLAY: information if v does not equal 0
4760EXAMPLE: example secondary_with_irreducible_ones_no_molien; shows an example
4761THEORY:  The secondary invariants are calculated by finding a basis (in terms of
4762         monomials) of the basering modulo the primary invariants, mapping those
4763         to invariants with the Reynolds operator and using these images or
4764         their power products such that they are linearly independent modulo the
4765         primary invariants (see paper "Some Algorithms in Invariant Theory of
4766         Finite Groups" by Kemper and Steel (1997)).
4767{ int i;
4768  degBound=0;
4769 //--------------------- checking input and setting verbose mode --------------
4770  if (size(#)==1 or size(#)==2)
4771  { if (typeof(#[size(#)])=="int")
4772    { if (size(#)==2)
4773      { if (typeof(#[size(#)-1])=="intvec")
4774        { intvec deg_vec=#[size(#)-1];
4775        }
4776        else
4777        { "ERROR:   the third parameter should be an <intvec>";
4778          return();
4779        }
4780      }
4781      int v=#[size(#)];
4782    }
4783    else
4784    { if (size(#)==1)
4785      { if (typeof(#[size(#)])=="intvec")
4786        { intvec deg_vec=#[size(#)];
4787          int v=0;
4788        }
4789        else
4790        { "ERROR:   the third parameter should be an <intvec>";
4791          return();
4792        }
4793      }
4794      else
4795      { "ERROR:   wrong list of parameters";
4796        return();
4797      }
4798    }
4799  }
4800  else
4801  { if (size(#)>2)
4802    { "ERROR:   there are too many parameters";
4803      return();
4804    }
4805    int v=0;
4806  }
4807  int n=nvars(basering);               // n is the number of variables, as well
4808                                       // as the size of the matrices, as well
4809                                       // as the number of primary invariants,
4810                                       // we should get
4811  if (ncols(P)<>n)
4812  { "ERROR:   The first parameter ought to be the matrix of the primary";
4813    "         invariants."
4814    return();
4815  }
4816  if (ncols(REY)<>n)
4817  { "ERROR:   The second parameter ought to be the Reynolds operator."
4818    return();
4819  }
4820  if (v && voice==2)
4821  { "";
4822  }
4823  int j, m, d;
4824  int max=1;
4825  for (j=1;j<=n;j=j+1)
4826  { max=max*deg(P[j]);
4827  }
4828  max=max/nrows(REY);
4829  if (v)
4830  { "  We need to find "+string(max)+" secondary invariants.";
4831    "";
4832    "  In degree 0 we have: 1";
4833    "";
4834  }
4835 //------------------------ initializing variables ----------------------------
4836  intmat PP;
4837  poly pp;
4838  int k;
4839  intvec irreducible_deg_vec;
4840  ideal sP=std(ideal(P));
4841  ideal B,TEST,IS;
4842  ideal S=1;                           // 1 is the first secondary invariant
4843  int counter=1;
4844  i=0;
4845  if (defined(deg_vec)<>voice)
4846  { intvec deg_vec;
4847  }
4848  int l=1;
4849 //------------------- generating secondary invariants ------------------------
4850  while (counter<>max)
4851  { i=i+1;
4852    if (deg_vec[l]<>i)
4853    { if (v)
4854      { "  Searching in degree "+string(i)+"...";
4855      }
4856      TEST=sP;
4857      if (IS[1]<>0)
4858      { PP=power_products(irreducible_deg_vec,i);  // generating all power
4859      }                                // products of irreducible secondary
4860                                       // invariants
4861      if (size(ideal(PP))<>0)
4862      { for (j=1;j<=ncols(PP);j=j+1)   // going through all those power products
4863        { pp=1;
4864          for (k=1;k<=nrows(PP);k=k+1)
4865          { pp=pp*IS[1,k]^PP[k,j];
4866          }
4867          if (reduce(pp,TEST)<>0)
4868          { S=S,pp;
4869            counter=counter+1;
4870            if (v)
4871            { "  We find: "+string(pp);
4872            }
4873            if (counter<>max)
4874            { TEST=std(TEST+ideal(NF(pp,TEST))); // should soon be replaced by
4875                                                 // next line
4876            // TEST=std(TEST,NF(pp,TEST));
4877            }
4878            else
4879            { break;
4880            }
4881          }
4882        }
4883      }
4884      if (max<>counter)
4885      { B=sort_of_invariant_basis(sP,REY,i,max); // B contains images of
4886                                       // kbase(sP,i) under the Reynolds
4887                                       // operator that are linearly independent
4888                                       // and that don't reduce to 0 modulo sP
4889        for (j=1;j<=ncols(B);j=j+1)
4890        { if (reduce(B[j],TEST)<>0)    // B[j] should be added
4891          { S=S,B[j];
4892            IS=IS+ideal(B[j]);
4893            if (irreducible_deg_vec[1]==0)
4894            { irreducible_deg_vec[1]=i;
4895            }
4896            else
4897            { irreducible_deg_vec=irreducible_deg_vec,i;
4898            }
4899            counter=counter+1;
4900            if (v)
4901            { "  We find: "+string(B[j]);
4902            }
4903            if (counter==max)
4904            { break;
4905            }
4906            else
4907            { if (j<>ncols(B))
4908              { TEST=std(TEST+ideal(NF(B[j],TEST))); // should soon be replaced
4909                                                     // by next line
4910              // TEST=std(TEST,NF(B[j],TEST));
4911              }
4912            }
4913          }
4914        }
4915      }
4916    }
4917    else
4918    { if (size(deg_vec)==l)
4919      { l=1;
4920      }
4921      else
4922      { l=l+1;
4923      }
4924    }
4925  }
4926  if (v)
4927  { "";
4928  }
4929  if (v)
4930  { "  We're done!";
4931    "";
4932  }
4933  return(matrix(S),matrix(IS));
4934}
4935example
4936{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
4937  echo=2;
4938         ring R=0,(x,y,z),dp;
4939         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
4940         list L=primary_invariants(A,intvec(1,1,0));
4941         matrix S,IS=secondary_with_irreducible_ones_no_molien(L[1..2]);
4942         print(S);
4943         print(IS);
4944}
4945
4946proc secondary_not_cohen_macaulay (matrix P, list #)
4947USAGE:   secondary_not_cohen_macaulay(P,G1,G2,...[,v]);
4948         P: a 1xn <matrix> with primary invariants, G1,G2,...: nxn <matrices>
4949         generating a finite matrix group, v: an optional <int>
4950ASSUME:  n is the number of variables of the basering
4951RETURN:  secondary invariants of the invariant ring (type <matrix>)
4952DISPLAY: information if v does not equal 0
4953EXAMPLE: example secondary_not_cohen_macaulay; shows an example
4954THEORY:  The secondary invariants are generated following "Generating Invariant
4955         Rings of Finite Groups over Arbitrary Fields" by Kemper (1996, to
4956         appear in JSC).
4957{ int i, j;
4958  degBound=0;
4959  def br=basering;
4960  int n=nvars(br);                     // n is the number of variables, as well
4961                                       // as the size of the matrices, as well
4962                                       // as the number of primary invariants,
4963                                       // we should get -
4964  if (size(#)>0)                       // checking input and setting verbose
4965  { if (typeof(#[size(#)])=="int")
4966    { int gen_num=size(#)-1;
4967      if (gen_num==0)
4968      { "ERROR:   There are no generators of the finite matrix group given.";
4969        return();
4970      }
4971      int v=#[size(#)];
4972      for (i=1;i<=gen_num;i=i+1)
4973      { if (typeof(#[i])<>"matrix")
4974        { "ERROR:   These parameters should be generators of the finite matrix group.";
4975          return();
4976        }
4977        if ((n<>nrows(#[i])) or (n<>ncols(#[i])))
4978        { "ERROR:   matrices need to be square and of the same dimensions";
4979          return();
4980        }
4981      }
4982    }
4983    else
4984    { int v=0;
4985      int gen_num=size(#);
4986      for (i=1;i<=gen_num;i=i+1)
4987      { if (typeof(#[i])<>"matrix")
4988        { "ERROR:   These parameters should be generators of the finite matrix group.";
4989          return();
4990        }
4991        if ((n<>nrows(#[i])) or (n<>ncols(#[i])))
4992        { "ERROR:   matrices need to be square and of the same dimensions";
4993          return();
4994        }
4995      }
4996    }
4997  }
4998  else
4999  { "ERROR:   There are no generators of the finite matrix group given.";
5000    return();
5001  }
5002  if (ncols(P)<>n)
5003  { "ERROR:   The first parameter ought to be the matrix of the primary";
5004    "         invariants."
5005    return();
5006  }
5007  if (v && voice==2)
5008  { "";
5009  }
5010  ring alskdfalkdsj=0,x,dp;
5011  matrix M[1][2]=1,(1-x)^n;            // we look at our primary invariants as
5012  export alskdfalkdsj;
5013  export M;
5014  setring br;                          // such of the subgroup that only
5015  matrix REY=matrix(maxideal(1));      // contains the identity, this means that
5016                                       // ch does not divide the order anymore,
5017                                       // this means that we can make use of the
5018                                       // Molien series again - M[1,1]/M[1,2] is
5019                                       // the Molien series of that group, we
5020                                       // now calculate the secondary invariants
5021                                       // of this subgroup in the usual fashion
5022                                       // where the primary invariants are the
5023                                       // ones from the bigger group
5024  if (v)
5025  { "  The procedure secondary_charp() is called to calculate secondary invariants";
5026    "  of the invariant ring of the trivial group with respect to the primary";
5027    "  invariants found previously.";
5028    "";
5029  }
5030  matrix trivialS=secondary_charp(P,REY,"alskdfalkdsj",v);
5031  kill alskdfalkdsj;
5032  // now we have those secondary invariants
5033  int k=ncols(trivialS);               // k is the number of the secondary
5034                                       // invariants, we just calculated
5035  if (v)
5036  { "  We calculate secondary invariants from the ones found for the trivial";
5037    "  subgroup.";
5038    "";
5039  }
5040  map f;                               // used to let generators act on
5041                                       // secondary invariants with respect to
5042                                       // the trivial group -
5043  matrix M(1)[gen_num][k];             // M(1) will contain a module
5044  ideal B;
5045  for (i=1;i<=gen_num;i=i+1)
5046  { B=ideal(matrix(maxideal(1))*transpose(#[i]));   // image of the various
5047                                       // variables under the i-th generator -
5048    f=br,B;                            // the corresponding mapping -
5049    B=f(trivialS)-trivialS;            // these relations should be 0 -
5050    M(1)[i,1..k]=B[1..k];              // we will look for the syzygies of M(1)
5051  }
5052  module M(2)=res(M(1),2)[2];
5053  int m=ncols(M(2));                   // number of generators of the module
5054                                       // M(2) -
5055  // the following steps calculates the intersection of the module M(2) with
5056  // the algebra A^k where A denote the subalgebra of the usual polynomial
5057  // ring, generated by the primary invariants
5058  string mp=string(minpoly);           // generating a ring where we can do
5059                                       // elimination
5060  execute "ring R=("+charstr(br)+"),(x(1..n),y(1..n),h),dp;";
5061  execute "minpoly=number("+mp+");";
5062  map f=br,maxideal(1);                // canonical mapping
5063  matrix M[k][m+k*n];
5064  M[1..k,1..m]=matrix(f(M(2)));        // will contain a module -
5065  matrix P=f(P);                       // primary invariants in the new ring
5066  for (i=1;i<=n;i=i+1)
5067  { for (j=1;j<=k;j=j+1)
5068    { M[j,m+(i-1)*k+j]=y(i)-P[1,i];
5069    }
5070  }
5071  M=elim(module(M),1,n);               // eliminating x(1..n), std-calculation
5072                                       // is done internally -
5073  M=homog(module(M),h);                // homogenize for 'minbase'
5074  M=minbase(module(M));
5075  setring br;
5076  ideal substitute=maxideal(1),ideal(P),1;
5077  f=R,substitute;                      // replacing y(1..n) by primary
5078                                       // invariants -
5079  M(2)=f(M);                           // M(2) is the new module
5080  m=ncols(M(2));
5081  matrix S[1][m];
5082  S=matrix(trivialS)*matrix(M(2));     // S now contains the secondary
5083                                       // invariants
5084  for (i=1; i<=m;i=i+1)
5085  { S[1,i]=S[1,i]/leadcoef(S[1,i]); // making elements nice
5086  }
5087  S=sort(ideal(S))[1];
5088  if (v)
5089  { "  These are the secondary invariants: ";
5090    for (i=1;i<=m;i=i+1)
5091    { "   "+string(S[1,i]);
5092    }
5093    "";
5094    "  We're done!";
5095    "";
5096  }
5097  if ((v or (voice==2)) && (m>1))
5098  { "  WARNING: The invariant ring might not have a Hironaka decomposition";
5099    "           if the characteristic of the coefficient field divides the";
5100    "           group order.";
5101  }
5102  return(S);
5103}
5104example
5105{ echo=2;
5106           ring R=2,(x,y,z),dp;
5107           matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
5108           list L=primary_invariants(A);
5109           matrix S=secondary_not_cohen_macaulay(L[1],A);
5110           print(S);
5111}
5112
5113proc invariant_ring (list #)
5114USAGE:   invariant_ring(G1,G2,...[,flags]);
5115         G1,G2,...: <matrices> generating a finite matrix group, flags: an
5116         optional <intvec> with three entries: if the first one equals 0, the
5117         program attempts to compute the Molien series and Reynolds operator,
5118         if it equals 1, the program is told that the Molien series should not
5119         be computed, if it equals -1 characteristic 0 is simulated, i.e. the
5120         Molien series is computed as if the base field were characteristic 0
5121         (the user must choose a field of large prime characteristic, e.g.
5122         32003) and if the first one is anything else, it means that the
5123         characteristic of the base field divides the group order (i.e. it will
5124         not even be attempted to compute the Reynolds operator or Molien
5125         series), the second component should give the size of intervals
5126         between canceling common factors in the expansion of the Molien series,
5127         0 (the default) means only once after generating all terms, in prime
5128         characteristic also a negative number can be given to indicate that
5129         common factors should always be canceled when the expansion is simple
5130         (the root of the extension field does not occur among the coefficients)
5131RETURN:  primary and secondary invariants (both of type <matrix>) generating the
5132         invariant ring with respect to the matrix group generated by the
5133         matrices in the input and irreducible secondary invariants (type
5134         <matrix>) if the Molien series was available
5135DISPLAY: information about the various stages of the program if the third flag
5136         does not equal 0
5137EXAMPLE: example invariant_ring; shows an example
5138THEORY:  Bases of homogeneous invariants are generated successively and those
5139         are chosen as primary invariants that lower the dimension of the ideal
5140         generated by the previously found invariants (see paper "Generating a
5141         Noetherian Normalization of the Invariant Ring of a Finite Group" by
5142         Decker, Heydtmann, Schreyer (1997) to appear in JSC). In the
5143         non-modular case secondary invariants are calculated by finding a
5144         basis (in terms of monomials) of the basering modulo the primary
5145         invariants, mapping to invariants with the Reynolds operator and using
5146         those or their power products such that they are linearly independent
5147         modulo the primary invariants (see paper "Some Algorithms in Invariant
5148         Theory of Finite Groups" by Kemper and Steel (1997)). In the modular
5149         case they are generated according to "Generating Invariant Rings of
5150         Finite Groups over Arbitrary Fields" by Kemper (1996, to appear in
5151         JSC).
5152{ if (size(#)==0)
5153  { "ERROR:   There are no generators given.";
5154    return();
5155  }
5156  int ch=char(basering);               // the algorithms depend very much on the
5157                                       // characteristic of the ground field -
5158  int n=nvars(basering);               // n is the number of variables, as well
5159                                       // as the size of the matrices, as well
5160                                       // as the number of primary invariants,
5161                                       // we should get
5162  int gen_num;
5163  int mol_flag, v;
5164 //------------------- checking input and setting flags -----------------------
5165  if (typeof(#[size(#)])=="intvec")
5166  { if (size(#[size(#)])<>3)
5167    { "ERROR:   The <intvec> should have three entries.";
5168      return();
5169    }
5170    gen_num=size(#)-1;
5171    mol_flag=#[size(#)][1];
5172    if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0)))
5173    { "ERROR:   the second component of <intvec> should be >=0";
5174      return();
5175    }
5176    int interval=#[size(#)][2];
5177    v=#[size(#)][3];
5178  }
5179  else
5180  { gen_num=size(#);
5181    mol_flag=0;
5182    int interval=0;
5183    v=0;
5184  }
5185 //----------------------------------------------------------------------------
5186  if (mol_flag==0)                     // calculation Molien series will be
5187  { if (ch==0)                         // attempted -
5188    { matrix REY,M=reynolds_molien(#[1..gen_num],intvec(0,interval,v)); // one
5189                                       // will contain Reynolds operator and the
5190                                       // other enumerator and denominator of
5191                                       // Molien series
5192      matrix P=primary_char0(REY,M,v);
5193      matrix S,IS=secondary_char0(P,REY,M,v);
5194      return(P,S,IS);
5195    }
5196    else
5197    { list L=group_reynolds(#[1..gen_num],v);
5198      if (L[1]<>0)                     // testing whether we are in the modular
5199      { string newring="aksldfalkdsflkj"; // case
5200        if (minpoly==0)
5201        { if (v)
5202          { "  We are dealing with the non-modular case.";
5203          }
5204          molien(L[3..size(L)],newring,L[2],intvec(0,interval,v));
5205          matrix P=primary_charp(L[1],newring,v);
5206          matrix S,IS=secondary_charp(P,L[1],newring,v];
5207          if (voice==2)
5208          { kill aksldfalkdsflkj;
5209          }
5210          return(P,S,IS);
5211        }
5212        else
5213        { if (v)
5214          { "  Since it is impossible for this programme to calculate the Molien
5215 series for";
5216            "  invariant rings over extension fields of prime characteristic, we
5217 have to";
5218            "  continue without it.";
5219            "";
5220
5221          }
5222          list l=primary_charp_no_molien(L[1],v);
5223          if (size(l)==2)
5224          { matix S=secondary_no_molien(l[1],L[1],l[2],v);
5225          }
5226          else
5227          { matix S=secondary_no_molien(l[1],L[1],v);
5228          }
5229          return(l[1],S);
5230        }
5231      }
5232      else                             // the modular case
5233      { if (v)
5234        { "  There is also no Molien series or Reynolds operator, we can make use of...";
5235          "";
5236          "  We can start looking for primary invariants...";
5237          "";
5238        }
5239        matrix P=primary_charp_without(#[1..gen_num],v);
5240        matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
5241        return(P,S);
5242      }
5243    }
5244  }
5245  if (mol_flag==1)                     // the user wants no calculation of the
5246  { list L=group_reynolds(#[1..gen_num],v); // Molien series
5247    if (ch==0)
5248    { list l=primary_char0_no_molien(L[1],v);
5249      if (size(l)==2)
5250      { matix S=secondary_no_molien(l[1],L[1],l[2],v);
5251      }
5252      else
5253      { matix S=secondary_no_molien(l[1],L[1],v);
5254      }
5255      return(l[1],S);
5256    }
5257    else
5258    { if (L[1]<>0)                     // testing whether we are in the modular
5259      { list l=primary_charp_no_molien(L[1],v); // case
5260        if (size(l)==2)
5261        { matix S=secondary_no_molien(l[1],L[1],l[2],v);
5262        }
5263        else
5264        { matix S=secondary_no_molien(l[1],L[1],v);
5265        }
5266        return(l[1],S);
5267      }
5268      else                             // the modular case
5269      { if (v)
5270        { "  We can start looking for primary invariants...";
5271          "";
5272        }
5273        matrix P=primary_charp_without(#[1..gen_num],v);
5274        matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
5275        return(L[1],S);
5276      }
5277    }
5278  }
5279  if (mol_flag==-1)
5280  { if (ch==0)
5281    { "ERROR:   Characteristic 0 can only be simulated in characteristic p>>0.
5282";
5283      return();
5284    }
5285    list L=group_reynolds(#[1..gen_num],v);
5286    string newring="aksldfalkdsflkj";
5287    molien(L[2..size(L)],newring,intvec(1,interval,v));
5288    matrix P=primary_charp(L[1],newring,v);
5289    matrix S,IS=secondary_charp(P,L[1],newring,v);
5290    kill aksldfalkdsflkj;
5291    return(P,S,IS);
5292  }
5293  else                                 // the user specified that the
5294  { if (ch==0)                         // characteristic divides the group order
5295    { "ERROR:   The characteristic cannot divide the group order when it is 0.
5296";
5297      return();
5298    }
5299    if (v)
5300    { "";
5301    }
5302    matrix P=primary_charp_without(#[1..gen_num],v);
5303    matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
5304    return(L[1],S);
5305  }
5306}
5307example
5308{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
5309  echo=2;
5310         ring R=0,(x,y,z),dp;
5311         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
5312         matrix P,S,IS=invariant_ring(A);
5313         print(P);
5314         print(S);
5315         print(IS);
5316}
5317
5318proc invariant_ring_random (list #)
5319USAGE:   invariant_ring_random(G1,G2,...,r[,flags]);
5320         G1,G2,...: <matrices> generating a finite matrix group, r: an <int>
5321         where -|r| to |r| is the range of coefficients of random
5322         combinations of bases elements that serve as primary invariants,
5323         flags: an optional <intvec> with three entries: if the first one equals
5324         0, the program attempts to compute the Molien series and Reynolds
5325         operator, if it equals 1, the program is told that the Molien series
5326         should not be computed, if it equals -1 characteristic 0 is simulated,
5327         i.e. the Molien series is computed as if the base field were
5328         characteristic 0 (the user must choose a field of large prime
5329         characteristic, e.g.  32003) and if the first one is anything else, it
5330         means that the characteristic of the base field divides the group order
5331         (i.e. it will not even be attempted to compute the Reynolds operator or
5332         Molien series), the second component should give the size of intervals
5333         between canceling common factors in the expansion of the Molien series,
5334         0 (the default) means only once after generating all terms, in prime
5335         characteristic also a negative number can be given to indicate that
5336         common factors should always be canceled when the expansion is simple
5337         (the root of the extension field does not occur among the coefficients)
5338RETURN:  primary and secondary invariants (both of type <matrix>) generating the
5339         invariant ring with respect to the matrix group generated by the
5340         matrices in the input and irreducible secondary invariants (type
5341         <matrix>) if the Molien series was available
5342DISPLAY: information about the various stages of the program if the third flag
5343         does not equal 0
5344EXAMPLE: example invariant_ring_random; shows an example
5345THEORY:  is the same as for invariant_ring except that random combinations of
5346         basis elements are chosen as candidates for primary invariants and
5347         hopefully they lower the dimension of the previously found primary
5348         invariants by the right amount.
5349{ if (size(#)<2)
5350  { "ERROR:   There are too few parameters.";
5351    return();
5352  }
5353  int ch=char(basering);               // the algorithms depend very much on the
5354                                       // characteristic of the ground field
5355  int n=nvars(basering);               // n is the number of variables, as well
5356                                       // as the size of the matrices, as well
5357                                       // as the number of primary invariants,
5358                                       // we should get
5359  int gen_num;
5360  int mol_flag, v;
5361 //------------------- checking input and setting flags -----------------------
5362  if (typeof(#[size(#)])=="intvec" && typeof(#[size(#)-1])=="int")
5363  { if (size(#[size(#)])<>3)
5364    { "ERROR:   <intvec> should have three entries.";
5365      return();
5366    }
5367    gen_num=size(#)-2;
5368    mol_flag=#[size(#)][1];
5369    if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0)))
5370    { "ERROR:   the second component of <intvec> should be >=0";
5371      return();
5372    }
5373    int interval=#[size(#)][2];
5374    v=#[size(#)][3];
5375    int max=#[size(#)-1];
5376    if (gen_num==0)
5377    { "ERROR:   There are no generators of a finite matrix group given.";
5378      return();
5379    }
5380  }
5381  else
5382  { if (typeof(#[size(#)])=="int")
5383    { gen_num=size(#)-1;
5384      mol_flag=0;
5385      int interval=0;
5386      v=0;
5387      int max=#[size(#)];
5388    }
5389   else
5390    { "ERROR:   If the two last parameters are not <int> and <intvec>, the last";
5391      "         parameter should be an <int>.";
5392      return();
5393    }
5394  }
5395  for (int i=1;i<=gen_num;i=i+1)
5396  { if (typeof(#[i])=="matrix")
5397    { if (nrows(#[i])<>n or ncols(#[i])<>n)
5398      { "ERROR:   The number of variables of the base ring needs to be the same";
5399        "         as the dimension of the square matrices";
5400        return();
5401      }
5402    }
5403    else
5404    { "ERROR:   The first parameters should be a list of matrices";
5405      return();
5406    }
5407  }
5408 //----------------------------------------------------------------------------
5409  if (mol_flag==0)
5410  { if (ch==0)
5411    { matrix REY,M=reynolds_molien(#[1..gen_num],intvec(0,interval,v)); // one
5412                                       // will contain Reynolds operator and the
5413                                       // other enumerator and denominator of
5414                                       // Molien series
5415      matrix P=primary_char0_random(REY,M,max,v);
5416      matrix S,IS=secondary_char0(P,REY,M,v);
5417      return(P,S,IS);
5418    }
5419    else
5420    { list L=group_reynolds(#[1..gen_num],v);
5421      if (L[1]<>0)                     // testing whether we are in the modular
5422      { string newring="aksldfalkdsflkj"; // case
5423        if (minpoly==0)
5424        { if (v)
5425          { "  We are dealing with the non-modular case.";
5426          }
5427          molien(L[3..size(L)],newring,L[2],intvec(0,interval,v));
5428          matrix P=primary_charp_random(L[1],newring,max,v);
5429          matrix S,IS=secondary_charp(P,L[1],newring,v);
5430          if (voice==2)
5431          { kill aksldfalkdsflkj;
5432          }
5433          return(P,S,IS);
5434        }
5435        else
5436        { if (v)
5437          { "  Since it is impossible for this programme to calculate the Molien
5438 series for";
5439            "  invariant rings over extension fields of prime characteristic, we
5440 have to";
5441            "  continue without it.";
5442            "";
5443
5444          }
5445          list l=primary_charp_no_molien_random(L[1],max,v);
5446          if (size(l)==2)
5447          { matix S=secondary_no_molien(l[1],L[1],l[2],v);
5448          }
5449          else
5450          { matix S=secondary_no_molien(l[1],L[1],v);
5451          }
5452          return(l[1],S);
5453        }
5454      }
5455      else                             // the modular case
5456      { if (v)
5457        { "  There is also no Molien series, we can make use of...";
5458          "";
5459          "  We can start looking for primary invariants...";
5460          "";
5461        }
5462        matrix P=primary_charp_without_random(#[1..gen_num],max,v);
5463        matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
5464        return(L[1],S);
5465      }
5466    }
5467  }
5468  if (mol_flag==1)                     // the user wants no calculation of the
5469  { list L=group_reynolds(#[1..gen_num],v); // Molien series
5470    if (ch==0)
5471    { list l=primary_char0_no_molien_random(L[1],max,v);
5472      if (size(l)==2)
5473      { matix S=secondary_no_molien(l[1],L[1],l[2],v);
5474      }
5475      else
5476      { matix S=secondary_no_molien(l[1],L[1],v);
5477      }
5478      return(l[1],S);
5479    }
5480    else
5481    { if (L[1]<>0)                     // testing whether we are in the modular
5482      { list l=primary_charp_no_molien_random(L[1],max,v); // case
5483        if (size(l)==2)
5484        { matix S=secondary_no_molien(l[1],L[1],l[2],v);
5485        }
5486        else
5487        { matix S=secondary_no_molien(l[1],L[1],v);
5488        }
5489        return(l[1],S);
5490      }
5491      else                             // the modular case
5492      { if (v)
5493        { "  We can start looking for primary invariants...";
5494          "";
5495        }
5496        matrix P=primary_charp_without_random(#[1..gen_num],max,v);
5497        matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
5498        return(L[1],S);
5499      }
5500    }
5501  }
5502  if (mol_flag==-1)
5503  { if (ch==0)
5504    { "ERROR:   Characteristic 0 can only be simulated in characteristic p>>0.
5505";
5506      return();
5507    }
5508    list L=group_reynolds(#[1..gen_num],v);
5509    string newring="aksldfalkdsflkj";
5510    molien(L[2..size(L)],newring,intvec(1,v));
5511    matrix P=primary_charp_random(L[1],newring,max,v);
5512    matrix S,IS=secondary_charp(P,L[1],newring,v);
5513    kill aksldfalkdsflkj;
5514    return(P,S,IS);
5515  }
5516  else                                 // the user specified that the
5517  { if (ch==0)                         // characteristic divides the group order
5518    { "ERROR:   The characteristic cannot divide the group order when it is 0.
5519";
5520      return();
5521    }
5522    if (v)
5523    { "";
5524    }
5525    matrix P=primary_charp_without_random(#[1..gen_num],max,v);
5526    matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
5527    return(L[1],S);
5528  }
5529}
5530example
5531{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
5532  echo=2;
5533         ring R=0,(x,y,z),dp;
5534         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
5535         matrix P,S,IS=invariant_ring_random(A,1);
5536         print(P);
5537         print(S);
5538         print(IS);
5539}
5540
5541proc algebra_containment (poly p, matrix A)
5542USAGE:   algebra_containment(p,A);
5543         p: arbitrary <poly>, A: a 1xm <matrix> giving generators of a
5544         subalgebra of the basering
5545RETURN:  1 (TRUE) (type <int>) if p is contained in the subalgebra
5546         0 (FALSE) (type <int>) if <poly> is not contained
5547DISPLAY: a representation of p in terms of algebra generators A[1,i]=y(i) if p
5548         is contained in the subalgebra
5549EXAMPLE: example algebra_containment; shows an example
5550THEORY:  The ideal of algebraic relations of the algebra generators f1,...,fm
5551         given by A is computed introducing new variables y(i) and the product
5552         order: x^a*y^b > y^d*y^e if x^a > x^d or else if y^b > y^e. p reduces
5553         to a polynomial only in the y(i) <=> p is contained in the subring
5554         generated by the polynomials in A.
5555{ degBound=0;
5556  if (nrows(A)==1)
5557  { def br=basering;
5558    int n=nvars(br);
5559    int m=ncols(A);
5560    string mp=string(minpoly);
5561    execute "ring R=("+charstr(br)+"),(x(1..n),y(1..m)),(dp(n),dp(m));";
5562    execute "minpoly=number("+mp+");";
5563    ring R=0,(x(1..n),y(1..m)),(dp(n),dp(m));
5564    ideal vars=x(1..n);
5565    map emb=br,vars;
5566    ideal A=ideal(emb(A));
5567    ideal check=emb(p);
5568    for (int i=1;i<=m;i=i+1)
5569    { A[i]=A[i]-y(i);
5570    }
5571    A=std(A);
5572    check[1]=reduce(check[1],A);
5573    A=elim(check,1,n);
5574    if (A[1]<>0)
5575    { "\/\/ "+string(check);
5576      return(1);
5577    }
5578    else
5579    { return(0);
5580    }
5581  }
5582  else
5583  { "ERROR:   <matrix> may only have one row";
5584    return();
5585  }
5586}
5587example
5588{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
5589  echo=2;
5590         ring R=0,(x,y,z),dp;
5591         matrix A[1][7]=x2+y2,z2,x4+y4,1,x2z-1y2z,xyz,x3y-1xy3;
5592         poly p1=x10z3-x8y2z3+2x6y4z3-2x4y6z3+x2y8z3-y10z3+x6z4+3x4y2z4+3x2y4z4+y6z4;
5593         algebra_containment(p1,A);
5594         poly p2=z;
5595         algebra_containment(p2,A);
5596}
5597
5598proc module_containment(poly p, matrix P, matrix S)
5599USAGE:   module_containment(p,P,S);
5600         p: arbitrary <poly>, P: a 1xn <matrix> giving generators of an algebra,
5601         S: a 1xt <matrix> giving generators of a module over the algebra
5602         generated by P
5603ASSUME:  n is the number of variables in the basering and the generators in P
5604         are algebraically independent
5605RETURNS: 1 (TRUE) (type <int>) if p is contained in the ring
5606         0 (FALSE) type <int>) if p is not contained
5607DISPLAY: the representation of p in terms of algebra generators P[1,i]=z(i) and
5608         module generators S[1,j]=y(j) if p is contained in the module
5609EXAMPLE: example module_containment; shows an example
5610THEORY:  The ideal of algebraic relations of all the generators p1,...,pn,
5611         s1,...,st given by P and S is computed introducing new variables y(j),
5612         z(i) and the product order: x^a*y^b*z^c > x^d*y^e*z^f if x^a > x^d
5613         with respect to the lp ordering or else if z^c > z^f with respect to
5614         the dp ordering or else if y^b > y^e with respect to the lp ordering
5615         again. p reduces to a polynomial only in the y(j) and z(i) linear in
5616         the z(i)) <=> p is contained in the module.
5617{ def br=basering;
5618  degBound=0;
5619  int n=nvars(br);
5620  if (ncols(P)==n and nrows(P)==1 and nrows(S)==1)
5621  { int m=ncols(S);
5622    string mp=string(minpoly);
5623    execute "ring R=("+charstr(br)+"),(x(1..n),y(1..m),z(1..n)),(lp(n),dp(m),lp(n));";
5624    execute "minpoly=number("+mp+");";
5625    ideal vars=x(1..n);
5626    map emb=br,vars;
5627    matrix P=emb(P);
5628    matrix S=emb(S);
5629    ideal check=emb(p);
5630    ideal I;
5631    for (int i=1;i<=m;i=i+1)
5632    { I[i]=S[1,i]-y(i);
5633    }
5634    for (i=1;i<=n;i=i+1)
5635    { I[m+i]=P[1,i]-z(i);
5636    }
5637    I=std(I);
5638    check[1]=reduce(check[1],I);
5639    I=elim(check,1,n);                 // checking whether all variables from
5640    if (I[1]<>0)                       // the former ring have disappeared
5641    { "\/\/ "+string(check);
5642      return(1);
5643    }
5644    else
5645    { return(0);
5646    }
5647  }
5648  else
5649  { "ERROR:   the first <matrix> must have the same number of columns as the";
5650    "         basering and both <matrices> may only have one row";
5651    return();
5652  }
5653}
5654example
5655{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
5656  echo=2;
5657         ring R=0,(x,y,z),dp;
5658         matrix P[1][3]=x2+y2,z2,x4+y4;
5659         matrix S[1][4]=1,x2z-1y2z,xyz,x3y-1xy3;
5660         poly p1=x10z3-x8y2z3+2x6y4z3-2x4y6z3+x2y8z3-y10z3+x6z4+3x4y2z4+3x2y4z4+y6z4;
5661         module_containment(p1,P,S);
5662         poly p2=z;
5663         module_containment(p2,P,S);
5664}
5665
5666proc orbit_variety (matrix F,string newring)
5667USAGE:   orbit_variety(F,s);
5668         F: a 1xm <matrix> defing an invariant ring, s: a <string> giving the
5669         name for a new ring
5670RETURN:  a Groebner basis (type <ideal>, named G) for the ideal defining the
5671         orbit variety (i.e. the syzygy ideal) in the new ring (named `s`)
5672EXAMPLE: example orbit_variety; shows an example
5673THEORY:  The ideal of algebraic relations of the invariant ring generators is
5674         calculated, then the variables of the original ring are eliminated and
5675         the polynomials that are left over define the orbit variety
5676{ if (newring=="")
5677  { "ERROR:   the second parameter may not be an empty <string>";
5678    return();
5679  }
5680  if (nrows(F)==1)
5681  { def br=basering;
5682    int n=nvars(br);
5683    int m=ncols(F);
5684    string mp=string(minpoly);
5685    execute "ring R=("+charstr(br)+"),("+varstr(br)+",y(1..m)),dp;";
5686    execute "minpoly=number("+mp+");";
5687    ideal I=ideal(imap(br,F));
5688    for (int i=1;i<=m;i=i+1)
5689    { I[i]=I[i]-y(i);
5690    }
5691    I=elim(I,1,n);
5692    execute "ring "+newring+"=("+charstr(br)+"),(y(1..m)),dp(m);";
5693    execute "minpoly=number("+mp+");";
5694    ideal vars;
5695    for (i=2;i<=n;i=i+1)
5696    { vars[i]=0;
5697    }
5698    vars=vars,y(1..m);
5699    map emb=R,vars;
5700    ideal G=emb(I);
5701    kill emb, vars, R;
5702    keepring `newring`;
5703    // execute "keepring "+newring+";";
5704    return();
5705  }
5706  else
5707  { "ERROR:   the <matrix> may only have one row";
5708    return();
5709  }
5710}
5711example
5712{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
5713  echo=2;
5714         ring R=0,(x,y,z),dp;
5715         matrix F[1][7]=x2+y2,z2,x4+y4,1,x2z-1y2z,xyz,x3y-1xy3;
5716         string newring="E";
5717         orbit_variety(F,newring);
5718         print(G);
5719         basering;
5720}
5721
5722proc relative_orbit_variety(ideal I,matrix F,string newring)
5723USAGE:   relative_orbit_variety(I,F,s);
5724         I: an <ideal> invariant under the action of a group, F: a 1xm
5725         <matrix> defining the invariant ring of this group, s: a <string>
5726         giving a name for a new ring
5727RETURN:  a Groebner basis (type <ideal>, named G) for the ideal defining the
5728         relative orbit variety with respect to I in the new ring (named s)
5729EXAMPLE: example relative_orbit_variety; shows an example
5730THEORY:  A Groebner basis of the ideal of algebraic relations of the invariant
5731         ring generators is calculated, then one of the basis elements plus the
5732         ideal generators. The variables of the original ring are eliminated and
5733         the polynomials that are left over define thecrelative orbit variety
5734         with respect to I.
5735{ if (newring=="")
5736  { "ERROR:   the third parameter may not be empty a <string>";
5737    return();
5738  }
5739  degBound=0;
5740  if (nrows(F)==1)
5741  { def br=basering;
5742    int n=nvars(br);
5743    int m=ncols(F);
5744    string mp=string(minpoly);
5745    execute "ring R=("+charstr(br)+"),("+varstr(br)+",y(1..m)),lp;";
5746    execute "minpoly=number("+mp+");";
5747    ideal J=ideal(imap(br,F));
5748    ideal I=imap(br,I);
5749    for (int i=1;i<=m;i=i+1)
5750    { J[i]=J[i]-y(i);
5751    }
5752    J=std(J);
5753    J=J,I;
5754    option(redSB);
5755    J=std(J);
5756    ideal vars;
5757    //for (i=1;i<=n;i=i+1)
5758    //{ vars[i]=0;
5759    //}
5760    vars[n]=0;
5761    vars=vars,y(1..m);
5762    map emb=R,vars;
5763    ideal G=emb(J);
5764    J=J-G;
5765    for (i=1;i<=ncols(G);i=i+1)
5766    { if (J[i]<>0)
5767      { G[i]=0;
5768      }
5769    }
5770    G=compress(G);
5771    execute "ring "+newring+"=("+charstr(br)+"),(y(1..m)),lp;";
5772    execute "minpoly=number("+mp+");";
5773    ideal vars;
5774    for (i=2;i<=n;i=i+1)
5775    { vars[i]=0;
5776    }
5777    vars=vars,y(1..m);
5778    map emb=R,vars;
5779    ideal G=emb(G);
5780    kill vars, emb;
5781    keepring `newring`;
5782    // execute "keepring "+newring+";";
5783    return();
5784  }
5785  else
5786  { "ERROR:   the <matrix> may only have one row";
5787    return();
5788  }
5789}
5790example
5791{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.6.3:";
5792  echo=2;
5793         ring R=0,(x,y,z),dp;
5794         matrix F[1][3]=x+y+z,xy+xz+yz,xyz;
5795         ideal I=x2+y2+z2-1,x2y+y2z+z2x-2x-2y-2z,xy2+yz2+zx2-2x-2y-2z;
5796         string newring="E";
5797         relative_orbit_variety(I,F,newring);
5798         print(G);
5799         basering;
5800}
5801
5802proc image_of_variety(ideal I,matrix F)
5803USAGE:   image_of_variety(I,F);
5804         I: an arbitray <ideal>, F: a 1xm <matrix> defining an invariant ring
5805         of a some matrix group
5806RETURN:  the <ideal> defining the image under that group of the variety defined
5807         by I
5808EXAMPLE: example image_of_variety; shows an example
5809THEORY:  relative_orbit_variety(I,F,s) is called and the newly introduced
5810         variables in the output are replaced by the generators of the
5811         invariant ring. This ideal in the original variables defines the image
5812         of the variety defined by I
5813{ if (nrows(F)==1)
5814  { def br=basering;
5815    int n=nvars(br);
5816    string newring="E";
5817    relative_orbit_variety(I,F,newring);
5818    execute "ring R=("+charstr(br)+"),("+varstr(br)+","+varstr(E)+"),lp;";
5819    ideal F=imap(br,F);
5820    for (int i=1;i<=n;i=i+1)
5821    { F=0,F;
5822    }
5823    setring br;
5824    map emb2=E,F;
5825    return(compress(emb2(G)));
5826  }
5827  else
5828  { "ERROR:   the <matrix> may only have one row";
5829    return();
5830  }
5831}
5832example
5833{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.6.8:";
5834  echo=2;
5835         ring R=0,(x,y,z),dp;
5836         matrix F[1][3]=x+y+z,xy+xz+yz,xyz;
5837         ideal I=xy;
5838         print(image_of_variety(I,F));
5839}
Note: See TracBrowser for help on using the repository browser.