source: git/Singular/LIB/rinvar.lib @ 900cea

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