1 | /////////////////////////////////////////////////////////////////////////////// |
---|
2 | version="$Id: qhmoduli.lib,v 1.15 2009-04-14 08:36:17 Singular Exp $"; |
---|
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 / 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 | poly mPoly = minpoly; |
---|
133 | kill RME3; |
---|
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); |
---|
154 | kill RME4; |
---|
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 | } |
---|
165 | dbprint(dbPrt, " |
---|
166 | // 'ModEqn' created a new ring. |
---|
167 | // To see the ring, type (if the name of the ring is R): |
---|
168 | show(R); |
---|
169 | // To access the ideal of the moduli space of semiquasihomogeneous singularities |
---|
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 | } |
---|
178 | example |
---|
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 | |
---|
190 | proc QuotientEquations(ideal G, ideal Gaction, ideal embedding, list#) |
---|
191 | "USAGE: QuotientEquations(G,action,emb [, opt]); ideal G,action,emb;int opt |
---|
192 | PURPOSE: compute the quotient of the variety given by the parameterization |
---|
193 | 'emb' by the linear action 'action' of the algebraic group G. |
---|
194 | ASSUME: 'action' is linear, G must be finite if the Reynolds operator is |
---|
195 | needed (i.e., NullCone(G,action) returns some non-invariant polys) |
---|
196 | RETURN: polynomial ring over a simple extension of the ground field of the |
---|
197 | basering, containing the ideals 'id' and 'embedid'. |
---|
198 | - 'id' contains the equations of the quotient, if opt = 1; |
---|
199 | if opt = 0, 2, 'id' contains generators of the coordinate ring R |
---|
200 | of the quotient (Spec(R) is the quotient) |
---|
201 | - 'embedid' = 0, if opt = 1; |
---|
202 | if opt = 0, 2, it is the ideal defining the equivariant embedding |
---|
203 | OPTIONS: 1 compute equations of the quotient, |
---|
204 | 2 use a primary decomposition when computing the Reynolds operator,@* |
---|
205 | to combine options, add their value, default: opt =3. |
---|
206 | EXAMPLE: 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 ! |
---|
226 | //for(i = 1; i <= ncols(Gaction); i ++) { wt[i] = 1;}; |
---|
227 | for(i = 1; i <= nvars(basering); i ++) { wt[i] = 1;}; |
---|
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 | } |
---|
253 | kill RQER; |
---|
254 | export(id); |
---|
255 | export(embedid); |
---|
256 | return(RQES); |
---|
257 | } |
---|
258 | |
---|
259 | /////////////////////////////////////////////////////////////////////////////// |
---|
260 | |
---|
261 | proc UpperMonomials(poly f, list #) |
---|
262 | "USAGE: UpperMonomials(poly f, [intvec w]) |
---|
263 | PURPOSE: compute the upper monomials of the milnor algebra of f. |
---|
264 | ASSUME: f is quasihomogeneous (w.r.t. w) |
---|
265 | RETURN: ideal |
---|
266 | EXAMPLE: 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 | } |
---|
280 | example |
---|
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 | |
---|
289 | proc 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') |
---|
292 | PURPOSE: compute the induced action of the stabilizer G of f on T_, where |
---|
293 | T_ is given by the upper monomials B of the Milnor algebra of f. |
---|
294 | ASSUME: f is quasihomogeneous |
---|
295 | RETURN: polynomial ring over the same ground field, containing the ideals |
---|
296 | 'actionid' and 'stabid'. |
---|
297 | - 'actionid' is the ideal defining the induced action of Gf on T_ @* |
---|
298 | - 'stabid' is the ideal of the stabilizer Gf in the new ring |
---|
299 | EXAMPLE: 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); |
---|
387 | kill RAAR; |
---|
388 | dbprint(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, |
---|
393 | // where f is the quasihomogeneous principal part, type |
---|
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 | } |
---|
400 | example |
---|
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 | |
---|
412 | proc StabOrder(list #) |
---|
413 | "USAGE: StabOrder(f); poly f |
---|
414 | PURPOSE: compute the order of the stabilizer group of f. |
---|
415 | ASSUME: f quasihomogeneous polynomial with an isolated singularity at 0 |
---|
416 | RETURN: int |
---|
417 | GLOBAL: 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 | |
---|
432 | proc StabEqn(poly f) |
---|
433 | "USAGE: StabEqn(f); f polynomial |
---|
434 | PURPOSE: compute the equations of the isometry group of f. |
---|
435 | ASSUME: f semiquasihomogeneous polynomial with an isolated singularity at 0 |
---|
436 | RETURN: list of two rings 'S1', 'S2' |
---|
437 | - 'S1' contians the equations of the stabilizer (ideal 'stabid') @* |
---|
438 | - 'S2' contains the action of the stabilizer (ideal 'actionid') |
---|
439 | EXAMPLE: example StabEqn; shows an example |
---|
440 | GLOBAL: varSubsList, contains the index j s.t. x(i) -> x(i)t(j) ... |
---|
441 | " |
---|
442 | { |
---|
443 | dbprint(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 | } |
---|
459 | example |
---|
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 | |
---|
470 | proc StabEqnId(ideal data, intvec wt) |
---|
471 | "USAGE: StabEqn(I, w); I ideal, w intvec |
---|
472 | PURPOSE: compute the equations of the isometry group of the ideal I, |
---|
473 | each generator of I is fixed by the stabilizer. |
---|
474 | ASSUME: I semiquasihomogeneous ideal wrt 'w' with an isolated singularity at 0 |
---|
475 | RETURN: list of two rings 'S1', 'S2' |
---|
476 | - 'S1' contians the equations of the stabilizer (ideal 'stabid') @* |
---|
477 | - 'S2' contains the action of the stabilizer (ideal 'actionid') |
---|
478 | EXAMPLE: example StabEqnId; shows an example |
---|
479 | GLOBAL: varSubsList, contains the index j s.t. t(i) -> t(i)t(j) ... |
---|
480 | " |
---|
481 | { |
---|
482 | int i, j, c, k, r, nrVars, offset, n, sln, dbPrt; |
---|
483 | list Variables, rd, temp, sList, varSubsList; |
---|
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); |
---|
492 | Variables = StabVar(wt); // possible quasihomogeneous substitutions |
---|
493 | nrVars = 0; |
---|
494 | for(i = 1; i <= size(wt); i++) |
---|
495 | { |
---|
496 | nrVars = nrVars + size(Variables[i]); |
---|
497 | } |
---|
498 | |
---|
499 | // set the new basering needed for the substitutions |
---|
500 | |
---|
501 | varString = "s(1.." + string(nrVars) + ")"; |
---|
502 | if(npars(basering) == 1) |
---|
503 | { |
---|
504 | parString = "(0, " + parstr(basering) + ")"; |
---|
505 | } |
---|
506 | else { parString = "0"; } |
---|
507 | |
---|
508 | def RSTB = basering; |
---|
509 | mPoly = minpoly; |
---|
510 | ringSTR = "ring RSTR = " + parString + ", (" + varstr(basering) + ", " + varString + "), dp;"; // dp |
---|
511 | ringSTR1 = "ring RSTT = " + parString + ", (" + varString + ", " + varstr(basering) + "), dp;"; |
---|
512 | |
---|
513 | if(defined(RSTR)) { kill RSTR;} |
---|
514 | if(defined(RSTT)) { kill RSTT;} |
---|
515 | execute(ringSTR1); // this ring is only used for the result, where the variables |
---|
516 | export(RSTT); // are s(1..m),t(1..n), as needed for Derksens algorithm (NullCone) |
---|
517 | minpoly = number(imap(RSTB, mPoly)); |
---|
518 | |
---|
519 | execute(ringSTR); |
---|
520 | export(RSTR); |
---|
521 | minpoly = number(imap(RSTB, mPoly)); |
---|
522 | poly f, f1, g, h, vars, pp; // f1 is the polynomial after subs, |
---|
523 | ideal allEqns, qhsubs, actionid, stabid, J; |
---|
524 | list ringList; // all t(i)`s which do not appear in f1 |
---|
525 | ideal data = simplify(imap(RSTB, data), 2); |
---|
526 | |
---|
527 | // generate the quasihomogeneous substitution map F |
---|
528 | |
---|
529 | nrVars = 0; |
---|
530 | offset = 0; |
---|
531 | for(i = 1; i <= size(wt); i++) |
---|
532 | { // build the substitution t(i) -> ... |
---|
533 | if(i > 1) { offset = offset + size(Variables[i - 1]); } |
---|
534 | g = 0; |
---|
535 | for(j = 1; j <= size(Variables[i]); j++) |
---|
536 | { |
---|
537 | pp = 1; |
---|
538 | for(k = 2; k <= size(Variables[i][j]); k++) |
---|
539 | { |
---|
540 | pp = pp * var(Variables[i][j][k]); |
---|
541 | if(Variables[i][j][k] == i) { varSubsList[i] = offset + j;} |
---|
542 | } |
---|
543 | g = g + s(offset + j) * pp; |
---|
544 | } |
---|
545 | qhsubs[i] = g; |
---|
546 | } |
---|
547 | dbprint(dbPrt, " qhasihomogenous substituion =" + string(qhsubs)); |
---|
548 | map F = RSTR, qhsubs; |
---|
549 | kill varSubsList; |
---|
550 | |
---|
551 | // get the equations of the stabilizer by comparing coefficients |
---|
552 | // in the equation f = F(f). |
---|
553 | |
---|
554 | vars = RingVarProduct(Table("i", "i", 1, size(wt))); |
---|
555 | |
---|
556 | allEqns = 0; |
---|
557 | |
---|
558 | matrix newcoMx, coMx; |
---|
559 | int d; |
---|
560 | for(r = 1; r <= ncols(data); r++) |
---|
561 | { |
---|
562 | |
---|
563 | f = data[r]; |
---|
564 | f1 = F(f); |
---|
565 | d = deg(f); |
---|
566 | newcoMx = coef(f1, vars); // coefficients of F(f) |
---|
567 | coMx = coef(f, vars); // coefficients of f |
---|
568 | |
---|
569 | for(i = 1; i <= ncols(newcoMx); i++) |
---|
570 | { // build the system of eqns via coeff. comp. |
---|
571 | j = 1; |
---|
572 | h = 0; |
---|
573 | while(j <= ncols(coMx)) |
---|
574 | { // all monomials in f |
---|
575 | if(coMx[j][1] == newcoMx[i][1]) { h = coMx[j][2]; j = ncols(coMx) + 1;} |
---|
576 | else {j = j + 1;} |
---|
577 | } |
---|
578 | J = J, newcoMx[i][2] - h; // add equation |
---|
579 | } |
---|
580 | allEqns = allEqns, J; |
---|
581 | |
---|
582 | } |
---|
583 | allEqns = std(allEqns); |
---|
584 | |
---|
585 | // simplify the equations, i.e., if s(i) in J then remove s(i) from J |
---|
586 | // and from the basering |
---|
587 | |
---|
588 | sList = SimplifyIdeal(allEqns, n, "s"); |
---|
589 | stabid = sList[1]; |
---|
590 | map phi = basering, sList[2]; |
---|
591 | sln = size(sList[3]) - n; |
---|
592 | |
---|
593 | // change the substitution |
---|
594 | |
---|
595 | actionid = phi(qhsubs); |
---|
596 | |
---|
597 | // change to new ring, auxillary construction |
---|
598 | |
---|
599 | setring(RSTT); |
---|
600 | ideal actionid, stabid; |
---|
601 | |
---|
602 | actionid = imap(RSTR, actionid); |
---|
603 | stabid = imap(RSTR, stabid); |
---|
604 | export(stabid); |
---|
605 | export(actionid); |
---|
606 | ringList[2] = RSTT; |
---|
607 | |
---|
608 | dbprint(dbPrt, " substitution = " + string(actionid)); |
---|
609 | dbprint(dbPrt, " equations of stabilizer = " + string(stabid)); |
---|
610 | |
---|
611 | varString = "s(1.." + string(sln) + ")"; |
---|
612 | ringSTR = "ring RSTS = " + parString + ", (" + varString + "), dp;"; |
---|
613 | execute(ringSTR); |
---|
614 | minpoly = number(imap(RSTB, mPoly)); |
---|
615 | ideal stabid = std(imap(RSTR, stabid)); |
---|
616 | export(stabid); |
---|
617 | ringList[1] = RSTS; |
---|
618 | dbprint(dbPrt, " |
---|
619 | // 'StabEqnId' created a list of 2 rings. |
---|
620 | // To see the rings, type (if the name of your list is stab): |
---|
621 | show(stab); |
---|
622 | // To access the 1-st ring and map (and similair for the others), type: |
---|
623 | def S1 = stab[1]; setring S1; stabid; |
---|
624 | // S1/stabid is the coordinate ring of the variety of the |
---|
625 | // stabilizer, say G. If G x K^n --> K^n is the action of G on |
---|
626 | // K^n, then the ideal 'actionid' in the second ring describes |
---|
627 | // the dual map on the ring level. |
---|
628 | // To access the 2-nd ring and map (and similair for the others), type: |
---|
629 | def S2 = stab[2]; setring S2; actionid; |
---|
630 | "); |
---|
631 | return(ringList); |
---|
632 | } |
---|
633 | example |
---|
634 | {"EXAMPLE:"; echo = 2; |
---|
635 | ring B = 0,(x,y,z), ls; |
---|
636 | ideal I = x2,y3,z6; |
---|
637 | intvec w = 3,2,1; |
---|
638 | list stab = StabEqnId(I, w); |
---|
639 | def S1 = stab[1]; setring S1; stabid; |
---|
640 | def S2 = stab[2]; setring S2; actionid; |
---|
641 | } |
---|
642 | |
---|
643 | /////////////////////////////////////////////////////////////////////////////// |
---|
644 | static |
---|
645 | proc ArnoldFormMain(poly f, B, poly Fs, ideal reduceIdeal, int nrs, int nrt) |
---|
646 | "USAGE: ArnoldFormMain(f, B, Fs, rI, nrs, nrt); |
---|
647 | poly f,Fs; ideal B, rI; int nrs, nrt |
---|
648 | PURPOSE: compute the induced action of 'G_f' on T_, where f is the principal |
---|
649 | part and 'Fs' is the semiuniversal unfolding of 'f' with x_i |
---|
650 | substituted by actionid[i], 'B' is a list of upper basis monomials |
---|
651 | for the milnor algebra of 'f', 'nrs' = number of variables for 'G_f' |
---|
652 | and 'nrt' = dimension of T_ |
---|
653 | ASSUME: f is quasihomogeneous with an isolated singularity at 0, |
---|
654 | s(1..r), t(1..m) are parameters of the basering |
---|
655 | RETURN: poly |
---|
656 | EXAMPLE: example ArnoldAction; shows an example |
---|
657 | " |
---|
658 | { |
---|
659 | int i, j, d, ub, dbPrt; |
---|
660 | list upperBasis, basisDegList, gmonos, common, parts; |
---|
661 | ideal jacobianId, jacobIdstd, mapId; // needed for phi |
---|
662 | intvec wt = weight(f); |
---|
663 | matrix gCoeffMx; // for lift command |
---|
664 | poly newFs, g, gred, tt; // g = sum of all monomials of degree d, gred is needed for lift |
---|
665 | map phi; // the map from Arnold's Theorem |
---|
666 | |
---|
667 | dbPrt = printlevel-voice+2; |
---|
668 | jacobianId = jacob(f); |
---|
669 | jacobIdstd = std(jacobianId); |
---|
670 | newFs = Fs; |
---|
671 | for(i = 1; i <= size(B); i++) |
---|
672 | { |
---|
673 | basisDegList[i] = deg(B[i], wt); |
---|
674 | } |
---|
675 | ub = Max(basisDegList) + 1; // max degree of an upper monomial |
---|
676 | |
---|
677 | parts = MonosAndTerms(newFs - f, wt, ub); |
---|
678 | gmonos = parts[1]; |
---|
679 | d = deg(f, wt); |
---|
680 | |
---|
681 | for(i = d + 1; i < ub; i++) |
---|
682 | { // base[1] = monomials of degree i |
---|
683 | upperBasis[i] = SelectMonos(list(B, B), wt, i); // B must not contain 0's |
---|
684 | } |
---|
685 | |
---|
686 | // test if each monomial of Fs is contained in B, if not, |
---|
687 | // compute a substitution via Arnold's theorem and substitutite |
---|
688 | // it into newFs |
---|
689 | |
---|
690 | for(i = d + 1; i < ub; i = i + 1) |
---|
691 | { // ub instead of @UB |
---|
692 | dbprint(dbPrt, "-- degree = " + string(i) + " of " + string(ub - 1) + " ---------------------------"); |
---|
693 | if(size(newFs) < 80) { dbprint(dbPrt, " polynomial = " + string(newFs - f));} |
---|
694 | else { dbprint(dbPrt, " poly has deg (not weighted) " + string(deg(newFs)) + " and contains " + string(size(newFs)) + " monos");} |
---|
695 | |
---|
696 | // select monomials of degree i and intersect them with upperBasis[i] |
---|
697 | |
---|
698 | gmonos = SelectMonos(parts, wt, i); |
---|
699 | common = IntersectionQHM(upperBasis[i][1], gmonos[1]); |
---|
700 | if(size(common) == size(gmonos[1])) |
---|
701 | { |
---|
702 | dbprint(dbPrt, " no additional monomials "); |
---|
703 | } |
---|
704 | |
---|
705 | // other monomials than those in upperBasis occur, compute |
---|
706 | // the map constructed in the proof of Arnold's theorem |
---|
707 | // write g = c[i] * jacobianId[i] |
---|
708 | |
---|
709 | else |
---|
710 | { |
---|
711 | dbprint(dbPrt, " additional Monomials found, compute the map "); |
---|
712 | g = PSum(gmonos[2]); // sum of all monomials in g of degree i |
---|
713 | dbprint(dbPrt, " sum of degree " + string(i) + " is " + string(g)); |
---|
714 | |
---|
715 | gred = reduce(g, jacobIdstd); |
---|
716 | gCoeffMx = lift(jacobianId, g - gred); // compute c[i] |
---|
717 | mapId = var(1) - gCoeffMx[1][1]; // generate the map |
---|
718 | for(j = 2; j <= size(gCoeffMx); j++) |
---|
719 | { |
---|
720 | mapId[j] = var(j) - gCoeffMx[1][j]; |
---|
721 | } |
---|
722 | dbprint(dbPrt, " map = " + string(mapId)); |
---|
723 | // apply the map to newFs |
---|
724 | newFs = APSubstitution(newFs, mapId, reduceIdeal, wt, ub, nrs, nrt); |
---|
725 | parts = MonosAndTerms(newFs - f, wt, ub); // monos and terms of deg < ub |
---|
726 | newFs = PSum(parts[2]) + f; // result of APS... is already reduced |
---|
727 | dbprint(dbPrt, " monomials of degree " + string(i)); |
---|
728 | } |
---|
729 | } |
---|
730 | return(newFs); |
---|
731 | } |
---|
732 | |
---|
733 | /////////////////////////////////////////////////////////////////////////////// |
---|
734 | |
---|
735 | static proc MonosAndTerms(poly f, wt, int ub) |
---|
736 | "USAGE: MonosAndTerms(f, w, ub); poly f, intvec w, int ub |
---|
737 | PURPOSE: returns a list of all monomials and terms occuring in f of |
---|
738 | weighted degree < ub |
---|
739 | RETURN: list |
---|
740 | _[1] list of monomials |
---|
741 | _[2] list of terms |
---|
742 | EXAMPLE: example MonosAndTerms shows an example |
---|
743 | " |
---|
744 | { |
---|
745 | int i, k; |
---|
746 | list monomials, terms; |
---|
747 | poly mono, lcInv, data; |
---|
748 | |
---|
749 | data = jet(f, ub - 1, wt); |
---|
750 | k = 0; |
---|
751 | for(i = 1; i <= size(data); i++) |
---|
752 | { |
---|
753 | mono = lead(data[i]); |
---|
754 | if(deg(mono, wt) < ub) |
---|
755 | { |
---|
756 | k = k + 1; |
---|
757 | lcInv = 1/leadcoef(mono); |
---|
758 | monomials[k] = mono * lcInv; |
---|
759 | terms[k] = mono; |
---|
760 | } |
---|
761 | } |
---|
762 | return(list(monomials, terms)); |
---|
763 | } |
---|
764 | example |
---|
765 | {"EXAMPLE:"; echo = 2; |
---|
766 | ring B = 0,(x,y,z), lp; |
---|
767 | poly f = 6*x2 + 2*x3 + 9*x*y2 + z*y + x*z6; |
---|
768 | MonosAndTerms(f, intvec(2,1,1), 5); |
---|
769 | } |
---|
770 | |
---|
771 | /////////////////////////////////////////////////////////////////////////////// |
---|
772 | |
---|
773 | static proc SelectMonos(parts, intvec wt, int d) |
---|
774 | "USAGE: SelectMonos(parts, w, d); list/ideal parts, intvec w, int d |
---|
775 | PURPOSE: returns a list of all monomials and terms occuring in f of |
---|
776 | weighted degree = d |
---|
777 | RETURN: list |
---|
778 | _[1] list of monomials |
---|
779 | _[2] list of terms |
---|
780 | EXAMPLE: example SelectMonos; shows an example |
---|
781 | " |
---|
782 | { |
---|
783 | int i, k; |
---|
784 | list monomials, terms; |
---|
785 | poly mono; |
---|
786 | |
---|
787 | k = 0; |
---|
788 | for(i = 1; i <= size(parts[1]); i++) |
---|
789 | { |
---|
790 | mono = parts[1][i]; |
---|
791 | if(deg(mono, wt) == d) |
---|
792 | { |
---|
793 | k++; |
---|
794 | monomials[k] = mono; |
---|
795 | terms[k] = parts[2][i]; |
---|
796 | } |
---|
797 | } |
---|
798 | return(list(monomials, terms)); |
---|
799 | } |
---|
800 | example |
---|
801 | {"EXAMPLE:"; echo = 2; |
---|
802 | ring B = 0,(x,y,z), lp; |
---|
803 | poly f = 6*x2 + 2*x3 + 9*x*y2 + z*y + x*z6; |
---|
804 | list mt = MonosAndTerms(f, intvec(2,1,1), 5); |
---|
805 | SelectMonos(mt, intvec(2,1,1), 4); |
---|
806 | } |
---|
807 | |
---|
808 | /////////////////////////////////////////////////////////////////////////////// |
---|
809 | |
---|
810 | static proc Expand(substitution, degVec, ideal reduceI, intvec w1, int ub, list truncated) |
---|
811 | "USAGE: Expand(substitution, degVec, reduceI, w, ub, truncated); |
---|
812 | ideal/list substitution, list/intvec degVec, ideal reduceI, intvec w, |
---|
813 | int ub, list truncated |
---|
814 | PURPOSE: substitute 'substitution' in the monomial given by the list of |
---|
815 | exponents 'degVec', omit all terms of weighted degree > ub and reduce |
---|
816 | the result w.r.t. 'reduceI'. If truncated[i] = 0 then the result is |
---|
817 | stored for later use. |
---|
818 | RETURN: poly |
---|
819 | NOTE: used by APSubstitution |
---|
820 | GLOBAL: computedPowers |
---|
821 | " |
---|
822 | { |
---|
823 | int i, minDeg; |
---|
824 | list powerList; |
---|
825 | poly g, h; |
---|
826 | |
---|
827 | // compute substitution[1]^degVec[1],...,subs[n]^degVec[n] |
---|
828 | |
---|
829 | for(i = 1; i <= ncols(substitution); i++) |
---|
830 | { |
---|
831 | if(size(substitution[i]) < 3 || degVec[i] < 4) |
---|
832 | { |
---|
833 | powerList[i] = reduce(substitution[i]^degVec[i], reduceI); // new |
---|
834 | } // directly for small exponents |
---|
835 | else |
---|
836 | { |
---|
837 | powerList[i] = PolyPower1(i, substitution[i], degVec[i], reduceI, w1, truncated[i], ub); |
---|
838 | } |
---|
839 | } |
---|
840 | // multiply the terms obtained by using PolyProduct(); |
---|
841 | g = powerList[1]; |
---|
842 | minDeg = w1[1] * degVec[1]; |
---|
843 | for(i = 2; i <= ncols(substitution); i++) |
---|
844 | { |
---|
845 | g = jet(g, ub - w1[i] * degVec[i] - 1, w1); |
---|
846 | h = jet(powerList[i], ub - minDeg - 1, w1); |
---|
847 | g = PolyProduct(g, h, reduceI, w1, ub); |
---|
848 | if(g == 0) { Print(" g = 0 "); break;} |
---|
849 | minDeg = minDeg + w1[i] * degVec[i]; |
---|
850 | } |
---|
851 | return(g); |
---|
852 | } |
---|
853 | |
---|
854 | /////////////////////////////////////////////////////////////////////////////// |
---|
855 | |
---|
856 | static proc PolyProduct(poly g1, poly h1, ideal reduceI, intvec wt, int ub) |
---|
857 | "USAGE: PolyProduct(g, h, reduceI, wt, ub); poly g, h; ideal reduceI, |
---|
858 | intvec wt, int ub. |
---|
859 | PURPOSE: compute g*h and reduce it w.r.t 'reduceI' and omit terms of weighted |
---|
860 | degree > ub. |
---|
861 | RETURN: poly |
---|
862 | NOTE: used by 'Expand' |
---|
863 | " |
---|
864 | { |
---|
865 | int SUBSMAXSIZE = 3000; |
---|
866 | int i, nrParts, sizeOfPart, currentPos, partSize, maxSIZE; |
---|
867 | poly g, h, gxh, prodComp, @g2; // replace @g2 by subst. |
---|
868 | |
---|
869 | g = g1; |
---|
870 | h = h1; |
---|
871 | |
---|
872 | if(size(g)*size(h) > SUBSMAXSIZE) |
---|
873 | { |
---|
874 | // divide the polynomials with more terms in parts s.t. |
---|
875 | // the product of each part with the other polynomial |
---|
876 | // has at most SUBMAXSIZE terms |
---|
877 | |
---|
878 | if(size(g) < size(h)) { poly @h = h; h = g; g = @h;@h = 0; } |
---|
879 | maxSIZE = SUBSMAXSIZE / size(h); |
---|
880 | //print(" SUBSMAXSIZE = "+string(SUBSMAXSIZE)+" exceeded by "+string(size(g)*size(h)) + ", maxSIZE = ", string(maxSIZE)); |
---|
881 | nrParts = size(g) / maxSIZE + 1; |
---|
882 | partSize = size(g) / nrParts; |
---|
883 | gxh = 0; // 'g times h' |
---|
884 | for(i = 1; i <= nrParts; i++) |
---|
885 | { |
---|
886 | //print(" loop #" + string(i) + " of " + string(nrParts)); |
---|
887 | currentPos = (i - 1) * partSize; |
---|
888 | if(i < nrParts) {sizeOfPart = partSize;} |
---|
889 | else { sizeOfPart = size(g) - (nrParts - 1) * partSize; print(" last #" + string(sizeOfPart) + " terms ");} |
---|
890 | prodComp = g[currentPos + 1..sizeOfPart + currentPos] * h; // multiply a part |
---|
891 | @g2 = jet(prodComp, ub - 1, wt); // eventual reduce ... |
---|
892 | if(size(@g2) < size(prodComp)) { print(" killed " + string(size(prodComp) - size(@g2)) + " terms ");} |
---|
893 | gxh = reduce(gxh + @g2, reduceI); |
---|
894 | } |
---|
895 | } |
---|
896 | else |
---|
897 | { |
---|
898 | gxh = reduce(jet(g * h,ub - 1, wt), reduceI); |
---|
899 | } // compute directly |
---|
900 | return(gxh); |
---|
901 | } |
---|
902 | |
---|
903 | /////////////////////////////////////////////////////////////////////////////// |
---|
904 | |
---|
905 | static proc PolyPower1(int varIndex, poly f, int e, ideal reduceI, intvec wt, |
---|
906 | int truncated, int ub) |
---|
907 | "USAGE: PolyPower1(i, f, e, reduceI, wt, truncated, ub);int i, e, ub;poly f; |
---|
908 | ideal reduceI; intvec wt; list truncated; |
---|
909 | PURPOSE: compute f^e, use previous computations if possible, and reduce it |
---|
910 | w.r.t reudecI and omit terms of weighted degree > ub. |
---|
911 | RETURN: poly |
---|
912 | NOTE: used by 'Expand' |
---|
913 | GLOBAL: 'computedPowers' |
---|
914 | " |
---|
915 | { |
---|
916 | int i, ordOfg, lb, maxPrecomputedPower; |
---|
917 | poly g, fn; |
---|
918 | |
---|
919 | if(e == 0) { return(1);} |
---|
920 | if(e == 1) { return(f);} |
---|
921 | if(f == 0) { return(1); } |
---|
922 | else |
---|
923 | { |
---|
924 | // test if f has been computed to some power |
---|
925 | if(computedPowers[varIndex][1] > 0) |
---|
926 | { |
---|
927 | maxPrecomputedPower = computedPowers[varIndex][1]; |
---|
928 | if(maxPrecomputedPower >= e) |
---|
929 | { |
---|
930 | // no computation necessary, f^e has already benn computed |
---|
931 | g = computedPowers[varIndex][2][e - 1]; |
---|
932 | //Print("No computation, from list : g = elem [", varIndex, ", 2, ", e - 1, "]"); |
---|
933 | lb = e + 1; |
---|
934 | } |
---|
935 | else { // f^d computed, where d < e |
---|
936 | g = computedPowers[varIndex][2][maxPrecomputedPower - 1]; |
---|
937 | ordOfg = maxPrecomputedPower * wt[varIndex]; |
---|
938 | lb = maxPrecomputedPower + 1; |
---|
939 | } |
---|
940 | } |
---|
941 | else |
---|
942 | { // no precomputed data |
---|
943 | lb = 2; |
---|
944 | ordOfg = wt[varIndex]; |
---|
945 | g = f; |
---|
946 | } |
---|
947 | for(i = lb; i <= e; i++) |
---|
948 | { |
---|
949 | fn = jet(f, ub - ordOfg - 1, wt); // reduce w.r.t. reduceI |
---|
950 | g = PolyProduct(g, fn, reduceI, wt, ub); |
---|
951 | ordOfg = ordOfg + wt[varIndex]; |
---|
952 | if(g == 0) { break; } |
---|
953 | if((i > maxPrecomputedPower) && !truncated) |
---|
954 | { |
---|
955 | if(maxPrecomputedPower == 0) |
---|
956 | { // init computedPowers |
---|
957 | computedPowers[varIndex] = list(i, list(g)); |
---|
958 | } |
---|
959 | computedPowers[varIndex][1] = i; // new degree |
---|
960 | computedPowers[varIndex][2][i - 1] = g; |
---|
961 | maxPrecomputedPower = i; |
---|
962 | } |
---|
963 | } |
---|
964 | } |
---|
965 | return(g); |
---|
966 | } |
---|
967 | |
---|
968 | /////////////////////////////////////////////////////////////////////////////// |
---|
969 | |
---|
970 | static proc RingVarsToList(list @index) |
---|
971 | { |
---|
972 | int i; |
---|
973 | list temp; |
---|
974 | |
---|
975 | for(i = 1; i <= size(@index); i++) { temp[i] = string(var(@index[i])); } |
---|
976 | return(temp); |
---|
977 | } |
---|
978 | |
---|
979 | /////////////////////////////////////////////////////////////////////////////// |
---|
980 | static |
---|
981 | proc APSubstitution(poly f, ideal substitution, ideal reduceIdeal, intvec wt, int ub, int nrs, int nrt) |
---|
982 | "USAGE: APSubstitution(f, subs, reduceI, w, ub, int nrs, int nrt); poly f |
---|
983 | ideal subs, reduceI, intvec w, int ub, nrs, nrt; |
---|
984 | nrs = number of parameters s(1..nrs), |
---|
985 | nrt = number of parameters t(1..nrt) |
---|
986 | PURPOSE: substitute 'subs' in f, omit all terms with weighted degree > ub and |
---|
987 | reduce the result w.r.t. 'reduceI'. |
---|
988 | RETURN: poly |
---|
989 | GLOBAL: 'computedPowers' |
---|
990 | " |
---|
991 | { |
---|
992 | int i, j, k, d, offset; |
---|
993 | int n = nvars(basering); |
---|
994 | list coeffList, parts, degVecList, degOfMonos; |
---|
995 | list computedPowers, truncatedQ, degOfSubs, @temp; |
---|
996 | string ringSTR, @ringVars; |
---|
997 | |
---|
998 | export(computedPowers); |
---|
999 | |
---|
1000 | // store arguments in strings |
---|
1001 | |
---|
1002 | def RASB = basering; |
---|
1003 | |
---|
1004 | parts = MonosAndTerms(f, wt, ub); |
---|
1005 | for(i = 1; i <= size(parts[1]); i = i + 1) |
---|
1006 | { |
---|
1007 | coeffList[i] = parts[2][i]/parts[1][i]; |
---|
1008 | degVecList[i] = leadexp(parts[1][i]); |
---|
1009 | degOfMonos[i] = deg(parts[1][i], wt); |
---|
1010 | } |
---|
1011 | |
---|
1012 | // built new basering with no parameters and order dp ! |
---|
1013 | // the parameters of the basering are appended to |
---|
1014 | // the variables of the basering ! |
---|
1015 | // set ideal mpoly = minpoly for reduction ! |
---|
1016 | |
---|
1017 | @ringVars = "(" + varstr(basering) + ", " + parstr(1) + ","; // precondition |
---|
1018 | if(nrs > 0) |
---|
1019 | { |
---|
1020 | @ringVars = @ringVars + "s(1.." + string(nrs) + "), "; |
---|
1021 | } |
---|
1022 | @ringVars = @ringVars + "t(1.." + string(nrt) + "))"; |
---|
1023 | ringSTR = "ring RASR = 0, " + @ringVars + ", dp;"; // new basering |
---|
1024 | |
---|
1025 | // built the "reduction" ring with the reduction ideal |
---|
1026 | |
---|
1027 | execute(ringSTR); |
---|
1028 | export(RASR); |
---|
1029 | ideal reduceIdeal, substitution, newSubs; |
---|
1030 | intvec w1, degVec; |
---|
1031 | list minDeg, coeffList, degList; |
---|
1032 | poly f, g, h, subsPoly; |
---|
1033 | |
---|
1034 | w1 = wt; // new weights |
---|
1035 | offset = nrs + nrt + 1; |
---|
1036 | for(i = n + 1; i <= offset + n; i = i + 1) { w1[i] = 0; } |
---|
1037 | |
---|
1038 | reduceIdeal = std(imap(RASB, reduceIdeal)); // omit later ! |
---|
1039 | coeffList = imap(RASB, coeffList); |
---|
1040 | substitution = imap(RASB, substitution); |
---|
1041 | |
---|
1042 | f = imap(RASB, f); |
---|
1043 | |
---|
1044 | for(i = 1; i <= n; i++) |
---|
1045 | { // all "base" variables |
---|
1046 | computedPowers[i] = list(0); |
---|
1047 | for(j = 1; j <= size(substitution[i]); j++) { degList[j] = deg(substitution[i][j], w1);} |
---|
1048 | degOfSubs[i] = degList; |
---|
1049 | } |
---|
1050 | |
---|
1051 | // substitute in each monomial seperately |
---|
1052 | |
---|
1053 | g = 0; |
---|
1054 | for(i = 1; i <= size(degVecList); i++) |
---|
1055 | { |
---|
1056 | truncatedQ = Table("0", "i", 1, n); |
---|
1057 | newSubs = 0; |
---|
1058 | degVec = degVecList[i]; |
---|
1059 | d = degOfMonos[i]; |
---|
1060 | |
---|
1061 | // check if some terms in the substitution can be omitted |
---|
1062 | // degVec = list of exponents of the monomial m |
---|
1063 | // minDeg[j] denotes the weighted degree of the monomial m' |
---|
1064 | // where m' is the monomial m without the j-th variable |
---|
1065 | |
---|
1066 | for(j = 1; j <= size(degVec); j++) { minDeg[j] = d - degVec[j] * wt[j]; } |
---|
1067 | |
---|
1068 | for(j = 1; j <= size(degVec); j++) |
---|
1069 | { |
---|
1070 | subsPoly = 0; // set substitution to 0 |
---|
1071 | if(degVec[j] > 0) |
---|
1072 | { |
---|
1073 | // if variable occurs then check if |
---|
1074 | // substitution[j][k] * (linear part)^(degVec[j]-1) + minDeg[j] < ub |
---|
1075 | // i.e. look for the smallest possible combination in subs[j]^degVec[j] |
---|
1076 | // which comes from the term substitution[j][k]. This term is multiplied |
---|
1077 | // with the rest of the monomial, which has at least degree minDeg[j]. |
---|
1078 | // If the degree of this product is < ub then substitution[j][k] contributes |
---|
1079 | // to the result and cannot be omitted |
---|
1080 | |
---|
1081 | for(k = 1; k <= size(substitution[j]); k++) |
---|
1082 | { |
---|
1083 | if(degOfSubs[j][k] + (degVec[j] - 1) * wt[j] + minDeg[j] < ub) |
---|
1084 | { |
---|
1085 | subsPoly = subsPoly + substitution[j][k]; |
---|
1086 | } |
---|
1087 | } |
---|
1088 | } |
---|
1089 | newSubs[j] = subsPoly; // set substitution |
---|
1090 | if(substitution[j] - subsPoly != 0) { truncatedQ[j] = 1;} // mark that substitution[j] is truncated |
---|
1091 | } |
---|
1092 | h = Expand(newSubs, degVec, reduceIdeal, w1, ub, truncatedQ) * coeffList[i]; // already reduced |
---|
1093 | g = reduce(g + h, reduceIdeal); |
---|
1094 | } |
---|
1095 | kill computedPowers; |
---|
1096 | setring RASB; |
---|
1097 | poly fnew = imap(RASR, g); |
---|
1098 | kill RASR; |
---|
1099 | return(fnew); |
---|
1100 | } |
---|
1101 | |
---|
1102 | /////////////////////////////////////////////////////////////////////////////// |
---|
1103 | |
---|
1104 | static proc StabVar(intvec wt) |
---|
1105 | "USAGE: StabVar(w); intvec w |
---|
1106 | PURPOSE: compute the indicies for quasihomogeneous substitutions of each |
---|
1107 | variable. |
---|
1108 | ASSUME: f semiquasihomogeneous polynomial with an isolated singularity at 0 |
---|
1109 | RETURN: list |
---|
1110 | _[i] list of combinations for var(i) (i must be appended |
---|
1111 | to each comb) |
---|
1112 | GLOBAL: 'varSubsList', contains the index j s.t. x(i) -> x(i)t(j) ... |
---|
1113 | " |
---|
1114 | { |
---|
1115 | int i, j, k, uw, ic; |
---|
1116 | list varList, Variables, subs; |
---|
1117 | string str, varString; |
---|
1118 | |
---|
1119 | varList = StabVarComb(wt); |
---|
1120 | for(i = 1; i <= size(wt); i = i + 1) |
---|
1121 | { |
---|
1122 | subs = 0; |
---|
1123 | // built linear substituitons |
---|
1124 | for(j = 1; j <= size(varList[1][i]); j++) |
---|
1125 | { |
---|
1126 | subs[j] = list(i) + list(varList[1][i][j]); |
---|
1127 | } |
---|
1128 | Variables[i] = subs; |
---|
1129 | if(size(varList[2][i]) > 0) |
---|
1130 | { |
---|
1131 | // built nonlinear substituitons |
---|
1132 | subs = 0; |
---|
1133 | for(j = 1; j <= size(varList[2][i]); j++) |
---|
1134 | { |
---|
1135 | subs[j] = list(i) + varList[2][i][j]; |
---|
1136 | } |
---|
1137 | Variables[i] = Variables[i] + subs; |
---|
1138 | } |
---|
1139 | } |
---|
1140 | return(Variables); |
---|
1141 | } |
---|
1142 | |
---|
1143 | /////////////////////////////////////////////////////////////////////////////// |
---|
1144 | |
---|
1145 | static proc StabVarComb(intvec wt) |
---|
1146 | "USAGE: StabVarComb(w); intvec w |
---|
1147 | PURPOSE: list all possible indices of indeterminates for a quasihom. subs. |
---|
1148 | RETURN: list |
---|
1149 | _[1] linear substitutions |
---|
1150 | _[2] nonlinear substiutions |
---|
1151 | GLOBAL: 'varSubsList', contains the index j s.t. x(i) -> x(i)t(j) ... |
---|
1152 | " |
---|
1153 | { |
---|
1154 | int mmi, mma, ii, j, k, uw, ic; |
---|
1155 | list index, indices, usedWeights, combList, combinations; |
---|
1156 | list linearSubs, nonlinearSubs; |
---|
1157 | list partitions, subs, temp; // subs[i] = substitution for var(i) |
---|
1158 | |
---|
1159 | linearSubs = Table("0", "i", 1, size(wt)); |
---|
1160 | nonlinearSubs = Table("0", "i", 1, size(wt)); |
---|
1161 | |
---|
1162 | uw = 0; |
---|
1163 | ic = 0; |
---|
1164 | mmi = Min(wt); |
---|
1165 | mma = Max(wt); |
---|
1166 | |
---|
1167 | for(ii = mmi; ii <= mma; ii++) |
---|
1168 | { |
---|
1169 | if(ContainedQ(wt, ii)) |
---|
1170 | { // find variables of weight ii |
---|
1171 | k = 0; |
---|
1172 | index = 0; |
---|
1173 | // collect the indices of all variables of weight i |
---|
1174 | for(j = 1; j <= size(wt); j++) |
---|
1175 | { |
---|
1176 | if(wt[j] == ii) |
---|
1177 | { |
---|
1178 | k++; |
---|
1179 | index[k] = j; |
---|
1180 | } |
---|
1181 | } |
---|
1182 | uw++; |
---|
1183 | usedWeights[uw] = ii; |
---|
1184 | ic++; |
---|
1185 | indices[ii] = index; |
---|
1186 | |
---|
1187 | // linear part of the substitution |
---|
1188 | |
---|
1189 | for(j = 1; j <= size(index); j++) |
---|
1190 | { |
---|
1191 | linearSubs[index[j]] = index; |
---|
1192 | } |
---|
1193 | |
---|
1194 | // nonlinear part of the substitution |
---|
1195 | |
---|
1196 | if(uw > 1) |
---|
1197 | { // variables of least weight do not allow nonlinear subs. |
---|
1198 | |
---|
1199 | partitions = Partitions(ii, delete(usedWeights, uw)); |
---|
1200 | for(j = 1; j <= size(partitions); j++) |
---|
1201 | { |
---|
1202 | combinations[j] = AllCombinations(partitions[j], indices); |
---|
1203 | } |
---|
1204 | for(j = 1; j <= size(index); j++) |
---|
1205 | { |
---|
1206 | nonlinearSubs[index[j]] = FlattenQHM(combinations); // flatten one level ! |
---|
1207 | } |
---|
1208 | } |
---|
1209 | } |
---|
1210 | } |
---|
1211 | combList[1] = linearSubs; |
---|
1212 | combList[2] = nonlinearSubs; |
---|
1213 | return(combList); |
---|
1214 | } |
---|
1215 | |
---|
1216 | /////////////////////////////////////////////////////////////////////////////// |
---|
1217 | |
---|
1218 | static proc AllCombinations(list partition, list indices) |
---|
1219 | "USAGE: AllCombinations(partition,indices); list partition, indices) |
---|
1220 | PURPOSE: all combinations for a given partititon |
---|
1221 | RETURN: list |
---|
1222 | GLOBAL: varSubsList, contains the index j s.t. x(i) -> x(i)t(j) ... |
---|
1223 | " |
---|
1224 | { |
---|
1225 | int i, k, m, ok, p, offset; |
---|
1226 | list nrList, indexList; |
---|
1227 | |
---|
1228 | k = 0; |
---|
1229 | offset = 0; |
---|
1230 | i = 1; |
---|
1231 | ok = 1; |
---|
1232 | m = partition[1]; |
---|
1233 | while(ok) |
---|
1234 | { |
---|
1235 | if(i > size(partition)) |
---|
1236 | { |
---|
1237 | ok = 0; |
---|
1238 | p = 0; |
---|
1239 | } |
---|
1240 | else { p = partition[i];} |
---|
1241 | if(p == m) { i = i + 1;} |
---|
1242 | else |
---|
1243 | { |
---|
1244 | k = k + 1; |
---|
1245 | nrList[k] = i - 1 - offset; |
---|
1246 | offset = offset + i - 1; |
---|
1247 | indexList[k] = indices[m]; |
---|
1248 | if(ok) { m = partition[i];} |
---|
1249 | } |
---|
1250 | } |
---|
1251 | return(AllCombinationsAux(nrList, indexList)); |
---|
1252 | } |
---|
1253 | |
---|
1254 | /////////////////////////////////////////////////////////////////////////////// |
---|
1255 | |
---|
1256 | static proc AllSingleCombinations(int n, list index) |
---|
1257 | "USAGE: AllSingleCombinations(n index); int n, list index |
---|
1258 | PURPOSE: all combinations for var(n) |
---|
1259 | RETURN: list |
---|
1260 | " |
---|
1261 | { |
---|
1262 | int i, j, k; |
---|
1263 | list comb, newC, temp, newIndex; |
---|
1264 | |
---|
1265 | if(n == 1) |
---|
1266 | { |
---|
1267 | for(i = 1; i <= size(index); i++) |
---|
1268 | { |
---|
1269 | temp = index[i]; |
---|
1270 | comb[i] = temp; |
---|
1271 | } |
---|
1272 | return(comb); |
---|
1273 | } |
---|
1274 | if(size(index) == 1) |
---|
1275 | { |
---|
1276 | temp = Table(string(index[1]), "i", 1, n); |
---|
1277 | comb[1] = temp; |
---|
1278 | return(comb); |
---|
1279 | } |
---|
1280 | newIndex = index; |
---|
1281 | for(i = 1; i <= size(index); i = i + 1) |
---|
1282 | { |
---|
1283 | if(i > 1) { newIndex = delete(newIndex, 1); } |
---|
1284 | newC = AllSingleCombinations(n - 1, newIndex); |
---|
1285 | k = size(comb); |
---|
1286 | temp = 0; |
---|
1287 | for(j = 1; j <= size(newC); j++) |
---|
1288 | { |
---|
1289 | temp[1] = index[i]; |
---|
1290 | temp = temp + newC[j]; |
---|
1291 | comb[k + j] = temp; |
---|
1292 | temp = 0; |
---|
1293 | } |
---|
1294 | } |
---|
1295 | return(comb); |
---|
1296 | } |
---|
1297 | |
---|
1298 | /////////////////////////////////////////////////////////////////////////////// |
---|
1299 | |
---|
1300 | static proc AllCombinationsAux(list parts, list index) |
---|
1301 | "USAGE: AllCombinationsAux(parts ,index); list parts, index |
---|
1302 | PURPOSE: all compbinations for nonlinear substituiton |
---|
1303 | RETURN: list |
---|
1304 | " |
---|
1305 | { |
---|
1306 | int i, j, k; |
---|
1307 | list comb, firstC, restC; |
---|
1308 | |
---|
1309 | if(size(parts) == 0 || size(index) == 0) { return(comb);} |
---|
1310 | |
---|
1311 | firstC = AllSingleCombinations(parts[1], index[1]); |
---|
1312 | restC = AllCombinationsAux(delete(parts, 1), delete(index, 1)); |
---|
1313 | |
---|
1314 | if(size(restC) == 0) { comb = firstC;} |
---|
1315 | else |
---|
1316 | { |
---|
1317 | for(i = 1; i <= size(firstC); i++) |
---|
1318 | { |
---|
1319 | k = size(comb); |
---|
1320 | for(j = 1; j <= size(restC); j++) |
---|
1321 | { |
---|
1322 | //elem = firstC[i] + restC[j]; |
---|
1323 | // comb[k + j] = elem; |
---|
1324 | comb[k + j] = firstC[i] + restC[j]; |
---|
1325 | } |
---|
1326 | } |
---|
1327 | } |
---|
1328 | return(comb); |
---|
1329 | } |
---|
1330 | |
---|
1331 | /////////////////////////////////////////////////////////////////////////////// |
---|
1332 | |
---|
1333 | static proc Partitions(int n, list nr) |
---|
1334 | "USAGE: Partitions(n, nr); int n, list nr |
---|
1335 | PURPOSE: partitions of n consisting of elements from nr |
---|
1336 | RETURN: list |
---|
1337 | " |
---|
1338 | { |
---|
1339 | int i, j, k; |
---|
1340 | list parts, temp, restP, newP, decP; |
---|
1341 | |
---|
1342 | if(size(nr) == 0) { return(list());} |
---|
1343 | if(size(nr) == 1) |
---|
1344 | { |
---|
1345 | if(NumFactor(nr[1], n) > 0) |
---|
1346 | { |
---|
1347 | parts[1] = Table(string(nr[1]), "i", 1, NumFactor(nr[1], n)); |
---|
1348 | } |
---|
1349 | return(parts); |
---|
1350 | } |
---|
1351 | else |
---|
1352 | { |
---|
1353 | parts = Partitions(n, nr[1]); |
---|
1354 | restP = Partitions(n, delete(nr, 1)); |
---|
1355 | |
---|
1356 | parts = parts + restP; |
---|
1357 | for(i = 1; i <= n / nr[1]; i = i + 1) |
---|
1358 | { |
---|
1359 | temp = Table(string(nr[1]), "i", 1, i); |
---|
1360 | decP = Partitions(n - i*nr[1], delete(nr, 1)); |
---|
1361 | k = size(parts); |
---|
1362 | for(j = 1; j <= size(decP); j++) |
---|
1363 | { |
---|
1364 | newP = temp + decP[j]; // new partition |
---|
1365 | if(!ContainedQ(parts, newP, 1)) |
---|
1366 | { |
---|
1367 | k = k + 1; |
---|
1368 | parts[k] = newP; |
---|
1369 | } |
---|
1370 | } |
---|
1371 | } |
---|
1372 | } |
---|
1373 | return(parts); |
---|
1374 | } |
---|
1375 | |
---|
1376 | /////////////////////////////////////////////////////////////////////////////// |
---|
1377 | |
---|
1378 | static proc NumFactor(int a, int b) |
---|
1379 | " USAGE: NumFactor(a, b); int a, b |
---|
1380 | PURPOSE: if b divides a then return b/a, else return 0 |
---|
1381 | RETURN: int |
---|
1382 | " |
---|
1383 | { |
---|
1384 | int c = int(b/a); |
---|
1385 | if(c*a == b) { return(c); } |
---|
1386 | else {return(0)} |
---|
1387 | } |
---|
1388 | |
---|
1389 | /////////////////////////////////////////////////////////////////////////////// |
---|
1390 | |
---|
1391 | static proc Table(string cmd, string iterator, int lb, int ub) |
---|
1392 | " USAGE: Table(cmd,i, lb, ub); string cmd, i; int lb, ub |
---|
1393 | PURPOSE: generate a list of size ub - lb + 1 s.t. _[i] = cmd(i) |
---|
1394 | RETURN: list |
---|
1395 | " |
---|
1396 | { |
---|
1397 | list data; |
---|
1398 | execute("int " + iterator + ";"); |
---|
1399 | |
---|
1400 | for(int @i = lb; @i <= ub; @i++) |
---|
1401 | { |
---|
1402 | execute(iterator + " = " + string(@i)); |
---|
1403 | execute("data[" + string(@i) + "] = " + cmd + ";"); |
---|
1404 | } |
---|
1405 | return(data); |
---|
1406 | } |
---|
1407 | |
---|
1408 | /////////////////////////////////////////////////////////////////////////////// |
---|
1409 | |
---|
1410 | static proc FlattenQHM(list data) |
---|
1411 | " USAGE: FlattenQHM(n, nr); list data |
---|
1412 | PURPOSE: flatten the list (one level) 'data', which is a list of lists |
---|
1413 | RETURN: list |
---|
1414 | " |
---|
1415 | { |
---|
1416 | int i, j, c; |
---|
1417 | list fList, temp; |
---|
1418 | |
---|
1419 | c = 1; |
---|
1420 | |
---|
1421 | for(i = 1; i <= size(data); i++) |
---|
1422 | { |
---|
1423 | for(j = 1; j <= size(data[i]); j++) |
---|
1424 | { |
---|
1425 | fList[c] = data[i][j]; |
---|
1426 | c = c + 1; |
---|
1427 | } |
---|
1428 | } |
---|
1429 | return(fList); |
---|
1430 | } |
---|
1431 | |
---|
1432 | /////////////////////////////////////////////////////////////////////////////// |
---|
1433 | |
---|
1434 | static proc IntersectionQHM(list l1, list l2) |
---|
1435 | // Type : list |
---|
1436 | // Purpose : Intersection of l1 and l2 |
---|
1437 | { |
---|
1438 | list l; |
---|
1439 | int b, c; |
---|
1440 | |
---|
1441 | c = 1; |
---|
1442 | |
---|
1443 | for(int i = 1; i <= size(l1); i++) |
---|
1444 | { |
---|
1445 | b = ContainedQ(l2, l1[i]); |
---|
1446 | if(b == 1) |
---|
1447 | { |
---|
1448 | l[c] = l1[i]; |
---|
1449 | c++; |
---|
1450 | } |
---|
1451 | } |
---|
1452 | return(l); |
---|
1453 | } |
---|
1454 | |
---|
1455 | /////////////////////////////////////////////////////////////////////////////// |
---|
1456 | |
---|
1457 | static proc FirstEntryQHM(data, elem) |
---|
1458 | // Type : int |
---|
1459 | // Purpose : position of first entry equal to elem in data (from left to right) |
---|
1460 | { |
---|
1461 | int i, pos; |
---|
1462 | |
---|
1463 | i = 0; |
---|
1464 | pos = 0; |
---|
1465 | while(i < size(data)) |
---|
1466 | { |
---|
1467 | i++; |
---|
1468 | if(data[i] == elem) { pos = i; break;} |
---|
1469 | } |
---|
1470 | return(pos); |
---|
1471 | } |
---|
1472 | |
---|
1473 | /////////////////////////////////////////////////////////////////////////////// |
---|
1474 | |
---|
1475 | static proc PSum(e) |
---|
1476 | { |
---|
1477 | poly f; |
---|
1478 | for(int i = size(e);i>=1;i--) |
---|
1479 | { |
---|
1480 | f = f + e[i]; |
---|
1481 | } |
---|
1482 | return(f); |
---|
1483 | } |
---|
1484 | |
---|
1485 | /////////////////////////////////////////////////////////////////////////////// |
---|
1486 | |
---|
1487 | proc Max(data) |
---|
1488 | "USAGE: Max(data); intvec/list of integers |
---|
1489 | PURPOSE: find the maximal integer contained in 'data' |
---|
1490 | RETURN: list |
---|
1491 | ASSUME: 'data' contians only integers and is not empty |
---|
1492 | " |
---|
1493 | { |
---|
1494 | int i; |
---|
1495 | int max = data[1]; |
---|
1496 | |
---|
1497 | for(i = size(data); i>1;i--) |
---|
1498 | { |
---|
1499 | if(data[i] > max) { max = data[i]; } |
---|
1500 | } |
---|
1501 | return(max); |
---|
1502 | } |
---|
1503 | example |
---|
1504 | {"EXAMPLE:"; echo = 2; |
---|
1505 | Max(list(1,2,3)); |
---|
1506 | } |
---|
1507 | |
---|
1508 | /////////////////////////////////////////////////////////////////////////////// |
---|
1509 | |
---|
1510 | proc Min(data) |
---|
1511 | "USAGE: Min(data); intvec/list of integers |
---|
1512 | PURPOSE: find the minimal integer contained in 'data' |
---|
1513 | RETURN: list |
---|
1514 | ASSUME: 'data' contians only integers and is not empty |
---|
1515 | " |
---|
1516 | { |
---|
1517 | int i; |
---|
1518 | int min = data[1]; |
---|
1519 | |
---|
1520 | for(i = size(data);i>1; i--) |
---|
1521 | { |
---|
1522 | if(data[i] < min) { min = data[i]; } |
---|
1523 | } |
---|
1524 | return(min); |
---|
1525 | } |
---|
1526 | example |
---|
1527 | {"EXAMPLE:"; echo = 2; |
---|
1528 | Min(intvec(1,2,3)); |
---|
1529 | } |
---|
1530 | |
---|
1531 | /////////////////////////////////////////////////////////////////////////////// |
---|