source: git/Singular/LIB/rinvar.lib @ 6fe3a0

fieker-DuValspielwiese
Last change on this file since 6fe3a0 was 6fe3a0, checked in by Hans Schönemann <hannes@…>, 19 years ago
*hannes/lossen: new libs, new tests git-svn-id: file:///usr/local/Singular/svn/trunk@7891 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 36.1 KB
Line 
1// Last change 10.12.2000 (TB)
2///////////////////////////////////////////////////////////////////////////////
3version="$Id: rinvar.lib,v 1.9 2005-04-25 10:13:07 Singular Exp $";
4category="Invariant theory";
5info="
6LIBRARY:  rinvar.lib      Invariant Rings of Reductive Groups
7AUTHOR:   Thomas Bayer,   tbayer@in.tum.de
8          http://wwwmayr.informatik.tu-muenchen.de/personen/bayert/
9          Current Adress: Institut fuer Informatik, TU Muenchen
10OVERVIEW:
11 Implementation based on Derksen's algorithm. Written in the frame of the
12 diploma thesis (advisor: Prof. Gert-Martin Greuel) 'Computations of moduli
13 spaces of semiquasihomogenous singularities and an implementation in Singular'
14
15PROCEDURES:
16 HilbertSeries(I, w);           Hilbert series of the ideal I w.r.t. weight w
17 HilbertWeights(I, w);          weighted degrees of the generators of I
18 ImageVariety(I, F);            ideal of the image variety F(variety(I))
19 ImageGroup(G, F);              ideal of G w.r.t. the induced representation
20 InvariantRing(G, Gaction);     generators of the invariant ring of G
21 InvariantQ(f, G, Gaction);     decide if f is invariant w.r.t. G
22 LinearizeAction(G, Gaction);   linearization of the action 'Gaction' of G
23 LinearActionQ(action,s,t);     decide if action is linear in var(s..nvars)
24 LinearCombinationQ(base, f);   decide if f is in the linear hull of 'base'
25 MinimalDecomposition(f,s,t);   minimal decomposition of f (like coef)
26 NullCone(G, act);              ideal of the nullcone of the action 'act' of G
27 ReynoldsImage(RO,f);           image of f under the Reynolds operator 'RO'
28 ReynoldsOperator(G, Gaction);  Reynolds operator of the group G
29 SimplifyIdeal(I[,m,s]);        simplify the ideal I (try to reduce variables)
30 TransferIdeal(R,name,nA);      transfer the ideal 'name' from R to basering
31
32SEE ALSO: qhmoduli_lib, zeroset_lib
33";
34
35LIB "presolve.lib";
36LIB "elim.lib";
37LIB "zeroset.lib";
38
39///////////////////////////////////////////////////////////////////////////////
40
41proc EquationsOfEmbedding(ideal embedding, int nrs)
42"USAGE:   EquationsOfEmbedding(embedding, s); ideal embedding; int s;
43PURPOSE: compute the ideal of the variety parameterized by 'embedding' by
44         implicitation and change the variables to the old ones.
45RETURN:  ideal
46ASSUME:  nvars(basering) = n, size(embedding) = r and s = n - r.
47         The polynomials of embedding contain only var(s + 1 .. n).
48NOTE:    the result is the Zariski closure of the parameterized variety
49EXAMPLE: example  EquationsOfEmbedding; shows an example
50"
51{
52  ideal tvars;
53
54  for(int i = nrs + 1; i <= nvars(basering); i++) { tvars[i - nrs] = var(i); }
55
56  def RE1 = ImageVariety(ideal(0), embedding);  // implicitation of the parameterization
57  // map F = RE1, tvars;
58        map F = RE1, maxideal(1);
59  return(F(imageid));
60}
61example
62{"EXAMPLE:";  echo = 2;
63  ring R   = 0,(s(1..5), t(1..4)),dp;
64  ideal emb = t(1), t(2), t(3), t(3)^2;
65  ideal I = EquationsOfEmbedding(emb, 5);
66  I;
67}
68
69///////////////////////////////////////////////////////////////////////////////
70
71proc ImageGroup(ideal Grp, ideal Gaction)
72"USAGE:   ImageGroup(G, action); ideal G, action;
73PURPOSE: compute the ideal of the image of G in GL(m,K) induced by the linear
74         action 'action', where G is an algebraic group and 'action' defines
75         an action of G on K^m (size(action) = m).
76RETURN:  ring, a polynomial ring over the same ground field as the basering,
77         containing the ideals 'groupid' and 'actionid'.
78         - 'groupid' is the ideal of the image of G (order <= order of G)
79         - 'actionid' defines the linear action of 'groupid' on K^m.
80NOTE:    'action' and 'actionid' have the same orbits
81         all variables which give only rise to 0's in the m x m matrices of G
82         have been omitted.
83ASSUME:  basering K[s(1..r),t(1..m)] has r + m variables, G is the ideal of an
84         algebraic group and F is an action of G on K^m. G contains only the
85         variables s(1)...s(r). The action 'action' is given by polynomials
86         f_1,...,f_m in basering, s.t. on the ring level we have
87         K[t_1,...,t_m] --> K[s_1,...,s_r,t_1,...,t_m]/G
88         t_i   -->  f_i(s_1,...,s_r,t_1,...,t_m)
89EXAMPLE: example ImageGroup; shows an example
90"
91{
92  int i, j, k, newVars, nrt, imageSize, dbPrt;
93  ideal matrixEntries;
94  matrix coMx;
95  poly tVars, mPoly;
96  string ringSTR1, ringSTR2, order;
97
98  dbPrt = printlevel-voice+2;
99  dbprint(dbPrt, "Image Group of " + string(Grp) + ", action =  "
100                                   + string(Gaction));
101  def RIGB = basering;
102  mPoly = minpoly;
103  tVars = 1;
104  k = 0;
105
106  // compute the representation of G induced by Gaction, i.e., a matrix
107  // of size(Gaction) x size(Gaction) and polynomials in s(1),...,s(r) as
108  // entries
109  // the matrix is represented as the list 'matrixEntries' where
110  // the entries which are always 0 are omittet.
111
112  for(i = 1; i <= ncols(Gaction); i++) {
113    tVars = tVars * var(i + nvars(basering) - ncols(Gaction));
114  }
115  for(i = 1; i <= ncols(Gaction); i++){
116    coMx = coef(Gaction[i], tVars);
117    for(j = 1; j <= ncols(coMx); j++){
118      k++;
119      matrixEntries[k] = coMx[2, j];
120    }
121  }
122  newVars = size(matrixEntries);
123  nrt = ncols(Gaction);
124
125  // this matrix defines an embedding of G into GL(m, K).
126  // in the next step the ideal of this image is computed
127  // note that we have omitted all variables which give give rise
128  // only to 0's. Note that z(1..newVars) are slack variables
129
130  order = "(dp(" + string(nvars(basering)) + "), dp);";
131  ringSTR1 = "ring RIGR = (" + charstr(basering) + "), (" + varstr(basering)
132                             + ", z(1.." + string(newVars) + "))," + order;
133  execute(ringSTR1);
134  minpoly = number(imap(RIGB, mPoly));
135  ideal I1, I2, Gn, G, F, mEntries, newGaction;
136  G = imap(RIGB, Grp);
137  F = imap(RIGB, Gaction);
138  mEntries = imap(RIGB, matrixEntries);
139
140  // prepare the ideals needed to compute the image
141  // and compute the new action of the image on K^m
142
143  for(i=1;i<=size(mEntries);i++){ I1[i] = var(i + nvars(RIGB))-mEntries[i]; }
144  I1 = std(I1);
145
146  for(i = 1; i <= ncols(F); i++) { newGaction[i] = reduce(F[i], I1); }
147  I2 = G, I1;
148  I2 = std(I2);
149  Gn = nselect(I2, 1, nvars(RIGB));
150  imageSize = ncols(Gn);
151
152  // create a new basering which might contain more variables
153  // s(1..newVars) as the original basering and map the ideal
154  // Gn (contians only z(1..newVars)) to this ring
155
156  ringSTR2 = "ring RIGS = (" + charstr(basering) + "), (s(1.."
157                     + string(newVars) + "), t(1.." + string(nrt) + ")), lp;";
158  execute(ringSTR2);
159  minpoly = number(imap(RIGB, mPoly));
160  ideal mapIdeal, groupid, actionid;
161  int offset;
162
163  // construct the map F : RIGB -> RIGS
164
165  for(i=1;i<=nvars(RIGB)-nrt;i++) { mapIdeal[i] = 0;}            // s(i)-> 0
166  offset = nvars(RIGB) - nrt;
167  for(i=1;i<=nrt;i++) { mapIdeal[i+offset] = var(newVars + i);}  // t(i)->t(i)
168  offset = offset + nrt;
169  for(i=1;i<=newVars;i++) { mapIdeal[i + offset] = var(i);}      // z(i)->s(i)
170
171  // map Gn and newGaction to RIGS
172
173  map F = RIGR, mapIdeal;
174  groupid = F(Gn);
175  actionid = F(newGaction);
176  export groupid, actionid;
177  dbprint(dbPrt+1, "
178// 'ImageGroup' created a new ring.
179// To see the ring, type (if the name 'R' was assigned to the return value):
180     show(R);
181// To access the ideal of the image of the input group and to access the new
182// action of the group, type
183     setring R;  groupid; actionid;
184");
185  setring RIGB;
186  return(RIGS);
187}
188example
189{"EXAMPLE:";  echo = 2;
190  ring B   = 0,(s(1..2), t(1..2)),dp;
191  ideal G = s(1)^3-1, s(2)^10-1;
192  ideal action = s(1)*s(2)^8*t(1), s(1)*s(2)^7*t(2);
193  def R = ImageGroup(G, action);
194  setring R;
195  groupid;
196  actionid;
197}
198
199///////////////////////////////////////////////////////////////////////////////
200
201proc HilbertWeights(ideal I, wt)
202"USAGE:   HilbertWeights(I, w); ideal I, intvec wt
203PURPOSE: compute the weights of the "slack" variables needed for the
204         computation of the algebraic relations of the generators of 'I' s.t.
205         the Hilbert driven 'std' can be used.
206RETURN: intvec
207ASSUME: basering = K[t_1,...,t_m,...], 'I' is quasihomogenous w.r.t. 'w' and
208        contains only polynomials in t_1,...,t_m
209"
210{
211  int offset = size(wt);
212  intvec wtn = wt;
213
214  for(int i = 1; i <= size(I); i++) { wtn[offset + i] = deg(I[i], wt); }
215  return(wtn);
216}
217
218///////////////////////////////////////////////////////////////////////////////
219
220proc HilbertSeries(ideal I, wt)
221"USAGE:   HilbertSeries(I, w); ideal I, intvec wt
222PURPOSE: compute the polynomial p of the Hilbert Series,represented by p/q, of
223         the ring K[t_1,...,t_m,y_1,...,y_r]/I1 where 'w' are the weights of
224         the variables, computed, e.g., by 'HilbertWeights', 'I1' is of the
225         form I[1] - y_1,...,I[r] - y_r and is quasihomogenous w.r.t. 'w'
226RETURN:  intvec
227NOTE:    the leading 0 of the result does not belong to p, but is needed in
228         the hilbert-driven 'std'.
229"
230{
231  int i;
232  intvec hs1;
233  matrix coMx;
234  poly f = 1;
235
236  for(i = 1; i <= ncols(I); i++) { f = f * (1 - var(1)^deg(I[i], wt));}
237  coMx = coeffs(f, var(1));
238  for(i = 1; i <= deg(f) + 1; i++) {
239    hs1[i] = int(coMx[i, 1]);
240  }
241  hs1[size(hs1) + 1] = 0;
242  return(hs1);
243}
244///////////////////////////////////////////////////////////////////////////////
245
246proc HilbertSeries1(wt)
247"USAGE:   HilbertSeries1(wt); ideal I, intvec wt
248PURPOSE: compute the polynomial p of the Hilbert Series represented by p/q of
249         the ring K[t_1,...,t_m,y_1,...,y_r]/I where I is a complete inter-
250         section and the generator I[i] has degree wt[i]
251RETURN:  poly
252"
253{
254  int i, j;
255  intvec hs1;
256  matrix ma;
257  poly f = 1;
258
259  for(i = 1; i <= size(wt); i++) { f = f * (1 - var(1)^wt[i]);}
260  ma = coef(f, var(1));
261  j = ncols(ma);
262  for(i = 0; i <= deg(f); i++) {
263    if(var(1)^i == ma[1, j]) {
264      hs1[i + 1] = int(ma[2, j]);
265      j--;
266    }
267    else { hs1[i + 1] = 0; }
268  }
269  hs1[size(hs1) + 1] = 0;
270  return(hs1);
271}
272
273///////////////////////////////////////////////////////////////////////////////
274
275proc ImageVariety(ideal I, F, list #)
276"USAGE:   ImageVariety(ideal I, F [, w]);ideal I; F is a list/ideal, intvec w.
277PURPOSE: compute the Zariski closure of the image of the variety of I under
278         the morphism F.
279NOTE:    if 'I' and 'F' are quasihomogenous w.r.t. 'w' then the Hilbert-driven
280         'std' is used.
281RETURN:  polynomial ring over the same groundfield, containing the ideal
282         'imageid'. The variables are Y(1),...,Y(k) where k = size(F)
283         - 'imageid' is the ideal of the Zariski closure of F(X) where
284           X is the variety of I.
285EXAMPLE: example ImageVariety; shows an example
286"
287{
288  int i, dbPrt, nrNewVars;
289  intvec wt, wth, hs1;
290  string ringSTR1, ringSTR2, order;
291
292  def RARB = basering;
293  nrNewVars = size(F);
294
295  dbPrt = printlevel-voice+2;
296  dbprint(dbPrt, "ImageVariety of " + string(I) + " under the map " + string(F));
297
298  if(size(#) > 0) { wt = #[1]; }
299
300  // create new ring for elimination, Y(1),...,Y(m) are slack variables.
301
302  poly mPoly = minpoly;
303  order = "(dp(" + string(nvars(basering)) + "), dp);";
304  ringSTR1 = "ring RAR1 = (" + charstr(basering) + "), (" + varstr(basering) + ", Y(1.." + string(nrNewVars) + ")), " + order;
305  ringSTR2 = "ring RAR2 = (" + charstr(basering) + "), Y(1.." + string(nrNewVars) + "), dp;";
306  execute(ringSTR1);
307  minpoly = number(imap(RARB, mPoly));
308
309  ideal I, J1, J2, Fm;
310
311  I = imap(RARB, I);
312  Fm = imap(RARB, F);
313
314  if(size(wt) > 1) {
315    wth = HilbertWeights(Fm, wt);
316    hs1 = HilbertSeries(Fm, wt);
317  }
318
319  // get the ideal of the graph of F : X -> Y and compute a standard basis
320
321  for(i = 1; i <= nrNewVars; i++) { J1[i] = var(i + nvars(RARB)) - Fm[i];}
322  J1 = J1, I;
323  if(size(wt) > 1) {
324    J1 = std(J1, hs1, wth);  // Hilbert-driven algorithm
325  }
326  else {
327    J1 = std(J1);
328  }
329
330  // forget all elements which contain other than the slack variables
331
332  J2 = nselect(J1, 1, nvars(RARB));
333
334  execute(ringSTR2);
335  minpoly = number(imap(RARB, mPoly));
336  ideal imageid = imap(RAR1, J2);
337  export(imageid);
338     dbprint(dbPrt+1, "
339// 'ImageVariety' created a new ring.
340// To see the ring, type (if the name 'R' was assigned to the return value):
341     show(R);
342// To access the ideal of the image variety, type
343     setring R;  imageid;
344");
345  setring RARB;
346  return(RAR2);
347}
348example
349{"EXAMPLE:";  echo = 2;
350  ring B   = 0,(x,y),dp;
351  ideal I  = x4 - y4;
352  ideal F  = x2, y2, x*y;
353  def R = ImageVariety(I, F);
354  setring R;
355  imageid;
356}
357
358///////////////////////////////////////////////////////////////////////////////
359
360proc LinearizeAction(ideal Grp, Gaction, int nrs)
361"USAGE:   LinearizeAction(G,action,r); ideal G, action; int r
362PURPOSE: linearize the group action 'action' and find an equivariant embedding
363         of K^m where m = size(action).
364ASSUME:  G contains only variables var(1..r) (r = nrs)
365basering = K[s(1..r),t(1..m)], K = Q or K = Q(a) and minpoly != 0.
366RETURN:  polynomial ring containing the ideals 'actionid', 'embedid', 'groupid'
367         - 'actionid' is the ideal defining the linearized action of G
368         - 'embedid' is a parameterization of an equivariant embedding (closed)
369         - 'groupid' is the ideal of G in the new ring
370NOTE:    set printlevel > 0 to see a trace
371EXAMPLE: example LinearizeAction; shows an example
372"
373{
374  def altring = basering;
375  int i, j, k, ok, loop, nrt, sizeOfDecomp, dbPrt;
376  intvec wt;
377  ideal action, basis, G, reduceIdeal;
378  matrix decompMx;
379  poly actCoeff;
380  string str, order, mPoly;
381
382  dbPrt = printlevel-voice+2;
383  dbprint(dbPrt, "LinearizeAction " + string(Gaction));
384  def RLAR = basering;
385  mPoly = "minpoly = " + string(minpoly) + ";";
386  order = ordstr(basering);
387  nrt = ncols(Gaction);
388  for(i = 1; i <= nrs; i++) { wt[i] = 0;}
389  for(i = nrs + 1; i <= nrs + nrt; i++) { basis[i - nrs] = var(i); wt[i] = 1;}
390  dbprint(dbPrt, "  basis = " + string(basis));
391  if(attrib(Grp, "isSB")) { G = Grp; }
392  else { G = std(Grp); }
393  reduceIdeal = G;
394  action = Gaction;
395  loop = 1;
396  i = 1;
397
398  // check if each component of 'action' is linear in t(1),...,t(nrt).
399
400  while(loop){
401    if(deg(action[i], wt) <= 1) {
402      sizeOfDecomp = 0;
403      dbprint(dbPrt, "  " + string(action[i]) + " is linear");
404    }
405    else {   // action[i] is not linear
406
407    // compute the minimal decomposition of action[i]
408    // action[i]=decompMx[1,1]*decompMx[2,1]+ ... +decompMx[1,k]*decompMx[2,k]
409    // decompMx[1,j] contains variables var(1)...var(nrs)
410    // decompMx[2,j] contains variables var(nrs + 1)...var(nvars(basering))
411
412    dbprint(dbPrt, "  " + string(action[i])
413                   + " is not linear, a minimal decomposition is :");
414    decompMx = MinimalDecomposition(action[i], nrs, nrt);
415    sizeOfDecomp = ncols(decompMx);
416    dbprint(dbPrt, decompMx);
417
418    for(j = 1; j <= sizeOfDecomp; j++) {   
419      // check if decompMx[2,j] is a linear combination of basis elements
420      actCoeff = decompMx[2, j];
421      ok = LinearCombinationQ(basis, actCoeff, nrt + nrs);
422      if(ok == 0) {
423
424        // nonlinear element, compute new component of the action
425
426        dbprint(dbPrt, "  the polynomial " + string(actCoeff)
427                   + " is not a linear combination of the elements of basis");
428        nrt++;
429        str = charstr(basering) + ", (" + varstr(basering)
430                   + ",t(" + string(nrt) + ")),";
431        if(defined(RLAB)) { kill(RLAB);}
432        def RLAB = basering;
433        if(defined(RLAR)) { kill(RLAR);}
434        execute("ring RLAR = " + str + "(" + order + ");");
435        execute(mPoly);
436
437        ideal basis, action, G, reduceIdeal;
438        matrix decompMx;
439        map F;
440        poly actCoeff;
441
442        wt[nrs + nrt] = 1;
443        basis = imap(RLAB, basis), imap(RLAB, actCoeff);
444        action = imap(RLAB, action);
445        decompMx = imap(RLAB, decompMx);
446        actCoeff = imap(RLAB, actCoeff);
447        G = imap(RLAB, G);
448        attrib(G, "isSB", 1);
449        reduceIdeal = imap(RLAB, reduceIdeal), actCoeff - var(nrs + nrt);
450
451        // compute action on the new basis element
452
453        for(k = 1; k <= nrs; k++) { F[k] = 0;}
454        for(k = nrs + 1; k < nrs + nrt; k++) { F[k] = action[k - nrs];}
455        actCoeff = reduce(F(actCoeff), G);
456        action[ncols(action) + 1] = actCoeff;
457        dbprint(dbPrt, "  extend basering by " + string(var(nrs + nrt)));
458        dbprint(dbPrt, "  new basis = " + string(basis));
459        dbprint(dbPrt, "  action of G on new basis element = "
460                    + string(actCoeff));
461        dbprint(dbPrt, "  decomp : " + string(decompMx[2, j]) + " -> "
462                    + string(var(nrs + nrt)));
463      } // end if
464      else {
465        dbprint(dbPrt, "  the polynomial " + string(actCoeff)
466                    + " is a linear combination of the elements of basis");
467      }
468    } // end for
469    reduceIdeal = std(reduceIdeal);
470    action[i] = reduce(action[i], reduceIdeal);
471    } // end else
472    if(i < ncols(action)) { i++;}
473    else {loop = 0;}
474  } // end while
475  if(defined(actionid)) { kill actionid; }
476  ideal actionid, embedid, groupid;
477  actionid = action;
478  embedid = basis;
479  groupid = G;
480  export actionid, embedid, groupid;
481dbprint(dbPrt+1, "
482// 'LinearizeAction' created a new ring.
483// To see the ring, type (if the name 'R' was assigned to the return value):
484     show(R);
485// To access the new action and the equivariant embedding, type
486     setring R; actionid; embedid; groupid
487");
488  setring altring;
489  return(RLAR);
490}
491example
492{"EXAMPLE:";  echo = 2;
493  ring B   = 0,(s(1..5), t(1..3)),dp;
494  ideal G =  s(3)-s(4), s(2)-s(5), s(4)*s(5), s(1)^2*s(4)+s(1)^2*s(5)-1, s(1)^2*s(5)^2-s(5), s(4)^4-s(5)^4+s(1)^2, s(1)^4+s(4)^3-s(5)^3, s(5)^5-s(1)^2*s(5);
495  ideal action = -s(4)*t(1)+s(5)*t(1), -s(4)^2*t(2)+2*s(4)^2*t(3)^2+s(5)^2*t(2), s(4)*t(3)+s(5)*t(3);
496  LinearActionQ(action, 5);
497  def R = LinearizeAction(G, action, 5);
498  setring R;
499  R;
500  actionid;
501  embedid;
502  groupid;
503  LinearActionQ(actionid, 5);
504}
505
506///////////////////////////////////////////////////////////////////////////////
507
508proc LinearActionQ(Gaction, int nrs)
509"USAGE:   LinearActionQ(action,nrs); ideal action, int nrs
510PURPOSE: check if the action defined by 'action' is linear w.r.t. the variables
511         var(nrs + 1...nvars(basering)).
512RETURN:  0 action not linear
513         1 action is linear
514EXAMPLE: example LinearActionQ; shows an example
515"
516{
517  int i, nrt, loop;
518  intvec wt;
519
520  nrt = ncols(Gaction);
521  for(i = 1; i <= nrs; i++) { wt[i] = 0;}
522  for(i = nrs + 1; i <= nrs + nrt; i++) { wt[i] = 1;}
523  loop = 1;
524  i = 1;
525  while(loop)
526  {
527    if(deg(Gaction[i], wt) > 1) { loop = 0; }
528    else
529    {
530      i++;
531      if(i > ncols(Gaction)) { loop = 0;}
532    }
533  }
534  return(i > ncols(Gaction));
535}
536example
537{"EXAMPLE:";  echo = 2;
538  ring R   = 0,(s(1..5), t(1..3)),dp;
539  ideal G =  s(3)-s(4), s(2)-s(5), s(4)*s(5), s(1)^2*s(4)+s(1)^2*s(5)-1,
540             s(1)^2*s(5)^2-s(5), s(4)^4-s(5)^4+s(1)^2, s(1)^4+s(4)^3-s(5)^3,
541             s(5)^5-s(1)^2*s(5);
542  ideal Gaction = -s(4)*t(1)+s(5)*t(1),
543                  -s(4)^2*t(2)+2*s(4)^2*t(3)^2+s(5)^2*t(2),
544                   s(4)*t(3)+s(5)*t(3);
545  LinearActionQ(Gaction, 5);
546  LinearActionQ(Gaction, 8);
547}
548
549///////////////////////////////////////////////////////////////////////////////
550
551proc LinearCombinationQ(ideal I, poly f)
552"USAGE:   LinearCombination(I, f); ideal I, poly f
553PURPOSE: test if f can be written as a linear combination of the generators of I.
554RETURN:  0 f is not a linear combination
555         1 f is a linear combination
556"
557{
558  int i, loop, sizeJ;
559  ideal J;
560
561  J = I, f;
562  sizeJ = size(J);
563
564  def RLC = ImageVariety(ideal(0), J);   // compute algebraic relations
565  setring RLC;
566  matrix coMx;
567  poly relation = 0;
568
569  loop = 1;
570  i = 1;
571  while(loop)
572  { // look for a linear relation containing Y(nr)
573    if(deg(imageid[i]) == 1)
574    {
575      coMx = coef(imageid[i], var(sizeJ));
576      if(coMx[1,1] == var(sizeJ))
577      {
578        relation = imageid[i];
579        loop = 0;
580      }
581    }
582    else
583    {
584      i++;
585      if(i > ncols(imageid)) { loop = 0;}
586    }
587  }
588  return(i <= ncols(imageid));
589}
590
591///////////////////////////////////////////////////////////////////////////////
592
593proc InvariantRing(ideal G, ideal action, list #)
594"USAGE:   InvariantRing(G, Gact [, opt]); ideal G, Gact; int opt
595PURPOSE: compute generators of the invariant ring of G w.r.t. the action 'Gact'
596ASSUME:  G is a finite group and 'Gact' is a linear action.
597RETURN:  ring R; this ring comes with the ideals 'invars' and 'groupid' and
598         with the poly 'newA':
599         - 'invars' contains the algebra generators of the invariant ring
600         - 'groupid' is the ideal of G in the new ring
601         - 'newA' is the new representation of the primitive root of the
602           minimal polynomial of the ring which was active when calling the
603           procedure (if the minpoly did not change, 'newA' is set to 'a').
604NOTE:    the minimal polynomial of the output ring depends on some random
605         choices
606EXAMPLE: example InvariantRing; shows an example
607"
608{
609  int i, ok, dbPrt, noReynolds, primaryDec;
610  ideal invarsGens, groupid;
611
612  dbPrt = printlevel-voice+2;
613
614  if(size(#) > 0) { primaryDec = #[1]; }
615  else { primaryDec = 0; }
616
617  dbprint(dbPrt, "InvariantRing of " + string(G));
618  dbprint(dbPrt, " action  = " + string(action));
619
620  if(!attrib(G, "isSB")) { groupid = std(G);}
621  else { groupid = G; }
622
623  // compute the nullcone of G by means of Derksen's algorithm
624
625  invarsGens = NullCone(groupid, action);  // compute nullcone of linear action
626  dbprint(dbPrt, " generators of zero-fibre ideal are " + string(invarsGens));
627
628  // make all generators of the nullcone invariant
629  // if necessary, compute the Reynolds Operator, i.e., find all elements
630  // of the variety defined by G. It might be necessary to extend the
631  // ground field.
632
633  def IRB = basering;
634  if(defined(RIRR)) { kill(RIRR);}
635  def RIRR = basering;
636  setring RIRR;
637//  export(RIRR);
638//  export(invarsGens);
639  noReynolds = 1;
640  dbprint(dbPrt, " nullcone is generated by " + string(size(invarsGens)));
641   dbprint(dbPrt, " degrees = " + string(maxdeg(invarsGens)));
642  for(i = 1; i <= ncols(invarsGens); i++){
643    ok = InvariantQ(invarsGens[i], groupid, action);
644    if(ok) { dbprint(dbPrt, string(i) + ": poly " + string(invarsGens[i])
645                                      + " is invariant");}
646    else {
647      if(noReynolds) {
648
649        // compute the Reynolds operator and change the ring !
650        noReynolds = 0;
651        def RORN = ReynoldsOperator(groupid, action, primaryDec);
652        setring RORN;
653        ideal groupid = std(id);
654        attrib(groupid, "isSB", 1);
655        ideal action = actionid;
656//        ideal invarsGens = TransferIdeal(RIRR, "invarsGens", newA);
657//        export(invarsGens);
658// lossen: ersetzt durch
659        setring RIRR;
660        string parName, minPoly;
661        if(npars(basering) == 0) {
662          parName = "a";
663          minPoly = "0";
664        }
665        else {
666          parName = parstr(basering);
667          minPoly = string(minpoly);
668        }
669        execute("ring RA1=0,(" + varstr(basering) + "," + parName + "), lp;");
670        if (minPoly!="0") { execute("ideal mpoly = std(" + minPoly + ");"); }
671        ideal I = imap(RIRR,invarsGens);
672        setring RORN;
673        map Phi = RA1, maxideal(1);
674        Phi[nvars(RORN) + 1] = newA;
675        ideal invarsGens = Phi(I);
676        kill Phi,RA1,RIRR;
677// end of ersetzt durch     
678
679      }
680      dbprint(dbPrt, string(i) + ": poly " + string(invarsGens[i])
681                               + " is NOT invariant");
682      invarsGens[i] = ReynoldsImage(ROelements, invarsGens[i]);
683      dbprint(dbPrt, " --> " + string(invarsGens[i]));
684    }
685  }
686  for(i = 1; i <= ncols(invarsGens); i++){
687    ok =  InvariantQ(invarsGens[i], groupid, action);
688    if(ok) { dbprint(dbPrt, string(i) + ": poly " + string(invarsGens[i])
689                                      + " is invariant"); }
690    else { print(string(i) + ": Fatal Error with Reynolds ");}
691  }
692  if(noReynolds == 0) {
693    def RIRS = RORN;
694    setring RIRS;
695    kill RORN;
696    export groupid;
697  }
698  else {
699    def RIRS = RIRR;
700    kill RIRR;
701    setring RIRS;
702    export groupid;
703  }
704  ideal invars = invarsGens;
705  kill invarsGens;
706  if (defined(ROelements)) { kill ROelements,actionid,zeroset,id; }
707  export invars;
708  dbprint(dbPrt+1, "
709// 'InvariantRing' created a new ring.
710// To see the ring, type (if the name 'R' was assigned to the return value):
711     show(R);
712// To access the generators of the invariant ring type
713     setring R; invars;
714// Note that the input group G is stored in R as the ideal 'groupid'; to
715// see it, type
716    groupid;
717// Note that 'InvariantRing' might change the minimal polynomial
718// The representation of the algebraic number is given by 'newA'
719");
720  setring IRB;
721  return(RIRS);
722}
723example
724{"EXAMPLE:";  echo = 2;
725  ring B = 0, (s(1..2), t(1..2)), dp;
726  ideal G = -s(1)+s(2)^3, s(1)^4-1;
727  ideal action = s(1)*t(1), s(2)*t(2);
728
729  def R = InvariantRing(std(G), action);
730  setring R;
731  invars;
732}
733
734///////////////////////////////////////////////////////////////////////////////
735
736proc InvariantQ(poly f, ideal G, action)
737"USAGE:   InvariantQ(f, G, action); poly f; ideal G, action
738PURPOSE: check if the polynomial f is invariant w.r.t. G where G acts via
739         'action' on K^m.
740ASSUME:  basering = K[s_1,...,s_m,t_1,...,t_m] where K = Q of K = Q(a) and
741         minpoly != 0, f contains only t_1,...,t_m, G is the ideal of an
742         algebraic group and a standardbasis.
743RETURN:  int;
744         0 if f is not invariant,
745         1 if f is invariant
746NOTE:    G need not be finite
747"
748{
749  def altring=basering;
750  map F;
751  if(deg(f) == 0) { return(1); }
752  for(int i = 1; i <= size(action); i++) {
753    F[nvars(basering) - size(action) + i] = action[i];
754  }
755  return(reduce(f - F(f), G) == 0);
756}
757
758///////////////////////////////////////////////////////////////////////////////
759
760proc MinimalDecomposition(poly f, int nrs, int nrt)
761"USAGE:   MinimalDecomposition(f,a,b); poly f; int a, b.
762PURPOSE: decompose f as a sum M[1,1]*M[2,1] + ... + M[1,r]*M[2,r] where M[1,i]
763         contains only s(1..a), M[2,i] contains only t(1...b) s.t. r is minimal
764ASSUME:  f polynomial in K[s(1..a),t(1..b)], K = Q or K = Q(a) and minpoly != 0
765RETURN:  2 x r matrix M s.t.  f = M[1,1]*M[2,1] + ... + M[1,r]*M[2,r]
766EXAMPLE: example MinimalDecomposition;
767"
768{
769  int i, sizeOfMx, changed, loop;
770  list initialTerms;
771  matrix coM1, coM2, coM, decompMx, auxM;
772  matrix m[2][2] = 0,1,1,0;
773  poly vars1, vars2;
774
775  if(f == 0) { return(decompMx); }
776
777  // first decompose f w.r.t. t(1..nrt)
778  // then  decompose f w.r.t. s(1..nrs)
779
780  vars1 = RingVarProduct(nrs+1..nrt+nrs);
781  vars2 = RingVarProduct(1..nrs);
782  coM1 = SimplifyCoefficientMatrix(m*coef(f, vars1));    // exchange rows of decomposition
783  coM2 = SimplifyCoefficientMatrix(coef(f, vars2));
784  if(ncols(coM2) < ncols(coM1)) {
785    auxM = coM1;
786    coM1 = coM2;
787    coM2 = auxM;
788  }
789  decompMx = coM1;        // decompMx is the smaller decomposition
790  if(ncols(decompMx) == 1) { return(decompMx);}      // n = 1 is minimal
791  changed = 0;
792  loop = 1;
793  i = 1;
794
795  // first loop, try coM1
796
797  while(loop) {
798    coM = MinimalDecomposition(f - coM1[1, i]*coM1[2, i], nrs, nrt);
799    if(size(coM) == 1) { sizeOfMx = 0; }    // coM = 0
800    else {sizeOfMx = ncols(coM); }      // number of columns
801    if(sizeOfMx + 1 < ncols(decompMx)) {    // shorter decomposition
802      changed = 1;
803      decompMx = coM;
804      initialTerms[1] =  coM1[1, i];
805      initialTerms[2] =  coM1[2, i];
806    }
807    if(sizeOfMx == 1) { loop = 0;}      // n = 2 is minimal
808    if(i < ncols(coM1)) {i++;}
809    else {loop = 0;}
810  }
811  if(sizeOfMx > 1) {            // n > 2
812    loop = 1;          // coM2 might yield
813    i = 1;            // a smaller decomposition
814  }
815
816  // first loop, try coM2
817
818  while(loop) {
819    coM = MinimalDecomposition(f - coM2[1, i]*coM2[2, i], nrs, nrt);
820    if(size(coM) == 1) { sizeOfMx = 0; }
821    else {sizeOfMx = ncols(coM); }
822    if(sizeOfMx + 1 < ncols(decompMx)) {
823      changed = 1;
824      decompMx = coM;
825      initialTerms[1] =  coM2[1, i];
826      initialTerms[2] =  coM2[2, i];
827    }
828    if(sizeOfMx == 1) { loop = 0;}
829    if(i < ncols(coM2)) {i++;}
830    else {loop = 0;}
831  }
832  if(!changed) { return(decompMx); }
833  if(size(decompMx) == 1) { matrix decompositionM[2][1];}
834  else  { matrix decompositionM[2][ncols(decompMx) + 1];}
835  decompositionM[1, 1] = initialTerms[1];
836  decompositionM[2, 1] = initialTerms[2];
837  if(size(decompMx) > 1) {
838    for(i = 1; i <= ncols(decompMx); i++) {
839      decompositionM[1, i + 1] = decompMx[1, i];
840      decompositionM[2, i + 1] = decompMx[2, i];
841    }
842  }
843  return(decompositionM);
844}
845example
846{"EXAMPLE:"; echo = 2;
847  ring R = 0, (s(1..2), t(1..2)), dp;
848  poly h = s(1)*(t(1) + t(1)^2) +  (t(2) + t(2)^2)*(s(1)^2 + s(2));
849  matrix M = MinimalDecomposition(h, 2, 2);
850  M;
851  M[1,1]*M[2,1] + M[1,2]*M[2,2] - h;
852}
853
854///////////////////////////////////////////////////////////////////////////////
855
856proc NullCone(ideal G, action)
857"USAGE:   NullCone(G, action); ideal G, action
858PURPOSE: compute the ideal of the nullcone of the linear action of G on K^n,
859         given by 'action', by means of Deksen's algorithm
860ASSUME:  basering = K[s(1..r),t(1..n)], K = Q or K = Q(a) and minpoly != 0,
861         G is an ideal of a reductive algebraic group in K[s(1..r)],
862         'action' is a linear group action of G on K^n (n = ncols(action))
863RETURN:  ideal of the nullcone of G.
864NOTE:    the generators of the nullcone are homogenous, but i.g. not invariant
865EXAMPLE: example NullCone; shows an example
866"
867{
868  int i, nt, dbPrt, offset, groupVars;
869  poly minPoly;
870  string ringSTR, vars, order;
871  def RNCB = basering;
872
873  // prepare the ring needed for the computation
874  // s(1...) variables of the group
875  // t(1...) variables of the affine space
876  // y(1...) additional 'slack' variables
877
878  nt = size(action);
879  order = "(dp(" + string(nvars(basering) - nt) + "), dp);";
880  vars = "(s(1.." + string(nvars(basering) - nt);
881  vars = vars +"),t(1.."+string(nt) + "), Y(1.." + string(nt) + "))," + order;
882  ringSTR = "ring RNCR = (" + charstr(basering) +  "),"  + vars;
883  // ring for the computation
884
885  minPoly = minpoly;
886  offset =  size(G) + nt;
887  execute(ringSTR);
888  def aaa=imap(RNCB, minPoly);
889  if (aaa!=0) { minpoly = number(aaa); }
890  ideal action, G, I, J, N, generators;
891  map F;
892  poly f;
893
894  // built the ideal of the graph of GxV -> V, (s,v) -> s(v), i.e.
895  // of the image of the map GxV -> GxVxV, (s,v) -> (s,v,s(v))
896
897  G = fetch(RNCB, G);
898  action = fetch(RNCB, action);
899  groupVars = nvars(basering) - 2*nt;
900  offset =  groupVars + nt;
901  I = G;
902  for(i = 1; i <= nt; i = i + 1) {
903    I = I, var(offset + i) - action[i];
904  }
905
906  J = std(I);  // takes long, try to improve
907
908  // eliminate
909
910  N = nselect(J, 1, groupVars);
911
912  // substitute
913  for(i = 1; i <= nvars(basering); i = i + 1) { F[i] = 0; }
914  for(i = groupVars + 1; i <= offset; i = i + 1) { F[i] = var(i); }
915
916  generators = mstd(F(N))[2];
917  setring RNCB;
918  return(fetch(RNCR, generators));
919}
920example
921{"EXAMPLE:";  echo = 2;
922  ring R = 0, (s(1..2), x, y), dp;
923  ideal G = -s(1)+s(2)^3, s(1)^4-1;
924  ideal action = s(1)*x, s(2)*y;
925
926  ideal inv = NullCone(G, action);
927  inv;
928}
929
930///////////////////////////////////////////////////////////////////////////////
931
932proc ReynoldsOperator(ideal Grp, ideal Gaction, list #)
933"USAGE:   ReynoldsOperator(G, action [, opt]); ideal G, action; int opt
934PURPOSE: compute the Reynolds operator of the group G which act via 'action'
935RETURN:  polynomial ring R over a simple extension of the groundfield of the
936         basering (the extension might be trivial), containing a list
937         'ROelements', the ideals 'id', 'actionid' and the polynomial 'newA'.
938         R = K(a)[s(1..r),t(1..n)].
939         - 'ROelements'  is a list of ideals, each ideal represents a
940           substitution map F : R -> R according to the zero-set of G
941         - 'id' is the ideal of G in the new ring
942         - 'newA' is the new representation of a' in terms of a. If the
943           basering does not contain a parameter then 'newA' = 'a'.
944ASSUME:  basering = K[s(1..r),t(1..n)], K = Q or K = Q(a') and minpoly != 0,
945         G is the ideal of a finite group in K[s(1..r)], 'action' is a linear
946         group action of G
947"
948{
949  def ROBR = basering;
950  int i, j, n, ns, primaryDec;
951  ideal G1 = Grp;
952  list solution, saction;
953  string str;
954
955  if(size(#) > 0) { primaryDec = #[1]; }
956  else { primaryDec = 0; }
957  kill #;
958
959  n = nvars(basering);
960  ns = n - size(Gaction);
961  for(i = ns + 1; i <= n; i++) { G1 = G1, var(i);}
962
963  def RORR = ZeroSet(G1, primaryDec);
964//  setring RORR;
965//  id = TransferIdeal(ROBR, "Grp", newA);  // defined in ZeroSet ...
966//  ideal actionid = TransferIdeal(ROBR, "Gaction", newA);
967// lossen: ersetzt durch
968  setring ROBR;
969  string parName, minPoly;
970  if(npars(basering) == 0) {
971    parName = "a";
972    minPoly = "0";
973  }
974  else {
975    parName = parstr(basering);
976    minPoly = string(minpoly);
977  }
978  execute("ring RA1=0,(" + varstr(basering) + "," + parName + "), lp;");
979  if (minPoly!="0") { execute("ideal mpoly = std(" + minPoly + ");"); }
980  ideal Grp = imap(ROBR,Grp);
981  ideal Gaction = imap(ROBR,Gaction);
982  setring RORR;
983  map Phi = RA1, maxideal(1);
984  Phi[nvars(RORR) + 1] = newA;
985  ideal id = Phi(Grp);
986  ideal actionid = Phi(Gaction);
987  kill parName,minPoly,Phi,RA1;
988// end of ersetzt durch   
989  list ROelements;
990  ideal Rf;
991  map groupElem;
992  poly h1, h2;
993
994  for(i = 1; i <= size(zeroset); i++) {
995    groupElem = zeroset[i];    // element of G
996    for(j = ns + 1; j<=n; j++) { groupElem[j] = var(j); } //do not change t's
997    for(j = 1; j <= n - ns; j++) {
998      h1 = actionid[j];
999      h2 = groupElem(h1);
1000      Rf[ns + j] = h2;
1001    }
1002    ROelements[i] = Rf;
1003  }
1004  export actionid, ROelements;
1005  setring ROBR;
1006  return(RORR);
1007}
1008
1009///////////////////////////////////////////////////////////////////////////////
1010
1011proc ReynoldsImage(list reynoldsOp, poly f)
1012"USAGE:   ReynoldsImage(RO, f); list RO, poly f
1013PURPOSE: compute the Reynolds image of the polynomial f where RO represents
1014         the Reynolds operator
1015RETURN:  poly
1016"
1017{
1018  def RIBR=basering;
1019  map F;
1020  poly h = 0;
1021
1022  for(int i = 1; i <= size(reynoldsOp); i++) {
1023    F = RIBR, reynoldsOp[i];
1024    h = h + F(f);
1025  }
1026  return(h/size(reynoldsOp));
1027}
1028
1029///////////////////////////////////////////////////////////////////////////////
1030
1031static proc SimplifyCoefficientMatrix(matrix coefMatrix)
1032"USAGE:   SimplifyCoefficientMatrix(M); M matrix coming from coef(...)
1033PURPOSE: simplify the matrix, i.e. find linear dependencies among the columns
1034RETURN:  matrix M, f = M[1,1]*M[2,1] + ... + M[1,n]*M[2,n]
1035"
1036{
1037  int i, j , loop;
1038  intvec columnList;
1039  matrix decompMx = coefMatrix;
1040
1041  loop = 1;
1042  i = 1;
1043  while(loop) {
1044    columnList = 1..i;            // current column
1045    for(j = i + 1; j <= ncols(decompMx); j++) {
1046      // test if decompMx[2, j] equals const * decompMx[2, i]
1047      if(LinearCombinationQ(ideal(decompMx[2, i]), decompMx[2, j])) {    // column not needed
1048        decompMx[1, i] = decompMx[1, i] +  decompMx[2, j] / decompMx[2, i] * decompMx[1, j];
1049      }
1050      else { columnList[size(columnList) + 1] = j; }
1051    }
1052    if(defined(auxM)) { kill(auxM);}
1053    matrix auxM[2][size(columnList)];      // built new matrix and omit
1054    for(j = 1; j <= size(columnList); j++) {    // the linear dependent colums
1055      auxM[1, j] = decompMx[1, columnList[j]];    // found above
1056      auxM[2, j] = decompMx[2, columnList[j]];
1057    }
1058    decompMx = auxM;
1059    if(i < ncols(decompMx) - 1) { i++;}
1060    else { loop = 0;}
1061  }
1062  return(decompMx);
1063}
1064
1065///////////////////////////////////////////////////////////////////////////////
1066
1067proc SimplifyIdeal(ideal I, list #)
1068"USAGE:   SimplifyIdeal(I [,m, name]); ideal I; int m, string name"
1069PURPOSE: simplify ideal I to the ideal I', do not change the names of the
1070         first m variables, new ideal I' might contain less variables.
1071         I' contains variables var(1..m)
1072RETURN: list
1073  _[1] ideal I'
1074  _[2] ideal representing a map phi to a ring with probably less vars. s.t.
1075       phi(I) = I'
1076  _[3] list of variables
1077  _[4] list from 'elimpart'
1078"
1079{
1080  int i, k, m;
1081  string nameCMD;
1082  ideal mId, In, mapId;  // ideal for the map
1083  list sList, result;
1084
1085  sList = elimpart(I);
1086  In = sList[1];
1087  mapId = sList[5];
1088
1089  if(size(#) > 0)
1090  {
1091    m = #[1];
1092    nameCMD = #[2];
1093  }
1094  else { m = 0;} // nvars(basering);
1095  k = 0;
1096  for(i = 1; i <= nvars(basering); i++)
1097  {
1098    if(sList[4][i] != 0)
1099    {
1100      k++;
1101      if(k <= m) { mId[i] = sList[4][i]; }
1102      else { execute("mId["+string(i) +"] = "+nameCMD+"("+string(k-m)+");");}
1103    }
1104    else { mId[i] = 0;}
1105  }
1106  map phi = basering, mId;
1107  result[1] = phi(In);
1108  result[2] = phi(mapId);
1109  result[3] = simplify(sList[4], 2);
1110  result[4] = sList;
1111  return(result);
1112}
1113
1114///////////////////////////////////////////////////////////////////////////////
1115
1116static proc RingVarProduct(index)
1117// list of indices
1118{
1119  poly f = 1;
1120  for(int i = 1; i <= size(index); i++)
1121  {
1122    f = f * var(index[i]);
1123  }
1124  return(f);
1125}
1126///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.