source: git/Singular/LIB/finvar.lib @ 5811fb

spielwiese
Last change on this file since 5811fb was 5811fb, checked in by Hans Schönemann <hannes@…>, 25 years ago
*** empty log message *** git-svn-id: file:///usr/local/Singular/svn/trunk@2702 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 221.4 KB
Line 
1// $Id: finvar.lib,v 1.15 1998-12-02 10:52:08 Singular Exp $
2// author: Agnes Eileen Heydtmann, email:agnes@math.uni-sb.de
3// last change: 98/11/05
4//////////////////////////////////////////////////////////////////////////////
5version="$Id: finvar.lib,v 1.15 1998-12-02 10:52:08 Singular Exp $"
6info="
7LIBRARY:  finvar.lib       LIBRARY TO CALCULATE INVARIANT RINGS & MORE
8
9 A library for computing polynomial invariants of finite matrix groups and
10 generators of related varieties. The algorithms are based on B. Sturmfels,
11 G. Kemper and Decker et al..
12 Author: Agnes E. Heydtmann, agnes@math.uni-sb.de
13
14 MAIN PROCEDURES:
15 invariant_ring()                  generators of the invariant ring (i.r.)
16 invariant_ring_random()           generators of the i.r., randomized alg.
17 primary_invariants()              primary invariants (p.i.)
18 primary_invariants_random()       primary invariants, randomized alg.
19
20 SUB-PROCEDURES:
21 cyclotomic()                      cyclotomic polynomial
22 group_reynolds()                  finite group and Reynolds operator (R.o.)
23 molien()                          Molien series (M.s.)
24 reynolds_molien()                 Reynolds operator and Molien series
25 partial_molien()                  partial expansion of Molien series
26 evaluate_reynolds()               image under the Reynolds operator
27 invariant_basis()                 basis of homogeneous invariants of a degree
28 invariant_basis_reynolds()        as invariant_basis(), with R.o.
29 primary_char0()                   primary invariants in char 0
30 primary_charp()                   primary invariant in char p
31 primary_char0_no_molien()         p.i., char 0, without Molien series
32 primary_charp_no_molien()         p.i., char p, without Molien series
33 primary_charp_without()           p.i., char p, without R.o. or Molien series
34 primary_char0_random()            primary invariants in char 0, randomized
35 primary_charp_random()            primary invariants in char p, randomized
36 primary_char0_no_molien_random()  p.i., char 0, without M.s., randomized
37 primary_charp_no_molien_random()  p.i., char p, without M.s., randomized
38 primary_charp_without_random()    p.i., char p, without R.o. or M.s., random.
39 power_products()                  exponents for power products
40 secondary_char0()                 secondary (s.i.) invariants in char 0
41 secondary_charp()                 secondary invariants in char p
42 secondary_no_molien()             secondary invariants, without Molien series
43 secondary_and_irreducibles_no_molien() s.i. & irreducible s.i., without M.s.
44 secondary_not_cohen_macaulay()    s.i. when invariant ring not Cohen-Macaulay
45 algebra_containment()             query of algebra containment
46 module_containment()              query of module containment
47 orbit_variety()                   ideal of the orbit variety
48 relative_orbit_variety()          ideal of a relative orbit variety
49 image_of_variety()                ideal of the image of a variety
50";
51////////////////////////////////////////////////////////////////////////////////
52// perhaps useful procedures (no help provided):
53// unique()                        is a matrix among other matrices?
54// exponent()                      gives the exponent of a number
55// sort_of_invariant_basis()       lin. ind. invariants of a degree mod p.i.
56// next_vector                     lists all of Z^n with first nonzero entry 1
57// int_number_map                  integers 1..q are maped to q field elements
58// search                          searches a number of p.i., char 0
59// p_search                        searches a number of p.i., char p
60// search_random                   searches a # of p.i., char 0, randomized
61// p_search_random                 searches a # of p.i., char p, randomized
62// concat_intmat                   concatenates two integer matrices
63////////////////////////////////////////////////////////////////////////////////
64LIB "matrix.lib";
65LIB "elim.lib";
66LIB "general.lib";
67////////////////////////////////////////////////////////////////////////////////
68
69////////////////////////////////////////////////////////////////////////////////
70// Checks whether the last parameter, being a matrix, is among the previous
71// parameters, also being matrices
72////////////////////////////////////////////////////////////////////////////////
73proc unique (list #)
74{ for (int i=1;i<size(#);i=i+1)
75  { if (#[i]==#[size(#)])
76    { return(0);
77    }
78  }
79  return(1);
80}
81
82proc cyclotomic (int i)
83"USAGE:   cyclotomic(i);
84         i: an <int> > 0
85RETURNS: the i-th cyclotomic polynomial (type <poly>) as one in the first ring
86         variable
87EXAMPLE: example cyclotomic; shows an example
88THEORY:  x^i-1 is divided by the j-th cyclotomic polynomial where j takes on the
89         value of proper divisors of i
90"
91{ if (i<=0)
92  { "ERROR:   the input should be > 0.";
93    return();
94  }
95  poly v1=var(1);
96  if (i==1)
97  { return(v1-1);                      // 1-st cyclotomic polynomial
98  }
99  poly min=v1^i-1;
100  matrix s[1][2];
101  min=min/(v1-1);                      // dividing by the 1-st cyclotomic
102                                       // polynomial
103  int j=2;
104  int n;
105  poly c;
106  int flag=1;
107  while(2*j<=i)                        // there are no proper divisors of i
108  { if ((i%j)==0)                      // greater than i/2
109    { if (flag==1)
110      { n=j;                           // n stores the first proper divisor of
111      }                                // i > 1
112      flag=0;
113      c=cyclotomic(j);                 // recursive computation
114      s=min,c;
115      s=matrix(syz(ideal(s)));         // dividing
116      min=s[2,1];
117    }
118    if (n*j==i)                        // the earliest possible point to break
119    { break;
120    }
121    j=j+1;
122  }
123  min=min/leadcoef(min);               // making sure that the leading
124  return(min);                         // coefficient is 1
125}
126example
127{ echo=2;
128          ring R=0,(x,y,z),dp;
129          print(cyclotomic(25));
130}
131
132proc group_reynolds (list #)
133"USAGE:   group_reynolds(G1,G2,...[,v]);
134         G1,G2,...: nxn <matrices> generating a finite matrix group, v: an
135         optional <int>
136ASSUME:  n is the number of variables of the basering, g the number of group
137         elements
138RETURN:  a <list>, the first list element will be a gxn <matrix> representing
139         the Reynolds operator if we are in the non-modular case; if the
140         characteristic is >0, minpoly==0 and the finite group non-cyclic the
141         second list element is an <int> giving the lowest common multiple of
142         the matrix group elements (used in molien); in general all other list
143         elements are nxn <matrices> listing all elements of the finite group
144DISPLAY: information if v does not equal 0
145EXAMPLE: example group_reynolds; shows an example
146THEORY:  The entire matrix group is generated by getting all left products of
147         the generators with the new elements from the last run through the loop
148         (or the generators themselves during the first run). All the ones that
149         have been generated before are thrown out and the program terminates
150         when there are no new elements found in one run. Additionally each time
151         a new group element is found the corresponding ring mapping of which
152         the Reynolds operator is made up is generated. They are stored in the
153         rows of the first return value.
154"
155{ int ch=char(basering);               // the existance of the Reynolds operator
156                                       // is dependent on the characteristic of
157                                       // the base field
158  int gen_num;                         // number of generators
159 //------------------------ making sure the input is okay ---------------------
160  if (typeof(#[size(#)])=="int")
161  { if (size(#)==1)
162    { "ERROR:   there are no matrices given among the parameters";
163      return();
164    }
165    int v=#[size(#)];
166    gen_num=size(#)-1;
167  }
168  else                                 // last parameter is not <int>
169  { int v=0;                           // no information is default
170    gen_num=size(#);
171  }
172  if (typeof(#[1])<>"matrix")
173  { "ERROR:   the parameters must be a list of matrices and maybe an <int>";
174    return();
175  }
176  int n=nrows(#[1]);
177  if (n<>nvars(basering))
178  { "ERROR:   the number of variables of the basering needs to be the same";
179    "         as the dimension of the matrices";
180    return();
181  }
182  if (n<>ncols(#[1]))
183  { "ERROR:   matrices need to be square and of the same dimensions";
184    return();
185  }
186  matrix vars=matrix(maxideal(1));     // creating an nx1-matrix containing the
187  vars=transpose(vars);                // variables of the ring -
188  matrix REY=#[1]*vars;                // calculating the first ring mapping -
189                                       // REY will contain the Reynolds
190                                       // operator -
191  matrix G(1)=#[1];                    // G(k) are elements of the group -
192  if (ch<>0 && minpoly==0 && gen_num<>1) // finding out of which order the group
193  { matrix I=diag(1,n);                // element is
194    matrix TEST=G(1);
195    int o1=1;
196    int o2;
197    while (TEST<>I)
198    { TEST=TEST*G(1);
199      o1=o1+1;
200    }
201  }
202  int i=1;
203 // -------------- doubles among the generators should be avoided -------------
204  for (int j=2;j<=gen_num;j=j+1)       // this loop adds the parameters to the
205  {                                    // group, leaving out doubles and
206                                       // checking whether the parameters are
207                                       // compatible with the task of the
208                                       // procedure
209    if (not(typeof(#[j])=="matrix"))
210    { "ERROR:   the parameters must be a list of matrices and maybe an <int>";
211      return();
212    }
213    if ((n!=nrows(#[j])) or (n!=ncols(#[j])))
214    { "ERROR:   matrices need to be square and of the same dimensions";
215       return();
216    }
217    if (unique(G(1..i),#[j]))
218    { i=i+1;
219      matrix G(i)=#[j];
220      if (ch<>0 && minpoly==0)         // finding out of which order the group
221      { TEST=G(i);                     // element is
222        o2=1;
223        while (TEST<>I)
224        { TEST=TEST*G(i);
225          o2=o2+1;
226        }
227        o1=o1*o2/gcd(o1,o2);           // lowest common multiple of the element
228      }                                // orders -
229      REY=concat(REY,#[j]*vars);       // adding ring homomorphisms to REY
230    }
231  }
232  int g=i;                             // G(1)..G(i) are generators without
233                                       // doubles - g generally is the number
234                                       // of elements in the group so far -
235  j=i;                                 // j is the number of new elements that
236                                       // we use as factors
237  int k, m, l;
238  if (v)
239  { "";
240    "  Generating the entire matrix group and the Reynolds operator...";
241    "";
242  }
243 // -------------- main loop that finds all the group elements ----------------
244  while (1)
245  { l=0;                               // l is the number of products we get in
246                                       // one going
247    for (m=g-j+1;m<=g;m=m+1)
248    { for (k=1;k<=i;k=k+1)
249      { l=l+1;
250        matrix P(l)=G(k)*G(m);         // possible new element
251      }
252    }
253    j=0;
254    for (k=1;k<=l;k=k+1)
255    { if (unique(G(1..g),P(k)))
256      { j=j+1;                         // a new factor for next run
257        g=g+1;
258        matrix G(g)=P(k);              // a new group element -
259        if (ch<>0 && minpoly==0 && i<>1) // finding out of which order the group
260        { TEST=G(g);                   // element is
261          o2=1;
262          while (TEST<>I)
263          { TEST=TEST*G(g);
264            o2=o2+1;
265          }
266          o1=o1*o2/gcd(o1,o2);         // lowest common multiple of the element
267        }                              // orders -
268        REY=concat(REY,P(k)*vars);     // adding new mapping to REY
269        if (v)
270        { "  Group element "+string(g)+" has been found.";
271        }
272      }
273      kill P(k);
274    }
275    if (j==0)                          // when we didn't add any new elements
276    { break;                           // in one run through the while loop
277    }                                  // we are done
278  }
279  if (v)
280  { if (g<=i)
281    { "  There are only "+string(g)+" group elements.";
282    }
283    "";
284  }
285  REY=transpose(REY);                  // when we evaluate the Reynolds operator
286                                       // later on, we actually want 1xn
287                                       // matrices
288  if (ch<>0)
289  { if ((g%ch)==0)
290    { if (voice==2)
291      { "WARNING: The characteristic of the coefficient field divides the group order.";
292        "         Proceed without the Reynolds operator!";
293      }
294      else
295      { if (v)
296        { "  The characteristic of the base field divides the group order.";
297          "  We have to continue without Reynolds operator...";
298          "";
299        }
300      }
301      kill REY;
302      matrix REY[1][1]=0;
303      return(REY,G(1..g));
304    }
305    if (minpoly==0)
306    { if (i>1)
307      { return(REY,o1,G(1..g));
308      }
309      return(REY,G(1..g));
310    }
311  }
312  if (v)
313  { "  Done generating the group and the Reynolds operator.";
314    "";
315  }
316  return(REY,G(1..g));
317}
318example
319{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
320  echo=2;
321         ring R=0,(x,y,z),dp;
322         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
323         list L=group_reynolds(A);
324         print(L[1]);
325         print(L[2..size(L)]);
326}
327
328////////////////////////////////////////////////////////////////////////////////
329// Returns i such that root^i==n, i.e. it heavily relies on the right input.
330////////////////////////////////////////////////////////////////////////////////
331proc exponent(number n, number root)
332{ int i=0;
333   while((n/root^i)<>1)
334   { i=i+1;
335   }
336   return(i);
337}
338
339proc molien (list #)
340"USAGE:   molien(G1,G2,...[,ringname,lcm,flags]);
341         G1,G2,...: nxn <matrices> generating a finite matrix group, ringname:
342         a <string> giving a name for a new ring of characteristic 0 for the
343         Molien series in case of prime characteristic, lcm: an <int> giving the
344         lowest common multiple of the elements' orders in case of prime
345         characteristic, minpoly==0 and a non-cyclic group, flags: an optional
346         <intvec> with three components: if the first element is not equal to 0
347         characteristic 0 is simulated, i.e. the Molien series is computed as
348         if the base field were characteristic 0 (the user must choose a field
349         of large prime characteristic, e.g. 32003), the second component should
350         give the size of intervals between canceling common factors in the
351         expansion of the Molien series, 0 (the default) means only once after
352         generating all terms, in prime characteristic also a negative number
353         can be given to indicate that common factors should always be canceled
354         when the expansion is simple (the root of the extension field does not
355         occur among the coefficients)
356ASSUME:  n is the number of variables of the basering, G1,G2... are the group
357         elements generated by group_reynolds(), lcm is the second return value
358         of group_reynolds()
359RETURN:  in case of characteristic 0 a 1x2 <matrix> giving enumerator and
360         denominator of Molien series; in case of prime characteristic a ring
361         with the name `ringname` of characteristic 0 is created where the same
362         Molien series (named M) is stored
363DISPLAY: information if the third component of flags does not equal 0
364EXAMPLE: example molien; shows an example
365THEORY:  In characteristic 0 the terms 1/det(1-xE) for all group elements of the
366         Molien series are computed in a straight forward way. In prime
367         characteristic a Brauer lift is involved. The returned matrix gives
368         enumerator and denominator of the expanded version where common factors
369         have been canceled.
370"
371{ def br=basering;                     // the Molien series depends on the
372  int ch=char(br);                     // characteristic of the coefficient
373                                       // field -
374  int g;                               // size of the group
375 //---------------------- making sure the input is okay -----------------------
376  if (typeof(#[size(#)])=="intvec")
377  { if (size(#[size(#)])==3)
378    { int mol_flag=#[size(#)][1];
379      if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0)))
380      { "ERROR:   the second component of <intvec> should be >=0"
381        return();
382      }
383      int interval=#[size(#)][2];
384      int v=#[size(#)][3];
385    }
386    else
387    { "ERROR:   <intvec> should have three components";
388      return();
389    }
390    if (ch<>0)
391    { if (typeof(#[size(#)-1])=="int")
392      { int r=#[size(#)-1];
393        if (typeof(#[size(#)-2])<>"string")
394        { "ERROR:   in characteristic p>0 a <string> must be given for the name of a new";
395          "         ring where the Molien series can be stored";
396          return();
397        }
398        else
399        { if (#[size(#)-2]=="")
400          { "ERROR:   <string> may not be empty";
401            return();
402          }
403          string newring=#[size(#)-2];
404          g=size(#)-3;
405        }
406      }
407      else
408      { if (typeof(#[size(#)-1])<>"string")
409        { "ERROR:   in characteristic p>0 a <string> must be given for the name of a new";
410          "         ring where the Molien series can be stored";
411          return();
412        }
413        else
414        { if (#[size(#)-1]=="")
415          { "ERROR:   <string> may not be empty";
416            return();
417          }
418          string newring=#[size(#)-1];
419          g=size(#)-2;
420          int r=g;
421        }
422      }
423    }
424    else                               // then <string> ist not needed
425    { g=size(#)-1;
426    }
427  }
428  else                                 // last parameter is not <intvec>
429  { int v=0;                           // no information is default
430    int mol_flag=0;                    // computing of Molien series is default
431    int interval=0;
432    if (ch<>0)
433    { if (typeof(#[size(#)])=="int")
434      { int r=#[size(#)];
435        if (typeof(#[size(#)-1])<>"string")
436        { "ERROR:   in characteristic p>0 a <string> must be given for the name of a new";
437          "         ring where the Molien series can be stored";
438            return();
439        }
440        else
441        { if (#[size(#)-1]=="")
442          { "ERROR:   <string> may not be empty";
443            return();
444          }
445          string newring=#[size(#)-1];
446          g=size(#)-2;
447        }
448      }
449      else
450      { if (typeof(#[size(#)])<>"string")
451        { "ERROR:   in characteristic p>0 a <string> must be given for the name of a new";
452          "         ring where the Molien series can be stored";
453          return();
454        }
455        else
456        { if (#[size(#)]=="")
457          { "ERROR:   <string> may not be empty";
458            return();
459          }
460          string newring=#[size(#)];
461          g=size(#)-1;
462          int r=g;
463        }
464      }
465    }
466    else
467    { g=size(#);
468    }
469  }
470  if (ch<>0)
471  { if ((g/r)*r<>g)
472   { "ERROR:   <int> should divide the group order."
473      return();
474    }
475  }
476  if (ch<>0)
477  { if ((g%ch)==0)
478    { if (voice==2)
479      { "WARNING: The characteristic of the coefficient field divides the group";
480        "         order. Proceed without the Molien series!";
481      }
482      else
483      { if (v)
484        { "  The characteristic of the base field divides the group order.";
485          "  We have to continue without Molien series...";
486          "";
487        }
488      }
489    }
490    if (minpoly<>0 && mol_flag==0)
491    { if (voice==2)
492      { "WARNING: It is impossible for this program to calculate the Molien series";
493        "         for finite groups over extension fields of prime characteristic.";
494      }
495      else
496      { if (v)
497        { "  Since it is impossible for this program to calculate the Molien series for";
498          "  invariant rings over extension fields of prime characteristic, we have to";
499          "  continue without it.";
500          "";
501        }
502      }
503      return();
504    }
505  }
506 //----------------------------------------------------------------------------
507  if (not(typeof(#[1])=="matrix"))
508  { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
509    return();
510  }
511  int n=nrows(#[1]);
512  if (n<>nvars(br))
513  { "ERROR:   the number of variables of the basering needs to be the same";
514    "         as the dimension of the square matrices";
515    return();
516  }
517  if (v && voice<>2)
518  { "";
519    "  Generating the Molien series...";
520    "";
521  }
522  if (v && voice==2)
523  { "";
524  }
525 //------------- calculating Molien series in characteristic 0 ----------------
526  if (ch==0)                           // when ch==0 we can calculate the Molien
527  { matrix I=diag(1,n);                // series in any case -
528    poly v1=maxideal(1)[1];            // the Molien series will be in terms of
529                                       // the first variable of the current
530                                       // ring -
531    matrix M[1][2];                    // M will contain the Molien series -
532    M[1,1]=0;                          // M[1,1] will be the numerator -
533    M[1,2]=1;                          // M[1,2] will be the denominator -
534    matrix s;                          // will help us canceling in the
535                                       // fraction
536    poly p;                            // will contain the denominator of the
537                                       // new term of the Molien series
538 //------------ computing 1/det(1+xE) for all E in the group ------------------
539    for (int j=1;j<=g;j=j+1)
540    { if (not(typeof(#[j])=="matrix"))
541      { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
542        return();
543      }
544      if ((n<>nrows(#[j])) or (n<>ncols(#[j])))
545      { "ERROR:   matrices need to be square and of the same dimensions";
546         return();
547      }
548      p=det(I-v1*#[j]);                // denominator of new term -
549      M[1,1]=M[1,1]*p+M[1,2];          // expanding M[1,1]/M[1,2] + 1/p
550      M[1,2]=M[1,2]*p;
551      if (interval<>0)                 // canceling common terms of denominator
552      { if ((j/interval)*interval==j or j==g) // and enumerator -
553        { s=matrix(syz(ideal(M)));     // once gcd() is faster than syz() these
554          M[1,1]=-s[2,1];              // three lines should be replaced by the
555          M[1,2]=s[1,1];               // following three
556          // p=gcd(M[1,1],M[1,2]);
557          // M[1,1]=M[1,1]/p;
558          // M[1,2]=M[1,2]/p;
559        }
560      }
561      if (v)
562      { "  Term "+string(j)+" of the Molien series has been computed.";
563      }
564    }
565    if (interval==0)                   // canceling common terms of denominator
566    {                                  // and enumerator -
567      s=matrix(syz(ideal(M)));         // once gcd() is faster than syz() these
568      M[1,1]=-s[2,1];                  // three lines should be replaced by the
569      M[1,2]=s[1,1];                   // following three
570      // p=gcd(M[1,1],M[1,2]);
571      // M[1,1]=M[1,1]/p;
572      // M[1,2]=M[1,2]/p;
573    }
574    map slead=br,ideal(0);
575    s=slead(M);
576    M[1,1]=1/s[1,1]*M[1,1];            // numerator and denominator have to have
577    M[1,2]=1/s[1,2]*M[1,2];            // a constant term of 1
578    if (v)
579    { "";
580      "  We are done calculating the Molien series.";
581      "";
582    }
583    return(M);
584  }
585 //---- calculating Molien series in prime characteristic with Brauer lift ----
586  if (ch<>0 && mol_flag==0)
587  { if (g<>1)
588    { matrix G(1..g)=#[1..g];
589      if (interval<0)
590      { string Mstring;
591      }
592 //------ preparing everything for Brauer lifts into characteristic 0 ---------
593      ring Q=0,x,dp;                   // we want to extend our ring as well as
594                                       // the ring of rational numbers Q to
595                                       // contain r-th primitive roots of unity
596                                       // in order to factor characteristic
597                                       // polynomials of group elements into
598                                       // linear factors and lift eigenvalues to
599                                       // characteristic 0 -
600      poly minq=cyclotomic(r);         // minq now contains the size-of-group-th
601                                       // cyclotomic polynomial of Q, it is
602                                       // irreducible there
603      ring `newring`=(0,e),x,dp;
604      map f=Q,ideal(e);
605      minpoly=number(f(minq));         // e is now a r-th primitive root of
606                                       // unity -
607      kill Q, f;                       // no longer needed -
608      poly p=1;                        // used to build the denominator of the
609                                       // new term in the Molien series
610      matrix s[1][2];                  // used for canceling -
611      matrix M[1][2]=0,1;              // will contain Molien series -
612      ring v1br=char(br),x,dp;         // we calculate the r-th cyclotomic
613      poly minp=cyclotomic(r);         // polynomial of the base field and pick
614      minp=factorize(minp)[1][2];      // an irreducible factor of it -
615      if (deg(minp)==1)                // in this case the base field contains
616      { ring bre=char(br),x,dp;        // r-th roots of unity already
617        map f1=v1br,ideal(0);
618        number e=-number((f1(minp)));  // e is a r-th primitive root of unity
619      }
620      else
621      { ring bre=(char(br),e),x,dp;
622        map f1=v1br,ideal(e);
623        minpoly=number(f1(minp));      // e is a r-th primitive root of unity
624      }
625      map f2=br,ideal(0);              // we need f2 to map our group elements
626                                       // to this new extension field bre
627      matrix xI=diag(x,n);
628      poly p;                          // used for the characteristic polynomial
629                                       // to factor -
630      list L;                          // will contain the linear factors of the
631      ideal F;                         // characteristic polynomial of the group
632      intvec C;                        // elements and their powers
633      int i, j, k;
634 // -------------- finding all the terms of the Molien series -----------------
635      for (i=1;i<=g;i=i+1)
636      { setring br;
637        if (not(typeof(#[i])=="matrix"))
638        { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
639          return();
640        }
641        if ((n<>nrows(#[i])) or (n<>ncols(#[i])))
642        { "ERROR:   matrices need to be square and of the same dimensions";
643           return();
644        }
645        setring bre;
646        p=det(xI-f2(G(i)));            // characteristic polynomial of G(i)
647        L=factorize(p);
648        F=L[1];
649        C=L[2];
650        for (j=2;j<=ncols(F);j=j+1)
651        { F[j]=-1*(F[j]-x);            // F[j] is now an eigenvalue of G(i),
652                                       // it is a power of a primitive r-th root
653                                       // of unity -
654          k=exponent(number(F[j]),e);  // F[j]==e^k
655          setring `newring`;
656          p=p*(1-x*(e^k))^C[j];        // building the denominator of the new
657          setring bre;                 // term
658        }
659//         -----------
660//         k=0;
661//         while(k<r)
662//         { map f3=basering,ideal(e^k);
663//           while (f3(p)==0)
664//           { p=p/(x-e^k);
665//             setring `newring`;
666//             p=p*(1-x*(e^k));        // building the denominator of the new
667//             setring bre;
668//           }
669//           kill f3;
670//           if (p==1)
671//           { break;
672//           }
673//           k=k+1;
674//         }
675        setring `newring`;
676        M[1,1]=M[1,1]*p+M[1,2];        // expanding M[1,1]/M[1,2] + 1/p
677        M[1,2]=M[1,2]*p;
678        if (interval<0)
679        { if (i<>g)
680          { Mstring=string(M);
681            for (j=1;j<=size(Mstring);j=j+1)
682            { if (Mstring[j]=="e")
683              { interval=0;
684                break;
685              }
686            }
687          }
688          if (interval<>0)
689          { s=matrix(syz(ideal(M)));   // once gcd() is faster than syz()
690            M[1,1]=-s[2,1];            // these three lines should be
691            M[1,2]=s[1,1];             // replaced by the following three
692            // p=gcd(M[1,1],M[1,2]);
693            // M[1,1]=M[1,1]/p;
694            // M[1,2]=M[1,2]/p;
695          }
696          else
697          { interval=-1;
698          }
699        }
700        else
701        { if (interval<>0)             // canceling common terms of denominator
702          { if ((i/interval)*interval==i or i==g) // and enumerator
703            { s=matrix(syz(ideal(M))); // once gcd() is faster than syz()
704              M[1,1]=-s[2,1];          // these three lines should be
705              M[1,2]=s[1,1];           // replaced by the following three
706              // p=gcd(M[1,1],M[1,2]);
707              // M[1,1]=M[1,1]/p;
708              // M[1,2]=M[1,2]/p;
709            }
710          }
711        }
712        p=1;
713        setring bre;
714        if (v)
715        { "  Term "+string(i)+" of the Molien series has been computed.";
716        }
717      }
718      if (v)
719      { "";
720      }
721      setring `newring`;
722      if (interval==0)                 // canceling common terms of denominator
723      {                                // and enumerator -
724        s=matrix(syz(ideal(M)));       // once gcd() is faster than syz() these
725        M[1,1]=-s[2,1];                // three lines should be replaced by the
726        M[1,2]=s[1,1];                 // following three
727        // p=gcd(M[1,1],M[1,2]);
728        // M[1,1]=M[1,1]/p;
729        // M[1,2]=M[1,2]/p;
730      }
731      map slead=`newring`,ideal(0);
732      s=slead(M);                      // forcing the constant term of numerator
733      M[1,1]=1/s[1,1]*M[1,1];          // and denominator to be 1
734      M[1,2]=1/s[1,2]*M[1,2];
735      kill slead;
736      kill s;
737      kill p;
738    }
739    else                               // if the group only contains an identity
740    { ring `newring`=0,x,dp;           // element, it is very easy to calculate
741      matrix M[1][2]=1,(1-x)^n;        // the Molien series
742    }
743    export `newring`;                  // we keep the ring where we computed the
744    export M;                          // Molien series in such that we can
745    setring br;                        // keep it
746    if (v)
747    { "  We are done calculating the Molien series.";
748      "";
749    }
750  }
751  else                                 // i.e. char<>0 and mol_flag<>0, the user
752  {                                    // has specified that we are dealing with
753                                       // a ring of large characteristic which
754                                       // can be treated like a ring of
755                                       // characteristic 0; we'll avoid the
756                                       // Brauer lifts
757 //----------------------- simulating characteristic 0 ------------------------
758    string chst=charstr(br);
759    for (int i=1;i<=size(chst);i=i+1)
760    { if (chst[i]==",")
761      { break;
762      }
763    }
764 //----------------- generating ring of characteristic 0 ----------------------
765    if (minpoly==0)
766    { if (i>size(chst))
767      { execute "ring "+newring+"=0,("+varstr(br)+"),("+ordstr(br)+")";
768      }
769      else
770      { chst=chst[i..size(chst)];
771        execute "ring "+newring+"=(0"+chst+"),("+varstr(br)+"),("+ordstr(br)+")";
772      }
773    }
774    else
775    { string minp=string(minpoly);
776      minp=minp[2..size(minp)-1];
777      chst=chst[i..size(chst)];
778      execute "ring "+newring+"=(0"+chst+"),("+varstr(br)+"),("+ordstr(br)+")";
779      execute "minpoly="+minp;
780    }
781    matrix I=diag(1,n);
782    poly v1=maxideal(1)[1];            // the Molien series will be in terms of
783                                       // the first variable of the current
784                                       // ring -
785    matrix M[1][2];                    // M will contain the Molien series -
786    M[1,1]=0;                          // M[1,1] will be the numerator -
787    M[1,2]=1;                          // M[1,2] will be the denominator -
788    matrix s;                          // will help us canceling in the
789                                       // fraction
790    poly p;                            // will contain the denominator of the
791                                       // new term of the Molien series
792    int j;
793    string links, rechts;
794 //----------------- finding all terms of the Molien series -------------------
795    for (i=1;i<=g;i=i+1)
796    { setring br;
797      if (not(typeof(#[i])=="matrix"))
798      { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
799        return();
800      }
801      if ((n<>nrows(#[i])) or (n<>ncols(#[i])))
802      { "ERROR:   matrices need to be square and of the same dimensions";
803         return();
804      }
805      string stM(i)=string(#[i]);
806      for (j=1;j<=size(stM(i));j=j+1)
807      { if (stM(i)[j]=="
808")
809        { links=stM(i)[1..j-1];
810          rechts=stM(i)[j+1..size(stM(i))];
811          stM(i)=links+rechts;
812        }
813      }
814      setring `newring`;
815      execute "matrix G(i)["+string(n)+"]["+string(n)+"]="+stM(i);
816      p=det(I-v1*G(i));                // denominator of new term -
817      M[1,1]=M[1,1]*p+M[1,2];          // expanding M[1,1]/M[1,2] + 1/p
818      M[1,2]=M[1,2]*p;
819      if (interval<>0)                 // canceling common terms of denominator
820      { if ((i/interval)*interval==i or i==g) // and enumerator
821        {
822          s=matrix(syz(ideal(M)));     // once gcd() is faster than syz() these
823          M[1,1]=-s[2,1];              // three lines should be replaced by the
824          M[1,2]=s[1,1];               // following three
825          // p=gcd(M[1,1],M[1,2]);
826          // M[1,1]=M[1,1]/p;
827          // M[1,2]=M[1,2]/p;
828        }
829      }
830      if (v)
831      { "  Term "+string(i)+" of the Molien series has been computed.";
832      }
833    }
834    if (interval==0)                   // canceling common terms of denominator
835    {                                  // and enumerator -
836      s=matrix(syz(ideal(M)));         // once gcd() is faster than syz() these
837      M[1,1]=-s[2,1];                  // three lines should be replaced by the
838      M[1,2]=s[1,1];                   // following three
839      // p=gcd(M[1,1],M[1,2]);
840      // M[1,1]=M[1,1]/p;
841      // M[1,2]=M[1,2]/p;
842    }
843    map slead=`newring`,ideal(0);
844    s=slead(M);
845    M[1,1]=1/s[1,1]*M[1,1];            // numerator and denominator have to have
846    M[1,2]=1/s[1,2]*M[1,2];            // a constant term of 1
847    if (v)
848    { "";
849      "  We are done calculating the Molien series.";
850      "";
851    }
852    kill G(1..g), s, slead, p, v1, I;
853    export `newring`;                  // we keep the ring where we computed the
854    export M;                          // the Molien series such that we can
855    setring br;                        // keep it
856  }
857}
858example
859{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
860  "         note the case of prime characteristic";
861  echo=2;
862         ring R=0,(x,y,z),dp;
863         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
864         list L=group_reynolds(A);
865         matrix M=molien(L[2..size(L)]);
866         print(M);
867         ring S=3,(x,y,z),dp;
868         string newring="alksdfjlaskdjf";
869         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
870         list L=group_reynolds(A);
871         molien(L[2..size(L)],newring);
872         setring alksdfjlaskdjf;
873         print(M);
874         setring S;
875         kill alksdfjlaskdjf;
876}
877
878proc reynolds_molien (list #)
879"USAGE:   reynolds_molien(G1,G2,...[,ringname,flags]);
880         G1,G2,...: nxn <matrices> generating a finite matrix group, ringname:
881         a <string> giving a name for a new ring of characteristic 0 for the
882         Molien series in case of prime characteristic, flags: an optional
883         <intvec> with three components: if the first element is not equal to 0
884         characteristic 0 is simulated, i.e. the Molien series is computed as
885         if the base field were characteristic 0 (the user must choose a field
886         of large prime characteristic, e.g. 32003) the second component should
887         give the size of intervals between canceling common factors in the
888         expansion of the Molien series, 0 (the default) means only once after
889         generating all terms, in prime characteristic also a negative number
890         can be given to indicate that common factors should always be canceled
891         when the expansion is simple (the root of the extension field does not
892         occur among the coefficients)
893ASSUME:  n is the number of variables of the basering, G1,G2... are the group
894         elements generated by group_reynolds(), g is the size of the group
895RETURN:  a gxn <matrix> representing the Reynolds operator is the first return
896         value and in case of characteristic 0 a 1x2 <matrix> giving enumerator
897         and denominator of Molien series is the second one; in case of prime
898         characteristic a ring with the name `ringname` of characteristic 0 is
899         created where the same Molien series (named M) is stored
900DISPLAY: information if the third component of flags does not equal 0
901EXAMPLE: example reynolds_molien; shows an example
902THEORY:  The entire matrix group is generated by getting all left products of
903         the generators with the new elements from the last run through the loop
904         (or the generators themselves during the first run). All the ones that
905         have been generated before are thrown out and the program terminates
906         when there are no new elements found in one run. Additionally each time
907         a new group element is found the corresponding ring mapping of which
908         the Reynolds operator is made up is generated. They are stored in the
909         rows of the first return value. In characteristic 0 the terms
910         1/det(1-xE) is computed whenever a new element E is found. In prime
911         characteristic a Brauer lift is involved and the terms are only
912         computed after the entire matrix group is generated (to avoid the
913         modular case). The returned matrix gives enumerator and denominator of
914         the expanded version where common factors have been canceled.
915"
916{ def br=basering;                     // the Molien series depends on the
917  int ch=char(br);                     // characteristic of the coefficient
918                                       // field
919  int gen_num;
920 //------------------- making sure the input is okay --------------------------
921  if (typeof(#[size(#)])=="intvec")
922  { if (size(#[size(#)])==3)
923    { int mol_flag=#[size(#)][1];
924      if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0)))
925      { "ERROR:   the second component of the <intvec> should be >=0";
926        return();
927      }
928      int interval=#[size(#)][2];
929      int v=#[size(#)][3];
930    }
931    else
932    { "ERROR:   <intvec> should have three components";
933      return();
934    }
935    if (ch<>0)
936    { if (typeof(#[size(#)-1])<>"string")
937      { "ERROR:   in characteristic p a <string> must be given for the name";
938        "         of a new ring where the Molien series can be stored";
939        return();
940      }
941      else
942      { if (#[size(#)-1]=="")
943        { "ERROR:   <string> may not be empty";
944          return();
945        }
946        string newring=#[size(#)-1];
947        gen_num=size(#)-2;
948      }
949    }
950    else                               // then <string> ist not needed
951    { gen_num=size(#)-1;
952    }
953  }
954  else                                 // last parameter is not <intvec>
955  { int v=0;                           // no information is default
956    int interval;
957    int mol_flag=0;                    // computing of Molien series is default
958    if (ch<>0)
959    { if (typeof(#[size(#)])<>"string")
960      { "ERROR:   in characteristic p a <string> must be given for the name";
961        "         of a new ring where the Molien series can be stored";
962        return();
963      }
964      else
965      { if (#[size(#)]=="")
966        { "ERROR:   <string> may not be empty";
967          return();
968        }
969        string newring=#[size(#)];
970        gen_num=size(#)-1;
971      }
972    }
973    else
974    { gen_num=size(#);
975    }
976  }
977 // ----------------- computing the terms with Brauer lift --------------------
978  if (ch<>0 && mol_flag==0)
979  { list L=group_reynolds(#[1..gen_num],v);
980    if (L[1]==0)
981    { if (voice==2)
982      { "WARNING: The characteristic of the coefficient field divides the group order.";
983        "         Proceed without the Reynolds operator or the Molien series!";
984        return();
985      }
986      if (v)
987      { "  The characteristic of the base field divides the group order.";
988        "  We have to continue without Reynolds operator or the Molien series...";
989        return();
990      }
991    }
992    if (minpoly<>0)
993    { if (voice==2)
994      { "WARNING: It is impossible for this program to calculate the Molien series";
995        "         for finite groups over extension fields of prime characteristic.";
996        return(L[1]);
997      }
998      else
999      { if (v)
1000        { "  Since it is impossible for this program to calculate the Molien series for";
1001          "  invariant rings over extension fields of prime characteristic, we have to";
1002          "  continue without it.";
1003          return(L[1]);
1004        }
1005      }
1006    }
1007    if (typeof(L[2])=="int")
1008    { molien(L[3..size(L)],newring,L[2],intvec(mol_flag,interval,v));
1009    }
1010    else
1011    { molien(L[2..size(L)],newring,intvec(mol_flag,interval,v));
1012    }
1013    return(L[1]);
1014  }
1015 //----------- computing Molien series in the straight forward way ------------
1016  if (ch==0)
1017  { if (typeof(#[1])<>"matrix")
1018    { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
1019      return();
1020    }
1021    int n=nrows(#[1]);
1022    if (n<>nvars(br))
1023    { "ERROR:   the number of variables of the basering needs to be the same";
1024      "         as the dimension of the matrices";
1025      return();
1026    }
1027    if (n<>ncols(#[1]))
1028    { "ERROR:   matrices need to be square and of the same dimensions";
1029      return();
1030    }
1031    matrix vars=matrix(maxideal(1));   // creating an nx1-matrix containing the
1032    vars=transpose(vars);              // variables of the ring -
1033    matrix A(1)=#[1]*vars;             // calculating the first ring mapping -
1034                                       // A(1) will contain the Reynolds
1035                                       // operator -
1036    poly v1=vars[1,1];                 // the Molien series will be in terms of
1037                                       // the first variable of the current
1038                                       // ring
1039    matrix I=diag(1,n);
1040    matrix A(2)[1][2];                 // A(2) will contain the Molien series -
1041    A(2)[1,1]=1;                       // A(2)[1,1] will be the numerator
1042    matrix G(1)=#[1];                  // G(k) are elements of the group -
1043    A(2)[1,2]=det(I-v1*(G(1)));        // A(2)[1,2] will be the denominator -
1044    matrix s;                          // will help us canceling in the
1045                                       // fraction
1046    poly p;                            // will contain the denominator of the
1047                                       // new term of the Molien series
1048    int i=1;
1049    for (int j=2;j<=gen_num;j=j+1)     // this loop adds the parameters to the
1050    {                                  // group, leaving out doubles and
1051                                       // checking whether the parameters are
1052                                       // compatible with the task of the
1053                                       // procedure
1054      if (not(typeof(#[j])=="matrix"))
1055      { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
1056        return();
1057      }
1058      if ((n!=nrows(#[j])) or (n!=ncols(#[j])))
1059      { "ERROR:   matrices need to be square and of the same dimensions";
1060         return();
1061      }
1062      if (unique(G(1..i),#[j]))
1063      { i=i+1;
1064        matrix G(i)=#[j];
1065        A(1)=concat(A(1),#[j]*vars);   // adding ring homomorphisms to A(1) -
1066        p=det(I-v1*#[j]);              // denominator of new term -
1067        A(2)[1,1]=A(2)[1,1]*p+A(2)[1,2]; // expanding A(2)[1,1]/A(2)[1,2] +1/p
1068        A(2)[1,2]=A(2)[1,2]*p;
1069        if (interval<>0)               // canceling common terms of denominator
1070        { if ((i/interval)*interval==i) // and enumerator
1071          {
1072            s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz() these
1073            A(2)[1,1]=-s[2,1];         // three lines should be replaced by the
1074            A(2)[1,2]=s[1,1];          // following three
1075            // p=gcd(A(2)[1,1],A(2)[1,2]);
1076            // A(2)[1,1]=A(2)[1,1]/p;
1077            // A(2)[1,2]=A(2)[1,2]/p;
1078          }
1079        }
1080      }
1081    }
1082    int g=i;                           // G(1)..G(i) are generators without
1083                                       // doubles - g generally is the number
1084                                       // of elements in the group so far -
1085    j=i;                               // j is the number of new elements that
1086                                       // we use as factors
1087    int k, m, l;
1088    if (v)
1089    { "";
1090      "  Generating the entire matrix group. Whenever a new group element is found,";
1091      "  the coressponding ring homomorphism of the Reynolds operator and the";
1092      "  corresponding term of the Molien series is generated.";
1093      "";
1094    }
1095 //----------- computing 1/det(I-xE) whenever a new element E is found --------
1096    while (1)
1097    { l=0;                             // l is the number of products we get in
1098                                       // one going
1099      for (m=g-j+1;m<=g;m=m+1)
1100      { for (k=1;k<=i;k=k+1)
1101        { l=l+1;
1102          matrix P(l)=G(k)*G(m);       // possible new element
1103        }
1104      }
1105      j=0;
1106      for (k=1;k<=l;k=k+1)
1107      { if (unique(G(1..g),P(k)))
1108        { j=j+1;                       // a new factor for next run
1109          g=g+1;
1110          matrix G(g)=P(k);            // a new group element -
1111          A(1)=concat(A(1),P(k)*vars); // adding new mapping to A(1)
1112          p=det(I-v1*P(k));            // denominator of new term
1113          A(2)[1,1]=A(2)[1,1]*p+A(2)[1,2];
1114          A(2)[1,2]=A(2)[1,2]*p;       // expanding A(2)[1,1]/A(2)[1,2] + 1/p -
1115          if (interval<>0)             // canceling common terms of denominator
1116          { if ((g/interval)*interval==g) // and enumerator
1117            {
1118              s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz()
1119              A(2)[1,1]=-s[2,1];       // these three lines should be replaced
1120              A(2)[1,2]=s[1,1];        // by the following three
1121              // p=gcd(A(2)[1,1],A(2)[1,2]);
1122              // A(2)[1,1]=A(2)[1,1]/p;
1123              // A(2)[1,2]=A(2)[1,2]/p;
1124            }
1125          }
1126          if (v)
1127          { "  Group element "+string(g)+" has been found.";
1128          }
1129        }
1130        kill P(k);
1131      }
1132      if (j==0)                        // when we didn't add any new elements
1133      { break;                         // in one run through the while loop
1134      }                                // we are done
1135    }
1136    if (v)
1137    { if (g<=i)
1138      { "  There are only "+string(g)+" group elements.";
1139      }
1140      "";
1141    }
1142    A(1)=transpose(A(1));              // when we evaluate the Reynolds operator
1143                                       // later on, we actually want 1xn
1144                                       // matrices
1145    if (interval==0)                   // canceling common terms of denominator
1146    {                                  // and enumerator -
1147      s=matrix(syz(ideal(A(2))));      // once gcd() is faster than syz()
1148      A(2)[1,1]=-s[2,1];               // these three lines should be replaced
1149      A(2)[1,2]=s[1,1];                // by the following three
1150      // p=gcd(A(2)[1,1],A(2)[1,2]);
1151      // A(2)[1,1]=A(2)[1,1]/p;
1152      // A(2)[1,2]=A(2)[1,2]/p;
1153    }
1154    if (interval<>0)                   // canceling common terms of denominator
1155    { if ((g/interval)*interval<>g)    // and enumerator
1156      {
1157        s=matrix(syz(ideal(A(2))));    // once gcd() is faster than syz()
1158        A(2)[1,1]=-s[2,1];             // these three lines should be replaced
1159        A(2)[1,2]=s[1,1];              // by the following three
1160        // p=gcd(A(2)[1,1],A(2)[1,2]);
1161        // A(2)[1,1]=A(2)[1,1]/p;
1162        // A(2)[1,2]=A(2)[1,2]/p;
1163      }
1164    }
1165    map slead=br,ideal(0);
1166    s=slead(A(2));
1167    A(2)[1,1]=1/s[1,1]*A(2)[1,1];      // numerator and denominator have to have
1168    A(2)[1,2]=1/s[1,2]*A(2)[1,2];      // a constant term of 1
1169    if (v)
1170    { "  Now we are done calculating Molien series and Reynolds operator.";
1171      "";
1172    }
1173    return(A(1..2));
1174  }
1175 //------------------------ simulating characteristic 0 -----------------------
1176  else                                 // if ch<>0 and mol_flag<>0
1177  { if (typeof(#[1])<>"matrix")
1178    { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
1179      return();
1180    }
1181    int n=nrows(#[1]);
1182    if (n<>nvars(br))
1183    { "ERROR:   the number of variables of the basering needs to be the same";
1184      "         as the dimension of the matrices";
1185      return();
1186    }
1187    if (n<>ncols(#[1]))
1188    { "ERROR:   matrices need to be square and of the same dimensions";
1189      return();
1190    }
1191    matrix vars=matrix(maxideal(1));   // creating an nx1-matrix containing the
1192    vars=transpose(vars);              // variables of the ring -
1193    matrix A(1)=#[1]*vars;             // calculating the first ring mapping -
1194                                       // A(1) will contain the Reynolds
1195                                       // operator
1196    string chst=charstr(br);
1197    for (int i=1;i<=size(chst);i=i+1)
1198    { if (chst[i]==",")
1199      { break;
1200      }
1201    }
1202    if (minpoly==0)
1203    { if (i>size(chst))
1204      { execute "ring "+newring+"=0,("+varstr(br)+"),("+ordstr(br)+")";
1205      }
1206      else
1207      { chst=chst[i..size(chst)];
1208        execute "ring "+newring+"=(0"+chst+"),("+varstr(br)+"),("+ordstr(br)+")";
1209      }
1210    }
1211    else
1212    { string minp=string(minpoly);
1213      minp=minp[2..size(minp)-1];
1214      chst=chst[i..size(chst)];
1215      execute "ring "+newring+"=(0"+chst+"),("+varstr(br)+"),("+ordstr(br)+")";
1216      execute "minpoly="+minp;
1217    }
1218    poly v1=var(1);                    // the Molien series will be in terms of
1219                                       // the first variable of the current
1220                                       // ring
1221    matrix I=diag(1,n);
1222    int o;
1223    setring br;
1224    matrix G(1)=#[1];
1225    string links, rechts;
1226    string stM(1)=string(#[1]);
1227    for (o=1;o<=size(stM(1));o=o+1)
1228    { if (stM(1)[o]=="
1229")
1230      { links=stM(1)[1..o-1];
1231        rechts=stM(1)[o+1..size(stM(1))];
1232        stM(1)=links+rechts;
1233      }
1234    }
1235    setring `newring`;
1236    execute "matrix G(1)["+string(n)+"]["+string(n)+"]="+stM(1);
1237    matrix A(2)[1][2];                 // A(2) will contain the Molien series -
1238    A(2)[1,1]=1;                       // A(2)[1,1] will be the numerator
1239    A(2)[1,2]=det(I-v1*(G(1)));        // A(2)[1,2] will be the denominator -
1240    matrix s;                          // will help us canceling in the
1241                                       // fraction
1242    poly p;                            // will contain the denominator of the
1243                                       // new term of the Molien series
1244    i=1;
1245    for (int j=2;j<=gen_num;j=j+1)     // this loop adds the parameters to the
1246    {                                  // group, leaving out doubles and
1247                                       // checking whether the parameters are
1248                                       // compatible with the task of the
1249                                       // procedure
1250      setring br;
1251      if (not(typeof(#[j])=="matrix"))
1252      { "ERROR:   the parameters must be a list of matrices and maybe an <intvec>";
1253        return();
1254      }
1255      if ((n!=nrows(#[j])) or (n!=ncols(#[j])))
1256      { "ERROR:   matrices need to be square and of the same dimensions";
1257         return();
1258      }
1259      if (unique(G(1..i),#[j]))
1260      { i=i+1;
1261        matrix G(i)=#[j];
1262        A(1)=concat(A(1),G(i)*vars);   // adding ring homomorphisms to A(1)
1263        string stM(i)=string(G(i));
1264        for (o=1;o<=size(stM(i));o=o+1)
1265        { if (stM(i)[o]=="
1266")
1267          { links=stM(i)[1..o-1];
1268            rechts=stM(i)[o+1..size(stM(i))];
1269            stM(i)=links+rechts;
1270          }
1271        }
1272        setring `newring`;
1273        execute "matrix G(i)["+string(n)+"]["+string(n)+"]="+stM(i);
1274        p=det(I-v1*G(i));              // denominator of new term -
1275        A(2)[1,1]=A(2)[1,1]*p+A(2)[1,2]; // expanding A(2)[1,1]/A(2)[1,2] +1/p
1276        A(2)[1,2]=A(2)[1,2]*p;
1277        if (interval<>0)               // canceling common terms of denominator
1278        { if ((i/interval)*interval==i) // and enumerator
1279          {
1280            s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz() these
1281            A(2)[1,1]=-s[2,1];         // three lines should be replaced by the
1282            A(2)[1,2]=s[1,1];          // following three
1283            // p=gcd(A(2)[1,1],A(2)[1,2]);
1284            // A(2)[1,1]=A(2)[1,1]/p;
1285            // A(2)[1,2]=A(2)[1,2]/p;
1286          }
1287        }
1288        setring br;
1289      }
1290    }
1291    int g=i;                           // G(1)..G(i) are generators without
1292                                       // doubles - g generally is the number
1293                                       // of elements in the group so far -
1294    j=i;                               // j is the number of new elements that
1295                                       // we use as factors
1296    int k, m, l;
1297    if (v)
1298    { "";
1299      "  Generating the entire matrix group. Whenever a new group element is found,";
1300      "  the coressponding ring homomorphism of the Reynolds operator and the";
1301      "  corresponding term of the Molien series is generated.";
1302      "";
1303    }
1304 // taking all elements in a ring of characteristic 0 and computing the terms
1305 // of the Molien series there
1306    while (1)
1307    { l=0;                             // l is the number of products we get in
1308                                       // one going
1309      for (m=g-j+1;m<=g;m=m+1)
1310      { for (k=1;k<=i;k=k+1)
1311        { l=l+1;
1312          matrix P(l)=G(k)*G(m);       // possible new element
1313        }
1314      }
1315      j=0;
1316      for (k=1;k<=l;k=k+1)
1317      { if (unique(G(1..g),P(k)))
1318        { j=j+1;                       // a new factor for next run
1319          g=g+1;
1320          matrix G(g)=P(k);            // a new group element -
1321          A(1)=concat(A(1),G(g)*vars); // adding new mapping to A(1)
1322          string stM(g)=string(G(g));
1323          for (o=1;o<=size(stM(g));o=o+1)
1324          { if (stM(g)[o]=="
1325")
1326            { links=stM(g)[1..o-1];
1327              rechts=stM(g)[o+1..size(stM(g))];
1328              stM(g)=links+rechts;
1329            }
1330          }
1331          setring `newring`;
1332          execute "matrix G(g)["+string(n)+"]["+string(n)+"]="+stM(g);
1333          p=det(I-v1*G(g));            // denominator of new term
1334          A(2)[1,1]=A(2)[1,1]*p+A(2)[1,2];
1335          A(2)[1,2]=A(2)[1,2]*p;       // expanding A(2)[1,1]/A(2)[1,2] + 1/p -
1336          if (interval<>0)             // canceling common terms of denominator
1337          { if ((g/interval)*interval==g) // and enumerator
1338            {
1339              s=matrix(syz(ideal(A(2)))); // once gcd() is faster than syz()
1340              A(2)[1,1]=-s[2,1];       // these three lines should be replaced
1341              A(2)[1,2]=s[1,1];        // by the following three
1342              // p=gcd(A(2)[1,1],A(2)[1,2]);
1343              // A(2)[1,1]=A(2)[1,1]/p;
1344              // A(2)[1,2]=A(2)[1,2]/p;
1345            }
1346          }
1347          if (v)
1348          { "  Group element "+string(g)+" has been found.";
1349          }
1350          setring br;
1351        }
1352        kill P(k);
1353      }
1354      if (j==0)                        // when we didn't add any new elements
1355      { break;                         // in one run through the while loop
1356      }                                // we are done
1357    }
1358    if (v)
1359    { if (g<=i)
1360      { "  There are only "+string(g)+" group elements.";
1361      }
1362      "";
1363    }
1364    A(1)=transpose(A(1));              // when we evaluate the Reynolds operator
1365                                       // later on, we actually want 1xn
1366                                       // matrices
1367    setring `newring`;
1368    if (interval==0)                   // canceling common terms of denominator
1369    {                                  // and enumerator -
1370      s=matrix(syz(ideal(A(2))));      // once gcd() is faster than syz()
1371      A(2)[1,1]=-s[2,1];               // these three lines should be replaced
1372      A(2)[1,2]=s[1,1];                // by the following three
1373      // p=gcd(A(2)[1,1],A(2)[1,2]);
1374      // A(2)[1,1]=A(2)[1,1]/p;
1375      // A(2)[1,2]=A(2)[1,2]/p;
1376    }
1377    if (interval<>0)                   // canceling common terms of denominator
1378    { if ((g/interval)*interval<>g)    // and enumerator
1379      {
1380        s=matrix(syz(ideal(A(2))));    // once gcd() is faster than syz()
1381        A(2)[1,1]=-s[2,1];             // these three lines should be replaced
1382        A(2)[1,2]=s[1,1];              // by the following three
1383        // p=gcd(A(2)[1,1],A(2)[1,2]);
1384        // A(2)[1,1]=A(2)[1,1]/p;
1385        // A(2)[1,2]=A(2)[1,2]/p;
1386      }
1387    }
1388    map slead=`newring`,ideal(0);
1389    s=slead(A(2));
1390    A(2)[1,1]=1/s[1,1]*A(2)[1,1];      // numerator and denominator have to have
1391    A(2)[1,2]=1/s[1,2]*A(2)[1,2];      // a constant term of 1
1392    if (v)
1393    { "  Now we are done calculating Molien series and Reynolds operator.";
1394      "";
1395    }
1396    matrix M=A(2);
1397    kill G(1..g), s, slead, p, v1, I, A(2);
1398    export `newring`;                  // we keep the ring where we computed the
1399    export M;                          // the Molien series such that we can
1400    setring br;                        // keep it
1401    return(A(1));
1402  }
1403}
1404example
1405{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
1406  "         note the case of prime characteristic";
1407  echo=2;
1408         ring R=0,(x,y,z),dp;
1409         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
1410         matrix REY,M=reynolds_molien(A);
1411         print(REY);
1412         print(M);
1413         ring S=3,(x,y,z),dp;
1414         string newring="Qadjoint";
1415         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
1416         matrix REY=reynolds_molien(A,newring);
1417         print(REY);
1418         setring Qadjoint;
1419         print(M);
1420         setring S;
1421         kill Qadjoint;
1422}
1423
1424proc partial_molien (matrix M, int n, list #)
1425"USAGE:   partial_molien(M,n[,p]);
1426         M: a 1x2 <matrix>, n: an <int> indicating  number of terms in the
1427         expansion, p: an optional <poly>
1428ASSUME:  M is the return value of molien or the second return value of
1429         reynolds_molien, p ought to be the second return value of a previous
1430         run of partial_molien and avoids recalculating known terms
1431RETURN:  n terms (type <poly>) of the partial expansion of the Molien series
1432         (first n if there is no third parameter given, otherwise the next n
1433         terms depending on a previous calculation) and an intermediate result
1434         (type <poly>) of the calculation to be used as third parameter in a next
1435         run of partial_molien
1436THEORY:  The following calculation is implemented:
1437         (1+a1x+a2x^2+...+anx^n)/(1+b1x+b2x^2+...+bmx^m)=(1+(a1-b1)x+...
1438         (1+b1x+b2x^2+...+bmx^m)
1439         -----------------------
1440            (a1-b1)x+(a2-b2)x^2+...
1441            (a1-b1)x+b1(a1-b1)x^2+...
1442EXAMPLE: example partial_molien; shows an example
1443"
1444{ poly A(2);                           // A(2) will contain the return value of
1445                                       // the intermediate result
1446  if (char(basering)<>0)
1447  { "ERROR:   you have to change to a basering of characteristic 0, one in";
1448    "         which the Molien series is defined";
1449  }
1450  if (ncols(M)==2 && nrows(M)==1 && n>0 && size(#)<2)
1451  { def br=basering;                   // keeping track of the old ring
1452    map slead=br,ideal(0);
1453    matrix s=slead(M);
1454    if (s[1,1]<>1 || s[1,2]<>1)
1455    { "ERROR:   the constant terms of enumerator and denominator are not 1";
1456      return();
1457    }
1458
1459    if (size(#)==0)
1460    { A(2)=M[1,1];                     // if a third parameter is not given, the
1461                                       // intermediate result from the last run
1462                                       // corresponds to the numerator - we need
1463    }                                  // its smallest term
1464    else
1465    { if (typeof(#[1])=="poly")
1466      { A(2)=#[1];                     // if a third term is given we 'start'
1467      }                                // with its smallest term
1468      else
1469      { "ERROR:   <poly> as third parameter expected";
1470        return();
1471      }
1472    }
1473    poly A(1)=M[1,2];                  // denominator of Molien series (for now)
1474    string mp=string(minpoly);
1475    execute "ring R=("+charstr(br)+"),("+varstr(br)+"),ds;";
1476    execute "minpoly=number("+mp+");";
1477    poly A(1)=0;                       // A(1) will contain the sum of n terms -
1478    poly min;                          // min will be our smallest term -
1479    poly A(2)=fetch(br,A(2));          // fetching A(2) and M[1,2] into R
1480    poly den=fetch(br,A(1));
1481    for (int i=1; i<=n; i=i+1)         // getting n terms and adding them up
1482    { min=lead(A(2));
1483      A(1)=A(1)+min;
1484      A(2)=A(2)-min*den;
1485    }
1486    setring br;                        // moving A(1) and A(2) back in the
1487    A(1)=fetch(R,A(1));                // actual ring for output
1488    A(2)=fetch(R,A(2));
1489    return(A(1..2));
1490  }
1491  else
1492  { "ERROR:   the first parameter has to be a 1x2-matrix, i.e. the matrix";
1493    "         returned by the procedure 'reynolds_molien', the second one";
1494    "         should be > 0 and there should be no more than 3 parameters;"
1495    return();
1496  }
1497}
1498example
1499{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
1500  echo=2;
1501         ring R=0,(x,y,z),dp;
1502         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
1503         matrix REY,M=reynolds_molien(A);
1504         poly p(1..2);
1505         p(1..2)=partial_molien(M,5);
1506         p(1);
1507         p(1..2)=partial_molien(M,5,p(2));
1508         p(1);
1509}
1510
1511proc evaluate_reynolds (matrix REY, ideal I)
1512"USAGE:   evaluate_reynolds(REY,I);
1513         REY: a <matrix> representing the Reynolds operator, I: an arbitrary
1514         <ideal>
1515ASSUME:  REY is the first return value of group_reynolds() or reynolds_molien()
1516RETURNS: image of the polynomials defining I under the Reynolds operator
1517         (type <ideal>)
1518NOTE:    the characteristic of the coefficient field of the polynomial ring
1519         should not divide the order of the finite matrix group
1520EXAMPLE: example evaluate_reynolds; shows an example
1521THEORY:  REY has been constructed in such a way that each row serves as a ring
1522         mapping of which the Reynolds operator is made up.
1523"
1524{ def br=basering;
1525  int n=nvars(br);
1526  if (ncols(REY)==n)
1527  { int m=nrows(REY);                  // we need m to 'cut' the ring
1528                                       // homomorphisms 'out' of REY and to
1529                                       // divide by the group order in the end
1530    int num_poly=ncols(I);
1531    matrix MI=matrix(I);
1532    matrix MiI[1][num_poly];
1533    map pREY;
1534    matrix rowREY[1][n];
1535    for (int i=1;i<=m;i=i+1)
1536    { rowREY=REY[i,1..n];
1537      pREY=br,ideal(rowREY);           // f is now the i-th ring homomorphism
1538      MiI=pREY(MI)+MiI;
1539    }
1540    MiI=(1/number(m))*MiI;
1541    return(ideal(MiI));
1542  }
1543  else
1544  { "ERROR:   the number of columns in the <matrix> should be the same as the";
1545    "         number of variables in the basering; in fact it should be first";
1546    "         return value of group_reynolds() or reynolds_molien().";
1547    return();
1548  }
1549}
1550example
1551{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
1552  echo=2;
1553         ring R=0,(x,y,z),dp;
1554         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
1555         list L=group_reynolds(A);
1556         ideal I=x2,y2,z2;
1557         print(evaluate_reynolds(L[1],I));
1558}
1559
1560proc invariant_basis (int g,list #)
1561"USAGE:   invariant_basis(g,G1,G2,...);
1562         g: an <int> indicating of which degree (>0) the homogeneous basis
1563         shoud be, G1,G2,...: <matrices> generating a finite matrix group
1564RETURNS: the basis (type <ideal>) of the space of invariants of degree g
1565EXAMPLE: example invariant_basis; shows an example
1566THEORY:  A general polynomial of degree g is generated and the generators of the
1567         matrix group applied. The difference ought to be 0 and this way a
1568         system of linear equations is created. It is solved by computing
1569         syzygies.
1570"
1571{ if (g<=0)
1572  { "ERROR:   the first parameter should be > 0";
1573    return();
1574  }
1575  def br=basering;
1576  ideal mon=sort(maxideal(g))[1];      // needed for constructing a general
1577  int m=ncols(mon);                    // homogeneous polynomial of degree g
1578  mon=sort(mon,intvec(m..1))[1];
1579  int a=size(#);
1580  int i;
1581  int n=nvars(br);
1582 //---------------------- checking that the input is ok -----------------------
1583  for (i=1;i<=a;i=i+1)
1584  { if (typeof(#[i])=="matrix")
1585    { if (nrows(#[i])==n && ncols(#[i])==n)
1586      { matrix G(i)=#[i];
1587      }
1588      else
1589      { "ERROR:   the number of variables of the base ring needs to be the same";
1590        "         as the dimension of the square matrices";
1591        return();
1592      }
1593    }
1594    else
1595    { "ERROR:   the last parameters should be a list of matrices";
1596      return();
1597    }
1598  }
1599 //----------------------------------------------------------------------------
1600  execute "ring T=("+charstr(br)+"),("+varstr(br)+",p(1..m)),lp;";
1601  // p(1..m) are the general coefficients of the general polynomial of degree g
1602  execute "ideal vars="+varstr(br)+";";
1603  map f;
1604  ideal mon=imap(br,mon);
1605  poly P=0;
1606  for (i=m;i>=1;i=i-1)
1607  { P=P+p(i)*mon[i];                   // P is the general polynomial
1608  }
1609  ideal I;                             // will help substituting variables in P
1610                                       // by linear combinations of variables -
1611  poly Pnew,temp;                      // Pnew is P with substitutions -
1612  matrix S[m*a][m];                    // will contain system of linear
1613                                       // equations
1614  int j,k;
1615 //------------------- building the system of linear equations ----------------
1616  for (i=1;i<=a;i=i+1)
1617  { I=ideal(matrix(vars)*transpose(imap(br,G(i))));
1618    I=I,p(1..m);
1619    f=T,I;
1620    Pnew=f(P);
1621    for (j=1;j<=m;j=j+1)
1622    { temp=P/mon[j]-Pnew/mon[j];
1623      for (k=1;k<=m;k=k+1)
1624      { S[m*(i-1)+j,k]=temp/p(k);
1625      }
1626    }
1627  }
1628 //----------------------------------------------------------------------------
1629  setring br;
1630  map f=T,ideal(0);
1631  matrix S=f(S);
1632  matrix s=matrix(syz(S));             // s contains a basis of the space of
1633                                       // solutions -
1634  ideal I=ideal(matrix(mon)*s);        // I contains a basis of homogeneous
1635  if (I[1]<>0)                         // invariants of degree d
1636  { for (i=1;i<=ncols(I);i=i+1)
1637    { I[i]=I[i]/leadcoef(I[i]);        // setting leading coefficients to 1
1638    }
1639  }
1640  return(I);
1641}
1642example
1643{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
1644  echo=2;
1645           ring R=0,(x,y,z),dp;
1646           matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
1647           print(invariant_basis(2,A));
1648}
1649
1650proc invariant_basis_reynolds (matrix REY,int d,list #)
1651"USAGE:   invariant_basis_reynolds(REY,d[,flags]);
1652         REY: a <matrix> representing the Reynolds operator, d: an <int>
1653         indicating of which degree (>0) the homogeneous basis shoud be, flags:
1654         an optional <intvec> with two entries: its first component gives the
1655         dimension of the space (default <0 meaning unknown) and its second
1656         component is used as the number of polynomials that should be mapped
1657         to invariants during one call of evaluate_reynolds if the dimension of
1658         the space is unknown or the number such that number x dimension
1659         polynomials are mapped to invariants during one call of
1660         evaluate_reynolds
1661ASSUME:  REY is the first return value of group_reynolds() or reynolds_molien()
1662         and flags[1] given by partial_molien
1663RETURN:  the basis (type <ideal>) of the space of invariants of degree d
1664EXAMPLE: example invariant_basis_reynolds; shows an example
1665THEORY:  Monomials of degree d are mapped to invariants with the Reynolds
1666         operator. A linearly independent set is generated with the help of
1667         minbase.
1668"
1669{
1670 //---------------------- checking that the input is ok -----------------------
1671  if (d<=0)
1672  { "  ERROR:   the second parameter should be > 0";
1673     return();
1674  }
1675  if (size(#)>1)
1676  { "  ERROR:   there should be at most three parameters";
1677    return();
1678  }
1679  if (size(#)==1)
1680  { if (typeof(#[1])<>"intvec")
1681    { "  ERROR: the third parameter should be of type <intvec>";
1682      return();
1683    }
1684    if (size(#[1])<>2)
1685    { "  ERROR: there should be two components in <intvec>";
1686      return();
1687    }
1688    else
1689    { int cd=#[1][1];
1690      int step_fac=#[1][2];
1691    }
1692    if (step_fac<=0)
1693    { "  ERROR: the second component of <intvec> should be > 0";
1694      return();
1695    }
1696    if (cd==0)
1697    { return(ideal(0));
1698    }
1699  }
1700  else
1701  { int step_fac=1;
1702    int cd=-1;
1703  }
1704  if (ncols(REY)<>nvars(basering))
1705  { "ERROR:   the number of columns in the <matrix> should be the same as the";
1706    "         number of variables in the basering; in fact it should be first";
1707    "         return value of group_reynolds() or reynolds_molien().";
1708    return();
1709  }
1710 //----------------------------------------------------------------------------
1711  ideal mon=sort(maxideal(d))[1];
1712  degBound=d;
1713  int j=ncols(mon);
1714  mon=sort(mon,intvec(j..1))[1];
1715  ideal B;                             // will contain the basis
1716  if (cd<0)
1717  { if (step_fac>j)                    // all of mon will be mapped to
1718    { B=evaluate_reynolds(REY,mon);    // invariants at once
1719      B=minbase(B);
1720      degBound=0;
1721      return(B);
1722    }
1723  }
1724  else
1725  { if (step_fac*cd>j)                 // all of mon will be mapped to
1726    { B=evaluate_reynolds(REY,mon);    // invariants at once
1727      B=minbase(B);
1728      degBound=0;
1729      return(B);
1730    }
1731  }
1732  int i,k;
1733  int upper_bound=0;
1734  int lower_bound=0;
1735  ideal part_mon;                      // a part of mon of size step_fac*cd
1736  while (1)
1737  { lower_bound=upper_bound+1;
1738    if (cd<0)
1739    { upper_bound=upper_bound+step_fac;
1740    }
1741    else
1742    { upper_bound=upper_bound+step_fac*cd;
1743    }
1744    if (upper_bound>j)
1745    { upper_bound=j;
1746    }
1747    part_mon=mon[lower_bound..upper_bound];
1748    B=minbase(B+evaluate_reynolds(REY,part_mon));
1749    if ((ncols(B)==cd and B[1]<>0) or upper_bound==j)
1750    { degBound=0;
1751      return(B);
1752    }
1753  }
1754}
1755example
1756{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
1757  echo=2;
1758           ring R=0,(x,y,z),dp;
1759           matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
1760           intvec flags=0,1,0;
1761           matrix REY,M=reynolds_molien(A,flags);
1762           flags=8,6;
1763           print(invariant_basis_reynolds(REY,6,flags));
1764}
1765
1766////////////////////////////////////////////////////////////////////////////////
1767// This procedure generates linearly independent invariant polynomials of degree
1768// d that do not reduce to 0 modulo the primary invariants. It does this by
1769// applying the Reynolds operator to the monomials returned by kbase(sP,d). The
1770// result is used when computing secondary invariants.
1771////////////////////////////////////////////////////////////////////////////////
1772proc sort_of_invariant_basis (ideal sP,matrix REY,int d,int step_fac)
1773{ ideal mon=kbase(sP,d);
1774  degBound=d;
1775  int j=ncols(mon);
1776  int i;
1777  mon=sort(mon,intvec(j..1))[1];
1778  ideal B;                             // will contain the "sort of basis"
1779  if (step_fac>j)
1780  { B=compress(evaluate_reynolds(REY,mon));
1781    for (i=1;i<=ncols(B);i=i+1)        // those are taken our that are o mod sP
1782    { if (reduce(B[i],sP)==0)
1783      { B[i]=0;
1784      }
1785    }
1786    B=minbase(B);                      // here are the linearly independent ones
1787    degBound=0;
1788    return(B);
1789  }
1790  int upper_bound=0;
1791  int lower_bound=0;
1792  ideal part_mon;                      // parts of mon
1793  while (1)
1794  { lower_bound=upper_bound+1;
1795    upper_bound=upper_bound+step_fac;
1796    if (upper_bound>j)
1797    { upper_bound=j;
1798    }
1799    part_mon=mon[lower_bound..upper_bound];
1800    part_mon=compress(evaluate_reynolds(REY,part_mon));
1801    for (i=1;i<=ncols(part_mon);i=i+1)
1802    { if (reduce(part_mon[i],sP)==0)
1803      { part_mon[i]=0;
1804      }
1805    }
1806    B=minbase(B+part_mon);             // here are the linearly independent ones
1807    if (upper_bound==j)
1808    { degBound=0;
1809      return(B);
1810    }
1811  }
1812}
1813
1814////////////////////////////////////////////////////////////////////////////////
1815// Procedure returning the succeeding vector after vec. It is used to list
1816// all the vectors of Z^n with first nonzero entry 1. They are listed by
1817// increasing sum of the absolute value of their entries.
1818////////////////////////////////////////////////////////////////////////////////
1819proc next_vector(intmat vec)
1820{ int n=ncols(vec);                    // p: >0, n: <0, p0: >=0, n0: <=0
1821  for (int i=1;i<=n;i=i+1)             // finding out which is the first
1822  { if (vec[1,i]<>0)                   // component <>0
1823    { break;
1824    }
1825  }
1826  intmat new[1][n];
1827  if (i>n)                             // 0,...,0 --> 1,0....,0
1828  { new[1,1]=1;
1829    return(new);
1830  }
1831  if (i==n)                            // 0,...,1 --> 1,1,0,...,0
1832  { new[1,1..2]=1,1;
1833    return(new);
1834  }
1835  if (i==n-1)
1836  { if (vec[1,n]==0)                   // 0,...,0,1,0 --> 0,...,0,1
1837    { new[1,n]=1;
1838      return(new);
1839    }
1840    if (vec[1,n]>0)                    // 0,..,0,1,p --> 0,...,0,1,-p
1841    { new[1,1..n]=vec[1,1..n-1],-vec[1,n];
1842      return(new);
1843    }
1844    new[1,1..2]=1,1-vec[1,n];          // 0,..,0,1,n --> 1,1-n,0,..,0
1845    return(new);
1846  }
1847  if (i>1)
1848  { intmat temp[1][n-i+1]=vec[1,i..n]; // 0,...,0,1,*,...,* --> 1,*,...,*
1849    temp=next_vector(temp);
1850    new[1,i..n]=temp[1,1..n-i+1];
1851    return(new);
1852  }                                    // case left: 1,*,...,*
1853  for (i=2;i<=n;i=i+1)
1854  { if (vec[1,i]>0)                    // make first positive negative and
1855    { vec[1,i]=-vec[1,i];              // return
1856      return(vec);
1857    }
1858    else
1859    { vec[1,i]=-vec[1,i];              // make all negatives before positives
1860    }                                  // positive
1861  }
1862  for (i=2;i<=n-1;i=i+1)               // case: 1,p,...,p after 1,n,...,n
1863  { if (vec[1,i]>0)
1864    { vec[1,2]=vec[1,i]-1;             // shuffleing things around...
1865      if (i>2)                         // same sum of absolute values of entries
1866      { vec[1,i]=0;
1867      }
1868      vec[1,i+1]=vec[1,i+1]+1;
1869      return(vec);
1870    }
1871  }                                    // case left: 1,0,...,0 --> 1,1,0,...,0
1872  new[1,2..3]=1,vec[1,n];              // and: 1,0,...,0,1 --> 0,1,1,0,...,0
1873  return(new);
1874}
1875
1876////////////////////////////////////////////////////////////////////////////////
1877// Maps integers to elements of the base field. It is only called if the base
1878// field is of prime characteristic. If the base field has q elements (depending
1879// on minpoly) 1..q is mapped to those q elements.
1880////////////////////////////////////////////////////////////////////////////////
1881proc int_number_map (int i)
1882{ int p=char(basering);
1883  if (minpoly==0)                      // if no minpoly is given, we have p
1884  { i=i%p;                             // elements in the field
1885    return(number(i));
1886  }
1887  int d=pardeg(minpoly);
1888  if (i<0)
1889  { int bool=1;
1890    i=(-1)*i;
1891  }
1892  i=i%p^d;                             // base field has p^d elements -
1893  number a=par(1);                     // a is the root of the minpoly - we have
1894  number out=0;                        // to construct a linear combination of
1895  int j=1;                             // a^k
1896  int k;
1897  while (1)
1898  { if (i<p^j)                         // finding an upper bound on i
1899    { for (k=0;k<j-1;k=k+1)
1900      { out=out+((i/p^k)%p)*a^k;       // finding how often p^k is contained in
1901      }                                // i
1902      out=out+(i/p^(j-1))*a^(j-1);
1903      if (defined(bool)==voice)
1904      { return((-1)*out);
1905      }
1906      return(out);
1907    }
1908    j=j+1;
1909  }
1910}
1911
1912////////////////////////////////////////////////////////////////////////////////
1913// This procedure finds dif primary invariants in degree d. It returns all
1914// primary invariants found so far. The coefficients lie in a field of
1915// characteristic 0.
1916////////////////////////////////////////////////////////////////////////////////
1917proc search (int n,int d,ideal B,int cd,ideal P,ideal sP,int i,int dif,int dB,ideal CI)
1918{ intmat vec[1][cd];                   // the coefficients for the next
1919                                       // combination -
1920  degBound=0;
1921  poly test_poly;                      // the linear combination to test
1922  int test_dim;
1923  intvec h;                            // Hilbert series
1924  int j=i+1;
1925  matrix tB=transpose(B);
1926  ideal TEST;
1927  while(j<=i+dif)
1928  { CI=CI+ideal(var(j)^d);             // homogeneous polynomial of the same
1929                                       // degree as the one we're looking for is
1930                                       // added
1931    // h=hilb(std(CI),1);
1932    dB=dB+d-1;                         // used as degBound
1933    while(1)
1934    { vec=next_vector(vec);            // next vector
1935      test_poly=(vec*tB)[1,1];
1936      // degBound=dB;
1937      TEST=sP+ideal(test_poly);
1938      attrib(TEST,"isSB",1);
1939      test_dim=dim(TEST);
1940      // degBound=0;
1941      if (n-test_dim==j)               // the dimension has been lowered by one
1942      { sP=TEST;
1943        break;
1944      }
1945      // degBound=dB;
1946      TEST=std(sP+ideal(test_poly));   // should soon be replaced by next line
1947      // TEST=std(sP,test_poly,h);        // Hilbert driven std-calculation
1948      test_dim=dim(TEST);
1949      // degBound=0;
1950      if (n-test_dim==j)               // the dimension has been lowered by one
1951      { sP=TEST;
1952        break;
1953      }
1954    }
1955    P[j]=test_poly;                    // test_poly ist added to primary
1956    j=j+1;                             // invariants
1957  }
1958  return(P,sP,CI,dB);
1959}
1960
1961////////////////////////////////////////////////////////////////////////////////
1962// This procedure finds at most dif primary invariants in degree d. It returns
1963// all primary invariants found so far. The coefficients lie in the field of
1964// characteristic p>0.
1965////////////////////////////////////////////////////////////////////////////////
1966proc p_search (int n,int d,ideal B,int cd,ideal P,ideal sP,int i,int dif,int dB,ideal CI)
1967{ def br=basering;
1968  degBound=0;
1969  matrix vec(1)[1][cd];                // starting with 0-vector -
1970  intmat new[1][cd];                   // the coefficients for the next
1971                                       // combination -
1972  matrix pnew[1][cd];                  // new needs to be mapped into br -
1973  int counter=1;                       // counts the vectors
1974  int j;
1975  int p=char(br);
1976  if (minpoly<>0)
1977  { int ext_deg=pardeg(minpoly);       // field has p^d elements
1978  }
1979  else
1980  { int ext_deg=1;                     // field has p^d elements
1981  }
1982  poly test_poly;                      // the linear combination to test
1983  int test_dim;
1984  ring R=0,x,dp;                       // just to calculate next variable
1985                                       // bound -
1986  number bound=(number(p)^(ext_deg*cd)-1)/(number(p)^ext_deg-1)+1; // this is
1987                                       // how many linearly independent vectors
1988                                       // of size cd exist having entries in the
1989                                       // base field of br
1990  setring br;
1991  intvec h;                            // Hilbert series
1992  int k=i+1;
1993  matrix tB=transpose(B);
1994  ideal TEST;
1995  while (k<=i+dif)
1996  { CI=CI+ideal(var(k)^d);             // homogeneous polynomial of the same
1997                                       // degree as the one we're looking for is
1998                                       // added
1999    // h=hilb(std(CI),1);
2000    dB=dB+d-1;                         // used as degBound
2001    setring R;
2002    while (number(counter)<>bound)     // otherwise, we are done
2003    { setring br;
2004      new=next_vector(new);
2005      for (j=1;j<=cd;j=j+1)
2006      { pnew[1,j]=int_number_map(new[1,j]); // mapping an integer into br
2007      }
2008      if (unique(vec(1..counter),pnew)) // checking whether we tried pnew before
2009      { counter=counter+1;
2010        matrix vec(counter)=pnew;      // keeping track of the ones we tried -
2011        test_poly=(vec(counter)*tB)[1,1]; // linear combination -
2012        // degBound=dB;
2013        TEST=sP+ideal(test_poly);
2014        attrib(TEST,"isSB",1);
2015        test_dim=dim(TEST);
2016        // degBound=0;
2017        if (n-test_dim==k)             // the dimension has been lowered by one
2018        { sP=TEST;
2019          setring R;
2020          break;
2021        }
2022        // degBound=dB;
2023        TEST=std(sP+ideal(test_poly)); // should soon to be replaced by next
2024                                       // line
2025        // TEST=std(sP,test_poly,h);      // Hilbert driven std-calculation
2026        test_dim=dim(TEST);
2027        // degBound=0;
2028        if (n-test_dim==k)             // the dimension has been lowered by one
2029        { sP=TEST;
2030          setring R;
2031          break;
2032        }
2033      }
2034      setring R;
2035    }
2036    if (number(counter)<=bound)
2037    { setring br;
2038      P[k]=test_poly;                  // test_poly ist added to primary
2039    }                                  // invariants
2040    else
2041    { setring br;
2042      CI=CI[1..size(CI)-1];
2043      return(P,sP,CI,dB-d+1);
2044    }
2045    k=k+1;
2046  }
2047  return(P,sP,CI,dB);
2048}
2049
2050proc primary_char0 (matrix REY,matrix M,list #)
2051"USAGE:   primary_char0(REY,M[,v]);
2052         REY: a <matrix> representing the Reynolds operator, M: a 1x2 <matrix>
2053         representing the Molien series, v: an optional <int>
2054ASSUME:  REY is the first return value of group_reynolds or reynolds_molien and
2055         M the one of molien or the second one of reynolds_molien
2056DISPLAY: information about the various stages of the programme if v does not
2057         equal 0
2058RETURN:  primary invariants (type <matrix>) of the invariant ring
2059EXAMPLE: example primary_char0; shows an example
2060THEORY:  Bases of homogeneous invariants are generated successively and those
2061         are chosen as primary invariants that lower the dimension of the ideal
2062         generated by the previously found invariants (see paper \"Generating a
2063         Noetherian Normalization of the Invariant Ring of a Finite Group\" by
2064         Decker, Heydtmann, Schreyer (1997) to appear in JSC).
2065"
2066{ degBound=0;
2067  if (char(basering)<>0)
2068  { "ERROR:   primary_char0 should only be used with rings of characteristic 0.";
2069    return();
2070  }
2071 //----------------- checking input and setting verbose mode ------------------
2072  if (size(#)>1)
2073  { "ERROR:   primary_char0 can only have three parameters.";
2074    return();
2075  }
2076  if (size(#)==1)
2077  { if (typeof(#[1])<>"int")
2078    { "ERROR:   The third parameter should be of type <int>.";
2079      return();
2080    }
2081    else
2082    { int v=#[1];
2083    }
2084  }
2085  else
2086  { int v=0;
2087  }
2088  int n=nvars(basering);               // n is the number of variables, as well
2089                                       // as the size of the matrices, as well
2090                                       // as the number of primary invariants,
2091                                       // we should get
2092  if (ncols(REY)<>n)
2093  { "ERROR:   First parameter ought to be the Reynolds operator."
2094    return();
2095  }
2096  if (ncols(M)<>2 or nrows(M)<>1)
2097  { "ERROR:   Second parameter ought to be the Molien series."
2098    return();
2099  }
2100 //----------------------------------------------------------------------------
2101  if (v && voice<>2)
2102  { "  We can start looking for primary invariants...";
2103    "";
2104  }
2105  if (v && voice==2)
2106  { "";
2107  }
2108 //------------------------- initializing variables ---------------------------
2109  int dB;
2110  poly p(1..2);                        // p(1) will be used for single terms of
2111                                       // the partial expansion, p(2) to store
2112  p(1..2)=partial_molien(M,1);         // the intermediate result -
2113  poly v1=var(1);                      // we need v1 to split off coefficients
2114                                       // in the partial expansion of M (which
2115                                       // is in terms of the first variable) -
2116  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
2117                                       // space of invariants of degree d,
2118                                       // newdim: dimension the ideal generated
2119                                       // the primary invariants plus basis
2120                                       // elements, dif=n-i-newdim, i.e. the
2121                                       // number of new primary invairants that
2122                                       // should be added in this degree -
2123  ideal P,Pplus,sPplus,CI,B;           // P: will contain primary invariants,
2124                                       // Pplus: P+B, CI: a complete
2125                                       // intersection with the same Hilbert
2126                                       // function as P
2127  ideal sP=std(P);
2128  dB=1;                                // used as degree bound
2129  int i=0;
2130 //-------------- loop that searches for primary invariants  ------------------
2131  while(1)                             // repeat until n primary invariants are
2132  {                                    // found -
2133    p(1..2)=partial_molien(M,1,p(2));  // next term of the partial expansion -
2134    d=deg(p(1));                       // degree where we'll search -
2135    cd=int(coef(p(1),v1)[2,1]);        // dimension of the homogeneous space of
2136                                       // inviarants of degree d
2137    if (v)
2138    { "  Computing primary invariants in degree "+string(d)+":";
2139    }
2140    B=invariant_basis_reynolds(REY,d,intvec(cd,6)); // basis of invariants of
2141                                       // degree d
2142    if (B[1]<>0)
2143    { Pplus=P+B;
2144      sPplus=std(Pplus);
2145      newdim=dim(sPplus);
2146      dif=n-i-newdim;
2147    }
2148    else
2149    { dif=0;
2150    }
2151    if (dif<>0)                        // we have to find dif new primary
2152    {                                  // invariants
2153      if (cd<>dif)
2154      { P,sP,CI,dB=search(n,d,B,cd,P,sP,i,dif,dB,CI); // searching for dif invariants
2155      }                                // i.e. we can take all of B
2156      else
2157      { for(j=i+1;j>i+dif;j=j+1)
2158        { CI=CI+ideal(var(j)^d);
2159        }
2160        dB=dB+dif*(d-1);
2161        P=Pplus;
2162        sP=sPplus;
2163      }
2164      if (v)
2165      { for (j=1;j<=dif;j=j+1)
2166        { "  We find: "+string(P[i+j]);
2167        }
2168      }
2169      i=i+dif;
2170      if (i==n)                        // found all primary invariants
2171      { if (v)
2172        { "";
2173          "  We found all primary invariants.";
2174          "";
2175        }
2176        return(matrix(P));
2177      }
2178    }                                  // done with degree d
2179  }
2180}
2181example
2182{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
2183  echo=2;
2184         ring R=0,(x,y,z),dp;
2185         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
2186         matrix REY,M=reynolds_molien(A);
2187         matrix P=primary_char0(REY,M);
2188         print(P);
2189}
2190
2191proc primary_charp (matrix REY,string ring_name,list #)
2192"USAGE:   primary_charp(REY,ringname[,v]);
2193         REY: a <matrix> representing the Reynolds operator, ringname: a
2194         <string> giving the name of a ring where the Molien series is stored,
2195         v: an optional <int>
2196ASSUME:  REY is the first return value of group_reynolds or reynolds_molien and
2197         ringname gives the name of a ring of characteristic 0 that has been
2198         created by molien or reynolds_molien
2199DISPLAY: information about the various stages of the programme if v does not
2200         equal 0
2201RETURN:  primary invariants (type <matrix>) of the invariant ring
2202EXAMPLE: example primary_charp; shows an example
2203THEORY:  Bases of homogeneous invariants are generated successively and those
2204         are chosen as primary invariants that lower the dimension of the ideal
2205         generated by the previously found invariants (see paper \"Generating a
2206         Noetherian Normalization of the Invariant Ring of a Finite Group\" by
2207         Decker, Heydtmann, Schreyer (1997) to appear in JSC).
2208"
2209{ degBound=0;
2210// ---------------- checking input and setting verbose mode -------------------
2211  if (char(basering)==0)
2212  { "ERROR:   primary_charp should only be used with rings of characteristic p>0.";
2213    return();
2214  }
2215  if (size(#)>1)
2216  { "ERROR:   primary_charp can only have three parameters.";
2217    return();
2218  }
2219  if (size(#)==1)
2220  { if (typeof(#[1])<>"int")
2221    { "ERROR:   The third parameter should be of type <int>.";
2222      return();
2223    }
2224    else
2225    { int v=#[1];
2226    }
2227  }
2228  else
2229  { int v=0;
2230  }
2231  def br=basering;
2232  int n=nvars(br);                     // n is the number of variables, as well
2233                                       // as the size of the matrices, as well
2234                                       // as the number of primary invariants,
2235                                       // we should get
2236  if (ncols(REY)<>n)
2237  { "ERROR:   First parameter ought to be the Reynolds operator."
2238    return();
2239  }
2240  if (typeof(`ring_name`)<>"ring")
2241  { "ERROR:   Second parameter ought to the name of a ring where the Molien";
2242    "         is stored.";
2243    return();
2244  }
2245 //----------------------------------------------------------------------------
2246  if (v && voice<>2)
2247  { "  We can start looking for primary invariants...";
2248    "";
2249  }
2250  if (v && voice==2)
2251  { "";
2252  }
2253 //----------------------- initializing variables -----------------------------
2254  int dB;
2255  setring `ring_name`;                 // the Molien series is stores here -
2256  poly p(1..2);                        // p(1) will be used for single terms of
2257                                       // the partial expansion, p(2) to store
2258  p(1..2)=partial_molien(M,1);         // the intermediate result -
2259  poly v1=var(1);                      // we need v1 to split off coefficients
2260                                       // in the partial expansion of M (which
2261                                       // is in terms of the first variable)
2262  setring br;
2263  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
2264                                       // space of invariants of degree d,
2265                                       // newdim: dimension the ideal generated
2266                                       // the primary invariants plus basis
2267                                       // elements, dif=n-i-newdim, i.e. the
2268                                       // number of new primary invairants that
2269                                       // should be added in this degree -
2270  ideal P,Pplus,sPplus,CI,B;           // P: will contain primary invariants,
2271                                       // Pplus: P+B, CI: a complete
2272                                       // intersection with the same Hilbert
2273                                       // function as P
2274  ideal sP=std(P);
2275  dB=1;                                // used as degree bound
2276  int i=0;
2277 //---------------- loop that searches for primary invariants -----------------
2278  while(1)                             // repeat until n primary invariants are
2279  {                                    // found
2280    setring `ring_name`;
2281    p(1..2)=partial_molien(M,1,p(2));  // next term of the partial expansion -
2282    d=deg(p(1));                       // degree where we'll search -
2283    cd=int(coef(p(1),v1)[2,1]);        // dimension of the homogeneous space of
2284                                       // inviarants of degree d
2285    setring br;
2286    if (v)
2287    { "  Computing primary invariants in degree "+string(d)+":";
2288    }
2289    B=invariant_basis_reynolds(REY,d,intvec(cd,6)); // basis of invariants of
2290                                       // degree d
2291    if (B[1]<>0)
2292    { Pplus=P+B;
2293      sPplus=std(Pplus);
2294      newdim=dim(sPplus);
2295      dif=n-i-newdim;
2296    }
2297    else
2298    { dif=0;
2299    }
2300    if (dif<>0)                        // we have to find dif new primary
2301    {                                  // invariants
2302      if (cd<>dif)
2303      { P,sP,CI,dB=p_search(n,d,B,cd,P,sP,i,dif,dB,CI);
2304      }
2305      else                             // i.e. we can take all of B
2306      { for(j=i+1;j>i+dif;j=j+1)
2307        { CI=CI+ideal(var(j)^d);
2308        }
2309        dB=dB+dif*(d-1);
2310        P=Pplus;
2311        sP=sPplus;
2312      }
2313      if (v)
2314      { for (j=1;j<=size(P)-i;j=j+1)
2315        { "  We find: "+string(P[i+j]);
2316        }
2317      }
2318      i=size(P);
2319      if (i==n)                        // found all primary invariants
2320      { if (v)
2321        { "";
2322          "  We found all primary invariants.";
2323          "";
2324        }
2325        return(matrix(P));
2326      }
2327    }                                  // done with degree d
2328  }
2329}
2330example
2331{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed into";
2332  "         characteristic 3)";
2333  echo=2;
2334         ring R=3,(x,y,z),dp;
2335         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
2336         list L=group_reynolds(A);
2337         string newring="alskdfj";
2338         molien(L[2..size(L)],newring);
2339         matrix P=primary_charp(L[1],newring);
2340         if(system("with","Namespaces")) { kill Top::`newring`; }
2341         kill `newring`;
2342         print(P);
2343}
2344
2345proc primary_char0_no_molien (matrix REY, list #)
2346"USAGE:   primary_char0_no_molien(REY[,v]);
2347         REY: a <matrix> representing the Reynolds operator, v: an optional
2348         <int>
2349ASSUME:  REY is the first return value of group_reynolds or reynolds_molien
2350DISPLAY: information about the various stages of the programme if v does not
2351         equal 0
2352RETURN:  primary invariants (type <matrix>) of the invariant ring and an
2353         <intvec> listing some of the degrees where no non-trivial homogeneous
2354         invariants are to be found
2355EXAMPLE: example primary_char0_no_molien; shows an example
2356THEORY:  Bases of homogeneous invariants are generated successively and those
2357         are chosen as primary invariants that lower the dimension of the ideal
2358         generated by the previously found invariants (see paper \"Generating a
2359         Noetherian Normalization of the Invariant Ring of a Finite Group\" by
2360         Decker, Heydtmann, Schreyer (1997) to appear in JSC).
2361"
2362{ degBound=0;
2363 //-------------- checking input and setting verbose mode ---------------------
2364  if (char(basering)<>0)
2365  { "ERROR:   primary_char0_no_molien should only be used with rings of";
2366    "         characteristic 0.";
2367    return();
2368  }
2369  if (size(#)>1)
2370  { "ERROR:   primary_char0_no_molien can only have two parameters.";
2371    return();
2372  }
2373  if (size(#)==1)
2374  { if (typeof(#[1])<>"int")
2375    { "ERROR:   The second parameter should be of type <int>.";
2376      return();
2377    }
2378    else
2379    { int v=#[1];
2380    }
2381  }
2382  else
2383  { int v=0;
2384  }
2385  int n=nvars(basering);               // n is the number of variables, as well
2386                                       // as the size of the matrices, as well
2387                                       // as the number of primary invariants,
2388                                       // we should get
2389  if (ncols(REY)<>n)
2390  { "ERROR:   First parameter ought to be the Reynolds operator."
2391    return();
2392  }
2393 //----------------------------------------------------------------------------
2394  if (v && voice<>2)
2395  { "  We can start looking for primary invariants...";
2396    "";
2397  }
2398  if (v && voice==2)
2399  { "";
2400  }
2401 //----------------------- initializing variables -----------------------------
2402  int dB;
2403  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
2404                                       // space of invariants of degree d,
2405                                       // newdim: dimension the ideal generated
2406                                       // the primary invariants plus basis
2407                                       // elements, dif=n-i-newdim, i.e. the
2408                                       // number of new primary invairants that
2409                                       // should be added in this degree -
2410  ideal P,Pplus,CI,B;                  // P: will contain primary invariants,
2411                                       // Pplus: P+B, CI: a complete
2412                                       // intersection with the same Hilbert
2413                                       // function as P
2414  ideal sP=std(P);
2415  dB=1;                                // used as degree bound -
2416  d=0;                                 // initializing
2417  int i=0;
2418  intvec deg_vector;
2419 //------------------ loop that searches for primary invariants ---------------
2420  while(1)                             // repeat until n primary invariants are
2421  {                                    // found -
2422    d=d+1;                             // degree where we'll search
2423    if (v)
2424    { "  Computing primary invariants in degree "+string(d)+":";
2425    }
2426    B=invariant_basis_reynolds(REY,d,intvec(-1,6)); // basis of invariants of
2427                                       // degree d
2428    if (B[1]<>0)
2429    { Pplus=P+B;
2430      newdim=dim(std(Pplus));
2431      dif=n-i-newdim;
2432    }
2433    else
2434    { dif=0;
2435      deg_vector=deg_vector,d;
2436    }
2437    if (dif<>0)                        // we have to find dif new primary
2438    {                                  // invariants
2439      cd=size(B);
2440      if (cd<>dif)
2441      { P,sP,CI,dB=search(n,d,B,cd,P,sP,i,dif,dB,CI);
2442      }
2443      else                             // i.e. we can take all of B
2444      { for(j=i+1;j<=i+dif;j=j+1)
2445        { CI=CI+ideal(var(j)^d);
2446        }
2447        dB=dB+dif*(d-1);
2448        P=Pplus;
2449        sP=std(P);
2450      }
2451      if (v)
2452      { for (j=1;j<=dif;j=j+1)
2453        { "  We find: "+string(P[i+j]);
2454        }
2455      }
2456      i=i+dif;
2457      if (i==n)                        // found all primary invariants
2458      { if (v)
2459        { "";
2460          "  We found all primary invariants.";
2461          "";
2462        }
2463        if (deg_vector==0)
2464        { return(matrix(P));
2465        }
2466        else
2467        { return(matrix(P),compress(deg_vector));
2468        }
2469      }
2470    }                                  // done with degree d
2471    else
2472    { if (v)
2473      { "  None here...";
2474      }
2475    }
2476  }
2477}
2478example
2479{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
2480  echo=2;
2481         ring R=0,(x,y,z),dp;
2482         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
2483         list L=group_reynolds(A);
2484         list l=primary_char0_no_molien(L[1]);
2485         print(l[1]);
2486}
2487
2488proc primary_charp_no_molien (matrix REY, list #)
2489"USAGE:   primary_charp_no_molien(REY[,v]);
2490         REY: a <matrix> representing the Reynolds operator, v: an optional
2491         <int>
2492ASSUME:  REY is the first return value of group_reynolds or reynolds_molien
2493DISPLAY: information about the various stages of the programme if v does not
2494         equal 0
2495RETURN:  primary invariants (type <matrix>) of the invariant ring  and an
2496         <intvec> listing some of the degrees where no non-trivial homogeneous
2497         invariants are to be found
2498EXAMPLE: example primary_charp_no_molien; shows an example
2499THEORY:  Bases of homogeneous invariants are generated successively and those
2500         are chosen as primary invariants that lower the dimension of the ideal
2501         generated by the previously found invariants (see paper \"Generating a
2502         Noetherian Normalization of the Invariant Ring of a Finite Group\" by
2503         Decker, Heydtmann, Schreyer (1997) to appear in JSC).
2504"
2505{ degBound=0;
2506 //----------------- checking input and setting verbose mode ------------------
2507  if (char(basering)==0)
2508  { "ERROR:   primary_charp_no_molien should only be used with rings of";
2509    "         characteristic p>0.";
2510    return();
2511  }
2512  if (size(#)>1)
2513  { "ERROR:   primary_charp_no_molien can only have two parameters.";
2514    return();
2515  }
2516  if (size(#)==1)
2517  { if (typeof(#[1])<>"int")
2518    { "ERROR:   The second parameter should be of type <int>.";
2519      return();
2520    }
2521    else
2522    { int v=#[1];
2523    }
2524  }
2525  else
2526  { int v=0;
2527  }
2528  int n=nvars(basering);               // n is the number of variables, as well
2529                                       // as the size of the matrices, as well
2530                                       // as the number of primary invariants,
2531                                       // we should get
2532  if (ncols(REY)<>n)
2533  { "ERROR:   First parameter ought to be the Reynolds operator."
2534    return();
2535  }
2536 //----------------------------------------------------------------------------
2537  if (v && voice<>2)
2538  { "  We can start looking for primary invariants...";
2539    "";
2540  }
2541  if (v && voice==2)
2542  { "";
2543  }
2544 //-------------------- initializing variables --------------------------------
2545  int dB;
2546  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
2547                                       // space of invariants of degree d,
2548                                       // newdim: dimension the ideal generated
2549                                       // the primary invariants plus basis
2550                                       // elements, dif=n-i-newdim, i.e. the
2551                                       // number of new primary invairants that
2552                                       // should be added in this degree -
2553  ideal P,Pplus,sPplus,CI,B;           // P: will contain primary invariants,
2554                                       // Pplus: P+B, CI: a complete
2555                                       // intersection with the same Hilbert
2556                                       // function as P
2557  ideal sP=std(P);
2558  dB=1;                                // used as degree bound -
2559  d=0;                                 // initializing
2560  int i=0;
2561  intvec deg_vector;
2562 //------------------ loop that searches for primary invariants ---------------
2563  while(1)                             // repeat until n primary invariants are
2564  {                                    // found -
2565    d=d+1;                             // degree where we'll search
2566    if (v)
2567    { "  Computing primary invariants in degree "+string(d)+":";
2568    }
2569    B=invariant_basis_reynolds(REY,d,intvec(-1,6)); // basis of invariants of
2570                                       // degree d
2571    if (B[1]<>0)
2572    { Pplus=P+B;
2573      sPplus=std(Pplus);
2574      newdim=dim(sPplus);
2575      dif=n-i-newdim;
2576    }
2577    else
2578    { dif=0;
2579      deg_vector=deg_vector,d;
2580    }
2581    if (dif<>0)                        // we have to find dif new primary
2582    {                                  // invariants
2583      cd=size(B);
2584      if (cd<>dif)
2585      { P,sP,CI,dB=p_search(n,d,B,cd,P,sP,i,dif,dB,CI);
2586      }
2587      else                             // i.e. we can take all of B
2588      { for(j=i+1;j<=i+dif;j=j+1)
2589        { CI=CI+ideal(var(j)^d);
2590        }
2591        dB=dB+dif*(d-1);
2592        P=Pplus;
2593        sP=sPplus;
2594      }
2595      if (v)
2596      { for (j=1;j<=size(P)-i;j=j+1)
2597        { "  We find: "+string(P[i+j]);
2598        }
2599      }
2600      i=size(P);
2601      if (i==n)                        // found all primary invariants
2602      { if (v)
2603        { "";
2604          "  We found all primary invariants.";
2605          "";
2606        }
2607        if (deg_vector==0)
2608        { return(matrix(P));
2609        }
2610        else
2611        { return(matrix(P),compress(deg_vector));
2612        }
2613      }
2614    }                                  // done with degree d
2615    else
2616    { if (v)
2617      { "  None here...";
2618      }
2619    }
2620  }
2621}
2622example
2623{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed into";
2624  "         characteristic 3)";
2625  echo=2;
2626         ring R=3,(x,y,z),dp;
2627         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
2628         list L=group_reynolds(A);
2629         list l=primary_charp_no_molien(L[1]);
2630         print(l[1]);
2631}
2632
2633proc primary_charp_without (list #)
2634"USAGE:   primary_charp_without(G1,G2,...[,v]);
2635         G1,G2,...: <matrices> generating a finite matrix group, v: an optional
2636         <int>
2637DISPLAY: information about the various stages of the programme if v does not
2638         equal 0
2639RETURN:  primary invariants (type <matrix>) of the invariant ring
2640EXAMPLE: example primary_charp_without; shows an example
2641THEORY:  Bases of homogeneous invariants are generated successively and those
2642         are chosen as primary invariants that lower the dimension of the ideal
2643         generated by the previously found invariants (see paper \"Generating a
2644         Noetherian Normalization of the Invariant Ring of a Finite Group\" by
2645         Decker, Heydtmann, Schreyer (1997) to appear in JSC). No Reynolds
2646         operator or Molien series is used.
2647"
2648{ degBound=0;
2649 //--------------------- checking input and setting verbose mode --------------
2650  if (char(basering)==0)
2651  { "ERROR:   primary_charp_without should only be used with rings of";
2652    "         characteristic 0.";
2653    return();
2654  }
2655  if (size(#)==0)
2656  { "ERROR:   There are no parameters.";
2657    return();
2658  }
2659  if (typeof(#[size(#)])=="int")
2660  { int v=#[size(#)];
2661    int gen_num=size(#)-1;
2662    if (gen_num==0)
2663    { "ERROR:   There are no generators of a finite matrix group given.";
2664      return();
2665    }
2666  }
2667  else
2668  { int v=0;
2669    int gen_num=size(#);
2670  }
2671  int n=nvars(basering);               // n is the number of variables, as well
2672                                       // as the size of the matrices, as well
2673                                       // as the number of primary invariants,
2674                                       // we should get
2675  for (int i=1;i<=gen_num;i=i+1)
2676  { if (typeof(#[i])=="matrix")
2677    { if (nrows(#[i])<>n or ncols(#[i])<>n)
2678      { "ERROR:   The number of variables of the base ring needs to be the same";
2679        "         as the dimension of the square matrices";
2680        return();
2681      }
2682    }
2683    else
2684    { "ERROR:   The first parameters should be a list of matrices";
2685      return();
2686    }
2687  }
2688 //----------------------------------------------------------------------------
2689  if (v && voice==2)
2690  { "";
2691  }
2692 //---------------------------- initializing variables ------------------------
2693  int dB;
2694  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
2695                                       // space of invariants of degree d,
2696                                       // newdim: dimension the ideal generated
2697                                       // the primary invariants plus basis
2698                                       // elements, dif=n-i-newdim, i.e. the
2699                                       // number of new primary invairants that
2700                                       // should be added in this degree -
2701  ideal P,Pplus,sPplus,CI,B;           // P: will contain primary invariants,
2702                                       // Pplus: P+B, CI: a complete
2703                                       // intersection with the same Hilbert
2704                                       // function as P
2705  ideal sP=std(P);
2706  dB=1;                                // used as degree bound -
2707  d=0;                                 // initializing
2708  i=0;
2709  intvec deg_vector;
2710 //-------------------- loop that searches for primary invariants -------------
2711  while(1)                             // repeat until n primary invariants are
2712  {                                    // found -
2713    d=d+1;                             // degree where we'll search
2714    if (v)
2715    { "  Computing primary invariants in degree "+string(d)+":";
2716    }
2717    B=invariant_basis(d,#[1..gen_num]); // basis of invariants of degree d
2718    if (B[1]<>0)
2719    { Pplus=P+B;
2720      sPplus=std(Pplus);
2721      newdim=dim(sPplus);
2722      dif=n-i-newdim;
2723    }
2724    else
2725    { dif=0;
2726      deg_vector=deg_vector,d;
2727    }
2728    if (dif<>0)                        // we have to find dif new primary
2729    {                                  // invariants
2730      cd=size(B);
2731      if (cd<>dif)
2732      { P,sP,CI,dB=p_search(n,d,B,cd,P,sP,i,dif,dB,CI);
2733      }
2734      else                             // i.e. we can take all of B
2735      { for(j=i+1;j<=i+dif;j=j+1)
2736        { CI=CI+ideal(var(j)^d);
2737        }
2738        dB=dB+dif*(d-1);
2739        P=Pplus;
2740        sP=sPplus;
2741      }
2742      if (v)
2743      { for (j=1;j<=size(P)-i;j=j+1)
2744        { "  We find: "+string(P[i+j]);
2745        }
2746      }
2747      i=size(P);
2748      if (i==n)                        // found all primary invariants
2749      { if (v)
2750        { "";
2751          "  We found all primary invariants.";
2752          "";
2753        }
2754        return(matrix(P));
2755      }
2756    }                                  // done with degree d
2757    else
2758    { if (v)
2759      { "  None here...";
2760      }
2761    }
2762  }
2763}
2764example
2765{ echo=2;
2766         ring R=2,(x,y,z),dp;
2767         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
2768         matrix P=primary_charp_without(A);
2769         print(P);
2770}
2771
2772proc primary_invariants (list #)
2773"USAGE:   primary_invariants(G1,G2,...[,flags]);
2774         G1,G2,...: <matrices> generating a finite matrix group, flags: an
2775         optional <intvec> with three entries, if the first one equals 0 (also
2776         the default), the programme attempts to compute the Molien series and
2777         Reynolds operator, if it equals 1, the programme is told that the
2778         Molien series should not be computed, if it equals -1 characteristic 0
2779         is simulated, i.e. the Molien series is computed as if the base field
2780         were characteristic 0 (the user must choose a field of large prime
2781         characteristic, e.g. 32003) and if the first one is anything else, it
2782         means that the characteristic of the base field divides the group
2783         order, the second component should give the size of intervals between
2784         canceling common factors in the expansion of the Molien series, 0 (the
2785         default) means only once after generating all terms, in prime
2786         characteristic also a negative number can be given to indicate that
2787         common factors should always be canceled when the expansion is simple
2788         (the root of the extension field does not occur among the coefficients)
2789DISPLAY: information about the various stages of the programme if the third
2790         flag does not equal 0
2791RETURN:  primary invariants (type <matrix>) of the invariant ring and if
2792         computable Reynolds operator (type <matrix>) and Molien series (type
2793         <matrix>), if the first flag is 1 and we are in the non-modular case
2794         then an <intvec> is returned giving some of the degrees where no
2795         non-trivial homogeneous invariants can be found
2796EXAMPLE: example primary_invariants; shows an example
2797THEORY:  Bases of homogeneous invariants are generated successively and those
2798         are chosen as primary invariants that lower the dimension of the ideal
2799         generated by the previously found invariants (see paper \"Generating a
2800         Noetherian Normalization of the Invariant Ring of a Finite Group\" by
2801         Decker, Heydtmann, Schreyer (1997) to appear in JSC).
2802"
2803{
2804 // ----------------- checking input and setting flags ------------------------
2805  if (size(#)==0)
2806  { "ERROR:   There are no parameters.";
2807    return();
2808  }
2809  int ch=char(basering);               // the algorithms depend very much on the
2810                                       // characteristic of the ground field
2811  int n=nvars(basering);               // n is the number of variables, as well
2812                                       // as the size of the matrices, as well
2813                                       // as the number of primary invariants,
2814                                       // we should get
2815  int gen_num;
2816  int mol_flag,v;
2817  if (typeof(#[size(#)])=="intvec")
2818  { if (size(#[size(#)])<>3)
2819    { "ERROR:   <intvec> should have three entries.";
2820      return();
2821    }
2822    gen_num=size(#)-1;
2823    mol_flag=#[size(#)][1];
2824    if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag==-1)))
2825    { "ERROR:   the second component of <intvec> should be >=0";
2826      return();
2827    }
2828    int interval=#[size(#)][2];
2829    v=#[size(#)][3];
2830    if (gen_num==0)
2831    { "ERROR:   There are no generators of a finite matrix group given.";
2832      return();
2833    }
2834  }
2835  else
2836  { gen_num=size(#);
2837    mol_flag=0;
2838    int interval=0;
2839    v=0;
2840  }
2841  for (int i=1;i<=gen_num;i=i+1)
2842  { if (typeof(#[i])=="matrix")
2843    { if (nrows(#[i])<>n or ncols(#[i])<>n)
2844      { "ERROR:   The number of variables of the base ring needs to be the same";
2845        "         as the dimension of the square matrices";
2846        return();
2847      }
2848    }
2849    else
2850    { "ERROR:   The first parameters should be a list of matrices";
2851      return();
2852    }
2853  }
2854 //----------------------------------------------------------------------------
2855  if (mol_flag==0)
2856  { if (ch==0)
2857    { matrix REY,M=reynolds_molien(#[1..gen_num],intvec(mol_flag,interval,v));
2858                                       // one will contain Reynolds operator and
2859                                       // the other enumerator and denominator
2860                                       // of Molien series
2861      matrix P=primary_char0(REY,M,v);
2862      return(P,REY,M);
2863    }
2864    else
2865    { list L=group_reynolds(#[1..gen_num],v);
2866      if (L[1]<>0)                     // testing whether we are in the modular
2867      { string newring="aksldfalkdsflkj"; // case
2868        if (minpoly==0)
2869        { if (v)
2870          { "  We are dealing with the non-modular case.";
2871          }
2872          if (typeof(L[2])=="int")
2873          { molien(L[3..size(L)],newring,L[2],intvec(mol_flag,interval,v));
2874          }
2875          else
2876          { molien(L[2..size(L)],newring,intvec(mol_flag,interval,v));
2877          }
2878          matrix P=primary_charp(L[1],newring,v);
2879          return(P,L[1],newring);
2880        }
2881        else
2882        { if (v)
2883          { "  Since it is impossible for this programme to calculate the Molien series for";
2884            "  invariant rings over extension fields of prime characteristic, we have to";
2885            "  continue without it.";
2886            "";
2887
2888          }
2889          list l=primary_charp_no_molien(L[1],v);
2890          if (size(l)==2)
2891          { return(l[1],L[1],l[2]);
2892          }
2893          else
2894          { return(l[1],L[1]);
2895          }
2896        }
2897      }
2898      else                             // the modular case
2899      { if (v)
2900        { "  There is also no Molien series, we can make use of...";
2901          "";
2902          "  We can start looking for primary invariants...";
2903          "";
2904        }
2905        return(primary_charp_without(#[1..gen_num],v));
2906      }
2907    }
2908  }
2909  if (mol_flag==1)                     // the user wants no calculation of the
2910  { list L=group_reynolds(#[1..gen_num],v); // Molien series
2911    if (ch==0)
2912    { list l=primary_char0_no_molien(L[1],v);
2913      if (size(l)==2)
2914      { return(l[1],L[1],l[2]);
2915      }
2916      else
2917      { return(l[1],L[1]);
2918      }
2919    }
2920    else
2921    { if (L[1]<>0)                     // testing whether we are in the modular
2922      { list l=primary_charp_no_molien(L[1],v); // case
2923        if (size(l)==2)
2924        { return(l[1],L[1],l[2]);
2925        }
2926        else
2927        { return(l[1],L[1]);
2928        }
2929      }
2930      else                             // the modular case
2931      { if (v)
2932        { "  We can start looking for primary invariants...";
2933          "";
2934        }
2935        return(primary_charp_without(#[1..gen_num],v));
2936      }
2937    }
2938  }
2939  if (mol_flag==-1)
2940  { if (ch==0)
2941    { "ERROR:   Characteristic 0 can only be simulated in characteristic p>>0.";
2942      return();
2943    }
2944    list L=group_reynolds(#[1..gen_num],v);
2945    string newring="aksldfalkdsflkj";
2946    molien(L[2..size(L)],newring,intvec(1,interval,v));
2947    matrix P=primary_charp(L[1],newring,v);
2948    return(P,L[1],newring);
2949  }
2950  else                                 // the user specified that the
2951  { if (ch==0)                         // characteristic divides the group order
2952    { "ERROR:   The characteristic cannot divide the group order when it is 0.";
2953      return();
2954    }
2955    if (v)
2956    { "";
2957    }
2958    return(primary_charp_without(#[1..gen_num],v));
2959  }
2960}
2961example
2962{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
2963  echo=2;
2964         ring R=0,(x,y,z),dp;
2965         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
2966         list L=primary_invariants(A);
2967         print(L[1]);
2968}
2969
2970////////////////////////////////////////////////////////////////////////////////
2971// This procedure finds dif primary invariants in degree d. It returns all
2972// primary invariants found so far. The coefficients lie in a field of
2973// characteristic 0.
2974////////////////////////////////////////////////////////////////////////////////
2975proc search_random (int n,int d,ideal B,int cd,ideal P,int i,int dif,int dB,ideal CI,int max)
2976{ string answer;
2977  degBound=0;
2978  int j,k,test_dim,flag;
2979  matrix test_matrix[1][dif];          // the linear combination to test
2980  intvec h;                            // Hilbert series
2981  for (j=i+1;j<=i+dif;j=j+1)
2982  { CI=CI+ideal(var(j)^d);             // homogeneous polynomial of the same
2983                                       // degree as the one we're looking for
2984                                       // is added
2985  }
2986  ideal TEST;
2987  // h=hilb(std(CI),1);
2988  dB=dB+dif*(d-1);                     // used as degBound
2989  while (1)
2990  { test_matrix=matrix(B)*random(max,cd,dif);
2991    // degBound=dB;
2992    TEST=P+ideal(test_matrix);
2993    attrib(TEST,"isSB",1);
2994    test_dim=dim(TEST);
2995    // degBound=0;
2996    if (n-test_dim==i+dif)
2997    { break;
2998    }
2999    // degBound=dB;
3000    test_dim=dim(std(TEST));
3001    // test_dim=dim(std(TEST,h));         // Hilbert driven std-calculation
3002    // degBound=0;
3003    if (n-test_dim==i+dif)
3004    { break;
3005    }
3006    else
3007    { "HELP:    The "+string(dif)+" random combination(s) of the "+string(cd)+" basis elements with";
3008      "         coefficients in the range from -"+string(max)+" to "+string(max)+" did not lower the";
3009      "         dimension by "+string(dif)+". You can abort, try again or give a new range:";
3010      answer="";
3011      while (answer<>"n
3012" && answer<>"y
3013")
3014      { "         Do you want to abort (y/n)?";
3015        answer=read("");
3016      }
3017      if (answer=="y
3018")
3019      { flag=1;
3020        break;
3021      }
3022      answer="";
3023      while (answer<>"n
3024" && answer<>"y
3025")
3026      { "         Do you want to try again (y/n)?";
3027        answer=read("");
3028      }
3029      if (answer=="n
3030")
3031      { flag=1;
3032        while (flag)
3033        { "         Give a new <int> > "+string(max)+" that bounds the range of coefficients:";
3034          answer=read("");
3035          for (j=1;j<=size(answer)-1;j=j+1)
3036          { for (k=0;k<=9;k=k+1)
3037            { if (answer[j]==string(k))
3038              { break;
3039              }
3040            }
3041            if (k>9)
3042            { flag=1;
3043              break;
3044            }
3045            flag=0;
3046          }
3047          if (not(flag))
3048          { execute "test_dim="+string(answer[1..size(answer)]);
3049            if (test_dim<=max)
3050            { flag=1;
3051            }
3052            else
3053            { max=test_dim;
3054            }
3055          }
3056        }
3057      }
3058    }
3059  }
3060  if (not(flag))
3061  { P[(i+1)..(i+dif)]=test_matrix[1,1..dif];
3062  }
3063  return(P,CI,dB);
3064}
3065
3066////////////////////////////////////////////////////////////////////////////////
3067// This procedure finds at most dif primary invariants in degree d. It returns
3068// all primary invariants found so far. The coefficients lie in the field of
3069// characteristic p>0.
3070////////////////////////////////////////////////////////////////////////////////
3071proc p_search_random (int n,int d,ideal B,int cd,ideal P,int i,int dif,int dB,ideal CI,int max)
3072{ string answer;
3073  degBound=0;
3074  int j,k,test_dim,flag;
3075  matrix test_matrix[1][dif];          // the linear combination to test
3076  intvec h;                            // Hilbert series
3077  ideal TEST;
3078  while (dif>0)
3079  { for (j=i+1;j<=i+dif;j=j+1)
3080    { CI=CI+ideal(var(j)^d);           // homogeneous polynomial of the same
3081                                       // degree as the one we're looking for
3082                                       // is added
3083    }
3084    // h=hilb(std(CI),1);
3085    dB=dB+dif*(d-1);                   // used as degBound
3086    test_matrix=matrix(B)*random(max,cd,dif);
3087    // degBound=dB;
3088    TEST=P+ideal(test_matrix);
3089    attrib(TEST,"isSB",1);
3090    test_dim=dim(TEST);
3091    // degBound=0;
3092    if (n-test_dim==i+dif)
3093    { break;
3094    }
3095    // degBound=dB;
3096    test_dim=dim(std(TEST));
3097    // test_dim=dim(std(TEST,h));         // Hilbert driven std-calculation
3098    // degBound=0;
3099    if (n-test_dim==i+dif)
3100    { break;
3101    }
3102    else
3103    { "HELP:    The "+string(dif)+" random combination(s) of the "+string(cd)+" basis elements with";
3104      "         coefficients in the range from -"+string(max)+" to "+string(max)+" did not lower the";
3105      "         dimension by "+string(dif)+". You can abort, try again, lower the number of";
3106      "         combinations searched for by 1 or give a larger coefficient range:";
3107      answer="";
3108      while (answer<>"n
3109" && answer<>"y
3110")
3111      { "         Do you want to abort (y/n)?";
3112        answer=read("");
3113      }
3114      if (answer=="y
3115")
3116      { flag=1;
3117        break;
3118      }
3119      answer="";
3120      while (answer<>"n
3121" && answer<>"y
3122")
3123      { "         Do you want to try again (y/n)?";
3124        answer=read("");
3125      }
3126      if (answer=="n
3127")
3128      { answer="";
3129        while (answer<>"n
3130" && answer<>"y
3131")
3132        { "         Do you want to lower the number of combinations by 1 (y/n)?";
3133          answer=read("");
3134        }
3135        if (answer=="y
3136")
3137        { dif=dif-1;
3138        }
3139        else
3140        { flag=1;
3141          while (flag)
3142          { "         Give a new <int> > "+string(max)+" that bounds the range of coefficients:";
3143            answer=read("");
3144            for (j=1;j<=size(answer)-1;j=j+1)
3145            { for (k=0;k<=9;k=k+1)
3146              { if (answer[j]==string(k))
3147                { break;
3148                }
3149              }
3150              if (k>9)
3151              { flag=1;
3152                break;
3153              }
3154              flag=0;
3155            }
3156            if (not(flag))
3157            { execute "test_dim="+string(answer[1..size(answer)]);
3158              if (test_dim<=max)
3159              { flag=1;
3160              }
3161              else
3162              { max=test_dim;
3163              }
3164            }
3165          }
3166        }
3167      }
3168    }
3169    CI=CI[1..i];
3170    dB=dB-dif*(d-1);
3171  }
3172  if (dif && not(flag))
3173  { P[(i+1)..(i+dif)]=test_matrix[1,1..dif];
3174  }
3175  if (dif && flag)
3176  { P[n+1]=0;
3177  }
3178  return(P,CI,dB);
3179}
3180
3181proc primary_char0_random (matrix REY,matrix M,int max,list #)
3182"USAGE:   primary_char0_random(REY,M,r[,v]);
3183         REY: a <matrix> representing the Reynolds operator, M: a 1x2 <matrix>
3184         representing the Molien series, r: an <int> where -|r| to |r| is the
3185         range of coefficients of the random combinations of bases elements,
3186         v: an optional <int>
3187ASSUME:  REY is the first return value of group_reynolds or reynolds_molien and
3188         M the one of molien or the second one of reynolds_molien
3189DISPLAY: information about the various stages of the programme if v does not
3190         equal 0
3191RETURN:  primary invariants (type <matrix>) of the invariant ring
3192EXAMPLE: example primary_char0_random; shows an example
3193THEORY:  Bases of homogeneous invariants are generated successively and random
3194         linear combinations are chosen as primary invariants that lower the
3195         dimension of the ideal generated by the previously found invariants
3196         (see paper \"Generating a Noetherian Normalization of the Invariant Ring
3197         of a Finite Group\" by Decker, Heydtmann, Schreyer (1997) to appear in
3198         JSC).
3199"
3200{ degBound=0;
3201  if (char(basering)<>0)
3202  { "ERROR:   primary_char0_random should only be used with rings of";
3203    "         characteristic 0.";
3204    return();
3205  }
3206 //----------------- checking input and setting verbose mode ------------------
3207  if (size(#)>1)
3208  { "ERROR:   primary_char0_random can only have four parameters.";
3209    return();
3210  }
3211  if (size(#)==1)
3212  { if (typeof(#[1])<>"int")
3213    { "ERROR:   The fourth parameter should be of type <int>.";
3214      return();
3215    }
3216    else
3217    { int v=#[1];
3218    }
3219  }
3220  else
3221  { int v=0;
3222  }
3223  int n=nvars(basering);               // n is the number of variables, as well
3224                                       // as the size of the matrices, as well
3225                                       // as the number of primary invariants,
3226                                       // we should get
3227  if (ncols(REY)<>n)
3228  { "ERROR:   First parameter ought to be the Reynolds operator."
3229    return();
3230  }
3231  if (ncols(M)<>2 or nrows(M)<>1)
3232  { "ERROR:   Second parameter ought to be the Molien series."
3233    return();
3234  }
3235 //----------------------------------------------------------------------------
3236  if (v && voice<>2)
3237  { "  We can start looking for primary invariants...";
3238    "";
3239  }
3240  if (v && voice==2)
3241  { "";
3242  }
3243 //------------------------- initializing variables ---------------------------
3244  int dB;
3245  poly p(1..2);                        // p(1) will be used for single terms of
3246                                       // the partial expansion, p(2) to store
3247  p(1..2)=partial_molien(M,1);         // the intermediate result -
3248  poly v1=var(1);                      // we need v1 to split off coefficients
3249                                       // in the partial expansion of M (which
3250                                       // is in terms of the first variable) -
3251  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
3252                                       // space of invariants of degree d,
3253                                       // newdim: dimension the ideal generated
3254                                       // the primary invariants plus basis
3255                                       // elements, dif=n-i-newdim, i.e. the
3256                                       // number of new primary invairants that
3257                                       // should be added in this degree -
3258  ideal P,Pplus,CI,B;                  // P: will contain primary invariants,
3259                                       // Pplus: P+B,CI: a complete
3260                                       // intersection with the same Hilbert
3261                                       // function as P -
3262  dB=1;                                // used as degree bound
3263  int i=0;
3264 //-------------- loop that searches for primary invariants  ------------------
3265  while(1)                             // repeat until n primary invariants are
3266  {                                    // found -
3267    p(1..2)=partial_molien(M,1,p(2));  // next term of the partial expansion -
3268    d=deg(p(1));                       // degree where we'll search -
3269    cd=int(coef(p(1),v1)[2,1]);        // dimension of the homogeneous space of
3270                                       // inviarants of degree d
3271    if (v)
3272    { "  Computing primary invariants in degree "+string(d)+":";
3273    }
3274    B=invariant_basis_reynolds(REY,d,intvec(cd,6)); // basis of invariants of
3275                                       // degree d
3276    if (B[1]<>0)
3277    { Pplus=P+B;
3278      newdim=dim(std(Pplus));
3279      dif=n-i-newdim;
3280    }
3281    else
3282    { dif=0;
3283    }
3284    if (dif<>0)                        // we have to find dif new primary
3285    {                                  // invariants
3286      if (cd<>dif)
3287      { P,CI,dB=search_random(n,d,B,cd,P,i,dif,dB,CI,max); // searching for
3288      }                                // dif invariants -
3289      else                             // i.e. we can take all of B
3290      { for(j=i+1;j>i+dif;j=j+1)
3291        { CI=CI+ideal(var(j)^d);
3292        }
3293        dB=dB+dif*(d-1);
3294        P=Pplus;
3295      }
3296      if (ncols(P)==i)
3297      { "WARNING: The return value is not a set of primary invariants, but";
3298        "         polynomials qualifying as the first "+string(i)+" primary invariants.";
3299        return(matrix(P));
3300      }
3301      if (v)
3302      { for (j=1;j<=dif;j=j+1)
3303        { "  We find: "+string(P[i+j]);
3304        }
3305      }
3306      i=i+dif;
3307      if (i==n)                        // found all primary invariants
3308      { if (v)
3309        { "";
3310          "  We found all primary invariants.";
3311          "";
3312        }
3313        return(matrix(P));
3314      }
3315    }                                  // done with degree d
3316  }
3317}
3318example
3319{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
3320  echo=2;
3321         ring R=0,(x,y,z),dp;
3322         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
3323         matrix REY,M=reynolds_molien(A);
3324         matrix P=primary_char0_random(REY,M,1);
3325         print(P);
3326}
3327
3328proc primary_charp_random (matrix REY,string ring_name,int max,list #)
3329"USAGE:   primary_charp_random(REY,ringname,r[,v]);
3330         REY: a <matrix> representing the Reynolds operator, ringname: a
3331         <string> giving the name of a ring where the Molien series is stored,
3332         r: an <int> where -|r| to |r| is the range of coefficients of the
3333         random combinations of bases elements, v: an optional <int>
3334ASSUME:  REY is the first return value of group_reynolds or reynolds_molien and
3335         ringname gives the name of a ring of characteristic 0 that has been
3336         created by molien or reynolds_molien
3337DISPLAY: information about the various stages of the programme if v does not
3338         equal 0
3339RETURN:  primary invariants (type <matrix>) of the invariant ring
3340EXAMPLE: example primary_charp_random; shows an example
3341THEORY:  Bases of homogeneous invariants are generated successively and random
3342         linear combinations are chosen as primary invariants that lower the
3343         dimension of the ideal generated by the previously found invariants
3344         (see paper \"Generating a Noetherian Normalization of the Invariant Ring
3345         of a Finite Group\" by Decker, Heydtmann, Schreyer (1997) to appear in
3346         JSC).
3347"
3348{ degBound=0;
3349 // ---------------- checking input and setting verbose mode ------------------
3350  if (char(basering)==0)
3351  { "ERROR:   primary_charp_random should only be used with rings of";
3352    "         characteristic p>0.";
3353    return();
3354  }
3355  if (size(#)>1)
3356  { "ERROR:   primary_charp_random can only have four parameters.";
3357    return();
3358  }
3359  if (size(#)==1)
3360  { if (typeof(#[1])<>"int")
3361    { "ERROR:   The fourth parameter should be of type <int>.";
3362      return();
3363    }
3364    else
3365    { int v=#[1];
3366    }
3367  }
3368  else
3369  { int v=0;
3370  }
3371  def br=basering;
3372  int n=nvars(br);                     // n is the number of variables, as well
3373                                       // as the size of the matrices, as well
3374                                       // as the number of primary invariants,
3375                                       // we should get
3376  if (ncols(REY)<>n)
3377  { "ERROR:   First parameter ought to be the Reynolds operator."
3378    return();
3379  }
3380  if (typeof(`ring_name`)<>"ring")
3381  { "ERROR:   Second parameter ought to the name of a ring where the Molien";
3382    "         is stored.";
3383    return();
3384  }
3385 //----------------------------------------------------------------------------
3386  if (v && voice<>2)
3387  { "  We can start looking for primary invariants...";
3388    "";
3389  }
3390  if (v && voice==2)
3391  { "";
3392  }
3393 //----------------------- initializing variables -----------------------------
3394  int dB;
3395  setring `ring_name`;                 // the Molien series is stores here -
3396  poly p(1..2);                        // p(1) will be used for single terms of
3397                                       // the partial expansion, p(2) to store
3398  p(1..2)=partial_molien(M,1);         // the intermediate result -
3399  poly v1=var(1);                      // we need v1 to split off coefficients
3400                                       // in the partial expansion of M (which
3401                                       // is in terms of the first variable)
3402  setring br;
3403  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
3404                                       // space of invariants of degree d,
3405                                       // newdim: dimension the ideal generated
3406                                       // the primary invariants plus basis
3407                                       // elements, dif=n-i-newdim, i.e. the
3408                                       // number of new primary invairants that
3409                                       // should be added in this degree -
3410  ideal P,Pplus,CI,B;                  // P: will contain primary invariants,
3411                                       // Pplus: P+B, CI: a complete
3412                                       // intersection with the same Hilbert
3413                                       // function as P -
3414  dB=1;                                // used as degree bound
3415  int i=0;
3416 //---------------- loop that searches for primary invariants -----------------
3417  while(1)                             // repeat until n primary invariants are
3418  {                                    // found
3419    setring `ring_name`;
3420    p(1..2)=partial_molien(M,1,p(2));  // next term of the partial expansion -
3421    d=deg(p(1));                       // degree where we'll search -
3422    cd=int(coef(p(1),v1)[2,1]);        // dimension of the homogeneous space of
3423                                       // inviarants of degree d
3424    setring br;
3425    if (v)
3426    { "  Computing primary invariants in degree "+string(d)+":";
3427    }
3428    B=invariant_basis_reynolds(REY,d,intvec(cd,6)); // basis of invariants of
3429                                       // degree d
3430    if (B[1]<>0)
3431    { Pplus=P+B;
3432      newdim=dim(std(Pplus));
3433      dif=n-i-newdim;
3434    }
3435    else
3436    { dif=0;
3437    }
3438    if (dif<>0)                        // we have to find dif new primary
3439    {                                  // invariants
3440      if (cd<>dif)
3441      { P,CI,dB=p_search_random(n,d,B,cd,P,i,dif,dB,CI,max);
3442      }
3443      else                             // i.e. we can take all of B
3444      { for(j=i+1;j>i+dif;j=j+1)
3445        { CI=CI+ideal(var(j)^d);
3446        }
3447        dB=dB+dif*(d-1);
3448        P=Pplus;
3449      }
3450      if (ncols(P)==n+1)
3451      { "WARNING: The first return value is not a set of primary invariants,";
3452        "         but polynomials qualifying as the first "+string(i)+" primary invariants.";
3453        return(matrix(P));
3454      }
3455      if (v)
3456      { for (j=1;j<=size(P)-i;j=j+1)
3457        { "  We find: "+string(P[i+j]);
3458        }
3459      }
3460      i=size(P);
3461      if (i==n)                        // found all primary invariants
3462      { if (v)
3463        { "";
3464          "  We found all primary invariants.";
3465          "";
3466        }
3467        return(matrix(P));
3468      }
3469    }                                  // done with degree d
3470  }
3471}
3472example
3473{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed into";
3474  "         characteristic 3)";
3475  echo=2;
3476         ring R=3,(x,y,z),dp;
3477         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
3478         list L=group_reynolds(A);
3479         string newring="alskdfj";
3480         molien(L[2..size(L)],newring);
3481         matrix P=primary_charp_random(L[1],newring,1);
3482         if(system("with","Namespaces")) { kill Top::`newring`; }
3483         kill `newring`;
3484         print(P);
3485}
3486
3487proc primary_char0_no_molien_random (matrix REY, int max, list #)
3488"USAGE:   primary_char0_no_molien_random(REY,r[,v]);
3489         REY: a <matrix> representing the Reynolds operator, r: an <int> where
3490         -|r| to |r| is the range of coefficients of the random combinations of
3491         bases elements, v: an optional <int>
3492ASSUME:  REY is the first return value of group_reynolds or reynolds_molien
3493DISPLAY: information about the various stages of the programme if v does not
3494         equal 0
3495RETURN:  primary invariants (type <matrix>) of the invariant ring  and an
3496         <intvec> listing some of the degrees where no non-trivial homogeneous
3497         invariants are to be found
3498EXAMPLE: example primary_char0_no_molien_random; shows an example
3499THEORY:  Bases of homogeneous invariants are generated successively and random
3500         linear combinations are chosen as primary invariants that lower the
3501         dimension of the ideal generated by the previously found invariants
3502         (see paper \"Generating a Noetherian Normalization of the Invariant Ring
3503         of a Finite Group\" by Decker, Heydtmann, Schreyer (1997) to appear in
3504         JSC).
3505"
3506{ degBound=0;
3507 //-------------- checking input and setting verbose mode ---------------------
3508  if (char(basering)<>0)
3509  { "ERROR:   primary_char0_no_molien_random should only be used with rings of";
3510    "         characteristic 0.";
3511    return();
3512  }
3513  if (size(#)>1)
3514  { "ERROR:   primary_char0_no_molien_random can only have three parameters.";
3515    return();
3516  }
3517  if (size(#)==1)
3518  { if (typeof(#[1])<>"int")
3519    { "ERROR:   The third parameter should be of type <int>.";
3520      return();
3521    }
3522    else
3523    { int v=#[1];
3524    }
3525  }
3526  else
3527  { int v=0;
3528  }
3529  int n=nvars(basering);               // n is the number of variables, as well
3530                                       // as the size of the matrices, as well
3531                                       // as the number of primary invariants,
3532                                       // we should get
3533  if (ncols(REY)<>n)
3534  { "ERROR:   First parameter ought to be the Reynolds operator."
3535    return();
3536  }
3537 //----------------------------------------------------------------------------
3538  if (v && voice<>2)
3539  { "  We can start looking for primary invariants...";
3540    "";
3541  }
3542  if (v && voice==2)
3543  { "";
3544  }
3545 //----------------------- initializing variables -----------------------------
3546  int dB;
3547  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
3548                                       // space of invariants of degree d,
3549                                       // newdim: dimension the ideal generated
3550                                       // the primary invariants plus basis
3551                                       // elements, dif=n-i-newdim, i.e. the
3552                                       // number of new primary invairants that
3553                                       // should be added in this degree -
3554  ideal P,Pplus,CI,B;                  // P: will contain primary invariants,
3555                                       // Pplus: P+B, CI: a complete
3556                                       // intersection with the same Hilbert
3557                                       // function as P -
3558  dB=1;                                // used as degree bound -
3559  d=0;                                 // initializing
3560  int i=0;
3561  intvec deg_vector;
3562 //------------------ loop that searches for primary invariants ---------------
3563  while(1)                             // repeat until n primary invariants are
3564  {                                    // found -
3565    d=d+1;                             // degree where we'll search
3566    if (v)
3567    { "  Computing primary invariants in degree "+string(d)+":";
3568    }
3569    B=invariant_basis_reynolds(REY,d,intvec(-1,6)); // basis of invariants of
3570                                       // degree d
3571    if (B[1]<>0)
3572    { Pplus=P+B;
3573      newdim=dim(std(Pplus));
3574      dif=n-i-newdim;
3575    }
3576    else
3577    { dif=0;
3578      deg_vector=deg_vector,d;
3579    }
3580    if (dif<>0)                        // we have to find dif new primary
3581    {                                  // invariants
3582      cd=size(B);
3583      if (cd<>dif)
3584      { P,CI,dB=search_random(n,d,B,cd,P,i,dif,dB,CI,max);
3585      }
3586      else                             // i.e. we can take all of B
3587      { for(j=i+1;j<=i+dif;j=j+1)
3588        { CI=CI+ideal(var(j)^d);
3589        }
3590        dB=dB+dif*(d-1);
3591        P=Pplus;
3592      }
3593      if (ncols(P)==i)
3594      { "WARNING: The first return value is not a set of primary invariants,";
3595        "         but polynomials qualifying as the first "+string(i)+" primary invariants.";
3596        return(matrix(P));
3597      }
3598      if (v)
3599      { for (j=1;j<=dif;j=j+1)
3600        { "  We find: "+string(P[i+j]);
3601        }
3602      }
3603      i=i+dif;
3604      if (i==n)                        // found all primary invariants
3605      { if (v)
3606        { "";
3607          "  We found all primary invariants.";
3608          "";
3609        }
3610        if (deg_vector==0)
3611        { return(matrix(P));
3612        }
3613        else
3614        { return(matrix(P),compress(deg_vector));
3615        }
3616      }
3617    }                                  // done with degree d
3618    else
3619    { if (v)
3620      { "  None here...";
3621      }
3622    }
3623  }
3624}
3625example
3626{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
3627  echo=2;
3628         ring R=0,(x,y,z),dp;
3629         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
3630         list L=group_reynolds(A);
3631         list l=primary_char0_no_molien_random(L[1],1);
3632         print(l[1]);
3633}
3634
3635proc primary_charp_no_molien_random (matrix REY, int max, list #)
3636"USAGE:   primary_charp_no_molien_random(REY,r[,v]);
3637         REY: a <matrix> representing the Reynolds operator, r: an <int> where
3638         -|r| to |r| is the range of coefficients of the random combinations of
3639         bases elements, v: an optional <int>
3640ASSUME:  REY is the first return value of group_reynolds or reynolds_molien
3641DISPLAY: information about the various stages of the programme if v does not
3642         equal 0
3643RETURN:  primary invariants (type <matrix>) of the invariant ring  and an
3644         <intvec> listing some of the degrees where no non-trivial homogeneous
3645         invariants are to be found
3646EXAMPLE: example primary_charp_no_molien_random; shows an example
3647THEORY:  Bases of homogeneous invariants are generated successively and random
3648         linear combinations are chosen as primary invariants that lower the
3649         dimension of the ideal generated by the previously found invariants
3650         (see paper \"Generating a Noetherian Normalization of the Invariant Ring
3651         of a Finite Group\" by Decker, Heydtmann, Schreyer (1997) to appear in
3652         JSC).
3653"
3654{ degBound=0;
3655 //----------------- checking input and setting verbose mode ------------------
3656  if (char(basering)==0)
3657  { "ERROR:   primary_charp_no_molien_random should only be used with rings of";
3658    "         characteristic p>0.";
3659    return();
3660  }
3661  if (size(#)>1)
3662  { "ERROR:   primary_charp_no_molien_random can only have three parameters.";
3663    return();
3664  }
3665  if (size(#)==1)
3666  { if (typeof(#[1])<>"int")
3667    { "ERROR:   The third parameter should be of type <int>.";
3668      return();
3669    }
3670    else
3671    { int v=#[1];
3672    }
3673  }
3674  else
3675  { int v=0;
3676  }
3677  int n=nvars(basering);               // n is the number of variables, as well
3678                                       // as the size of the matrices, as well
3679                                       // as the number of primary invariants,
3680                                       // we should get
3681  if (ncols(REY)<>n)
3682  { "ERROR:   First parameter ought to be the Reynolds operator."
3683    return();
3684  }
3685 //----------------------------------------------------------------------------
3686  if (v && voice<>2)
3687  { "  We can start looking for primary invariants...";
3688    "";
3689  }
3690  if (v && voice==2)
3691  { "";
3692  }
3693 //-------------------- initializing variables --------------------------------
3694  int dB;
3695  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
3696                                       // space of invariants of degree d,
3697                                       // newdim: dimension the ideal generated
3698                                       // the primary invariants plus basis
3699                                       // elements, dif=n-i-newdim, i.e. the
3700                                       // number of new primary invairants that
3701                                       // should be added in this degree -
3702  ideal P,Pplus,CI,B;                  // P: will contain primary invariants,
3703                                       // Pplus: P+B, CI: a complete
3704                                       // intersection with the same Hilbert
3705                                       // function as P -
3706  dB=1;                                // used as degree bound -
3707  d=0;                                 // initializing
3708  int i=0;
3709  intvec deg_vector;
3710 //------------------ loop that searches for primary invariants ---------------
3711  while(1)                             // repeat until n primary invariants are
3712  {                                    // found -
3713    d=d+1;                             // degree where we'll search
3714    if (v)
3715    { "  Computing primary invariants in degree "+string(d)+":";
3716    }
3717    B=invariant_basis_reynolds(REY,d,intvec(-1,6)); // basis of invariants of
3718                                       // degree d
3719    if (B[1]<>0)
3720    { Pplus=P+B;
3721      newdim=dim(std(Pplus));
3722      dif=n-i-newdim;
3723    }
3724    else
3725    { dif=0;
3726      deg_vector=deg_vector,d;
3727    }
3728    if (dif<>0)                        // we have to find dif new primary
3729    {                                  // invariants
3730      cd=size(B);
3731      if (cd<>dif)
3732      { P,CI,dB=p_search_random(n,d,B,cd,P,i,dif,dB,CI,max);
3733      }
3734      else                             // i.e. we can take all of B
3735      { for(j=i+1;j<=i+dif;j=j+1)
3736        { CI=CI+ideal(var(j)^d);
3737        }
3738        dB=dB+dif*(d-1);
3739        P=Pplus;
3740      }
3741      if (ncols(P)==n+1)
3742      { "WARNING: The first return value is not a set of primary invariants,";
3743        "         but polynomials qualifying as the first "+string(i)+" primary invariants.";
3744        return(matrix(P));
3745      }
3746      if (v)
3747      { for (j=1;j<=size(P)-i;j=j+1)
3748        { "  We find: "+string(P[i+j]);
3749        }
3750      }
3751      i=size(P);
3752      if (i==n)                        // found all primary invariants
3753      { if (v)
3754        { "";
3755          "  We found all primary invariants.";
3756          "";
3757        }
3758        if (deg_vector==0)
3759        { return(matrix(P));
3760        }
3761        else
3762        { return(matrix(P),compress(deg_vector));
3763        }
3764      }
3765    }                                  // done with degree d
3766    else
3767    { if (v)
3768      { "  None here...";
3769      }
3770    }
3771  }
3772}
3773example
3774{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed into";
3775  "         characteristic 3)";
3776  echo=2;
3777         ring R=3,(x,y,z),dp;
3778         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
3779         list L=group_reynolds(A);
3780         list l=primary_charp_no_molien_random(L[1],1);
3781         print(l[1]);
3782}
3783
3784proc primary_charp_without_random (list #)
3785"USAGE:   primary_charp_without_random(G1,G2,...,r[,v]);
3786         G1,G2,...: <matrices> generating a finite matrix group, r: an <int>
3787         where -|r| to |r| is the range of coefficients of the random
3788         combinations of bases elements, v: an optional <int>
3789DISPLAY: information about the various stages of the programme if v does not
3790         equal 0
3791RETURN:  primary invariants (type <matrix>) of the invariant ring
3792EXAMPLE: example primary_charp_without_random; shows an example
3793THEORY:  Bases of homogeneous invariants are generated successively and random
3794         linear combinations are chosen as primary invariants that lower the
3795         dimension of the ideal generated by the previously found invariants
3796         (see paper \"Generating a Noetherian Normalization of the Invariant Ring
3797         of a Finite Group\" by Decker, Heydtmann, Schreyer (1997) to appear in
3798         JSC). No Reynolds operator or Molien series is used.
3799"
3800{ degBound=0;
3801 //--------------------- checking input and setting verbose mode --------------
3802  if (char(basering)==0)
3803  { "ERROR:   primary_charp_without_random should only be used with rings of";
3804    "         characteristic 0.";
3805    return();
3806  }
3807  if (size(#)<2)
3808  { "ERROR:   There are too few parameters.";
3809    return();
3810  }
3811  if (typeof(#[size(#)])=="int" && typeof(#[size(#)-1])=="int")
3812  { int v=#[size(#)];
3813    int max=#[size(#)-1];
3814    int gen_num=size(#)-2;
3815    if (gen_num==0)
3816    { "ERROR:   There are no generators of a finite matrix group given.";
3817      return();
3818    }
3819  }
3820  else
3821  { if (typeof(#[size(#)])=="int")
3822    { int max=#[size(#)];
3823      int v=0;
3824      int gen_num=size(#)-1;
3825    }
3826    else
3827    { "ERROR:   The last parameter should be an <int>.";
3828      return();
3829    }
3830  }
3831  int n=nvars(basering);               // n is the number of variables, as well
3832                                       // as the size of the matrices, as well
3833                                       // as the number of primary invariants,
3834                                       // we should get
3835  for (int i=1;i<=gen_num;i=i+1)
3836  { if (typeof(#[i])=="matrix")
3837    { if (nrows(#[i])<>n or ncols(#[i])<>n)
3838      { "ERROR:   The number of variables of the base ring needs to be the same";
3839        "         as the dimension of the square matrices";
3840        return();
3841      }
3842    }
3843    else
3844    { "ERROR:   The first parameters should be a list of matrices";
3845      return();
3846    }
3847  }
3848 //----------------------------------------------------------------------------
3849  if (v && voice==2)
3850  { "";
3851  }
3852 //---------------------------- initializing variables ------------------------
3853  int dB;
3854  int j,d,cd,newdim,dif;               // d: current degree, cd: dimension of
3855                                       // space of invariants of degree d,
3856                                       // newdim: dimension the ideal generated
3857                                       // the primary invariants plus basis
3858                                       // elements, dif=n-i-newdim, i.e. the
3859                                       // number of new primary invairants that
3860                                       // should be added in this degree -
3861  ideal P,Pplus,CI,B;                  // P: will contain primary invariants,
3862                                       // Pplus: P+B, CI: a complete
3863                                       // intersection with the same Hilbert
3864                                       // function as P -
3865  dB=1;                                // used as degree bound -
3866  d=0;                                 // initializing
3867  i=0;
3868  intvec deg_vector;
3869 //-------------------- loop that searches for primary invariants -------------
3870  while(1)                             // repeat until n primary invariants are
3871  {                                    // found -
3872    d=d+1;                             // degree where we'll search
3873    if (v)
3874    { "  Computing primary invariants in degree "+string(d)+":";
3875    }
3876    B=invariant_basis(d,#[1..gen_num]); // basis of invariants of degree d
3877    if (B[1]<>0)
3878    { Pplus=P+B;
3879      newdim=dim(std(Pplus));
3880      dif=n-i-newdim;
3881    }
3882    else
3883    { dif=0;
3884      deg_vector=deg_vector,d;
3885    }
3886    if (dif<>0)                        // we have to find dif new primary
3887    {                                  // invariants
3888      cd=size(B);
3889      if (cd<>dif)
3890      { P,CI,dB=p_search_random(n,d,B,cd,P,i,dif,dB,CI,max);
3891      }
3892      else                             // i.e. we can take all of B
3893      { for(j=i+1;j<=i+dif;j=j+1)
3894        { CI=CI+ideal(var(j)^d);
3895        }
3896        dB=dB+dif*(d-1);
3897        P=Pplus;
3898      }
3899      if (ncols(P)==n+1)
3900      { "WARNING: The first return value is not a set of primary invariants,";
3901        "         but polynomials qualifying as the first "+string(i)+" primary invariants.";
3902        return(matrix(P));
3903      }
3904      if (v)
3905      { for (j=1;j<=size(P)-i;j=j+1)
3906        { "  We find: "+string(P[i+j]);
3907        }
3908      }
3909      i=size(P);
3910      if (i==n)                        // found all primary invariants
3911      { if (v)
3912        { "";
3913          "  We found all primary invariants.";
3914          "";
3915        }
3916        return(matrix(P));
3917      }
3918    }                                  // done with degree d
3919    else
3920    { if (v)
3921      { "  None here...";
3922      }
3923    }
3924  }
3925}
3926example
3927{ echo=2;
3928         ring R=2,(x,y,z),dp;
3929         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
3930         matrix P=primary_charp_without_random(A,1);
3931         print(P);
3932}
3933
3934proc primary_invariants_random (list #)
3935"USAGE:   primary_invariants_random(G1,G2,...,r[,flags]);
3936         G1,G2,...: <matrices> generating a finite matrix group, r: an <int>
3937         where -|r| to |r| is the range of coefficients of the random
3938         combinations of bases elements, flags: an optional <intvec> with three
3939         entries, if the first one equals 0 (also the default), the programme
3940         attempts to compute the Molien series and Reynolds operator, if it
3941         equals 1, the programme is told that the Molien series should not be
3942         computed, if it equals -1 characteristic 0 is simulated, i.e. the
3943         Molien series is computed as if the base field were characteristic 0
3944         (the user must choose a field of large prime characteristic, e.g.
3945         32003) and if the first one is anything else, it means that the
3946         characteristic of the base field divides the group order, the second
3947         component should give the size of intervals between canceling common
3948         factors in the expansion of the Molien series, 0 (the default) means
3949         only once after generating all terms, in prime characteristic also a
3950         negative number can be given to indicate that common factors should
3951         always be canceled when the expansion is simple (the root of the
3952         extension field does not occur among the coefficients)
3953DISPLAY: information about the various stages of the programme if the third
3954         flag does not equal 0
3955RETURN:  primary invariants (type <matrix>) of the invariant ring and if
3956         computable Reynolds operator (type <matrix>) and Molien series (type
3957         <matrix>), if the first flag is 1 and we are in the non-modular case
3958         then an <intvec> is returned giving some of the degrees where no
3959         non-trivial homogeneous invariants can be found
3960EXAMPLE: example primary_invariants_random; shows an example
3961THEORY:  Bases of homogeneous invariants are generated successively and random
3962         linear combinations are chosen as primary invariants that lower the
3963         dimension of the ideal generated by the previously found invariants
3964         (see paper \"Generating a Noetherian Normalization of the Invariant Ring
3965         of a Finite Group\" by Decker, Heydtmann, Schreyer (1997) to appear in
3966         JSC).
3967"
3968{
3969 // ----------------- checking input and setting flags ------------------------
3970  if (size(#)<2)
3971  { "ERROR:   There are too few parameters.";
3972    return();
3973  }
3974  int ch=char(basering);               // the algorithms depend very much on the
3975                                       // characteristic of the ground field
3976  int n=nvars(basering);               // n is the number of variables, as well
3977                                       // as the size of the matrices, as well
3978                                       // as the number of primary invariants,
3979                                       // we should get
3980  int gen_num;
3981  int mol_flag,v;
3982  if (typeof(#[size(#)])=="intvec" && typeof(#[size(#)-1])=="int")
3983  { if (size(#[size(#)])<>3)
3984    { "ERROR:   <intvec> should have three entries.";
3985      return();
3986    }
3987    gen_num=size(#)-2;
3988    mol_flag=#[size(#)][1];
3989    if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0)))
3990    { "ERROR:   the second component of <intvec> should be >=0";
3991      return();
3992    }
3993    int interval=#[size(#)][2];
3994    v=#[size(#)][3];
3995    int max=#[size(#)-1];
3996    if (gen_num==0)
3997    { "ERROR:   There are no generators of a finite matrix group given.";
3998      return();
3999    }
4000  }
4001  else
4002  { if (typeof(#[size(#)])=="int")
4003    { gen_num=size(#)-1;
4004      mol_flag=0;
4005      int interval=0;
4006      v=0;
4007      int max=#[size(#)];
4008    }
4009    else
4010    { "ERROR:   If the two last parameters are not <int> and <intvec>, the last";
4011      "         parameter should be an <int>.";
4012      return();
4013    }
4014  }
4015  for (int i=1;i<=gen_num;i=i+1)
4016  { if (typeof(#[i])=="matrix")
4017    { if (nrows(#[i])<>n or ncols(#[i])<>n)
4018      { "ERROR:   The number of variables of the base ring needs to be the same";
4019        "         as the dimension of the square matrices";
4020        return();
4021      }
4022    }
4023    else
4024    { "ERROR:   The first parameters should be a list of matrices";
4025      return();
4026    }
4027  }
4028 //----------------------------------------------------------------------------
4029  if (mol_flag==0)
4030  { if (ch==0)
4031    { matrix REY,M=reynolds_molien(#[1..gen_num],intvec(0,interval,v));
4032                                       // one will contain Reynolds operator and
4033                                       // the other enumerator and denominator
4034                                       // of Molien series
4035      matrix P=primary_char0_random(REY,M,max,v);
4036      return(P,REY,M);
4037    }
4038    else
4039    { list L=group_reynolds(#[1..gen_num],v);
4040      if (L[1]<>0)                     // testing whether we are in the modular
4041      { string newring="aksldfalkdsflkj"; // case
4042        if (minpoly==0)
4043        { if (v)
4044          { "  We are dealing with the non-modular case.";
4045          }
4046          molien(L[2..size(L)],newring,intvec(0,interval,v));
4047          matrix P=primary_charp_random(L[1],newring,max,v);
4048          return(P,L[1],newring);
4049        }
4050        else
4051        { if (v)
4052          { "  Since it is impossible for this programme to calculate the Molien series for";
4053            "  invariant rings over extension fields of prime characteristic, we have to";
4054            "  continue without it.";
4055            "";
4056
4057          }
4058          list l=primary_charp_no_molien_random(L[1],max,v);
4059          if (size(l)==2)
4060          { return(l[1],L[1],l[2]);
4061          }
4062          else
4063          { return(l[1],L[1]);
4064          }
4065        }
4066      }
4067      else                             // the modular case
4068      { if (v)
4069        { "  There is also no Molien series, we can make use of...";
4070          "";
4071          "  We can start looking for primary invariants...";
4072          "";
4073        }
4074        return(primary_charp_without_random(#[1..gen_num],max,v));
4075      }
4076    }
4077  }
4078  if (mol_flag==1)                     // the user wants no calculation of the
4079  { list L=group_reynolds(#[1..gen_num],v); // Molien series
4080    if (ch==0)
4081    { list l=primary_char0_no_molien_random(L[1],max,v);
4082      if (size(l)==2)
4083      { return(l[1],L[1],l[2]);
4084      }
4085      else
4086      { return(l[1],L[1]);
4087      }
4088    }
4089    else
4090    { if (L[1]<>0)                     // testing whether we are in the modular
4091      { list l=primary_charp_no_molien_random(L[1],max,v); // case
4092        if (size(l)==2)
4093        { return(l[1],L[1],l[2]);
4094        }
4095        else
4096        { return(l[1],L[1]);
4097        }
4098      }
4099      else                             // the modular case
4100      { if (v)
4101        { "  We can start looking for primary invariants...";
4102          "";
4103        }
4104        return(primary_charp_without_random(#[1..gen_num],max,v));
4105      }
4106    }
4107  }
4108  if (mol_flag==-1)
4109  { if (ch==0)
4110    { "ERROR:   Characteristic 0 can only be simulated in characteristic p>>0.";
4111      return();
4112    }
4113    list L=group_reynolds(#[1..gen_num],v);
4114    string newring="aksldfalkdsflkj";
4115    molien(L[2..size(L)],newring,intvec(1,interval,v));
4116    matrix P=primary_charp_random(L[1],newring,max,v);
4117    return(P,L[1],newring);
4118  }
4119  else                                 // the user specified that the
4120  { if (ch==0)                         // characteristic divides the group order
4121    { "ERROR:   The characteristic cannot divide the group order when it is 0.";
4122      return();
4123    }
4124    if (v)
4125    { "";
4126    }
4127    return(primary_charp_without_random(#[1..gen_num],max,v));
4128  }
4129}
4130example
4131{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
4132  echo=2;
4133         ring R=0,(x,y,z),dp;
4134         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
4135         list L=primary_invariants_random(A,1);
4136         print(L[1]);
4137}
4138
4139proc concat_intmat(intmat A,intmat B)
4140{ int n=nrows(A);
4141  int m1=ncols(A);
4142  int m2=ncols(B);
4143  intmat C[n][m1+m2];
4144  C[1..n,1..m1]=A[1..n,1..m1];
4145  C[1..n,m1+1..m1+m2]=B[1..n,1..m2];
4146  return(C);
4147}
4148
4149proc power_products(intvec deg_vec,int d)
4150"USAGE:   power_products(dv,d);
4151         dv: an <intvec> giving the degrees of homogeneous polynomials, d: the
4152         degree of the desired power products
4153RETURN:  a size(dv)*m <intmat> where each column ought to be interpreted as
4154         containing the exponents of the corresponding polynomials. The product
4155         of the powers is then homogeneous of degree d.
4156EXAMPLE: example power_products; gives an example
4157"
4158{ ring R=0,x,dp;
4159  if (d<=0)
4160  { "ERROR:   The <int> may not be <= 0";
4161    return();
4162  }
4163  int d_neu,j,nc;
4164  int s=size(deg_vec);
4165  intmat PP[s][1];
4166  intmat TEST[s][1];
4167  for (int i=1;i<=s;i=i+1)
4168  { if (i<0)
4169    { "ERROR:   The entries of <intvec> may not be <= 0";
4170      return();
4171    }
4172    d_neu=d-deg_vec[i];
4173    if (d_neu>0)
4174    { intmat PPd_neu=power_products(intvec(deg_vec[i..s]),d_neu);
4175      if (size(ideal(PPd_neu))<>0)
4176      { nc=ncols(PPd_neu);
4177        intmat PPd_neu_gross[s][nc];
4178        PPd_neu_gross[i..s,1..nc]=PPd_neu[1..s-i+1,1..nc];
4179        for (j=1;j<=nc;j=j+1)
4180        { PPd_neu_gross[i,j]=PPd_neu_gross[i,j]+1;
4181        }
4182        PP=concat_intmat(PP,PPd_neu_gross);
4183        kill PPd_neu_gross;
4184      }
4185      kill PPd_neu;
4186    }
4187    if (d_neu==0)
4188    { intmat PPd_neu[s][1];
4189      PPd_neu[i,1]=1;
4190      PP=concat_intmat(PP,PPd_neu);
4191      kill PPd_neu;
4192    }
4193  }
4194  if (matrix(PP)<>matrix(TEST))
4195  { PP=compress(PP);
4196  }
4197  return(PP);
4198}
4199example
4200{ echo=2;
4201         intvec dv=5,5,5,10,10;
4202         print(power_products(dv,10));
4203         print(power_products(dv,7));
4204}
4205
4206proc secondary_char0 (matrix P, matrix REY, matrix M, list #)
4207"USAGE:   secondary_char0(P,REY,M[,v]);
4208         P: a 1xn <matrix> with primary invariants, REY: a gxn <matrix>
4209         representing the Reynolds operator, M: a 1x2 <matrix> giving enumerator
4210         and denominator of the Molien series, v: an optional <int>
4211ASSUME:  n is the number of variables of the basering, g the size of the group,
4212         REY is the 1st return value of group_reynolds(), reynolds_molien() or
4213         the second one of primary_invariants(), M the return value of molien()
4214         or the second one of reynolds_molien() or the third one of
4215         primary_invariants()
4216RETURN:  secondary invariants of the invariant ring (type <matrix>) and
4217         irreducible secondary invariants (type <matrix>)
4218DISPLAY: information if v does not equal 0
4219EXAMPLE: example secondary_char0; shows an example
4220THEORY:  The secondary invariants are calculated by finding a basis (in terms of
4221         monomials) of the basering modulo the primary invariants, mapping those
4222         to invariants with the Reynolds operator and using these images or
4223         their power products such that they are linearly independent modulo the
4224         primary invariants (see paper \"Some Algorithms in Invariant Theory of
4225         Finite Groups\" by Kemper and Steel (1997)).
4226"
4227{ def br=basering;
4228  degBound=0;
4229 //----------------- checking input and setting verbose mode ------------------
4230  if (char(br)<>0)
4231  { "ERROR:   secondary_char0 should only be used with rings of characteristic 0.";
4232    return();
4233  }
4234  int i;
4235  if (size(#)>0)
4236  { if (typeof(#[size(#)])=="int")
4237    { int v=#[size(#)];
4238    }
4239    else
4240    { int v=0;
4241    }
4242  }
4243  else
4244  { int v=0;
4245  }
4246  int n=nvars(br);                     // n is the number of variables, as well
4247                                       // as the size of the matrices, as well
4248                                       // as the number of primary invariants,
4249                                       // we should get
4250  if (ncols(P)<>n)
4251  { "ERROR:   The first parameter ought to be the matrix of the primary";
4252    "         invariants."
4253    return();
4254  }
4255  if (ncols(REY)<>n)
4256  { "ERROR:   The second parameter ought to be the Reynolds operator."
4257    return();
4258  }
4259  if (ncols(M)<>2 or nrows(M)<>1)
4260  { "ERROR:   The third parameter ought to be the Molien series."
4261    return();
4262  }
4263  if (v && voice==2)
4264  { "";
4265  }
4266  int j, m, counter;
4267 //- finding the polynomial giving number and degrees of secondary invariants -
4268  poly p=1;
4269  for (j=1;j<=n;j=j+1)                 // calculating the denominator of the
4270  { p=p*(1-var(1)^deg(P[j]));          // Hilbert series of the ring generated
4271  }                                    // by the primary invariants -
4272  matrix s[1][2]=M[1,1]*p,M[1,2];      // s is used for canceling
4273  s=matrix(syz(ideal(s)));
4274  p=s[2,1];                            // the polynomial telling us where to
4275                                       // search for secondary invariants
4276  map slead=br,ideal(0);
4277  p=1/slead(p)*p;                      // smallest term of p needs to be 1
4278  if (v)
4279  { "  Polynomial telling us where to look for secondary invariants:";
4280    "   "+string(p);
4281    "";
4282  }
4283  matrix dimmat=coeffs(p,var(1));      // dimmat will contain the number of
4284                                       // secondary invariants, we need to find
4285                                       // of a certain degree -
4286  m=nrows(dimmat);                     // m-1 is the highest degree
4287  if (v)
4288  { "  In degree 0 we have: 1";
4289    "";
4290  }
4291 //-------------------------- initializing variables --------------------------
4292  intmat PP;
4293  poly pp;
4294  int k;
4295  intvec deg_vec;
4296  ideal sP=std(ideal(P));
4297  ideal TEST,B,IS;
4298  ideal S=1;                           // 1 is the first secondary invariant -
4299 //--------------------- generating secondary invariants ----------------------
4300  for (i=2;i<=m;i=i+1)                 // going through dimmat -
4301  { if (int(dimmat[i,1])<>0)           // when it is == 0 we need to find 0
4302    {                                  // elements in the current degree (i-1)
4303      if (v)
4304      { "  Searching in degree "+string(i-1)+", we need to find "+string(int(dimmat[i,1]))+" invariant(s)...";
4305      }
4306      TEST=sP;
4307      counter=0;                       // we'll count up to degvec[i]
4308      if (IS[1]<>0)
4309      { PP=power_products(deg_vec,i-1); // finding power products of irreducible
4310      }                                // secondary invariants
4311      if (size(ideal(PP))<>0)
4312      { for (j=1;j<=ncols(PP);j=j+1)   // going through all the power products
4313        { pp=1;
4314          for (k=1;k<=nrows(PP);k=k+1)
4315          { pp=pp*IS[1,k]^PP[k,j];
4316          }
4317          if (reduce(pp,TEST)<>0)
4318          { S=S,pp;
4319            counter=counter+1;
4320            if (v)
4321            { "  We find: "+string(pp);
4322            }
4323            if (int(dimmat[i,1])<>counter)
4324            { TEST=std(TEST+ideal(NF(pp,TEST))); // should be replaced by next
4325                                                 // line soon
4326            // TEST=std(TEST,NF(pp,TEST));
4327            }
4328            else
4329            { break;
4330            }
4331          }
4332        }
4333      }
4334      if (int(dimmat[i,1])<>counter)
4335      { B=sort_of_invariant_basis(sP,REY,i-1,int(dimmat[i,1])*6); // B contains
4336                                       // images of kbase(sP,i-1) under the
4337                                       // Reynolds operator that are linearly
4338                                       // independent and that don't reduce to
4339                                       // 0 modulo sP -
4340        if (counter==0 && ncols(B)==int(dimmat[i,1])) // then we can take all of
4341        { S=S,B;                       // B
4342          IS=IS+B;
4343          if (deg_vec[1]==0)
4344          { deg_vec=i-1;
4345            if (v)
4346            { "  We find: "+string(B[1]);
4347            }
4348            for (j=2;j<=int(dimmat[i,1]);j=j+1)
4349            { deg_vec=deg_vec,i-1;
4350              if (v)
4351              { "  We find: "+string(B[j]);
4352              }
4353            }
4354          }
4355          else
4356          { for (j=1;j<=int(dimmat[i,1]);j=j+1)
4357            { deg_vec=deg_vec,i-1;
4358              if (v)
4359              { "  We find: "+string(B[j]);
4360              }
4361            }
4362          }
4363        }
4364        else
4365        { j=0;                         // j goes through all of B -
4366          while (int(dimmat[i,1])<>counter) // need to find dimmat[i,1]
4367          {                            // invariants that are linearly
4368                                       // independent modulo TEST
4369            j=j+1;
4370            if (reduce(B[j],TEST)<>0)  // B[j] should be added
4371            { S=S,B[j];
4372              IS=IS+ideal(B[j]);
4373              if (deg_vec[1]==0)
4374              { deg_vec[1]=i-1;
4375              }
4376              else
4377              { deg_vec=deg_vec,i-1;
4378              }
4379              counter=counter+1;
4380              if (v)
4381              { "  We find: "+string(B[j]);
4382              }
4383              if (int(dimmat[i,1])<>counter)
4384              { TEST=std(TEST+ideal(NF(B[j],TEST))); // should be replaced by
4385                                                     // next line
4386              // TEST=std(TEST,NF(B[j],TEST));
4387              }
4388            }
4389          }
4390        }
4391      }
4392      if (v)
4393      { "";
4394      }
4395    }
4396  }
4397  if (v)
4398  { "  We're done!";
4399    "";
4400  }
4401  return(matrix(S),matrix(IS));
4402}
4403example
4404{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
4405  echo=2;
4406         ring R=0,(x,y,z),dp;
4407         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
4408         list L=primary_invariants(A);
4409         matrix S,IS=secondary_char0(L[1..3]);
4410         print(S);
4411         print(IS);
4412}
4413
4414proc secondary_charp (matrix P, matrix REY, string ring_name, list #)
4415"USAGE:   secondary_charp(P,REY,ringname[,v]);
4416         P: a 1xn <matrix> with primary invariants, REY: a gxn <matrix>
4417         representing the Reynolds operator, ringname: a <string> giving the
4418         name of a ring of characteristic 0 where the Molien series is stored,
4419         v: an optional <int>
4420ASSUME:  n is the number of variables of the basering, g the size of the group,
4421         REY is the 1st return value of group_reynolds(), reynolds_molien() or
4422         the second one of primary_invariants(), `ringname` is ring of
4423         characteristic 0 that has been created by molien() or reynolds_molien()
4424         or primary_invariants()
4425RETURN:  secondary invariants of the invariant ring (type <matrix>) and
4426         irreducible secondary invariants (type <matrix>)
4427DISPLAY: information if v does not equal 0
4428EXAMPLE: example secondary_charp; shows an example
4429THEORY:  The secondary invariants are calculated by finding a basis (in terms of
4430         monomials) of the basering modulo the primary invariants, mapping those
4431         to invariants with the Reynolds operator and using these images or
4432         their power products such that they are linearly independent modulo the
4433         primary invariants (see paper \"Some Algorithms in Invariant Theory of
4434         Finite Groups\" by Kemper and Steel (1997)).
4435"
4436{ def br=basering;
4437  degBound=0;
4438 //---------------- checking input and setting verbose mode -------------------
4439  if (char(br)==0)
4440  { "ERROR:   secondary_charp should only be used with rings of characteristic p>0.";
4441    return();
4442  }
4443  int i;
4444  if (size(#)>0)
4445  { if (typeof(#[size(#)])=="int")
4446    { int v=#[size(#)];
4447    }
4448    else
4449    { int v=0;
4450    }
4451  }
4452  else
4453  { int v=0;
4454  }
4455  int n=nvars(br);                     // n is the number of variables, as well
4456                                       // as the size of the matrices, as well
4457                                       // as the number of primary invariants,
4458                                       // we should get
4459  if (ncols(P)<>n)
4460  { "ERROR:   The first parameter ought to be the matrix of the primary";
4461    "         invariants."
4462    return();
4463  }
4464  if (ncols(REY)<>n)
4465  { "ERROR:   The second parameter ought to be the Reynolds operator."
4466    return();
4467  }
4468  if (typeof(`ring_name`)<>"ring")
4469  { "ERROR:   The <string> should give the name of the ring where the Molien."
4470    "         series is stored.";
4471    return();
4472  }
4473  if (v && voice==2)
4474  { "";
4475  }
4476  int j, m, counter, d;
4477  intvec deg_dim_vec;
4478 //- finding the polynomial giving number and degrees of secondary invariants -
4479  for (j=1;j<=n;j=j+1)
4480  { deg_dim_vec[j]=deg(P[j]);
4481  }
4482  setring `ring_name`;
4483  poly p=1;
4484  for (j=1;j<=n;j=j+1)                 // calculating the denominator of the
4485  { p=p*(1-var(1)^deg_dim_vec[j]);     // Hilbert series of the ring generated
4486  }                                    // by the primary invariants -
4487  matrix s[1][2]=M[1,1]*p,M[1,2];      // s is used for canceling
4488  s=matrix(syz(ideal(s)));
4489  p=s[2,1];                            // the polynomial telling us where to
4490                                       // search for secondary invariants
4491  map slead=basering,ideal(0);
4492  p=1/slead(p)*p;                      // smallest term of p needs to be 1
4493  if (v)
4494  { "  Polynomial telling us where to look for secondary invariants:";
4495    "   "+string(p);
4496    "";
4497  }
4498  matrix dimmat=coeffs(p,var(1));      // dimmat will contain the number of
4499                                       // secondary invariants, we need to find
4500                                       // of a certain degree -
4501  m=nrows(dimmat);                     // m-1 is the highest degree
4502  deg_dim_vec=1;
4503  for (j=2;j<=m;j=j+1)
4504  { deg_dim_vec=deg_dim_vec,int(dimmat[j,1]);
4505  }
4506  if (v)
4507  { "  In degree 0 we have: 1";
4508    "";
4509  }
4510 //------------------------ initializing variables ----------------------------
4511  setring br;
4512  intmat PP;
4513  poly pp;
4514  int k;
4515  intvec deg_vec;
4516  ideal sP=std(ideal(P));
4517  ideal TEST,B,IS;
4518  ideal S=1;                           // 1 is the first secondary invariant
4519 //------------------- generating secondary invariants ------------------------
4520  for (i=2;i<=m;i=i+1)                 // going through deg_dim_vec -
4521  { if (deg_dim_vec[i]<>0)             // when it is == 0 we need to find 0
4522    {                                  // elements in the current degree (i-1)
4523      if (v)
4524      { "  Searching in degree "+string(i-1)+", we need to find "+string(deg_dim_vec[i])+" invariant(s)...";
4525      }
4526      TEST=sP;
4527      counter=0;                       // we'll count up to degvec[i]
4528      if (IS[1]<>0)
4529      { PP=power_products(deg_vec,i-1); // generating power products of
4530      }                                // irreducible secondary invariants
4531      if (size(ideal(PP))<>0)
4532      { for (j=1;j<=ncols(PP);j=j+1)   // going through all of those
4533        { pp=1;
4534          for (k=1;k<=nrows(PP);k=k+1)
4535          { pp=pp*IS[1,k]^PP[k,j];
4536          }
4537          if (reduce(pp,TEST)<>0)
4538          { S=S,pp;
4539            counter=counter+1;
4540            if (v)
4541            { "  We find: "+string(pp);
4542            }
4543            if (deg_dim_vec[i]<>counter)
4544            { TEST=std(TEST+ideal(NF(pp,TEST))); // should be soon replaced by
4545                                                 // next line
4546            // TEST=std(TEST,NF(pp,TEST));
4547            }
4548            else
4549            { break;
4550            }
4551          }
4552        }
4553      }
4554      if (deg_dim_vec[i]<>counter)
4555      { B=sort_of_invariant_basis(sP,REY,i-1,deg_dim_vec[i]*6); // B contains
4556                                       // images of kbase(sP,i-1) under the
4557                                       // Reynolds operator that are linearly
4558                                       // independent and that don't reduce to
4559                                       // 0 modulo sP -
4560        if (counter==0 && ncols(B)==deg_dim_vec[i]) // then we can add all of B
4561        { S=S,B;
4562          IS=IS+B;
4563          if (deg_vec[1]==0)
4564          { deg_vec=i-1;
4565            if (v)
4566            { "  We find: "+string(B[1]);
4567            }
4568            for (j=2;j<=deg_dim_vec[i];j=j+1)
4569            { deg_vec=deg_vec,i-1;
4570              if (v)
4571              { "  We find: "+string(B[j]);
4572              }
4573            }
4574          }
4575          else
4576          { for (j=1;j<=deg_dim_vec[i];j=j+1)
4577            { deg_vec=deg_vec,i-1;
4578              if (v)
4579              { "  We find: "+string(B[j]);
4580              }
4581            }
4582          }
4583        }
4584        else
4585        { j=0;                         // j goes through all of B -
4586          while (deg_dim_vec[i]<>counter) // need to find deg_dim_vec[i]
4587          {                            // invariants that are linearly
4588                                       // independent modulo TEST
4589            j=j+1;
4590            if (reduce(B[j],TEST)<>0)   // B[j] should be added
4591            { S=S,B[j];
4592              IS=IS+ideal(B[j]);
4593              if (deg_vec[1]==0)
4594              { deg_vec[1]=i-1;
4595              }
4596              else
4597              { deg_vec=deg_vec,i-1;
4598              }
4599              counter=counter+1;
4600              if (v)
4601              { "  We find: "+string(B[j]);
4602              }
4603              if (deg_dim_vec[i]<>counter)
4604              { TEST=std(TEST+ideal(NF(B[j],TEST))); // should be soon replaced
4605                                                     // by next line
4606              // TEST=std(TEST,NF(B[j],TEST));
4607              }
4608            }
4609          }
4610        }
4611      }
4612      if (v)
4613      { "";
4614      }
4615    }
4616  }
4617  if (v)
4618  { "  We're done!";
4619    "";
4620  }
4621  if (ring_name=="aksldfalkdsflkj")
4622  { kill `ring_name`;
4623  }
4624  return(matrix(S),matrix(IS));
4625}
4626example
4627{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed into";
4628  "         characteristic 3)";
4629  echo=2;
4630         ring R=3,(x,y,z),dp;
4631         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
4632         list L=primary_invariants(A);
4633         matrix S,IS=secondary_charp(L[1..size(L)]);
4634         print(S);
4635         print(IS);
4636}
4637
4638proc secondary_no_molien (matrix P, matrix REY, list #)
4639"USAGE:   secondary_no_molien(P,REY[,deg_vec,v]);
4640         P: a 1xn <matrix> with primary invariants, REY: a gxn <matrix>
4641         representing the Reynolds operator, deg_vec: an optional <intvec>
4642         listing some degrees where no non-trivial homogeneous invariants can be
4643         found, v: an optional <int>
4644ASSUME:  n is the number of variables of the basering, g the size of the group,
4645         REY is the 1st return value of group_reynolds(), reynolds_molien() or
4646         the second one of primary_invariants(), deg_vec is the second return
4647         value of primary_char0_no_molien(), primary_charp_no_molien(),
4648         primary_char0_no_molien_random() or primary_charp_no_molien_random()
4649RETURN:  secondary invariants of the invariant ring (type <matrix>)
4650DISPLAY: information if v does not equal 0
4651EXAMPLE: example secondary_no_molien; shows an example
4652THEORY:  The secondary invariants are calculated by finding a basis (in terms of
4653         monomials) of the basering modulo the primary invariants, mapping those
4654         to invariants with the Reynolds operator and using these images as
4655         candidates for secondary invariants.
4656"
4657{ int i;
4658  degBound=0;
4659 //------------------ checking input and setting verbose ----------------------
4660  if (size(#)==1 or size(#)==2)
4661  { if (typeof(#[size(#)])=="int")
4662    { if (size(#)==2)
4663      { if (typeof(#[size(#)-1])=="intvec")
4664        { intvec deg_vec=#[size(#)-1];
4665        }
4666        else
4667        { "ERROR:   the third parameter should be an <intvec>";
4668          return();
4669        }
4670      }
4671      int v=#[size(#)];
4672    }
4673    else
4674    { if (size(#)==1)
4675      { if (typeof(#[size(#)])=="intvec")
4676        { intvec deg_vec=#[size(#)];
4677          int v=0;
4678        }
4679        else
4680        { "ERROR:   the third parameter should be an <intvec>";
4681          return();
4682        }
4683      }
4684      else
4685      { "ERROR:   wrong list of parameters";
4686        return();
4687      }
4688    }
4689  }
4690  else
4691  { if (size(#)>2)
4692    { "ERROR:   there are too many parameters";
4693      return();
4694    }
4695    int v=0;
4696  }
4697  int n=nvars(basering);               // n is the number of variables, as well
4698                                       // as the size of the matrices, as well
4699                                       // as the number of primary invariants,
4700                                       // we should get
4701  if (ncols(P)<>n)
4702  { "ERROR:   The first parameter ought to be the matrix of the primary";
4703    "         invariants."
4704    return();
4705  }
4706  if (ncols(REY)<>n)
4707  { "ERROR:   The second parameter ought to be the Reynolds operator."
4708    return();
4709  }
4710  if (v && voice==2)
4711  { "";
4712  }
4713  int j, m, d;
4714  int max=1;
4715  for (j=1;j<=n;j=j+1)
4716  { max=max*deg(P[j]);
4717  }
4718  max=max/nrows(REY);
4719  if (v)
4720  { "  We need to find "+string(max)+" secondary invariants.";
4721    "";
4722    "  In degree 0 we have: 1";
4723    "";
4724  }
4725 //------------------------- initializing variables ---------------------------
4726  ideal sP=std(ideal(P));
4727  ideal B, TEST;
4728  ideal S=1;                           // 1 is the first secondary invariant
4729  int counter=1;
4730  i=0;
4731  if (defined(deg_vec)<>voice)
4732  { intvec deg_vec;
4733  }
4734  int k=1;
4735 //--------------------- generating secondary invariants ----------------------
4736  while (counter<>max)
4737  { i=i+1;
4738    if (deg_vec[k]<>i)
4739    { if (v)
4740      { "  Searching in degree "+string(i)+"...";
4741      }
4742      B=sort_of_invariant_basis(sP,REY,i,max); // B contains images of
4743                                       // kbase(sP,i) under the Reynolds
4744                                       // operator that are linearly independent
4745                                       // and that don't reduce to 0 modulo sP
4746      TEST=sP;
4747      for (j=1;j<=ncols(B);j=j+1)
4748      {                                // that are linearly independent modulo
4749                                       // TEST
4750        if (reduce(B[j],TEST)<>0)      // B[j] should be added
4751        { S=S,B[j];
4752          counter=counter+1;
4753          if (v)
4754          { "  We find: "+string(B[j]);
4755          }
4756          if (counter==max)
4757          { break;
4758          }
4759          else
4760          { if (j<>ncols(B))
4761            { TEST=std(TEST+ideal(NF(B[j],TEST))); // should soon be replaced by
4762                                                   // next line
4763            // TEST=std(TEST,NF(B[j],TEST));
4764            }
4765          }
4766        }
4767      }
4768    }
4769    else
4770    { if (size(deg_vec)==k)
4771      { k=1;
4772      }
4773      else
4774      { k=k+1;
4775      }
4776    }
4777  }
4778  if (v)
4779  { "";
4780  }
4781  if (v)
4782  { "  We're done!";
4783    "";
4784  }
4785  return(matrix(S));
4786}
4787example
4788{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
4789  echo=2;
4790         ring R=0,(x,y,z),dp;
4791         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
4792         list L=primary_invariants(A,intvec(1,1,0));
4793         matrix S=secondary_no_molien(L[1..3]);
4794         print(S);
4795}
4796
4797proc secondary_and_irreducibles_no_molien (matrix P, matrix REY, list #)
4798"USAGE:   secondary_and_irreducibles_no_molien(P,REY[,v]);
4799         P: a 1xn <matrix> with primary invariants, REY: a gxn <matrix>
4800         representing the Reynolds operator, v: an optional <int>
4801ASSUME:  n is the number of variables of the basering, g the size of the group,
4802         REY is the 1st return value of group_reynolds(), reynolds_molien() or
4803         the second one of primary_invariants()
4804RETURN:  secondary invariants of the invariant ring (type <matrix>) and
4805         irreducible secondary invariants (type <matrix>)
4806DISPLAY: information if v does not equal 0
4807EXAMPLE: example secondary_and_irreducibles_no_molien; shows an example
4808THEORY:  The secondary invariants are calculated by finding a basis (in terms of
4809         monomials) of the basering modulo the primary invariants, mapping those
4810         to invariants with the Reynolds operator and using these images or
4811         their power products such that they are linearly independent modulo the
4812         primary invariants (see paper \"Some Algorithms in Invariant Theory of
4813         Finite Groups\" by Kemper and Steel (1997)).
4814"
4815{ int i;
4816  degBound=0;
4817 //--------------------- checking input and setting verbose mode --------------
4818  if (size(#)==1 or size(#)==2)
4819  { if (typeof(#[size(#)])=="int")
4820    { if (size(#)==2)
4821      { if (typeof(#[size(#)-1])=="intvec")
4822        { intvec deg_vec=#[size(#)-1];
4823        }
4824        else
4825        { "ERROR:   the third parameter should be an <intvec>";
4826          return();
4827        }
4828      }
4829      int v=#[size(#)];
4830    }
4831    else
4832    { if (size(#)==1)
4833      { if (typeof(#[size(#)])=="intvec")
4834        { intvec deg_vec=#[size(#)];
4835          int v=0;
4836        }
4837        else
4838        { "ERROR:   the third parameter should be an <intvec>";
4839          return();
4840        }
4841      }
4842      else
4843      { "ERROR:   wrong list of parameters";
4844        return();
4845      }
4846    }
4847  }
4848  else
4849  { if (size(#)>2)
4850    { "ERROR:   there are too many parameters";
4851      return();
4852    }
4853    int v=0;
4854  }
4855  int n=nvars(basering);               // n is the number of variables, as well
4856                                       // as the size of the matrices, as well
4857                                       // as the number of primary invariants,
4858                                       // we should get
4859  if (ncols(P)<>n)
4860  { "ERROR:   The first parameter ought to be the matrix of the primary";
4861    "         invariants."
4862    return();
4863  }
4864  if (ncols(REY)<>n)
4865  { "ERROR:   The second parameter ought to be the Reynolds operator."
4866    return();
4867  }
4868  if (v && voice==2)
4869  { "";
4870  }
4871  int j, m, d;
4872  int max=1;
4873  for (j=1;j<=n;j=j+1)
4874  { max=max*deg(P[j]);
4875  }
4876  max=max/nrows(REY);
4877  if (v)
4878  { "  We need to find "+string(max)+" secondary invariants.";
4879    "";
4880    "  In degree 0 we have: 1";
4881    "";
4882  }
4883 //------------------------ initializing variables ----------------------------
4884  intmat PP;
4885  poly pp;
4886  int k;
4887  intvec irreducible_deg_vec;
4888  ideal sP=std(ideal(P));
4889  ideal B,TEST,IS;
4890  ideal S=1;                           // 1 is the first secondary invariant
4891  int counter=1;
4892  i=0;
4893  if (defined(deg_vec)<>voice)
4894  { intvec deg_vec;
4895  }
4896  int l=1;
4897 //------------------- generating secondary invariants ------------------------
4898  while (counter<>max)
4899  { i=i+1;
4900    if (deg_vec[l]<>i)
4901    { if (v)
4902      { "  Searching in degree "+string(i)+"...";
4903      }
4904      TEST=sP;
4905      if (IS[1]<>0)
4906      { PP=power_products(irreducible_deg_vec,i);  // generating all power
4907      }                                // products of irreducible secondary
4908                                       // invariants
4909      if (size(ideal(PP))<>0)
4910      { for (j=1;j<=ncols(PP);j=j+1)   // going through all those power products
4911        { pp=1;
4912          for (k=1;k<=nrows(PP);k=k+1)
4913          { pp=pp*IS[1,k]^PP[k,j];
4914          }
4915          if (reduce(pp,TEST)<>0)
4916          { S=S,pp;
4917            counter=counter+1;
4918            if (v)
4919            { "  We find: "+string(pp);
4920            }
4921            if (counter<>max)
4922            { TEST=std(TEST+ideal(NF(pp,TEST))); // should soon be replaced by
4923                                                 // next line
4924            // TEST=std(TEST,NF(pp,TEST));
4925            }
4926            else
4927            { break;
4928            }
4929          }
4930        }
4931      }
4932      if (max<>counter)
4933      { B=sort_of_invariant_basis(sP,REY,i,max); // B contains images of
4934                                       // kbase(sP,i) under the Reynolds
4935                                       // operator that are linearly independent
4936                                       // and that don't reduce to 0 modulo sP
4937        for (j=1;j<=ncols(B);j=j+1)
4938        { if (reduce(B[j],TEST)<>0)    // B[j] should be added
4939          { S=S,B[j];
4940            IS=IS+ideal(B[j]);
4941            if (irreducible_deg_vec[1]==0)
4942            { irreducible_deg_vec[1]=i;
4943            }
4944            else
4945            { irreducible_deg_vec=irreducible_deg_vec,i;
4946            }
4947            counter=counter+1;
4948            if (v)
4949            { "  We find: "+string(B[j]);
4950            }
4951            if (counter==max)
4952            { break;
4953            }
4954            else
4955            { if (j<>ncols(B))
4956              { TEST=std(TEST+ideal(NF(B[j],TEST))); // should soon be replaced
4957                                                     // by next line
4958              // TEST=std(TEST,NF(B[j],TEST));
4959              }
4960            }
4961          }
4962        }
4963      }
4964    }
4965    else
4966    { if (size(deg_vec)==l)
4967      { l=1;
4968      }
4969      else
4970      { l=l+1;
4971      }
4972    }
4973  }
4974  if (v)
4975  { "";
4976  }
4977  if (v)
4978  { "  We're done!";
4979    "";
4980  }
4981  return(matrix(S),matrix(IS));
4982}
4983example
4984{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
4985  echo=2;
4986         ring R=0,(x,y,z),dp;
4987         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
4988         list L=primary_invariants(A,intvec(1,1,0));
4989         matrix S,IS=secondary_and_irreducibles_no_molien(L[1..2]);
4990         print(S);
4991         print(IS);
4992}
4993
4994proc secondary_not_cohen_macaulay (matrix P, list #)
4995"USAGE:   secondary_not_cohen_macaulay(P,G1,G2,...[,v]);
4996         P: a 1xn <matrix> with primary invariants, G1,G2,...: nxn <matrices>
4997         generating a finite matrix group, v: an optional <int>
4998ASSUME:  n is the number of variables of the basering
4999RETURN:  secondary invariants of the invariant ring (type <matrix>)
5000DISPLAY: information if v does not equal 0
5001EXAMPLE: example secondary_not_cohen_macaulay; shows an example
5002THEORY:  The secondary invariants are generated following \"Generating Invariant
5003         Rings of Finite Groups over Arbitrary Fields\" by Kemper (1996, to
5004         appear in JSC).
5005"
5006{ int i, j;
5007  degBound=0;
5008  def br=basering;
5009  int n=nvars(br);                     // n is the number of variables, as well
5010                                       // as the size of the matrices, as well
5011                                       // as the number of primary invariants,
5012                                       // we should get -
5013  if (size(#)>0)                       // checking input and setting verbose
5014  { if (typeof(#[size(#)])=="int")
5015    { int gen_num=size(#)-1;
5016      if (gen_num==0)
5017      { "ERROR:   There are no generators of the finite matrix group given.";
5018        return();
5019      }
5020      int v=#[size(#)];
5021      for (i=1;i<=gen_num;i=i+1)
5022      { if (typeof(#[i])<>"matrix")
5023        { "ERROR:   These parameters should be generators of the finite matrix group.";
5024          return();
5025        }
5026        if ((n<>nrows(#[i])) or (n<>ncols(#[i])))
5027        { "ERROR:   matrices need to be square and of the same dimensions";
5028          return();
5029        }
5030      }
5031    }
5032    else
5033    { int v=0;
5034      int gen_num=size(#);
5035      for (i=1;i<=gen_num;i=i+1)
5036      { if (typeof(#[i])<>"matrix")
5037        { "ERROR:   These parameters should be generators of the finite matrix group.";
5038          return();
5039        }
5040        if ((n<>nrows(#[i])) or (n<>ncols(#[i])))
5041        { "ERROR:   matrices need to be square and of the same dimensions";
5042          return();
5043        }
5044      }
5045    }
5046  }
5047  else
5048  { "ERROR:   There are no generators of the finite matrix group given.";
5049    return();
5050  }
5051  if (ncols(P)<>n)
5052  { "ERROR:   The first parameter ought to be the matrix of the primary";
5053    "         invariants."
5054    return();
5055  }
5056  if (v && voice==2)
5057  { "";
5058  }
5059  ring alskdfalkdsj=0,x,dp;
5060  matrix M[1][2]=1,(1-x)^n;            // we look at our primary invariants as
5061  export alskdfalkdsj;
5062  export M;
5063  setring br;                          // such of the subgroup that only
5064  matrix REY=matrix(maxideal(1));      // contains the identity, this means that
5065                                       // ch does not divide the order anymore,
5066                                       // this means that we can make use of the
5067                                       // Molien series again - M[1,1]/M[1,2] is
5068                                       // the Molien series of that group, we
5069                                       // now calculate the secondary invariants
5070                                       // of this subgroup in the usual fashion
5071                                       // where the primary invariants are the
5072                                       // ones from the bigger group
5073  if (v)
5074  { "  The procedure secondary_charp() is called to calculate secondary invariants";
5075    "  of the invariant ring of the trivial group with respect to the primary";
5076    "  invariants found previously.";
5077    "";
5078  }
5079  matrix trivialS=secondary_charp(P,REY,"alskdfalkdsj",v);
5080  kill alskdfalkdsj;
5081  // now we have those secondary invariants
5082  int k=ncols(trivialS);               // k is the number of the secondary
5083                                       // invariants, we just calculated
5084  if (v)
5085  { "  We calculate secondary invariants from the ones found for the trivial";
5086    "  subgroup.";
5087    "";
5088  }
5089  map f;                               // used to let generators act on
5090                                       // secondary invariants with respect to
5091                                       // the trivial group -
5092  matrix M(1)[gen_num][k];             // M(1) will contain a module
5093  ideal B;
5094  for (i=1;i<=gen_num;i=i+1)
5095  { B=ideal(matrix(maxideal(1))*transpose(#[i]));   // image of the various
5096                                       // variables under the i-th generator -
5097    f=br,B;                            // the corresponding mapping -
5098    B=f(trivialS)-trivialS;            // these relations should be 0 -
5099    M(1)[i,1..k]=B[1..k];              // we will look for the syzygies of M(1)
5100  }
5101  module M(2)=nres(M(1),2)[2];
5102  int m=ncols(M(2));                   // number of generators of the module
5103                                       // M(2) -
5104  // the following steps calculates the intersection of the module M(2) with
5105  // the algebra A^k where A denote the subalgebra of the usual polynomial
5106  // ring, generated by the primary invariants
5107  string mp=string(minpoly);           // generating a ring where we can do
5108                                       // elimination
5109  execute "ring R=("+charstr(br)+"),(x(1..n),y(1..n),h),dp;";
5110  execute "minpoly=number("+mp+");";
5111  map f=br,maxideal(1);                // canonical mapping
5112  matrix M[k][m+k*n];
5113  M[1..k,1..m]=matrix(f(M(2)));        // will contain a module -
5114  matrix P=f(P);                       // primary invariants in the new ring
5115  for (i=1;i<=n;i=i+1)
5116  { for (j=1;j<=k;j=j+1)
5117    { M[j,m+(i-1)*k+j]=y(i)-P[1,i];
5118    }
5119  }
5120  M=elim(module(M),1,n);               // eliminating x(1..n), std-calculation
5121                                       // is done internally -
5122  M=homog(module(M),h);                // homogenize for 'minbase'
5123  M=minbase(module(M));
5124  setring br;
5125  ideal substitute=maxideal(1),ideal(P),1;
5126  f=R,substitute;                      // replacing y(1..n) by primary
5127                                       // invariants -
5128  M(2)=f(M);                           // M(2) is the new module
5129  m=ncols(M(2));
5130  matrix S[1][m];
5131  S=matrix(trivialS)*matrix(M(2));     // S now contains the secondary
5132                                       // invariants
5133  for (i=1; i<=m;i=i+1)
5134  { S[1,i]=S[1,i]/leadcoef(S[1,i]); // making elements nice
5135  }
5136  S=sort(ideal(S))[1];
5137  if (v)
5138  { "  These are the secondary invariants: ";
5139    for (i=1;i<=m;i=i+1)
5140    { "   "+string(S[1,i]);
5141    }
5142    "";
5143    "  We're done!";
5144    "";
5145  }
5146  if ((v or (voice==2)) && (m>1))
5147  { "  WARNING: The invariant ring might not have a Hironaka decomposition";
5148    "           if the characteristic of the coefficient field divides the";
5149    "           group order.";
5150  }
5151  return(S);
5152}
5153example
5154{ echo=2;
5155           ring R=2,(x,y,z),dp;
5156           matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
5157           list L=primary_invariants(A);
5158           matrix S=secondary_not_cohen_macaulay(L[1],A);
5159           print(S);
5160}
5161
5162proc invariant_ring (list #)
5163"USAGE:   invariant_ring(G1,G2,...[,flags]);
5164         G1,G2,...: <matrices> generating a finite matrix group, flags: an
5165         optional <intvec> with three entries: if the first one equals 0, the
5166         program attempts to compute the Molien series and Reynolds operator,
5167         if it equals 1, the program is told that the Molien series should not
5168         be computed, if it equals -1 characteristic 0 is simulated, i.e. the
5169         Molien series is computed as if the base field were characteristic 0
5170         (the user must choose a field of large prime characteristic, e.g.
5171         32003) and if the first one is anything else, it means that the
5172         characteristic of the base field divides the group order (i.e. it will
5173         not even be attempted to compute the Reynolds operator or Molien
5174         series), the second component should give the size of intervals
5175         between canceling common factors in the expansion of the Molien series,
5176         0 (the default) means only once after generating all terms, in prime
5177         characteristic also a negative number can be given to indicate that
5178         common factors should always be canceled when the expansion is simple
5179         (the root of the extension field does not occur among the coefficients)
5180RETURN:  primary and secondary invariants (both of type <matrix>) generating the
5181         invariant ring with respect to the matrix group generated by the
5182         matrices in the input and irreducible secondary invariants (type
5183         <matrix>) if the Molien series was available
5184DISPLAY: information about the various stages of the program if the third flag
5185         does not equal 0
5186EXAMPLE: example invariant_ring; shows an example
5187THEORY:  Bases of homogeneous invariants are generated successively and those
5188         are chosen as primary invariants that lower the dimension of the ideal
5189         generated by the previously found invariants (see paper \"Generating a
5190         Noetherian Normalization of the Invariant Ring of a Finite Group\" by
5191         Decker, Heydtmann, Schreyer (1997) to appear in JSC). In the
5192         non-modular case secondary invariants are calculated by finding a
5193         basis (in terms of monomials) of the basering modulo the primary
5194         invariants, mapping to invariants with the Reynolds operator and using
5195         those or their power products such that they are linearly independent
5196         modulo the primary invariants (see paper \"Some Algorithms in Invariant
5197         Theory of Finite Groups\" by Kemper and Steel (1997)). In the modular
5198         case they are generated according to \"Generating Invariant Rings of
5199         Finite Groups over Arbitrary Fields\" by Kemper (1996, to appear in
5200         JSC).
5201"
5202{ if (size(#)==0)
5203  { "ERROR:   There are no generators given.";
5204    return();
5205  }
5206  int ch=char(basering);               // the algorithms depend very much on the
5207                                       // characteristic of the ground field -
5208  int n=nvars(basering);               // n is the number of variables, as well
5209                                       // as the size of the matrices, as well
5210                                       // as the number of primary invariants,
5211                                       // we should get
5212  int gen_num;
5213  int mol_flag, v;
5214 //------------------- checking input and setting flags -----------------------
5215  if (typeof(#[size(#)])=="intvec")
5216  { if (size(#[size(#)])<>3)
5217    { "ERROR:   The <intvec> should have three entries.";
5218      return();
5219    }
5220    gen_num=size(#)-1;
5221    mol_flag=#[size(#)][1];
5222    if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0)))
5223    { "ERROR:   the second component of <intvec> should be >=0";
5224      return();
5225    }
5226    int interval=#[size(#)][2];
5227    v=#[size(#)][3];
5228  }
5229  else
5230  { gen_num=size(#);
5231    mol_flag=0;
5232    int interval=0;
5233    v=0;
5234  }
5235 //----------------------------------------------------------------------------
5236  if (mol_flag==0)                     // calculation Molien series will be
5237  { if (ch==0)                         // attempted -
5238    { matrix REY,M=reynolds_molien(#[1..gen_num],intvec(0,interval,v)); // one
5239                                       // will contain Reynolds operator and the
5240                                       // other enumerator and denominator of
5241                                       // Molien series
5242      matrix P=primary_char0(REY,M,v);
5243      matrix S,IS=secondary_char0(P,REY,M,v);
5244      return(P,S,IS);
5245    }
5246    else
5247    { list L=group_reynolds(#[1..gen_num],v);
5248      if (L[1]<>0)                     // testing whether we are in the modular
5249      { string newring="aksldfalkdsflkj"; // case
5250        if (minpoly==0)
5251        { if (v)
5252          { "  We are dealing with the non-modular case.";
5253          }
5254          molien(L[3..size(L)],newring,L[2],intvec(0,interval,v));
5255          matrix P=primary_charp(L[1],newring,v);
5256          matrix S,IS=secondary_charp(P,L[1],newring,v);
5257          if (defined(aksldfalkdsflkj)==2)
5258          { kill aksldfalkdsflkj;
5259          }
5260          return(P,S,IS);
5261        }
5262        else
5263        { if (v)
5264          { "  Since it is impossible for this programme to calculate the Molien
5265 series for";
5266            "  invariant rings over extension fields of prime characteristic, we
5267 have to";
5268            "  continue without it.";
5269            "";
5270
5271          }
5272          list l=primary_charp_no_molien(L[1],v);
5273          if (size(l)==2)
5274          { matrix S=secondary_no_molien(l[1],L[1],l[2],v);
5275          }
5276          else
5277          { matix S=secondary_no_molien(l[1],L[1],v);
5278          }
5279          return(l[1],S);
5280        }
5281      }
5282      else                             // the modular case
5283      { if (v)
5284        { "  There is also no Molien series or Reynolds operator, we can make use of...";
5285          "";
5286          "  We can start looking for primary invariants...";
5287          "";
5288        }
5289        matrix P=primary_charp_without(#[1..gen_num],v);
5290        matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
5291        return(P,S);
5292      }
5293    }
5294  }
5295  if (mol_flag==1)                     // the user wants no calculation of the
5296  { list L=group_reynolds(#[1..gen_num],v); // Molien series
5297    if (ch==0)
5298    { list l=primary_char0_no_molien(L[1],v);
5299      if (size(l)==2)
5300      { matix S=secondary_no_molien(l[1],L[1],l[2],v);
5301      }
5302      else
5303      { matix S=secondary_no_molien(l[1],L[1],v);
5304      }
5305      return(l[1],S);
5306    }
5307    else
5308    { if (L[1]<>0)                     // testing whether we are in the modular
5309      { list l=primary_charp_no_molien(L[1],v); // case
5310        if (size(l)==2)
5311        { matix S=secondary_no_molien(l[1],L[1],l[2],v);
5312        }
5313        else
5314        { matix S=secondary_no_molien(l[1],L[1],v);
5315        }
5316        return(l[1],S);
5317      }
5318      else                             // the modular case
5319      { if (v)
5320        { "  We can start looking for primary invariants...";
5321          "";
5322        }
5323        matrix P=primary_charp_without(#[1..gen_num],v);
5324        matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
5325        return(L[1],S);
5326      }
5327    }
5328  }
5329  if (mol_flag==-1)
5330  { if (ch==0)
5331    { "ERROR:   Characteristic 0 can only be simulated in characteristic p>>0.
5332";
5333      return();
5334    }
5335    list L=group_reynolds(#[1..gen_num],v);
5336    string newring="aksldfalkdsflkj";
5337    molien(L[2..size(L)],newring,intvec(1,interval,v));
5338    matrix P=primary_charp(L[1],newring,v);
5339    matrix S,IS=secondary_charp(P,L[1],newring,v);
5340    kill aksldfalkdsflkj;
5341    return(P,S,IS);
5342  }
5343  else                                 // the user specified that the
5344  { if (ch==0)                         // characteristic divides the group order
5345    { "ERROR:   The characteristic cannot divide the group order when it is 0.
5346";
5347      return();
5348    }
5349    if (v)
5350    { "";
5351    }
5352    matrix P=primary_charp_without(#[1..gen_num],v);
5353    matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
5354    return(L[1],S);
5355  }
5356}
5357example
5358{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
5359  echo=2;
5360         ring R=0,(x,y,z),dp;
5361         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
5362         matrix P,S,IS=invariant_ring(A);
5363         print(P);
5364         print(S);
5365         print(IS);
5366}
5367
5368proc invariant_ring_random (list #)
5369"USAGE:   invariant_ring_random(G1,G2,...,r[,flags]);
5370         G1,G2,...: <matrices> generating a finite matrix group, r: an <int>
5371         where -|r| to |r| is the range of coefficients of random
5372         combinations of bases elements that serve as primary invariants,
5373         flags: an optional <intvec> with three entries: if the first one equals
5374         0, the program attempts to compute the Molien series and Reynolds
5375         operator, if it equals 1, the program is told that the Molien series
5376         should not be computed, if it equals -1 characteristic 0 is simulated,
5377         i.e. the Molien series is computed as if the base field were
5378         characteristic 0 (the user must choose a field of large prime
5379         characteristic, e.g.  32003) and if the first one is anything else, it
5380         means that the characteristic of the base field divides the group order
5381         (i.e. it will not even be attempted to compute the Reynolds operator or
5382         Molien series), the second component should give the size of intervals
5383         between canceling common factors in the expansion of the Molien series,
5384         0 (the default) means only once after generating all terms, in prime
5385         characteristic also a negative number can be given to indicate that
5386         common factors should always be canceled when the expansion is simple
5387         (the root of the extension field does not occur among the coefficients)
5388RETURN:  primary and secondary invariants (both of type <matrix>) generating the
5389         invariant ring with respect to the matrix group generated by the
5390         matrices in the input and irreducible secondary invariants (type
5391         <matrix>) if the Molien series was available
5392DISPLAY: information about the various stages of the program if the third flag
5393         does not equal 0
5394EXAMPLE: example invariant_ring_random; shows an example
5395THEORY:  is the same as for invariant_ring except that random combinations of
5396         basis elements are chosen as candidates for primary invariants and
5397         hopefully they lower the dimension of the previously found primary
5398         invariants by the right amount.
5399"
5400{ if (size(#)<2)
5401  { "ERROR:   There are too few parameters.";
5402    return();
5403  }
5404  int ch=char(basering);               // the algorithms depend very much on the
5405                                       // characteristic of the ground field
5406  int n=nvars(basering);               // n is the number of variables, as well
5407                                       // as the size of the matrices, as well
5408                                       // as the number of primary invariants,
5409                                       // we should get
5410  int gen_num;
5411  int mol_flag, v;
5412 //------------------- checking input and setting flags -----------------------
5413  if (typeof(#[size(#)])=="intvec" && typeof(#[size(#)-1])=="int")
5414  { if (size(#[size(#)])<>3)
5415    { "ERROR:   <intvec> should have three entries.";
5416      return();
5417    }
5418    gen_num=size(#)-2;
5419    mol_flag=#[size(#)][1];
5420    if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0)))
5421    { "ERROR:   the second component of <intvec> should be >=0";
5422      return();
5423    }
5424    int interval=#[size(#)][2];
5425    v=#[size(#)][3];
5426    int max=#[size(#)-1];
5427    if (gen_num==0)
5428    { "ERROR:   There are no generators of a finite matrix group given.";
5429      return();
5430    }
5431  }
5432  else
5433  { if (typeof(#[size(#)])=="int")
5434    { gen_num=size(#)-1;
5435      mol_flag=0;
5436      int interval=0;
5437      v=0;
5438      int max=#[size(#)];
5439    }
5440   else
5441    { "ERROR:   If the two last parameters are not <int> and <intvec>, the last";
5442      "         parameter should be an <int>.";
5443      return();
5444    }
5445  }
5446  for (int i=1;i<=gen_num;i=i+1)
5447  { if (typeof(#[i])=="matrix")
5448    { if (nrows(#[i])<>n or ncols(#[i])<>n)
5449      { "ERROR:   The number of variables of the base ring needs to be the same";
5450        "         as the dimension of the square matrices";
5451        return();
5452      }
5453    }
5454    else
5455    { "ERROR:   The first parameters should be a list of matrices";
5456      return();
5457    }
5458  }
5459 //----------------------------------------------------------------------------
5460  if (mol_flag==0)
5461  { if (ch==0)
5462    { matrix REY,M=reynolds_molien(#[1..gen_num],intvec(0,interval,v)); // one
5463                                       // will contain Reynolds operator and the
5464                                       // other enumerator and denominator of
5465                                       // Molien series
5466      matrix P=primary_char0_random(REY,M,max,v);
5467      matrix S,IS=secondary_char0(P,REY,M,v);
5468      return(P,S,IS);
5469    }
5470    else
5471    { list L=group_reynolds(#[1..gen_num],v);
5472      if (L[1]<>0)                     // testing whether we are in the modular
5473      { string newring="aksldfalkdsflkj"; // case
5474        if (minpoly==0)
5475        { if (v)
5476          { "  We are dealing with the non-modular case.";
5477          }
5478          molien(L[3..size(L)],newring,L[2],intvec(0,interval,v));
5479          matrix P=primary_charp_random(L[1],newring,max,v);
5480          matrix S,IS=secondary_charp(P,L[1],newring,v);
5481          if (voice==2)
5482          { kill aksldfalkdsflkj;
5483          }
5484          return(P,S,IS);
5485        }
5486        else
5487        { if (v)
5488          { "  Since it is impossible for this programme to calculate the Molien
5489 series for";
5490            "  invariant rings over extension fields of prime characteristic, we
5491 have to";
5492            "  continue without it.";
5493            "";
5494
5495          }
5496          list l=primary_charp_no_molien_random(L[1],max,v);
5497          if (size(l)==2)
5498          { matix S=secondary_no_molien(l[1],L[1],l[2],v);
5499          }
5500          else
5501          { matix S=secondary_no_molien(l[1],L[1],v);
5502          }
5503          return(l[1],S);
5504        }
5505      }
5506      else                             // the modular case
5507      { if (v)
5508        { "  There is also no Molien series, we can make use of...";
5509          "";
5510          "  We can start looking for primary invariants...";
5511          "";
5512        }
5513        matrix P=primary_charp_without_random(#[1..gen_num],max,v);
5514        matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
5515        return(L[1],S);
5516      }
5517    }
5518  }
5519  if (mol_flag==1)                     // the user wants no calculation of the
5520  { list L=group_reynolds(#[1..gen_num],v); // Molien series
5521    if (ch==0)
5522    { list l=primary_char0_no_molien_random(L[1],max,v);
5523      if (size(l)==2)
5524      { matix S=secondary_no_molien(l[1],L[1],l[2],v);
5525      }
5526      else
5527      { matix S=secondary_no_molien(l[1],L[1],v);
5528      }
5529      return(l[1],S);
5530    }
5531    else
5532    { if (L[1]<>0)                     // testing whether we are in the modular
5533      { list l=primary_charp_no_molien_random(L[1],max,v); // case
5534        if (size(l)==2)
5535        { matix S=secondary_no_molien(l[1],L[1],l[2],v);
5536        }
5537        else
5538        { matix S=secondary_no_molien(l[1],L[1],v);
5539        }
5540        return(l[1],S);
5541      }
5542      else                             // the modular case
5543      { if (v)
5544        { "  We can start looking for primary invariants...";
5545          "";
5546        }
5547        matrix P=primary_charp_without_random(#[1..gen_num],max,v);
5548        matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
5549        return(L[1],S);
5550      }
5551    }
5552  }
5553  if (mol_flag==-1)
5554  { if (ch==0)
5555    { "ERROR:   Characteristic 0 can only be simulated in characteristic p>>0.
5556";
5557      return();
5558    }
5559    list L=group_reynolds(#[1..gen_num],v);
5560    string newring="aksldfalkdsflkj";
5561    molien(L[2..size(L)],newring,intvec(1,v));
5562    matrix P=primary_charp_random(L[1],newring,max,v);
5563    matrix S,IS=secondary_charp(P,L[1],newring,v);
5564    kill aksldfalkdsflkj;
5565    return(P,S,IS);
5566  }
5567  else                                 // the user specified that the
5568  { if (ch==0)                         // characteristic divides the group order
5569    { "ERROR:   The characteristic cannot divide the group order when it is 0.
5570";
5571      return();
5572    }
5573    if (v)
5574    { "";
5575    }
5576    matrix P=primary_charp_without_random(#[1..gen_num],max,v);
5577    matrix S=secondary_not_cohen_macaulay(P,#[1..gen_num],v);
5578    return(L[1],S);
5579  }
5580}
5581example
5582{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
5583  echo=2;
5584         ring R=0,(x,y,z),dp;
5585         matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
5586         matrix P,S,IS=invariant_ring_random(A,1);
5587         print(P);
5588         print(S);
5589         print(IS);
5590}
5591
5592proc algebra_containment (poly p, matrix A)
5593"USAGE:   algebra_containment(p,A);
5594         p: arbitrary <poly>, A: a 1xm <matrix> giving generators of a
5595         subalgebra of the basering
5596RETURN:  1 (TRUE) (type <int>) if p is contained in the subalgebra
5597         0 (FALSE) (type <int>) if <poly> is not contained
5598DISPLAY: a representation of p in terms of algebra generators A[1,i]=y(i) if p
5599         is contained in the subalgebra
5600EXAMPLE: example algebra_containment; shows an example
5601THEORY:  The ideal of algebraic relations of the algebra generators f1,...,fm
5602         given by A is computed introducing new variables y(i) and the product
5603         order: x^a*y^b > y^d*y^e if x^a > x^d or else if y^b > y^e. p reduces
5604         to a polynomial only in the y(i) <=> p is contained in the subring
5605         generated by the polynomials in A.
5606"
5607{ degBound=0;
5608  if (nrows(A)==1)
5609  { def br=basering;
5610    int n=nvars(br);
5611    int m=ncols(A);
5612    string mp=string(minpoly);
5613    execute "ring R=("+charstr(br)+"),(x(1..n),y(1..m)),(dp(n),dp(m));";
5614    execute "minpoly=number("+mp+");";
5615    ideal vars=x(1..n);
5616    map emb=br,vars;
5617    ideal A=ideal(emb(A));
5618    ideal check=emb(p);
5619    for (int i=1;i<=m;i=i+1)
5620    { A[i]=A[i]-y(i);
5621    }
5622    A=std(A);
5623    check[1]=reduce(check[1],A);
5624    A=elim(check,1,n);
5625    if (A[1]<>0)
5626    { "\/\/ "+string(check);
5627      return(1);
5628    }
5629    else
5630    { return(0);
5631    }
5632  }
5633  else
5634  { "ERROR:   <matrix> may only have one row";
5635    return();
5636  }
5637}
5638example
5639{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
5640  echo=2;
5641         ring R=0,(x,y,z),dp;
5642         matrix A[1][7]=x2+y2,z2,x4+y4,1,x2z-1y2z,xyz,x3y-1xy3;
5643         poly p1=x10z3-x8y2z3+2x6y4z3-2x4y6z3+x2y8z3-y10z3+x6z4+3x4y2z4+3x2y4z4+y6z4;
5644         algebra_containment(p1,A);
5645         poly p2=z;
5646         algebra_containment(p2,A);
5647}
5648
5649proc module_containment(poly p, matrix P, matrix S)
5650"USAGE:   module_containment(p,P,S);
5651         p: arbitrary <poly>, P: a 1xn <matrix> giving generators of an algebra,
5652         S: a 1xt <matrix> giving generators of a module over the algebra
5653         generated by P
5654ASSUME:  n is the number of variables in the basering and the generators in P
5655         are algebraically independent
5656RETURNS: 1 (TRUE) (type <int>) if p is contained in the ring
5657         0 (FALSE) type <int>) if p is not contained
5658DISPLAY: the representation of p in terms of algebra generators P[1,i]=z(i) and
5659         module generators S[1,j]=y(j) if p is contained in the module
5660EXAMPLE: example module_containment; shows an example
5661THEORY:  The ideal of algebraic relations of all the generators p1,...,pn,
5662         s1,...,st given by P and S is computed introducing new variables y(j),
5663         z(i) and the product order: x^a*y^b*z^c > x^d*y^e*z^f if x^a > x^d
5664         with respect to the lp ordering or else if z^c > z^f with respect to
5665         the dp ordering or else if y^b > y^e with respect to the lp ordering
5666         again. p reduces to a polynomial only in the y(j) and z(i) linear in
5667         the z(i)) <=> p is contained in the module.
5668"
5669{ def br=basering;
5670  degBound=0;
5671  int n=nvars(br);
5672  if (ncols(P)==n and nrows(P)==1 and nrows(S)==1)
5673  { int m=ncols(S);
5674    string mp=string(minpoly);
5675    execute "ring R=("+charstr(br)+"),(x(1..n),y(1..m),z(1..n)),(lp(n),dp(m),lp(n));";
5676    execute "minpoly=number("+mp+");";
5677    ideal vars=x(1..n);
5678    map emb=br,vars;
5679    matrix P=emb(P);
5680    matrix S=emb(S);
5681    ideal check=emb(p);
5682    ideal I;
5683    for (int i=1;i<=m;i=i+1)
5684    { I[i]=S[1,i]-y(i);
5685    }
5686    for (i=1;i<=n;i=i+1)
5687    { I[m+i]=P[1,i]-z(i);
5688    }
5689    I=std(I);
5690    check[1]=reduce(check[1],I);
5691    I=elim(check,1,n);                 // checking whether all variables from
5692    if (I[1]<>0)                       // the former ring have disappeared
5693    { "\/\/ "+string(check);
5694      return(1);
5695    }
5696    else
5697    { return(0);
5698    }
5699  }
5700  else
5701  { "ERROR:   the first <matrix> must have the same number of columns as the";
5702    "         basering and both <matrices> may only have one row";
5703    return();
5704  }
5705}
5706example
5707{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
5708  echo=2;
5709         ring R=0,(x,y,z),dp;
5710         matrix P[1][3]=x2+y2,z2,x4+y4;
5711         matrix S[1][4]=1,x2z-1y2z,xyz,x3y-1xy3;
5712         poly p1=x10z3-x8y2z3+2x6y4z3-2x4y6z3+x2y8z3-y10z3+x6z4+3x4y2z4+3x2y4z4+y6z4;
5713         module_containment(p1,P,S);
5714         poly p2=z;
5715         module_containment(p2,P,S);
5716}
5717
5718proc orbit_variety (matrix F,string newring)
5719"USAGE:   orbit_variety(F,s);
5720         F: a 1xm <matrix> defing an invariant ring, s: a <string> giving the
5721         name for a new ring
5722RETURN:  a Groebner basis (type <ideal>, named G) for the ideal defining the
5723         orbit variety (i.e. the syzygy ideal) in the new ring (named `s`)
5724EXAMPLE: example orbit_variety; shows an example
5725THEORY:  The ideal of algebraic relations of the invariant ring generators is
5726         calculated, then the variables of the original ring are eliminated and
5727         the polynomials that are left over define the orbit variety
5728"
5729{ if (newring=="")
5730  { "ERROR:   the second parameter may not be an empty <string>";
5731    return();
5732  }
5733  if (nrows(F)==1)
5734  { def br=basering;
5735    int n=nvars(br);
5736    int m=ncols(F);
5737    string mp=string(minpoly);
5738    execute "ring R=("+charstr(br)+"),("+varstr(br)+",y(1..m)),dp;";
5739    execute "minpoly=number("+mp+");";
5740    ideal I=ideal(imap(br,F));
5741    for (int i=1;i<=m;i=i+1)
5742    { I[i]=I[i]-y(i);
5743    }
5744    I=elim(I,1,n);
5745    execute "ring "+newring+"=("+charstr(br)+"),(y(1..m)),dp(m);";
5746    execute "minpoly=number("+mp+");";
5747    ideal vars;
5748    for (i=2;i<=n;i=i+1)
5749    { vars[i]=0;
5750    }
5751    vars=vars,y(1..m);
5752    map emb=R,vars;
5753    ideal G=emb(I);
5754    kill emb, vars, R;
5755    keepring `newring`;
5756    // execute "keepring "+newring+";";
5757    return();
5758  }
5759  else
5760  { "ERROR:   the <matrix> may only have one row";
5761    return();
5762  }
5763}
5764example
5765{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:";
5766  echo=2;
5767         ring R=0,(x,y,z),dp;
5768         matrix F[1][7]=x2+y2,z2,x4+y4,1,x2z-1y2z,xyz,x3y-1xy3;
5769         string newring="E";
5770         orbit_variety(F,newring);
5771         print(G);
5772         basering;
5773}
5774
5775proc relative_orbit_variety(ideal I,matrix F,string newring)
5776"USAGE:   relative_orbit_variety(I,F,s);
5777         I: an <ideal> invariant under the action of a group, F: a 1xm
5778         <matrix> defining the invariant ring of this group, s: a <string>
5779         giving a name for a new ring
5780RETURN:  a Groebner basis (type <ideal>, named G) for the ideal defining the
5781         relative orbit variety with respect to I in the new ring (named s)
5782EXAMPLE: example relative_orbit_variety; shows an example
5783THEORY:  A Groebner basis of the ideal of algebraic relations of the invariant
5784         ring generators is calculated, then one of the basis elements plus the
5785         ideal generators. The variables of the original ring are eliminated and
5786         the polynomials that are left over define thecrelative orbit variety
5787         with respect to I.
5788"
5789{ if (newring=="")
5790  { "ERROR:   the third parameter may not be empty a <string>";
5791    return();
5792  }
5793  degBound=0;
5794  if (nrows(F)==1)
5795  { def br=basering;
5796    int n=nvars(br);
5797    int m=ncols(F);
5798    string mp=string(minpoly);
5799    execute "ring R=("+charstr(br)+"),("+varstr(br)+",y(1..m)),lp;";
5800    execute "minpoly=number("+mp+");";
5801    ideal J=ideal(imap(br,F));
5802    ideal I=imap(br,I);
5803    for (int i=1;i<=m;i=i+1)
5804    { J[i]=J[i]-y(i);
5805    }
5806    J=std(J);
5807    J=J,I;
5808    option(redSB);
5809    J=std(J);
5810    ideal vars;
5811    //for (i=1;i<=n;i=i+1)
5812    //{ vars[i]=0;
5813    //}
5814    vars[n]=0;
5815    vars=vars,y(1..m);
5816    map emb=R,vars;
5817    ideal G=emb(J);
5818    J=J-G;
5819    for (i=1;i<=ncols(G);i=i+1)
5820    { if (J[i]<>0)
5821      { G[i]=0;
5822      }
5823    }
5824    G=compress(G);
5825    execute "ring "+newring+"=("+charstr(br)+"),(y(1..m)),lp;";
5826    execute "minpoly=number("+mp+");";
5827    ideal vars;
5828    for (i=2;i<=n;i=i+1)
5829    { vars[i]=0;
5830    }
5831    vars=vars,y(1..m);
5832    map emb=R,vars;
5833    ideal G=emb(G);
5834    kill vars, emb;
5835    keepring `newring`;
5836    // execute "keepring "+newring+";";
5837    return();
5838  }
5839  else
5840  { "ERROR:   the <matrix> may only have one row";
5841    return();
5842  }
5843}
5844example
5845{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.6.3:";
5846  echo=2;
5847         ring R=0,(x,y,z),dp;
5848         matrix F[1][3]=x+y+z,xy+xz+yz,xyz;
5849         ideal I=x2+y2+z2-1,x2y+y2z+z2x-2x-2y-2z,xy2+yz2+zx2-2x-2y-2z;
5850         string newring="E";
5851         relative_orbit_variety(I,F,newring);
5852         print(G);
5853         basering;
5854}
5855
5856proc image_of_variety(ideal I,matrix F)
5857"USAGE:   image_of_variety(I,F);
5858         I: an arbitray <ideal>, F: a 1xm <matrix> defining an invariant ring
5859         of a some matrix group
5860RETURN:  the <ideal> defining the image under that group of the variety defined
5861         by I
5862EXAMPLE: example image_of_variety; shows an example
5863THEORY:  relative_orbit_variety(I,F,s) is called and the newly introduced
5864         variables in the output are replaced by the generators of the
5865         invariant ring. This ideal in the original variables defines the image
5866         of the variety defined by I
5867"
5868{ if (nrows(F)==1)
5869  { def br=basering;
5870    int n=nvars(br);
5871    string newring="E";
5872    relative_orbit_variety(I,F,newring);
5873    execute "ring R=("+charstr(br)+"),("+varstr(br)+","+varstr(E)+"),lp;";
5874    ideal F=imap(br,F);
5875    for (int i=1;i<=n;i=i+1)
5876    { F=0,F;
5877    }
5878    setring br;
5879    map emb2=E,F;
5880    return(compress(emb2(G)));
5881  }
5882  else
5883  { "ERROR:   the <matrix> may only have one row";
5884    return();
5885  }
5886}
5887example
5888{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.6.8:";
5889  echo=2;
5890         ring R=0,(x,y,z),dp;
5891         matrix F[1][3]=x+y+z,xy+xz+yz,xyz;
5892         ideal I=xy;
5893         print(image_of_variety(I,F));
5894}
Note: See TracBrowser for help on using the repository browser.