source: git/Singular/LIB/rinvar.lib @ 535a5a7

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