source: git/Singular/LIB/qhmoduli.lib @ e7778a

spielwiese
Last change on this file since e7778a was e7778a, checked in by Hans Schönemann <hannes@…>, 15 years ago
*hannes: fix git-svn-id: file:///usr/local/Singular/svn/trunk@11648 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 44.6 KB
RevLine 
[35f23d]1///////////////////////////////////////////////////////////////////////////////
[e7778a]2version="$Id: qhmoduli.lib,v 1.14 2009-04-08 14:08:25 Singular Exp $";
[fd3fb7]3category="Singularities";
[35f23d]4info="
[8942a5]5LIBRARY:  qhmoduli.lib    Moduli Spaces of Semi-Quasihomogeneous Singularities
[35f23d]6AUTHOR:   Thomas Bayer, email: bayert@in.tum.de
7
8PROCEDURES:
9 ArnoldAction(f, [G, w])  Induced action of G_f on T_.
10 ModEqn(f)                Equations of the moduli space for principal part f
11 QuotientEquations(G,A,I) Equations of Variety(I)/G w.r.t. action 'A'
12 StabEqn(f)               Equations of the stabilizer of f.
13 StabEqnId(I, w)          Equations of the stabilizer of the qhom. ideal I.
14 StabOrder(f)             Order of the stabilizer of f.
15 UpperMonomials(f, [w])   Upper basis of the Milnor algebra of f.
16
17 Max(data)                maximal integer contained in 'data'
18 Min(data)                minimal integer contained in  'data'
19";
20
21// NOTE: This library has been written in the frame of the diploma thesis
22// 'Computing moduli spaces of semiquasihomogeneous singularities and an
23//  implementation in Singular', Arbeitsgruppe Algebraische Geometrie,
24// Fachbereich Mathematik, University Kaiserslautern,
25// Advisor: Prof. Gert-Martin Greuel
26
27LIB "rinvar.lib";
28
29///////////////////////////////////////////////////////////////////////////////
30
31proc ModEqn(poly f, list #)
32"USAGE:   ModEqn(f [, opt]); poly f; int opt;
[1f92589]33PURPOSE: compute equations of the moduli space of semiquasihomogenos hypersurface         singularity with principal part f w.r.t. right equivalence
34ASSUME:  f quasihomogeneous polynomial with an isolated singularity at 0
[f04aaf]35RETURN:  polynomial ring, possibly a simple extension of the ground field of
[35f23d]36         the basering, containing the ideal 'modid'
37         - 'modid' is the ideal of the moduli space  if opt is even (> 0).
38           otherwise it contains generators of the coordinate ring R of the
39           moduli space (note : Spec(R) is the moduli space)
40OPTIONS: 1 compute equations of the mod. space,
[fd5013]41         2 use a primary decomposition,
42         4 compute E_f0, i.e., the image of G_f0,
43         to combine options, add their value, default: opt =7
[35f23d]44EXAMPLE: example ModEqn; shows an example
45"
46{
47  int sizeOfAction, i, dimT, nonLinearQ, milnorNr, dbPrt;
48  int imageQ, opt;
49  intvec wt;
50  ideal B;
51  list Gf, tIndex, sList;
52  string ringSTR;
53
54  dbPrt = printlevel-voice+2;
55  if(size(#) > 0) { opt = #[1]; }
56  else { opt = 7; }
57  if(opt / 4 > 0) { imageQ = 1; opt = opt - 4;}
58  else { imageQ = 0; }
59
60  wt = weight(f);
61  milnorNr = vdim(std(jacob(f)));
62  if(milnorNr == -1) {
63                ERROR("the polynomial " + string(f) + " has a nonisolated singularity at 0");
64        }       // singularity not isolated
65
66  // 1st step : compute a basis of T_
67
68  B = UpperMonomials(f, wt);
69  dimT = size(B);
70  dbprint(dbPrt, "moduli equations of f = " + string(f) + ", f has Milnor number = " + string(milnorNr));
71  dbprint(dbPrt, " upper basis = " + string(B));
72  if(size(B) > 1) {
73
74    // 2nd step : compute the stabilizer G_f of f
75
76    dbprint(dbPrt, " compute equations of the stabilizer of f, called G_f");
77    Gf = StabEqn(f);
78    dbprint(dbPrt, " order of the stabilizer = " + string(StabOrder(Gf)));
79
80    // 3rd step : compute the induced action of G_f on T_ by means of a theorem of Arnold
81
82    dbprint(dbPrt, " compute the induced action");
83    def RME1 = ArnoldAction(f, Gf, B);
84    setring(RME1);
85    export(RME1);
86    dbprint(dbPrt, " G_f = " + string(stabid));
87    dbprint(dbPrt, " action of G_f : " + string(actionid));
88
89    // 4th step : linearize the action of G_f
90
91    sizeOfAction = size(actionid);
92    def RME2 = LinearizeAction(stabid, actionid, nvars(Gf[1]));
93    setring RME2;
94    export(RME2);
[3b77465]95    kill RME1;
[35f23d]96
97    if(size(actionid) == sizeOfAction) { nonLinearQ = 0;}
98    else  {
99      nonLinearQ = 1;
100      dbprint(dbPrt, " linearized action = " + string(actionid));
101      dbprint(dbPrt, " embedding of T_ = " + string(embedid));
102    }
103
104
105
106    if(!imageQ) {        // do not compute the image of Gf
107      // 5th step : set E_f = G_f,
108      dbprint(dbPrt, " compute equations of the quotient T_/G_f");
109      def RME3 = basering;
110    }
111    else {
112
113      // 5th step : compute the ideal and the action of E_f
114
115      dbprint(dbPrt, " compute E_f");
116      def RME3 = ImageGroup(groupid, actionid);
117      setring(RME3);
118      ideal embedid = imap(RME2, embedid);
119      dbprint(dbPrt, " E_f  = (" + string(groupid) + ")");
120      dbprint(dbPrt, " action of E'f = " + string(actionid));
121      dbprint(dbPrt, " compute equations of the quotient T_/E_f");
122    }
123    export(RME3);
[3b77465]124    kill RME2;
[35f23d]125
126    // 6th step : compute the equations of the quotient T_/E_f
127
128                ideal G = groupid; ideal variety = embedid;
[3b77465]129                kill groupid,embedid;
[35f23d]130    def RME4 = QuotientEquations(G, actionid, variety, opt);
131    setring RME4;
132    poly mPoly = minpoly;
[3b77465]133    kill RME3;
[35f23d]134    export(RME4);
135
136    // simplify the ideal and create a new ring with propably less variables
137
138    if(opt == 1 || opt == 3) {      // equations computed ?
139      sList = SimplifyIdeal(id, 0, "Y");
140      ideal newid = sList[1];
141      dbprint(dbPrt, " number of equations = " + string(size(sList[1])));
142      dbprint(dbPrt, " number of variables = " + string(size(sList[3])));
143      ringSTR = "ring RME5 = (" + charstr(basering) + "), (Y(1.." + string(size(sList[3])) + ")),dp;";
144      execute(ringSTR);
145      minpoly = number(imap(RME4, mPoly));
146      ideal modid = imap(RME4, newid);
147    }
148    else {
149      def RME5 = RME4;
150      setring(RME5);
151      ideal modid = imap(RME4, id);
152    }
153    export(modid);
[3b77465]154    kill RME4;
[35f23d]155  }
156  else {
157                def RME5 = basering;
158                ideal modid = maxideal(1);
159                if(size(B) == 1) {                      // 1-dimensional
160                        modid[size(modid)] = 0;
161                        modid = simplify(modid,2);
162                }
163                export(modid);
164        }
165dbprint(dbPrt, "
166// 'ModEqn' created a new ring.
167// To see the ring, type (if the name of the ring is R):
168     show(R);
[1f92589]169// To access the ideal of the moduli space of semiquasihomogeneous singularities
[35f23d]170// with principal part f, type
171     def R = ModEqn(f); setring R;  modid;
172// 'modid' is the ideal of the moduli space.
173// if 'opt' = 0 or even, then 'modid' contains algebra generators of S s.t.
174// spec(S) = moduli space of f.
175");
176  return(RME5);
177}
178example
179{"EXAMPLE:";  echo = 2;
180  ring B   = 0,(x,y), ls;
181  poly f = -x4 + xy5;
182  def R = ModEqn(f);
183  setring R;
184  modid;
185}
186
187
188///////////////////////////////////////////////////////////////////////////////
189
190proc QuotientEquations(ideal G, ideal Gaction, ideal embedding, list#)
191"USAGE:   QuotientEquations(G,action,emb [, opt]); ideal G,action,emb;int opt
192PURPOSE: compute the quotient of the variety given by the parameterization
193         'emb'  by the linear action 'action' of the algebraic group G.
194ASSUME:  'action' is linear, G must be finite if the Reynolds operator is
195         needed (i.e., NullCone(G,action) returns some non-invariant polys)
[f04aaf]196RETURN:   polynomial ring over a simple extension of the ground field of the
[35f23d]197          basering, containing the ideals 'id' and 'embedid'.
[1f92589]198          - 'id' contains the equations of the quotient, if opt = 1;
199            if opt = 0, 2, 'id' contains generators of the coordinate ring R
[8942a5]200            of the quotient (Spec(R) is the quotient)
[1f92589]201          - 'embedid' = 0, if opt = 1;
202            if opt = 0, 2, it is the ideal defining the equivariant embedding
[35f23d]203OPTIONS: 1 compute equations of the quotient,
[fd5013]204         2 use a primary decomposition when computing the Reynolds operator,@*
205         to combine options, add their value, default: opt =3.
[35f23d]206EXAMPLE: example QuotientEquations; shows an example
207"
208{
209  int i, opt, primaryDec, relationsQ, dbPrt;
210  ideal Gf, variety;
211  intvec wt;
212
213  dbPrt = printlevel-voice+3;
214  if(size(#) > 0) { opt = #[1]; }
215  else { opt = 3; }
216
217  if(opt / 2 > 0) { primaryDec = 1; opt = opt - 2; }
218  else { primaryDec = 0; }
219  if(opt > 0) { relationsQ = 1;}
220  else { relationsQ = 0; }
221
222  Gf = std(G);
223  variety = EquationsOfEmbedding(embedding, nvars(basering) - size(Gaction));
224
225  if(size(variety) == 0) {    // use Hilbert function !
[9e9ec3]226    //for(i = 1; i <= ncols(Gaction); i ++) { wt[i] = 1;};
227    for(i = 1; i <= nvars(basering); i ++) { wt[i] = 1;};
[35f23d]228  }
229  def RQER = InvariantRing(Gf, Gaction, primaryDec);    // compute the nullcone of the linear action
230
231  def RQEB = basering;
232  setring(RQER);
233  export(RQER);
234
235  if(relationsQ > 0) {
236    dbprint(dbPrt, " compute equations of the variety (" + string(size(imap(RQER, invars))) + " invariants) ");
237    if(!defined(variety)) { ideal variety = imap(RQEB, variety); }
238    if(wt[1] > 0) {
239      def RQES = ImageVariety(variety, imap(RQER, invars), wt);
240    }
241    else {
242      def RQES = ImageVariety(variety, imap(RQER, invars));  // forget imap
243    }
244    setring(RQES);
245    ideal id = imageid;
246    ideal embedid = 0;
247  }
248  else {
249    def RQES = basering;
250    ideal id =  imap(RQER, invars);
251    ideal embedid = imap(RQEB, variety);
252  }
[3b77465]253  kill RQER;
[35f23d]254  export(id);
255  export(embedid);
256  return(RQES);
257}
258
259///////////////////////////////////////////////////////////////////////////////
260
261proc UpperMonomials(poly f, list #)
262"USAGE:   UpperMonomials(poly f, [intvec w])
263PURPOSE: compute the upper monomials of the milnor algebra of f.
[1f92589]264ASSUME:  f is quasihomogeneous (w.r.t. w)
[35f23d]265RETURN:  ideal
266EXAMPLE: example UpperMonomials; shows an example
267"
268{
269  int i,d;
270  intvec wt;
271  ideal I, J;
272
273  if(size(#) == 0) { wt = weight(f);}
274  else { wt = #[1];}
275   J = kbase(std(jacob(f)));
276  d = deg(f, wt);
277  for(i = 1; i <= size(J); i++) { if(deg(J[i], wt) > d) {I = I, J[i];} }
278  return(simplify(I, 2));
279}
280example
281{"EXAMPLE:";  echo = 2;
282  ring B   = 0,(x,y,z), ls;
283  poly f = -z5+y5+x2z+x2y;
284  UpperMonomials(f);
285}
286
287///////////////////////////////////////////////////////////////////////////////
288
289proc ArnoldAction(poly f, list #)
290"USAGE:   ArnoldAction(f, [Gf, B]); poly f; list Gf, B;
291         'Gf' is a list of two rings (coming from 'StabEqn')
[1f92589]292PURPOSE: compute the induced action of the stabilizer G of f on T_, where
[35f23d]293         T_ is given by the upper monomials B of the Milnor algebra of f.
[1f92589]294ASSUME:  f is quasihomogeneous
[f04aaf]295RETURN:  polynomial ring over the same ground field, containing the ideals
[35f23d]296         'actionid' and 'stabid'.
[1f92589]297         - 'actionid' is the ideal defining the induced action of Gf on T_ @*
[35f23d]298         - 'stabid' is the ideal of the stabilizer Gf in the new ring
299EXAMPLE: example ArnoldAction; shows an example
300"
301{
302  int i, offset, ub, pos, nrStabVars, dbPrt;
303  intvec wt = weight(f);
304  ideal B;
305  list Gf, parts, baseDeg;
306  string ringSTR1, ringSTR2, parName, namesSTR, varSTR;
307
308  dbPrt = printlevel-voice+2;
309  if(size(#) == 0) {
310    Gf = StabEqn(f);
311    B = UpperMonomials(f, wt);
312  }
313  else {
314    Gf = #[1];
315    if(size(#) > 1) { B = #[2];}
316    else {B = UpperMonomials(f, wt);}
317  }
318  if(size(B) == 0) { ERROR("the principal part " + string(f) + " has no upper monomials");}
319  for(i = 1; i <= size(B); i = i + 1) {
320    baseDeg[i] = deg(B[i], wt);
321  }
322  ub = Max(baseDeg) + 1;          // max degree of an upper mono.
323  def RAAB = basering;
324  def STR1 = Gf[1];
325  def STR2 = Gf[2];
326  nrStabVars = nvars(STR1);
327
328  dbprint(dbPrt, "ArnoldAction of f = ", f, ", upper base = " + string(B));
329
330  setring STR1;
331  poly mPoly = minpoly;
332  setring RAAB;
333
334  // setup new ring with s(..) and t(..) as parameters
335
336  varSTR = string(maxideal(1));
337  ringSTR2 = "ring RAAS = ";
338  if(npars(basering) == 1) {
339    parName = parstr(basering);
340    ringSTR2 = ringSTR2 + "(0, " + parstr(1) + "), ";
341  }
342  else {
343    parName = "a";
344    ringSTR2 = ringSTR2 + "0, ";
345  }
346  offset = 1 + nrStabVars;
347  namesSTR = "s(1.." + string(nrStabVars) + "), t(1.." + string(size(B)) + ")";
348  ringSTR2 = ringSTR2 + "(" + namesSTR + "), lp;";
349  ringSTR1 = "ring RAAR = (0, " + parName + "," + namesSTR + "), (" + varSTR + "), ls;";  // lp ?
350
351  execute(ringSTR1);
352  export(RAAR);
353  ideal upperBasis, stabaction, action, mPoly, reduceIdeal;
354  poly f, F, monos, h;
355
356  reduceIdeal = imap(STR1, mPoly), imap(STR1, stabid);
357  f = imap(RAAB, f);
358  F = f;
359  upperBasis = imap(RAAB, B);
360  for(i = 1; i <= size(upperBasis); i = i + 1) {
361    F = F + par(i + offset)*upperBasis[i];
362  }
363  monos = F - f;
364  stabaction = imap(STR2, actionid);
365
366  // action of the stabilizer on the semiuniversal unfolding of f
367
368  F = f + APSubstitution(monos, stabaction, reduceIdeal, wt, ub, nrStabVars, size(upperBasis));
369
370  // apply the theorem of Arnold
371
372  h = ArnoldFormMain(f, upperBasis, F, reduceIdeal, nrStabVars, size(upperBasis)) - f;
373
374  // extract the polynomials of the action of the stabilizer on T_
375
376  parts = MonosAndTerms(h, wt, ub);
377  for(i = 1; i <= size(parts[1]); i = i + 1) {
378    pos = FirstEntryQHM(upperBasis, parts[1][i]);
379    action[pos] = parts[2][i]/parts[1][i];
380  }
381  execute(ringSTR2);
382  minpoly = number(imap(STR1, mPoly));
383  ideal actionid = imap(RAAR, action);
384  ideal stabid = imap(STR1, stabid);
385  export(actionid);
386  export(stabid);
[3b77465]387  kill RAAR;
[35f23d]388dbprint(dbPrt, "
389// 'ArnoldAction' created a new ring.
390// To see the ring, type (if the name of the ring is R):
391     show(R);
392// To access the ideal of the stabilizer G of f and its group action,
[1f92589]393// where f is the quasihomogeneous principal part, type
[35f23d]394     def R = ArnoldAction(f); setring R;  stabid; actionid;
395// 'stabid' is the ideal of the group G and 'actionid' is the ideal defining
396// the group action of the group G on T_. Note: this action might be nonlinear
397");
398  return(RAAS);
399}
400example
401{"EXAMPLE:";  echo = 2;
402  ring B   = 0,(x,y,z), ls;
403  poly f = -z5+y5+x2z+x2y;
404  def R = ArnoldAction(f);
405  setring R;
406  actionid;
407  stabid;
408}
409
410///////////////////////////////////////////////////////////////////////////////
411
412proc StabOrder(list #)
[fd5013]413"USAGE:   StabOrder(f); poly f
[35f23d]414PURPOSE: compute the order of the stabilizer group of f.
[1f92589]415ASSUME:  f quasihomogeneous polynomial with an isolated singularity at 0
[35f23d]416RETURN:  int
417GLOBAL: varSubsList
418"
419{
420  list stab;
421
422  if(size(#) == 1) { stab = StabEqn(#[1]); }
423  else {  stab = #;}
424
425  def RSTO = stab[1];
426  setring(RSTO);
427  return(vdim(std(stabid)));
428}
429
430///////////////////////////////////////////////////////////////////////////////
431
432proc StabEqn(poly f)
433"USAGE:   StabEqn(f); f polynomial
434PURPOSE: compute the equations of the isometry group of f.
[1f92589]435ASSUME:  f semiquasihomogeneous polynomial with an isolated singularity at 0
[fd5013]436RETURN:  list of two rings 'S1', 'S2'
[1f92589]437         - 'S1' contians the equations of the stabilizer (ideal 'stabid') @*
[35f23d]438         - 'S2' contains the action of the stabilizer (ideal 'actionid')
439EXAMPLE: example StabEqn; shows an example
440GLOBAL: varSubsList, contains the index j s.t. x(i) -> x(i)t(j) ...
441"
442{
443dbprint(dbPrt, "
444// 'StabEqn' created a list of 2 rings.
445// To see the rings, type (if the name of your list is stab):
446     show(stab);
447// To access the 1-st ring and map (and similair for the others), type:
448     def S1 = stab[1]; setring S1;  stabid;
449// S1/stabid is the coordinate ring of the variety of the
450// stabilizer, say G. If G x K^n --> K^n is the action of G on
451// K^n, then the ideal 'actionid' in the second ring describes
452// the dual map on the ring level.
453// To access the 2-nd ring and map (and similair for the others), type:
454     def S2 = stab[2]; setring S2;  actionid;
455");
456
457        return(StabEqnId(ideal(f), qhweight(f)));
458}
459example
460{"EXAMPLE:";  echo = 2;
461  ring B = 0,(x,y,z), ls;
462  poly f = -z5+y5+x2z+x2y;
463  list stab = StabEqn(f);
464  def S1 = stab[1]; setring S1;  stabid;
465  def S2 = stab[2]; setring S2;  actionid;
466}
467
468///////////////////////////////////////////////////////////////////////////////
469
470proc StabEqnId(ideal data, intvec wt)
471"USAGE:   StabEqn(I, w); I ideal, w intvec
[fd5013]472PURPOSE: compute the equations of the isometry group of the ideal I,
[35f23d]473         each generator of I is fixed by the stabilizer.
[1f92589]474ASSUME:  I semiquasihomogeneous ideal wrt 'w' with an isolated singularity at 0
[fd5013]475RETURN:  list of two rings 'S1', 'S2'
[1f92589]476         - 'S1' contians the equations of the stabilizer (ideal 'stabid') @*
[35f23d]477         - 'S2' contains the action of the stabilizer (ideal 'actionid')
478EXAMPLE: example StabEqnId; shows an example
479GLOBAL: varSubsList, contains the index j s.t. t(i) -> t(i)t(j) ...
480"
481{
[e7778a]482  int i, j, c, k, r, nrVars, offset, n, sln, dbPrt;
[afd77b2]483  list Variables, rd, temp, sList, varSubsList;
[35f23d]484  poly mPoly;
485  string ringSTR, ringSTR1, varString, parString;
486
487  dbPrt = printlevel-voice+2;
488  dbprint(dbPrt, "StabilizerEquations of " + string(data));
489
490  export(varSubsList);
491  n = nvars(basering);
[afd77b2]492  Variables = StabVar(wt);    // possible quasihomogeneous substitutions
[35f23d]493  nrVars = 0;
494  for(i = 1; i <= size(wt); i = i + 1) {
[afd77b2]495    nrVars = nrVars + size(Variables[i]);
[35f23d]496  }
497
498  // set the new basering needed for the substitutions
499
500  varString = "s(1.." + string(nrVars) + ")";
501  if(npars(basering) == 1) {
502    parString = "(0, " + parstr(basering) + ")";
503  }
504  else { parString = "0"; }
505
506  def RSTB = basering;
507  mPoly = minpoly;
508  ringSTR = "ring RSTR = " + parString + ", (" + varstr(basering) + ", " + varString + "), dp;";  // dp
509        ringSTR1 = "ring RSTT = " + parString + ", (" + varString + ", " + varstr(basering) + "), dp;";
510
[3b77465]511  if(defined(RSTR)) { kill RSTR;}
512        if(defined(RSTT)) { kill RSTT;}
[35f23d]513        execute(ringSTR1);      // this ring is only used for the result, where the variables
514  export(RSTT);           // are s(1..m),t(1..n), as needed for Derksens algorithm (NullCone)
515  minpoly = number(imap(RSTB, mPoly));
516
517  execute(ringSTR);
518  export(RSTR);
519  minpoly = number(imap(RSTB, mPoly));
520  poly f, f1, g, h, vars, pp;      // f1 is the polynomial after subs,
521  ideal allEqns, qhsubs, actionid, stabid, J;
522  list ringList;          // all t(i)`s which do not appear in f1
523  ideal data = simplify(imap(RSTB, data), 2);
524
[1f92589]525  // generate the quasihomogeneous substitution map F
[35f23d]526
527  nrVars = 0;
528  offset = 0;
529  for(i = 1; i <= size(wt); i = i + 1) {    // build the substitution t(i) -> ...
[afd77b2]530    if(i > 1) { offset = offset + size(Variables[i - 1]); }
[35f23d]531    g = 0;
[afd77b2]532    for(j = 1; j <= size(Variables[i]); j = j + 1) {
[35f23d]533      pp = 1;
[afd77b2]534      for(k = 2; k <= size(Variables[i][j]); k = k + 1) {
535        pp = pp * var(Variables[i][j][k]);
536        if(Variables[i][j][k] == i) { varSubsList[i] = offset + j;}
[35f23d]537      }
538      g = g + s(offset + j) * pp;
539    }
540    qhsubs[i] = g;
541  }
542  dbprint(dbPrt, "  qhasihomogenous substituion =" + string(qhsubs));
543  map F = RSTR, qhsubs;
[3b77465]544  kill varSubsList;
[35f23d]545
546  // get the equations of the stabilizer by comparing coefficients
547  // in the equation f = F(f).
548
549  vars = RingVarProduct(Table("i", "i", 1, size(wt)));
550
551  allEqns = 0;
552
[e7778a]553  matrix newcoMx, coMx;
554  int d;
555  for(r = 1; r <= ncols(data); r++)
556  {
[35f23d]557
558  f = data[r];
559  f1 = F(f);
[7de8e4]560  d = deg(f);
561  newcoMx = coef(f1, vars);        // coefficients of F(f)
562  coMx = coef(f, vars);          // coefficients of f
[35f23d]563
564  for(i = 1; i <= ncols(newcoMx); i = i + 1) {      // build the system of eqns via coeff. comp.
565    j = 1;
566    h = 0;
567    while(j <= ncols(coMx)) {        // all monomials in f
568      if(coMx[j][1] == newcoMx[i][1]) { h = coMx[j][2]; j = ncols(coMx) + 1;}
569      else {j = j + 1;}
570    }
571    J = J, newcoMx[i][2] - h;        // add equation
572  }
573  allEqns =  allEqns, J;
574
575  }
576  allEqns = std(allEqns);
577
578  // simplify the equations, i.e., if s(i) in J then remove s(i) from J
579  // and from the basering
580
581  sList = SimplifyIdeal(allEqns, n, "s");
582  stabid = sList[1];
583  map phi = basering, sList[2];
584        sln = size(sList[3]) - n;
585
586  // change the substitution
587
588  actionid = phi(qhsubs);
589
590        // change to new ring, auxillary construction
591
592        setring(RSTT);
593        ideal actionid, stabid;
594
595        actionid = imap(RSTR, actionid);
596        stabid = imap(RSTR, stabid);
597        export(stabid);
598  export(actionid);
599  ringList[2] = RSTT;
600
601  dbprint(dbPrt, "  substitution = " + string(actionid));
602  dbprint(dbPrt, "  equations of stabilizer = " + string(stabid));
603
604  varString = "s(1.." + string(sln) + ")";
605  ringSTR = "ring RSTS = " + parString + ", (" + varString + "), dp;";
606  execute(ringSTR);
607  minpoly = number(imap(RSTB, mPoly));
608  ideal stabid = std(imap(RSTR, stabid));
609  export(stabid);
610  ringList[1] = RSTS;
611dbprint(dbPrt, "
612// 'StabEqnId' created a list of 2 rings.
613// To see the rings, type (if the name of your list is stab):
614     show(stab);
615// To access the 1-st ring and map (and similair for the others), type:
616     def S1 = stab[1]; setring S1;  stabid;
617// S1/stabid is the coordinate ring of the variety of the
618// stabilizer, say G. If G x K^n --> K^n is the action of G on
619// K^n, then the ideal 'actionid' in the second ring describes
620// the dual map on the ring level.
621// To access the 2-nd ring and map (and similair for the others), type:
622     def S2 = stab[2]; setring S2;  actionid;
623");
624  return(ringList);
625}
626example
627{"EXAMPLE:";  echo = 2;
628  ring B   = 0,(x,y,z), ls;
629  ideal I = x2,y3,z6;
630  intvec w = 3,2,1;
631  list stab = StabEqnId(I, w);
632  def S1 = stab[1]; setring S1;  stabid;
633  def S2 = stab[2]; setring S2;  actionid;
634}
635
636///////////////////////////////////////////////////////////////////////////////
637static
638proc ArnoldFormMain(poly f, B, poly Fs, ideal reduceIdeal, int nrs, int nrt)
639"USAGE:   ArnoldFormMain(f, B, Fs, rI, nrs, nrt);
640   poly f,Fs; ideal B, rI; int nrs, nrt
641PURPOSE: compute the induced action of 'G_f' on T_, where f is the principal
642         part and 'Fs' is the semiuniversal unfolding of 'f' with x_i
643         substituted by actionid[i], 'B' is a list of upper basis monomials
644         for the milnor algebra of 'f', 'nrs' = number of variables for 'G_f'
645         and 'nrt' = dimension of T_
[1f92589]646ASSUME:  f is quasihomogeneous with an isolated singularity at 0,
[35f23d]647         s(1..r), t(1..m) are parameters of the basering
648RETURN:  poly
649EXAMPLE: example ArnoldAction; shows an example
650"
651{
652  int i, j, d, ub, dbPrt;
653  list upperBasis, basisDegList, gmonos, common, parts;
654  ideal jacobianId, jacobIdstd, mapId;    // needed for phi
655  intvec wt = weight(f);
656  matrix gCoeffMx;        // for lift command
657  poly newFs, g, gred, tt;        // g = sum of all monomials of degree d, gred is needed for lift
658  map phi;          // the map from Arnold's Theorem
659
660  dbPrt = printlevel-voice+2;
661  jacobianId = jacob(f);
662  jacobIdstd = std(jacobianId);
663  newFs = Fs;
664  for(i = 1; i <= size(B); i = i + 1) {
665    basisDegList[i] = deg(B[i], wt);
666  }
667  ub = Max(basisDegList) + 1;          // max degree of an upper monomial
668
669  parts = MonosAndTerms(newFs - f, wt, ub);
670  gmonos = parts[1];
671  d = deg(f, wt);
672
673  for(i = d + 1; i < ub; i = i + 1) {    // base[1] = monomials of degree i
674    upperBasis[i] = SelectMonos(list(B, B), wt, i);    // B must not contain 0's
675  }
676
677  // test if each monomial of Fs is contained in B, if not,
678  // compute a substitution via Arnold's theorem and substitutite
679  // it into newFs
680
681  for(i = d + 1; i < ub; i = i + 1) {  // ub instead of @UB
682    dbprint(dbPrt, "-- degree = " + string(i) + " of " + string(ub - 1) + " ---------------------------");
683    if(size(newFs) < 80) { dbprint(dbPrt, "  polynomial = " + string(newFs - f));}
684    else {  dbprint(dbPrt, "  poly has deg (not weighted) " + string(deg(newFs)) + " and contains " + string(size(newFs)) + " monos");}
685
686    // select monomials of degree i and intersect them with upperBasis[i]
687
688    gmonos = SelectMonos(parts, wt, i);
689    common = IntersectionQHM(upperBasis[i][1], gmonos[1]);
690    if(size(common) == size(gmonos[1])) {
691      dbprint(dbPrt, " no additional monomials ");
692    }
693
694    // other monomials than those in upperBasis occur, compute
695    // the map constructed in the proof of Arnold's theorem
696    // write g = c[i] * jacobianId[i]
697
698    else {
699      dbprint(dbPrt, "  additional Monomials found, compute the map ");
700      g = PSum(gmonos[2]);      // sum of all monomials in g of degree i
701      dbprint(dbPrt, "  sum of degree " + string(i) + " is " + string(g));
702
703      gred = reduce(g, jacobIdstd);
704      gCoeffMx = lift(jacobianId, g - gred);    // compute c[i]
705      mapId = var(1) - gCoeffMx[1][1];    // generate the map
706      for(j = 2; j <= size(gCoeffMx); j = j + 1) {
707        mapId[j] = var(j) - gCoeffMx[1][j];
708      }
709      dbprint(dbPrt, "  map = " + string(mapId));
710      // apply the map to newFs
711      newFs = APSubstitution(newFs, mapId, reduceIdeal, wt, ub, nrs, nrt);
712      parts = MonosAndTerms(newFs - f, wt, ub);  // monos and terms of deg < ub
713      newFs = PSum(parts[2]) + f;      // result of APS... is already reduced
714      dbprint(dbPrt, "  monomials of degree " + string(i));
715    }
716  }
717  return(newFs);
718}
719
720///////////////////////////////////////////////////////////////////////////////
721
[70597c]722static proc MonosAndTerms(poly f, wt, int ub)
[35f23d]723"USAGE:   MonosAndTerms(f, w, ub); poly f, intvec w, int ub
724PURPOSE: returns a list of all monomials and terms occuring in f of
725         weighted degree < ub
726RETURN:  list
727         _[1]  list of monomials
[8942a5]728         _[2]  list of terms
[35f23d]729EXAMPLE: example MonosAndTerms shows an example
730"
731{
732  int i, k;
733  list monomials, terms;
734  poly mono, lcInv, data;
735
736  data = jet(f, ub - 1, wt);
737  k = 0;
738  for(i = 1; i <= size(data); i = i + 1) {
739    mono = lead(data[i]);
740    if(deg(mono, wt) < ub) {
741      k = k + 1;
742      lcInv = 1/leadcoef(mono);
743      monomials[k] = mono * lcInv;
744      terms[k] = mono;
745    }
746  }
747  return(list(monomials, terms));
748}
749example
750{"EXAMPLE:";  echo = 2;
751  ring B = 0,(x,y,z), lp;
752  poly f = 6*x2 + 2*x3 + 9*x*y2 + z*y + x*z6;
753  MonosAndTerms(f, intvec(2,1,1), 5);
754}
755
756///////////////////////////////////////////////////////////////////////////////
757
[70597c]758static proc SelectMonos(parts, intvec wt, int d)
[35f23d]759"USAGE:   SelectMonos(parts, w, d); list/ideal parts, intvec w, int d
760PURPOSE: returns a list of all monomials and terms occuring in f of
761         weighted degree = d
762RETURN:  list
[8942a5]763         _[1]  list of monomials
764         _[2]  list of terms
[35f23d]765EXAMPLE: example SelectMonos; shows an example
766"
767{
768  int i, k;
769  list monomials, terms;
770  poly mono;
771
772  k = 0;
773  for(i = 1; i <= size(parts[1]); i = i + 1) {
774    mono = parts[1][i];
775    if(deg(mono, wt) == d) {
776      k = k + 1;
777      monomials[k] = mono;
778      terms[k] = parts[2][i];
779    }
780  }
781  return(list(monomials, terms));
782}
783example
784{"EXAMPLE:";  echo = 2;
785  ring B = 0,(x,y,z), lp;
786  poly f = 6*x2 + 2*x3 + 9*x*y2 + z*y + x*z6;
787  list mt =  MonosAndTerms(f, intvec(2,1,1), 5);
788  SelectMonos(mt, intvec(2,1,1), 4);
789}
790
791///////////////////////////////////////////////////////////////////////////////
792
[70597c]793static proc Expand(substitution, degVec, ideal reduceI, intvec w1, int ub, list truncated)
[35f23d]794"USAGE:   Expand(substitution, degVec, reduceI, w, ub, truncated);
795         ideal/list substitution, list/intvec degVec, ideal reduceI, intvec w,
796         int ub, list truncated
797PURPOSE: substitute 'substitution' in the monomial given by the list of
798         exponents 'degVec', omit all terms of weighted degree > ub and reduce
799         the result w.r.t. 'reduceI'. If truncated[i] = 0 then the result is
800         stored for later use.
801RETURN:  poly
802NOTE:    used by APSubstitution
803GLOBAL:  computedPowers
804"
805{
806  int i, minDeg;
807  list powerList;
808  poly g, h;
809
810  // compute substitution[1]^degVec[1],...,subs[n]^degVec[n]
811
812  for(i = 1; i <= ncols(substitution); i = i + 1) {
813    if(size(substitution[i]) < 3 || degVec[i] < 4) {
814      powerList[i] = reduce(substitution[i]^degVec[i], reduceI); // new
815    }  // directly for small exponents
816    else {
817      powerList[i] = PolyPower1(i, substitution[i], degVec[i], reduceI, w1, truncated[i], ub);
818    }
819  }
820  // multiply the terms obtained by using PolyProduct();
821  g = powerList[1];
822  minDeg = w1[1] * degVec[1];
823  for(i = 2; i <= ncols(substitution); i = i + 1) {
824    g = jet(g, ub - w1[i] * degVec[i] - 1, w1);
825    h = jet(powerList[i], ub - minDeg - 1, w1);
826    g = PolyProduct(g, h, reduceI, w1, ub);
827    if(g == 0) { Print(" g = 0 "); break;}
828    minDeg = minDeg + w1[i] * degVec[i];
829  }
830  return(g);
831}
832
833///////////////////////////////////////////////////////////////////////////////
834
[70597c]835static proc PolyProduct(poly g1, poly h1, ideal reduceI, intvec wt, int ub)
[35f23d]836"USAGE:   PolyProduct(g, h, reduceI, wt, ub); poly g, h; ideal reduceI,
837          intvec wt, int ub.
838PURPOSE: compute g*h and reduce it w.r.t 'reduceI' and omit terms of weighted
839         degree > ub.
840RETURN:  poly
841NOTE:    used by 'Expand'
842"
843{
844  int SUBSMAXSIZE = 3000;    //
845  int i, nrParts, sizeOfPart, currentPos, partSize, maxSIZE;
846  poly g, h, gxh, prodComp, @g2;    // replace @g2 by subst.
847
848  g = g1;
849  h = h1;
850
851  if(size(g)*size(h) > SUBSMAXSIZE) {
852
853    // divide the polynomials with more terms in parts s.t.
854    // the product of each part with the other polynomial
855    // has at most SUBMAXSIZE terms
856
857    if(size(g) < size(h)) { poly @h = h; h = g; g = @h;@h = 0; }
858    maxSIZE = SUBSMAXSIZE / size(h);
859    //print(" SUBSMAXSIZE = "+string(SUBSMAXSIZE)+" exceeded by "+string(size(g)*size(h)) + ", maxSIZE = ", string(maxSIZE));
860    nrParts = size(g) / maxSIZE + 1;
861    partSize = size(g) / nrParts;
862    gxh = 0;  // 'g times h'
863    for(i = 1; i <= nrParts; i = i + 1){
864      //print(" loop #" + string(i) + " of " + string(nrParts));
865      currentPos = (i - 1) * partSize;
866      if(i < nrParts) {sizeOfPart = partSize;}
867      else { sizeOfPart = size(g) - (nrParts - 1) * partSize; print(" last #" + string(sizeOfPart) + " terms ");}
868      prodComp = g[currentPos + 1..sizeOfPart + currentPos] * h;  // multiply a part
869      @g2 = jet(prodComp, ub - 1, wt);  // eventual reduce ...
870      if(size(@g2) < size(prodComp)) { print(" killed " + string(size(prodComp) - size(@g2)) + " terms ");}
871      gxh =  reduce(gxh + @g2, reduceI);
872
873    }
874  }
875  else {
876    gxh = reduce(jet(g * h,ub - 1, wt), reduceI);
877  }  // compute directly
878  return(gxh);
879}
880
881///////////////////////////////////////////////////////////////////////////////
882
[3c4dcc]883static proc PolyPower1(int varIndex, poly f, int e, ideal reduceI, intvec wt,
[70597c]884                       int truncated, int ub)
[35f23d]885"USAGE:   PolyPower1(i, f, e, reduceI, wt, truncated, ub);int i, e, ub;poly f;
886         ideal reduceI; intvec wt; list truncated;
887PURPOSE: compute f^e, use previous computations if possible, and reduce it
888         w.r.t reudecI and omit terms of weighted degree > ub.
889RETURN:  poly
890NOTE:    used by 'Expand'
891GLOBAL:  'computedPowers'
892"
893{
894  int i, ordOfg, lb, maxPrecomputedPower;
895  poly g, fn;
896
897  if(e == 0) { return(1);}
898  if(e == 1) { return(f);}
899  if(f == 0) { return(1); }
900  else {
901
902    // test if f has been computed to some power
903
904    if(computedPowers[varIndex][1] > 0) {
905      maxPrecomputedPower = computedPowers[varIndex][1];
906      if(maxPrecomputedPower >= e) {
907        // no computation necessary, f^e has already benn computed
908        g = computedPowers[varIndex][2][e - 1];
909        //Print("No computation, from list : g = elem [", varIndex, ", 2, ", e - 1, "]");
910        lb = e + 1;
911      }
912      else {  // f^d computed, where d < e
913        g = computedPowers[varIndex][2][maxPrecomputedPower - 1];
914        ordOfg = maxPrecomputedPower * wt[varIndex];
915        lb = maxPrecomputedPower + 1;
916      }
917    }
918    else {    // no precomputed data
919      lb = 2;
920      ordOfg = wt[varIndex];
921      g = f;
922    }
923    for(i = lb; i <= e; i = i + 1) {
924      fn = jet(f, ub - ordOfg - 1, wt); // reduce w.r.t. reduceI
925      g = PolyProduct(g, fn, reduceI, wt, ub);
926      ordOfg = ordOfg + wt[varIndex];
927      if(g == 0) { break; }
928      if((i > maxPrecomputedPower) && !truncated) {
929        if(maxPrecomputedPower == 0) {  // init computedPowers
930          computedPowers[varIndex] = list(i, list(g));
931        }
932        computedPowers[varIndex][1] = i;  // new degree
933        computedPowers[varIndex][2][i - 1] = g;
934        maxPrecomputedPower = i;
935      }
936    }
937  }
938  return(g);
939}
940
941///////////////////////////////////////////////////////////////////////////////
942
[70597c]943static proc RingVarsToList(list @index)
[35f23d]944{
945  int i;
946  list temp;
947
948  for(i = 1; i <= size(@index); i = i + 1) { temp[i] = string(var(@index[i])); }
949  return(temp);
950}
951
952///////////////////////////////////////////////////////////////////////////////
[70597c]953static
[35f23d]954proc APSubstitution(poly f, ideal substitution, ideal reduceIdeal, intvec wt, int ub, int nrs, int nrt)
955"USAGE:   APSubstitution(f, subs, reduceI, w, ub, int nrs, int nrt); poly f
956         ideal subs, reduceI, intvec w, int ub, nrs, nrt;
957         nrs = number of parameters s(1..nrs),
958         nrt = number of parameters t(1..nrt)
959PURPOSE: substitute 'subs' in f, omit all terms with weighted degree > ub and
960         reduce the result w.r.t. 'reduceI'.
961RETURN:  poly
962GLOBAL:  'computedPowers'
963"
964{
965  int i, j, k, d, offset;
966  int n = nvars(basering);
967  list  coeffList, parts, degVecList, degOfMonos;
968  list computedPowers, truncatedQ, degOfSubs, @temp;
969  string ringSTR, @ringVars;
970
971  export(computedPowers);
972
973  // store arguments in strings
974
975  def RASB = basering;
976
977  parts = MonosAndTerms(f, wt, ub);
978  for(i = 1; i <= size(parts[1]); i = i + 1) {
979    coeffList[i] = parts[2][i]/parts[1][i];
980    degVecList[i] = leadexp(parts[1][i]);
981    degOfMonos[i] = deg(parts[1][i], wt);
982  }
983
984  // built new basering with no parameters and order dp !
985  // the parameters of the basering are appended to
986  // the variables of the basering !
987  // set ideal mpoly = minpoly for reduction !
988
989  @ringVars = "(" + varstr(basering) + ", " + parstr(1) + ",";  // precondition
990  if(nrs > 0) {
991    @ringVars = @ringVars + "s(1.." + string(nrs) + "), ";
992  }
993  @ringVars = @ringVars + "t(1.." + string(nrt) + "))";
994  ringSTR = "ring RASR = 0, " + @ringVars + ", dp;";    // new basering
995
996  // built the "reduction" ring with the reduction ideal
997
998  execute(ringSTR);
999  export(RASR);
1000  ideal reduceIdeal, substitution, newSubs;
1001  intvec w1, degVec;
1002  list minDeg, coeffList, degList;
1003  poly f, g, h, subsPoly;
1004
1005  w1 = wt;            // new weights
1006  offset = nrs + nrt + 1;
1007  for(i = n + 1; i <= offset + n; i = i + 1) { w1[i] = 0; }
1008
1009  reduceIdeal = std(imap(RASB, reduceIdeal)); // omit later !
1010  coeffList = imap(RASB, coeffList);
1011  substitution = imap(RASB, substitution);
1012
1013  f = imap(RASB, f);
1014
1015  for(i = 1; i <= n; i = i + 1) {      // all "base" variables
1016    computedPowers[i] = list(0);
1017    for(j = 1; j <= size(substitution[i]); j = j + 1) { degList[j] = deg(substitution[i][j], w1);}
1018    degOfSubs[i] = degList;
1019  }
1020
1021  // substitute in each monomial seperately
1022
1023  g = 0;
1024  for(i = 1; i <= size(degVecList); i = i + 1) {
1025    truncatedQ = Table("0", "i", 1, n);
1026    newSubs = 0;
1027    degVec = degVecList[i];
1028    d = degOfMonos[i];
1029
1030    // check if some terms in the substitution can be omitted
1031    // degVec = list of exponents of the monomial m
1032    // minDeg[j] denotes the weighted degree of the monomial m'
1033    // where m' is the monomial m without the j-th variable
1034
1035    for(j = 1; j <= size(degVec); j = j + 1) { minDeg[j] = d - degVec[j] * wt[j]; }
1036
1037    for(j = 1; j <= size(degVec); j = j + 1) {
1038      subsPoly = 0;      // set substitution to 0
1039      if(degVec[j] > 0) {
1040
1041        // if variable occurs then check if
1042        // substitution[j][k] * (linear part)^(degVec[j]-1) + minDeg[j] < ub
1043        // i.e. look for the smallest possible combination in subs[j]^degVec[j]
1044        // which comes from the term substitution[j][k]. This term is multiplied
1045        // with the rest of the monomial, which has at least degree minDeg[j].
1046        // If the degree of this product is < ub then substitution[j][k] contributes
1047        // to the result and cannot be omitted
1048
1049        for(k = 1; k <= size(substitution[j]); k = k + 1) {
1050          if(degOfSubs[j][k] + (degVec[j] - 1) * wt[j] + minDeg[j]  < ub) {
1051            subsPoly = subsPoly + substitution[j][k];
1052          }
1053        }
1054      }
1055      newSubs[j] = subsPoly;        // set substitution
1056      if(substitution[j] - subsPoly != 0) { truncatedQ[j] = 1;}  // mark that substitution[j] is truncated
1057    }
1058    h = Expand(newSubs, degVec, reduceIdeal, w1, ub, truncatedQ) * coeffList[i];  // already reduced
1059    g = reduce(g + h, reduceIdeal);
1060  }
[3b77465]1061  kill computedPowers;
[35f23d]1062  setring RASB;
1063  poly fnew = imap(RASR, g);
[3b77465]1064  kill RASR;
[35f23d]1065  return(fnew);
1066}
1067
1068///////////////////////////////////////////////////////////////////////////////
1069
1070static proc StabVar(intvec wt)
1071"USAGE:   StabVar(w); intvec w
[1f92589]1072PURPOSE: compute the indicies for quasihomogeneous substitutions of each
[35f23d]1073         variable.
[1f92589]1074ASSUME:  f semiquasihomogeneous polynomial with an isolated singularity at 0
[35f23d]1075RETURN:  list
1076         _[i]  list of combinations for var(i) (i must be appended
1077         to each comb)
1078GLOBAL:  'varSubsList', contains the index j s.t. x(i) -> x(i)t(j) ...
1079"
1080{
1081  int i, j, k, uw, ic;
[afd77b2]1082  list varList, Variables, subs;
[35f23d]1083  string str, varString;
1084
1085  varList = StabVarComb(wt);
1086  for(i = 1; i <= size(wt); i = i + 1) {
1087    subs = 0;
1088
1089    // built linear substituitons
1090    for(j = 1; j <= size(varList[1][i]); j = j + 1) {
1091      subs[j] = list(i) + list(varList[1][i][j]);
1092    }
[afd77b2]1093    Variables[i] = subs;
[35f23d]1094    if(size(varList[2][i]) > 0) {
1095
1096      // built nonlinear substituitons
1097      subs = 0;
1098      for(j = 1; j <= size(varList[2][i]); j = j + 1) {
1099        subs[j] = list(i) + varList[2][i][j];
1100      }
[afd77b2]1101      Variables[i] = Variables[i] + subs;
[35f23d]1102    }
1103  }
[afd77b2]1104  return(Variables);
[35f23d]1105}
1106
1107///////////////////////////////////////////////////////////////////////////////
1108
1109static proc StabVarComb(intvec wt)
1110"USAGE:   StabVarComb(w); intvec w
1111PURPOSE: list all possible indices of indeterminates for a quasihom. subs.
1112RETURN:  list
1113         _[1] linear substitutions
1114         _[2] nonlinear substiutions
1115GLOBAL: 'varSubsList', contains the index j s.t. x(i) -> x(i)t(j) ...
1116"
1117{
1118  int mi, ma, i, j, k, uw, ic;
1119  list index, indices, usedWeights, combList, combinations;
1120  list linearSubs, nonlinearSubs;
1121  list partitions, subs, temp;        // subs[i] = substitution for var(i)
1122
1123  linearSubs = Table("0", "i", 1, size(wt));
1124  nonlinearSubs = Table("0", "i", 1, size(wt));
1125
1126  uw = 0;
1127  ic = 0;
1128  mi = Min(wt);
1129  ma = Max(wt);
1130
1131  for(i = mi; i <= ma; i = i + 1) {
1132    if(ContainedQ(wt, i)) {    // find variables of weight i
1133      k = 0;
1134      index = 0;
1135      // collect the indices of all variables of weight i
1136      for(j = 1; j <= size(wt); j = j + 1) {
1137        if(wt[j] == i) {
1138          k = k + 1;
1139          index[k] = j;
1140        }
1141      }
1142      uw = uw + 1;
1143      usedWeights[uw] = i;
1144      ic = ic + 1;
1145      indices[i] = index;
1146
1147      // linear part of the substitution
1148
1149      for(j = 1; j <= size(index); j = j + 1) {
1150        linearSubs[index[j]] = index;
1151      }
1152
1153      // nonlinear part of the substitution
1154
1155      if(uw > 1) {    // variables of least weight do not allow nonlinear subs.
1156
1157      partitions = Partitions(i, delete(usedWeights, uw));
1158      for(j = 1; j <= size(partitions); j = j + 1) {
1159        combinations[j] = AllCombinations(partitions[j], indices);
1160      }
1161      for(j = 1; j <= size(index); j = j + 1) {
1162        nonlinearSubs[index[j]] = FlattenQHM(combinations);   // flatten one level !
1163      }
1164
1165      } // end if
1166
1167    }
1168  }
1169  combList[1] = linearSubs;
1170  combList[2] = nonlinearSubs;
1171  return(combList);
1172}
1173
1174///////////////////////////////////////////////////////////////////////////////
1175
[70597c]1176static proc AllCombinations(list partition, list indices)
[35f23d]1177"USAGE:   AllCombinations(partition,indices); list partition, indices)
[8942a5]1178PURPOSE: all combinations for a given partititon
[35f23d]1179RETURN:  list
1180GLOBAL: varSubsList, contains the index j s.t. x(i) -> x(i)t(j) ...
1181"
1182{
1183  int i, k, m, ok, p, offset;
1184  list nrList, indexList;
1185
1186  k = 0;
1187  offset = 0;
1188  i = 1;
1189  ok = 1;
1190  m = partition[1];
1191  while(ok) {
1192    if(i > size(partition)) {
1193      ok = 0;
1194      p = 0;
1195    }
1196    else { p = partition[i];}
1197    if(p == m) { i = i + 1;}
1198    else {
1199      k = k + 1;
1200      nrList[k] = i - 1 - offset;
1201      offset = offset + i - 1;
1202      indexList[k] = indices[m];
1203      if(ok) { m = partition[i];}
1204    }
1205  }
1206  return(AllCombinationsAux(nrList, indexList));
1207}
1208
1209///////////////////////////////////////////////////////////////////////////////
1210
[70597c]1211static proc AllSingleCombinations(int n, list index)
[35f23d]1212"USAGE:   AllSingleCombinations(n index); int n, list index
1213PURPOSE: all combinations for var(n)
1214RETURN:  list
1215"
1216{
1217  int i, j, k;
1218  list comb, newC, temp, newIndex;
1219
1220  if(n == 1) {
1221    for(i = 1; i <= size(index); i = i + 1) {
1222      temp = index[i];
1223      comb[i] = temp;
1224    }
1225    return(comb);
1226  }
1227  if(size(index) == 1) {
1228    temp = Table(string(index[1]), "i", 1, n);
1229    comb[1] = temp;
1230    return(comb);
1231  }
1232  newIndex = index;
1233  for(i = 1; i <= size(index); i = i + 1) {
1234    if(i > 1) { newIndex = delete(newIndex, 1); }
1235    newC = AllSingleCombinations(n - 1, newIndex);
1236    k = size(comb);
1237    temp = 0;
1238    for(j = 1; j <= size(newC); j = j + 1) {
1239      temp[1] = index[i];
1240      temp = temp + newC[j];
1241      comb[k + j] = temp;
1242      temp = 0;
1243    }
1244  }
1245  return(comb);
1246}
1247
1248///////////////////////////////////////////////////////////////////////////////
1249
[70597c]1250static proc AllCombinationsAux(list parts, list index)
[35f23d]1251"USAGE:  AllCombinationsAux(parts ,index); list parts, index
1252PURPOSE: all compbinations for nonlinear substituiton
1253RETURN:  list
1254"
1255{
1256  int i, j, k;
1257  list comb, firstC, restC;
1258
1259  if(size(parts) == 0 || size(index) == 0) { return(comb);}
1260
1261  firstC = AllSingleCombinations(parts[1], index[1]);
1262  restC = AllCombinationsAux(delete(parts, 1), delete(index, 1));
1263
1264  if(size(restC) == 0) { comb = firstC;}
1265  else {
1266    for(i = 1; i <= size(firstC); i = i + 1) {
1267      k = size(comb);
1268      for(j = 1; j <= size(restC); j = j + 1) {
1269        //elem = firstC[i] + restC[j];
1270        // comb[k + j] = elem;
1271        comb[k + j] = firstC[i] + restC[j];
1272      }
1273    }
1274  }
1275  return(comb);
1276}
1277
1278///////////////////////////////////////////////////////////////////////////////
1279
[70597c]1280static proc Partitions(int n, list nr)
[35f23d]1281"USAGE:   Partitions(n, nr); int n, list nr
1282PURPOSE: partitions of n consisting of elements from nr
1283RETURN:  list
1284"
1285{
1286  int i, j, k;
1287  list parts, temp, restP, newP, decP;
1288
1289  if(size(nr) == 0) { return(list());}
1290  if(size(nr) == 1) {
1291    if(NumFactor(nr[1], n) > 0) {
1292      parts[1] = Table(string(nr[1]), "i", 1, NumFactor(nr[1], n));
1293    }
1294    return(parts);
1295  }
1296  else {
1297    parts =  Partitions(n, nr[1]);
1298    restP = Partitions(n, delete(nr, 1));
1299
1300    parts = parts + restP;
1301    for(i = 1; i <= n / nr[1]; i = i + 1) {
1302      temp = Table(string(nr[1]), "i", 1, i);
1303      decP = Partitions(n - i*nr[1], delete(nr, 1));
1304      k = size(parts);
1305      for(j = 1; j <= size(decP); j = j + 1) {
1306        newP = temp + decP[j];        // new partition
1307        if(!ContainedQ(parts, newP, 1)) {
1308          k = k + 1;
1309          parts[k] = newP;
1310        }
1311      }
1312    }
1313  }
1314  return(parts);
1315}
1316
1317///////////////////////////////////////////////////////////////////////////////
1318
[70597c]1319static proc NumFactor(int a, int b)
[35f23d]1320" USAGE: NumFactor(a, b); int a, b
1321PURPOSE: if b divides a then return b/a, else return 0
1322RETURN:  int
1323"
1324{
1325  int c = int(b/a);
1326  if(c*a == b) { return(c); }
1327  else {return(0)}
1328}
1329
1330///////////////////////////////////////////////////////////////////////////////
1331
[70597c]1332static proc Table(string cmd, string iterator, int lb, int ub)
[35f23d]1333" USAGE: Table(cmd,i, lb, ub); string cmd, i; int lb, ub
1334PURPOSE: generate a list of size ub - lb + 1 s.t. _[i] = cmd(i)
1335RETURN:  list
1336"
1337{
1338  list data;
1339  execute("int " + iterator + ";");
1340
1341  for(int @i = lb; @i <= ub; @i++) {
1342    execute(iterator + " = " + string(@i));
1343    execute("data[" + string(@i) + "] = " + cmd + ";");
1344  }
1345  return(data);
1346}
1347
1348///////////////////////////////////////////////////////////////////////////////
1349
1350static proc FlattenQHM(list data)
1351" USAGE: FlattenQHM(n, nr); list data
1352PURPOSE: flatten the list (one level) 'data', which is a list of lists
1353RETURN:  list
1354"
1355{
1356  int i, j, c;
1357  list fList, temp;
1358
1359  c = 1;
1360
1361  for(i = 1; i <= size(data); i = i + 1) {
1362    for(j = 1; j <= size(data[i]); j = j + 1) {
1363      fList[c] = data[i][j];
1364      c = c + 1;
1365    }
1366  }
1367  return(fList);
1368}
1369
1370///////////////////////////////////////////////////////////////////////////////
1371
1372static proc IntersectionQHM(list l1, list l2)
1373// Type : list
1374// Purpose : Intersection of l1 and l2
1375{
1376  list l;
1377  int b, c;
1378
1379  c = 1;
1380
1381  for(int i = 1; i <= size(l1); i = i + 1) {
1382    b = ContainedQ(l2, l1[i]);
1383    if(b == 1) {
1384      l[c] = l1[i];
1385      c = c + 1;
1386    }
1387  }
1388  return(l);
1389}
1390
1391///////////////////////////////////////////////////////////////////////////////
1392
1393static proc FirstEntryQHM(data, elem)
1394// Type : int
1395// Purpose : position of first entry equal to elem in data (from left to right)
1396{
1397  int i, pos;
1398
1399  i = 0;
1400  pos = 0;
1401  while(!pos && i < size(data)) {
1402    i = i + 1;
1403    if(data[i] == elem) { pos = i;}
1404  }
1405  return(pos);
1406}
1407
1408///////////////////////////////////////////////////////////////////////////////
1409
1410static proc PSum(e)
1411{
1412  poly f;
1413  for(int i = 1; i <= size(e); i = i + 1) {
1414    f = f + e[i];
1415  }
1416  return(f);
1417
1418}
1419
1420///////////////////////////////////////////////////////////////////////////////
1421
1422proc Max(data)
[fd5013]1423"USAGE:   Max(data); intvec/list of integers
[35f23d]1424PURPOSE: find the maximal integer contained in 'data'
1425RETURN:  list
1426ASSUME:  'data' contians only integers and is not empty
1427"
1428{
1429  int i;
1430  int max = data[1];
1431
1432  for(i = 1; i <= size(data); i = i + 1) {
1433    if(data[i] > max) { max = data[i]; }
1434  }
1435  return(max);
1436}
[fd5013]1437example
1438{"EXAMPLE:";  echo = 2;
1439  Max(list(1,2,3));
1440}
[35f23d]1441
1442///////////////////////////////////////////////////////////////////////////////
1443
1444proc Min(data)
[fd5013]1445"USAGE:   Min(data); intvec/list of integers
[35f23d]1446PURPOSE: find the minimal integer contained in 'data'
1447RETURN:  list
1448ASSUME:  'data' contians only integers and is not empty
1449"
1450{
1451  int i;
1452  int min = data[1];
1453
1454  for(i = 1; i <= size(data); i = i + 1) {
1455    if(data[i] < min) { min = data[i]; }
1456  }
1457  return(min);
1458}
[fd5013]1459example
1460{"EXAMPLE:";  echo = 2;
1461  Min(intvec(1,2,3));
1462}
[35f23d]1463
1464///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.