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