1 | //////////////////////////////////////////////////////////////////////////// |
---|
2 | version="version realclassify.lib 4.0.0.0 Jun_2013 "; // $Id$ |
---|
3 | category="Singularities"; |
---|
4 | info=" |
---|
5 | LIBRARY: realclassify.lib Classification of real singularities |
---|
6 | AUTHOR: Magdaleen Marais, magdaleen@aims.ac.za |
---|
7 | Andreas Steenpass, steenpass@mathematik.uni-kl.de |
---|
8 | |
---|
9 | OVERVIEW: |
---|
10 | A library for classifying isolated hypersurface singularities over the reals |
---|
11 | w.r.t. right equivalence, based on the determinator of singularities by |
---|
12 | V.I. Arnold. This library is based on classify.lib by Kai Krueger, but |
---|
13 | handles the real case, while classify.lib does the complex classification. |
---|
14 | |
---|
15 | REFERENCES: |
---|
16 | Arnold, Varchenko, Gusein-Zade: Singularities of Differentiable Maps. |
---|
17 | Vol. 1: The classification of critical points caustics and wave fronts. |
---|
18 | Birkh\"auser, Boston 1985 |
---|
19 | |
---|
20 | Greuel, Lossen, Shustin: Introduction to singularities and deformations. |
---|
21 | Springer, Berlin 2007 |
---|
22 | |
---|
23 | PROCEDURES: |
---|
24 | realclassify(f); real classification of singularities of modality 0 and 1 |
---|
25 | realmorsesplit(f); splitting lemma in the real case |
---|
26 | milnornumber(f); Milnor number |
---|
27 | determinacy(f); an upper bound for the determinacy |
---|
28 | "; |
---|
29 | |
---|
30 | LIB "linalg.lib"; |
---|
31 | LIB "elim.lib"; |
---|
32 | LIB "primdec.lib"; |
---|
33 | LIB "classify.lib"; |
---|
34 | LIB "rootsur.lib"; |
---|
35 | LIB "rootsmr.lib"; |
---|
36 | LIB "atkins.lib"; |
---|
37 | LIB "solve.lib"; |
---|
38 | /////////////////////////////////////////////////////////////////////////////// |
---|
39 | proc realclassify(poly f, list #) |
---|
40 | " |
---|
41 | USAGE: realclassify(f[, format]); f poly, format string |
---|
42 | RETURN: A list containing (in this order) |
---|
43 | @* - the type of the singularity as a string, |
---|
44 | @* - the normal form, |
---|
45 | @* - the corank, the Milnor number, the inertia index and |
---|
46 | a bound for the determinacy as integers. |
---|
47 | @* The normal form involves parameters for singularities of modality |
---|
48 | greater than 0. The actual value of the parameters is not computed |
---|
49 | in most of the cases. If the value of the parameter is unknown, |
---|
50 | the normal form is given as a string with an \"a\" as the |
---|
51 | parameter. Otherwise, it is given as a polynomial. |
---|
52 | @* An optional string @code{format} can be provided. Its default |
---|
53 | value is \"short\" in which case the return value is the list |
---|
54 | described above. If set to \"nice\", a string is added at the end |
---|
55 | of this list, containing the result in a more readable form. |
---|
56 | NOTE: The classification is done over the real numbers, so in contrast to |
---|
57 | classify.lib, the signs of coefficients of monomials where even |
---|
58 | exponents occur matter. |
---|
59 | @* The ground field must be Q (the rational numbers). No field |
---|
60 | extensions of any kind nor floating point numbers are allowed. |
---|
61 | @* The monomial order must be local. |
---|
62 | @* The input polynomial must be contained in maxideal(2) and must be |
---|
63 | an isolated singularity of modality 0 or 1. The Milnor number is |
---|
64 | checked for being finite. |
---|
65 | SEE ALSO: classify |
---|
66 | KEYWORDS: Classification of singularities |
---|
67 | EXAMPLE: example realclassify; shows an example" |
---|
68 | { |
---|
69 | /* auxiliary variables */ |
---|
70 | int i, j; |
---|
71 | |
---|
72 | /* name for the basering */ |
---|
73 | def br = basering; |
---|
74 | |
---|
75 | /* read optional parameters */ |
---|
76 | int printcomments; |
---|
77 | if(size(#) > 0) |
---|
78 | { |
---|
79 | if(size(#) > 1 || typeof(#[1]) != "string") |
---|
80 | { |
---|
81 | ERROR("Wrong optional parameters."); |
---|
82 | } |
---|
83 | if(#[1] != "short" && #[1] != "nice") |
---|
84 | { |
---|
85 | ERROR("Wrong optional parameters."); |
---|
86 | } |
---|
87 | if(#[1] == "nice") |
---|
88 | { |
---|
89 | printcomments = 1; |
---|
90 | } |
---|
91 | } |
---|
92 | |
---|
93 | /* error check */ |
---|
94 | if(charstr(br) != "0") |
---|
95 | { |
---|
96 | ERROR("The ground field must be Q (the rational numbers)."); |
---|
97 | } |
---|
98 | int n = nvars(br); |
---|
99 | for(i = 1; i <= n; i++) |
---|
100 | { |
---|
101 | if(var(i) > 1) |
---|
102 | { |
---|
103 | ERROR("The monomial order must be local."); |
---|
104 | } |
---|
105 | } |
---|
106 | if(jet(f, 1) != 0) |
---|
107 | { |
---|
108 | ERROR("The input polynomial must be contained in maxideal(2)."); |
---|
109 | } |
---|
110 | |
---|
111 | /* compute Milnor number before continuing the error check */ |
---|
112 | int mu = milnornumber(f); |
---|
113 | |
---|
114 | /* continue error check */ |
---|
115 | if(mu < 1) |
---|
116 | { |
---|
117 | ERROR("The Milnor number of the input polynomial must be"+newline |
---|
118 | +"positive and finite."); |
---|
119 | } |
---|
120 | |
---|
121 | /* call classify before continuing the error check */ |
---|
122 | list dataFromClassify = prepRealclassify(f); |
---|
123 | int m = dataFromClassify[1]; // the modality of f |
---|
124 | string complextype = dataFromClassify[2]; // the complex type of f |
---|
125 | |
---|
126 | /* continue error check */ |
---|
127 | if(m > 1) |
---|
128 | { |
---|
129 | ERROR("The input polynomial must be a singularity of modality 0 or 1."); |
---|
130 | } |
---|
131 | |
---|
132 | /* apply splitting lemma */ |
---|
133 | list morse = realmorsesplit(f, mu); |
---|
134 | int cr = morse[1]; |
---|
135 | int lambda = morse[2]; |
---|
136 | int d = morse[3]; |
---|
137 | poly rf = morse[4]; |
---|
138 | |
---|
139 | /* determine the type */ |
---|
140 | string typeofsing; |
---|
141 | poly nf; |
---|
142 | poly monparam; // the monomial whose coefficient is the parameter |
---|
143 | // in the modality 1 cases, 0 otherwise |
---|
144 | string morecomments = newline; |
---|
145 | |
---|
146 | if(cr == 0) // case A[1] |
---|
147 | { |
---|
148 | typeofsing, nf = caseA1(rf, lambda, n); |
---|
149 | } |
---|
150 | if(cr == 1) // case A[k], k > 1 |
---|
151 | { |
---|
152 | typeofsing, nf = caseAk(rf, n); |
---|
153 | } |
---|
154 | if(cr == 2) |
---|
155 | { |
---|
156 | if(complextype[1,2] == "D[") // case D[k] |
---|
157 | { |
---|
158 | if(mu == 4) // case D[4] |
---|
159 | { |
---|
160 | typeofsing, nf = caseD4(rf); |
---|
161 | } |
---|
162 | else // case D[k], k > 4 |
---|
163 | { |
---|
164 | typeofsing, nf = caseDk(rf, mu); |
---|
165 | } |
---|
166 | } |
---|
167 | if(complextype == "E[6]") // case E[6] |
---|
168 | { |
---|
169 | typeofsing, nf = caseE6(rf); |
---|
170 | } |
---|
171 | if(complextype == "E[7]") // case E[7] |
---|
172 | { |
---|
173 | typeofsing, nf = caseE7(); |
---|
174 | } |
---|
175 | if(complextype == "E[8]") // case E[8] |
---|
176 | { |
---|
177 | typeofsing, nf = caseE8(); |
---|
178 | } |
---|
179 | if(typeofsing == "") |
---|
180 | { |
---|
181 | ERROR("This case is not yet implemented."); |
---|
182 | } |
---|
183 | } |
---|
184 | if(cr > 2) |
---|
185 | { |
---|
186 | ERROR("This case is not yet implemented."); |
---|
187 | } |
---|
188 | |
---|
189 | /* add the non-corank variables to the normal forms */ |
---|
190 | nf = addnondegeneratevariables(nf, lambda, cr); |
---|
191 | |
---|
192 | /* write normal form as a string in the cases with modality greater than 0 */ |
---|
193 | if(monparam != 0) |
---|
194 | { |
---|
195 | poly nf_tmp = nf; |
---|
196 | kill nf; |
---|
197 | def nf = modality1NF(nf_tmp, monparam); |
---|
198 | } |
---|
199 | |
---|
200 | /* write comments */ |
---|
201 | if(printcomments) |
---|
202 | { |
---|
203 | string comments = newline; |
---|
204 | comments = comments+"Type of singularity: " +typeofsing +newline |
---|
205 | +"Normal form: " +string(nf) +newline |
---|
206 | +"Corank: " +string(cr) +newline |
---|
207 | +"Milnor number: " +string(mu) +newline |
---|
208 | +"Inertia index: " +string(lambda)+newline |
---|
209 | +"Determinacy: <= "+string(d) +newline; |
---|
210 | if(morecomments != newline) |
---|
211 | { |
---|
212 | comments = comments+morecomments; |
---|
213 | } |
---|
214 | } |
---|
215 | |
---|
216 | /* return results */ |
---|
217 | if(printcomments) |
---|
218 | { |
---|
219 | return(list(typeofsing, nf, cr, mu, lambda, d, comments)); |
---|
220 | } |
---|
221 | else |
---|
222 | { |
---|
223 | return(list(typeofsing, nf, cr, mu, lambda, d)); |
---|
224 | } |
---|
225 | } |
---|
226 | example |
---|
227 | { |
---|
228 | "EXAMPLE:"; |
---|
229 | echo = 2; |
---|
230 | ring r = 0, (x,y,z), ds; |
---|
231 | poly f = (x2+3y-2z)^2+xyz-(x-y3+x2z3)^3; |
---|
232 | realclassify(f, "nice"); |
---|
233 | } |
---|
234 | |
---|
235 | /////////////////////////////////////////////////////////////////////////////// |
---|
236 | static proc caseA1(poly rf, int lambda, int n) |
---|
237 | { |
---|
238 | string typeofsing = "A[1]"; |
---|
239 | poly nf = 0; |
---|
240 | return(typeofsing, nf); |
---|
241 | } |
---|
242 | |
---|
243 | /////////////////////////////////////////////////////////////////////////////// |
---|
244 | static proc caseAk(poly rf, int n) |
---|
245 | { |
---|
246 | /* preliminaries */ |
---|
247 | string typeofsing; |
---|
248 | poly nf; |
---|
249 | |
---|
250 | int k = deg(lead(rf), 1:n)-1; |
---|
251 | if(k%2 == 0) |
---|
252 | { |
---|
253 | nf = var(1)^(k+1); |
---|
254 | typeofsing = "A["+string(k)+"]"; |
---|
255 | } |
---|
256 | else |
---|
257 | { |
---|
258 | if(leadcoef(rf) > 0) |
---|
259 | { |
---|
260 | nf = var(1)^(k+1); |
---|
261 | typeofsing = "A["+string(k)+"]+"; |
---|
262 | } |
---|
263 | else |
---|
264 | { |
---|
265 | nf = -var(1)^(k+1); |
---|
266 | typeofsing = "A["+string(k)+"]-"; |
---|
267 | } |
---|
268 | } |
---|
269 | return(typeofsing, nf); |
---|
270 | } |
---|
271 | |
---|
272 | /////////////////////////////////////////////////////////////////////////////// |
---|
273 | static proc caseD4(poly rf) |
---|
274 | { |
---|
275 | /* preliminaries */ |
---|
276 | string typeofsing; |
---|
277 | poly nf; |
---|
278 | def br = basering; |
---|
279 | map phi; |
---|
280 | |
---|
281 | rf = jet(rf, 3); |
---|
282 | number s1 = number(rf/(var(1)^3)); |
---|
283 | number s2 = number(rf/(var(2)^3)); |
---|
284 | if(s2 == 0 && s1 != 0) |
---|
285 | { |
---|
286 | phi = br, var(2), var(1); |
---|
287 | rf = phi(rf); |
---|
288 | } |
---|
289 | if(s1 == 0 && s2 == 0) |
---|
290 | { |
---|
291 | number t1 = number(rf/(var(1)^2*var(2))); |
---|
292 | number t2 = number(rf/(var(2)^2*var(1))); |
---|
293 | if(t1+t2 == 0) |
---|
294 | { |
---|
295 | phi = br, var(1)+2*var(2), var(2); |
---|
296 | rf = phi(rf); |
---|
297 | } |
---|
298 | else |
---|
299 | { |
---|
300 | phi = br, var(1)+var(2), var(2); |
---|
301 | rf = phi(rf); |
---|
302 | } |
---|
303 | } |
---|
304 | ring R = 0, y, dp; |
---|
305 | map phi = br, 1, y; |
---|
306 | poly rf = phi(rf); |
---|
307 | int k = nrroots(rf); |
---|
308 | setring(br); |
---|
309 | if(k == 3) |
---|
310 | { |
---|
311 | nf = var(1)^2*var(2)-var(2)^3; |
---|
312 | typeofsing = "D[4]-"; |
---|
313 | } |
---|
314 | else |
---|
315 | { |
---|
316 | nf = var(1)^2*var(2)+var(2)^3; |
---|
317 | typeofsing = "D[4]+"; |
---|
318 | } |
---|
319 | return(typeofsing, nf); |
---|
320 | } |
---|
321 | |
---|
322 | /////////////////////////////////////////////////////////////////////////////// |
---|
323 | static proc caseDk(poly rf, int mu) |
---|
324 | { |
---|
325 | /* preliminaries */ |
---|
326 | string typeofsing; |
---|
327 | poly nf; |
---|
328 | def br = basering; |
---|
329 | map phi; |
---|
330 | |
---|
331 | rf = jet(rf, mu-1); |
---|
332 | list factorization = factorize(jet(rf, 3)); |
---|
333 | list factors = factorization[1][2]; |
---|
334 | if(factorization[2][2] == 2) |
---|
335 | { |
---|
336 | factors = insert(factors, factorization[1][3], 1); |
---|
337 | } |
---|
338 | else |
---|
339 | { |
---|
340 | factors = insert(factors, factorization[1][3]); |
---|
341 | } |
---|
342 | factors[2] = factorization[1][1]*factors[2]; |
---|
343 | matrix T[2][2] = factors[1]/var(1), factors[1]/var(2), |
---|
344 | factors[2]/var(1), factors[2]/var(2); |
---|
345 | phi = br, luinverse(T)[2]*matrix(ideal(var(1), var(2)), 2, 1); |
---|
346 | rf = phi(rf); |
---|
347 | rf = jet(rf, mu-1); |
---|
348 | poly g; |
---|
349 | int i; |
---|
350 | for(i = 4; i < mu; i++) |
---|
351 | { |
---|
352 | g = jet(rf, i) - var(1)^2*var(2); |
---|
353 | if(g != 0) |
---|
354 | { |
---|
355 | phi = br, var(1)-(g/(var(1)*var(2)))/2, |
---|
356 | var(2)-(g/var(1)^i)*var(1)^(i-2); |
---|
357 | rf = phi(rf); |
---|
358 | rf = jet(rf, mu-1); |
---|
359 | } |
---|
360 | } |
---|
361 | number a = number(rf/var(2)^(mu-1)); |
---|
362 | if(a > 0) |
---|
363 | { |
---|
364 | typeofsing = "D["+string(mu)+"]+"; |
---|
365 | nf = var(1)^2*var(2)+var(2)^(mu-1); |
---|
366 | } |
---|
367 | else |
---|
368 | { |
---|
369 | typeofsing = "D["+string(mu)+"]-"; |
---|
370 | nf = var(1)^2*var(2)-var(2)^(mu-1); |
---|
371 | } |
---|
372 | return(typeofsing, nf); |
---|
373 | } |
---|
374 | |
---|
375 | /////////////////////////////////////////////////////////////////////////////// |
---|
376 | static proc caseE6(poly rf) |
---|
377 | { |
---|
378 | /* preliminaries */ |
---|
379 | string typeofsing; |
---|
380 | poly nf; |
---|
381 | def br = basering; |
---|
382 | map phi; |
---|
383 | |
---|
384 | poly g = jet(rf,3); |
---|
385 | number s = number(g/(var(1)^3)); |
---|
386 | if(s == 0) |
---|
387 | { |
---|
388 | phi = br, var(2), var(1); |
---|
389 | rf = phi(rf); |
---|
390 | g = jet(rf,3); |
---|
391 | } |
---|
392 | list Factors = factorize(g); |
---|
393 | poly g1 = Factors[1][2]; |
---|
394 | phi = br, (var(1)-(g1/var(2))*var(2))/(g1/var(1)), var(2); |
---|
395 | rf = phi(rf); |
---|
396 | rf = jet(rf,4); |
---|
397 | number w = number(rf/(var(2)^4)); |
---|
398 | if(w > 0) |
---|
399 | { |
---|
400 | typeofsing = "E[6]+"; |
---|
401 | nf = var(1)^3+var(2)^4; |
---|
402 | } |
---|
403 | else |
---|
404 | { |
---|
405 | typeofsing = "E[6]-"; |
---|
406 | nf = var(1)^3-var(2)^4; |
---|
407 | } |
---|
408 | return(typeofsing, nf); |
---|
409 | } |
---|
410 | |
---|
411 | /////////////////////////////////////////////////////////////////////////////// |
---|
412 | static proc caseE7() |
---|
413 | { |
---|
414 | string typeofsing = "E[7]"; |
---|
415 | poly nf = var(1)^3+var(1)*var(2)^3; |
---|
416 | return(typeofsing, nf); |
---|
417 | } |
---|
418 | |
---|
419 | /////////////////////////////////////////////////////////////////////////////// |
---|
420 | static proc caseE8() |
---|
421 | { |
---|
422 | string typeofsing = "E[8]"; |
---|
423 | poly nf = var(1)^3+var(2)^5; |
---|
424 | return(typeofsing, nf); |
---|
425 | } |
---|
426 | |
---|
427 | /////////////////////////////////////////////////////////////////////////////// |
---|
428 | /* |
---|
429 | print the normal form as a string for the modality 1 cases. |
---|
430 | The first argument is the normalform with parameter = 1, |
---|
431 | the second argument is the monomial whose coefficient is the parameter. |
---|
432 | */ |
---|
433 | static proc modality1NF(poly nf, poly monparam) |
---|
434 | { |
---|
435 | def br = basering; |
---|
436 | list lbr = ringlist(br); |
---|
437 | ring r = (0,a), x, dp; |
---|
438 | list lr = ringlist(r); |
---|
439 | setring(br); |
---|
440 | list lr = fetch(r, lr); |
---|
441 | lbr[1] = lr[1]; |
---|
442 | def s = ring(lbr); |
---|
443 | setring(s); |
---|
444 | poly nf = fetch(br, nf); |
---|
445 | poly monparam = fetch(br, monparam); |
---|
446 | nf = nf+(a-1)*monparam; |
---|
447 | string result = string(nf); |
---|
448 | setring(br); |
---|
449 | return(result); |
---|
450 | } |
---|
451 | |
---|
452 | /////////////////////////////////////////////////////////////////////////////// |
---|
453 | /* |
---|
454 | add squares of the non-degenerate variables (i.e. var(cr+1), ..., |
---|
455 | var(nvars(basering)) for corank cr) to the normalform nf, |
---|
456 | with signs according to the inertia index lambda |
---|
457 | */ |
---|
458 | static proc addnondegeneratevariables(poly nf, int lambda, int cr) |
---|
459 | { |
---|
460 | int n = nvars(basering); |
---|
461 | int i; |
---|
462 | for(i = cr+1; i <= n-lambda; i++) |
---|
463 | { |
---|
464 | nf = nf+var(i)^2; |
---|
465 | } |
---|
466 | for(i = n-lambda+1; i <= n ; i++) |
---|
467 | { |
---|
468 | nf = nf-var(i)^2; |
---|
469 | } |
---|
470 | return(nf); |
---|
471 | } |
---|
472 | |
---|
473 | /////////////////////////////////////////////////////////////////////////////// |
---|
474 | proc realmorsesplit(poly f, list #) |
---|
475 | " |
---|
476 | USAGE: realmorsesplit(f[, mu]); f poly, mu int |
---|
477 | RETURN: a list consisting of the corank of f, the inertia index, an upper |
---|
478 | bound for the determinacy, and the residual form of f |
---|
479 | NOTE: The characteristic of the basering must be zero, the monomial order |
---|
480 | must be local, f must be contained in maxideal(2) and the Milnor |
---|
481 | number of f must be finite. |
---|
482 | @* The Milnor number of f can be provided as an optional parameter in |
---|
483 | order to avoid that it is computed again. |
---|
484 | SEE ALSO: morsesplit |
---|
485 | KEYWORDS: Morse lemma; Splitting lemma |
---|
486 | EXAMPLE: example morsesplit; shows an example" |
---|
487 | { |
---|
488 | int i; |
---|
489 | |
---|
490 | /* error check */ |
---|
491 | if(char(basering) != 0) |
---|
492 | { |
---|
493 | ERROR("The characteristic must be zero."); |
---|
494 | } |
---|
495 | int n = nvars(basering); |
---|
496 | for(i = 1; i <= n; i++) |
---|
497 | { |
---|
498 | if(var(i) > 1) |
---|
499 | { |
---|
500 | ERROR("The monomial order must be local."); |
---|
501 | } |
---|
502 | } |
---|
503 | if(jet(f, 1) != 0) |
---|
504 | { |
---|
505 | ERROR("The input polynomial must be contained in maxideal(2)."); |
---|
506 | } |
---|
507 | |
---|
508 | /* get Milnor number before continuing error check */ |
---|
509 | int mu; |
---|
510 | if(size(#) > 0) // read optional parameter |
---|
511 | { |
---|
512 | if(size(#) > 1 || typeof(#[1]) != "int") |
---|
513 | { |
---|
514 | ERROR("Wrong optional parameters."); |
---|
515 | } |
---|
516 | else |
---|
517 | { |
---|
518 | mu = #[1]; |
---|
519 | } |
---|
520 | } |
---|
521 | else // compute Milnor number |
---|
522 | { |
---|
523 | mu = milnornumber(f); |
---|
524 | } |
---|
525 | |
---|
526 | /* continue error check */ |
---|
527 | if(mu < 0) |
---|
528 | { |
---|
529 | ERROR("The Milnor number of the input polynomial must be"+newline |
---|
530 | +"non-negative and finite."); |
---|
531 | } |
---|
532 | |
---|
533 | /* compute the determinacy */ |
---|
534 | int k = determinacy(f, mu); |
---|
535 | f = jet(f, k); |
---|
536 | |
---|
537 | /* get jet(f, 2) right */ |
---|
538 | matrix H = concat(jet(jacob(jacob(f)), 0)/2, unitmat(n)); |
---|
539 | H = sym_reduce(H); |
---|
540 | intvec perm_zero; |
---|
541 | intvec perm_neg; |
---|
542 | intvec perm_pos; |
---|
543 | int c; |
---|
544 | int lambda; |
---|
545 | for(i = 1; i <= n; i++) |
---|
546 | { |
---|
547 | if(H[i, i] == 0) |
---|
548 | { |
---|
549 | perm_zero = perm_zero, i; |
---|
550 | c++; |
---|
551 | } |
---|
552 | if(H[i, i] < 0) |
---|
553 | { |
---|
554 | perm_neg = perm_neg, i; |
---|
555 | lambda++; |
---|
556 | } |
---|
557 | if(H[i, i] > 0) |
---|
558 | { |
---|
559 | perm_pos = perm_pos, i; |
---|
560 | } |
---|
561 | } |
---|
562 | intvec perm; |
---|
563 | if(size(perm_zero) > 1) |
---|
564 | { |
---|
565 | perm = perm, perm_zero[2..size(perm_zero)]; |
---|
566 | } |
---|
567 | if(size(perm_neg) > 1) |
---|
568 | { |
---|
569 | perm = perm, perm_neg[2..size(perm_neg)]; |
---|
570 | } |
---|
571 | if(size(perm_pos) > 1) |
---|
572 | { |
---|
573 | perm = perm, perm_pos[2..size(perm_pos)]; |
---|
574 | } |
---|
575 | perm = perm[2..size(perm)]; |
---|
576 | matrix T[n][n]; |
---|
577 | matrix D[1][n]; |
---|
578 | for(i = 1; i <= n; i++) |
---|
579 | { |
---|
580 | T[1..n, i] = H[perm[i], (n+1)..(2*n)]; |
---|
581 | D[1, i] = H[perm[i], perm[i]]; |
---|
582 | } |
---|
583 | map phi = basering, matrix(maxideal(1))*transpose(T); |
---|
584 | f = phi(f); |
---|
585 | f = jet(f, k); |
---|
586 | |
---|
587 | /* separate the variables */ |
---|
588 | phi = basering, maxideal(1); |
---|
589 | map corank_part = basering, maxideal(1); |
---|
590 | for(i = c+1; i <= n; i++) |
---|
591 | { |
---|
592 | corank_part[i] = 0; |
---|
593 | } |
---|
594 | poly h = f-jet(f, 2)-corank_part(f); |
---|
595 | poly hi; |
---|
596 | while(h != 0) |
---|
597 | { |
---|
598 | for(i = c+1; i <= n; i++) |
---|
599 | { |
---|
600 | hi = h/var(i); |
---|
601 | phi[i] = var(i)-hi/(2*D[1, i]); |
---|
602 | h = h-hi*var(i); |
---|
603 | } |
---|
604 | f = phi(f); |
---|
605 | f = jet(f, k); |
---|
606 | h = f-jet(f, 2)-corank_part(f); |
---|
607 | } |
---|
608 | poly g = cleardenom(f-jet(f, 2)); |
---|
609 | |
---|
610 | return(list(c, lambda, k, g)); |
---|
611 | } |
---|
612 | example |
---|
613 | { |
---|
614 | "EXAMPLE:"; |
---|
615 | echo = 2; |
---|
616 | ring r = 0, (x,y,z), ds; |
---|
617 | poly f = (x2+3y-2z)^2+xyz-(x-y3+x2z3)^3; |
---|
618 | realmorsesplit(f); |
---|
619 | } |
---|
620 | |
---|
621 | /////////////////////////////////////////////////////////////////////////////// |
---|
622 | /* |
---|
623 | symmetric Gauss algorithm |
---|
624 | |
---|
625 | If A is not a square matrix, then the largest upper or left submatrix |
---|
626 | is assumed to be symmetric. |
---|
627 | */ |
---|
628 | proc sym_reduce(matrix A) |
---|
629 | { |
---|
630 | int r = nrows(A); |
---|
631 | int c = ncols(A); |
---|
632 | int n = r; |
---|
633 | if(n > c) |
---|
634 | { |
---|
635 | n = c; |
---|
636 | } |
---|
637 | poly q; |
---|
638 | int i, j; |
---|
639 | for(i = 1; i <= n; i++) |
---|
640 | { |
---|
641 | for(j = i+1; j <= n; j++) |
---|
642 | { |
---|
643 | if(A[i, j] != 0) |
---|
644 | { |
---|
645 | while(A[i, i] == 0) |
---|
646 | { |
---|
647 | A[1..r, i] = A[1..r, i]+A[1..r, j]; |
---|
648 | A[i, 1..c] = A[i, 1..c]+A[j, 1..c]; |
---|
649 | } |
---|
650 | q = A[i, j]/A[i, i]; |
---|
651 | A[1..r, j] = A[1..r, j]-q*A[1..r, i]; |
---|
652 | A[j, 1..c] = A[j, 1..c]-q*A[i, 1..c]; |
---|
653 | } |
---|
654 | } |
---|
655 | } |
---|
656 | return(A); |
---|
657 | } |
---|
658 | |
---|
659 | /////////////////////////////////////////////////////////////////////////////// |
---|
660 | /* |
---|
661 | - apply jet(f, k) |
---|
662 | - rewrite f as f = a*var(i)^2+p*var(i)+r with |
---|
663 | var(i)-free p and r |
---|
664 | */ |
---|
665 | static proc rewriteformorsesplit(poly f, int k, int i) |
---|
666 | { |
---|
667 | f = jet(f, k); |
---|
668 | matrix C = coeffs(f, var(i)); |
---|
669 | poly r = C[1,1]; |
---|
670 | poly p = C[2,1]; |
---|
671 | poly a = (f-r-p*var(i))/var(i)^2; |
---|
672 | return(f, a, p, r); |
---|
673 | } |
---|
674 | |
---|
675 | /////////////////////////////////////////////////////////////////////////////// |
---|
676 | proc milnornumber(poly f) |
---|
677 | " |
---|
678 | USAGE: milnornumber(f); f poly |
---|
679 | RETURN: Milnor number of f, or -1 if the Milnor number is not finite |
---|
680 | KEYWORDS: Milnor number |
---|
681 | NOTE: The monomial order must be local. |
---|
682 | EXAMPLE: example milnornumber; shows an example" |
---|
683 | { |
---|
684 | /* error check */ |
---|
685 | int i; |
---|
686 | for(i = nvars(basering); i > 0; i--) |
---|
687 | { |
---|
688 | if(var(i) > 1) |
---|
689 | { |
---|
690 | ERROR("The monomial order must be local."); |
---|
691 | } |
---|
692 | } |
---|
693 | |
---|
694 | return(vdim(std(jacob(f)))); |
---|
695 | } |
---|
696 | example |
---|
697 | { |
---|
698 | "EXAMPLE:"; |
---|
699 | echo = 2; |
---|
700 | ring r = 0, (x,y), ds; |
---|
701 | poly f = x3+y4; |
---|
702 | milnornumber(f); |
---|
703 | } |
---|
704 | |
---|
705 | /////////////////////////////////////////////////////////////////////////////// |
---|
706 | proc determinacy(poly f, list #) |
---|
707 | " |
---|
708 | USAGE: determinacy(f[, mu]); f poly, mu int |
---|
709 | RETURN: an upper bound for the determinacy of f |
---|
710 | NOTE: The characteristic of the basering must be zero, the monomial order |
---|
711 | must be local, f must be contained in maxideal(1) and the Milnor |
---|
712 | number of f must be finite. |
---|
713 | @* The Milnor number of f can be provided as an optional parameter in |
---|
714 | order to avoid that it is computed again. |
---|
715 | SEE ALSO: milnornumber, highcorner |
---|
716 | KEYWORDS: Determinacy |
---|
717 | EXAMPLE: example determinacy; shows an example" |
---|
718 | { |
---|
719 | /* auxiliary variables */ |
---|
720 | int i; |
---|
721 | |
---|
722 | /* error check */ |
---|
723 | if(char(basering) != 0) |
---|
724 | { |
---|
725 | ERROR("The characteristic must be zero."); |
---|
726 | } |
---|
727 | int n = nvars(basering); |
---|
728 | for(i = 1; i <= n; i++) |
---|
729 | { |
---|
730 | if(var(i) > 1) |
---|
731 | { |
---|
732 | ERROR("The monomial order must be local."); |
---|
733 | } |
---|
734 | } |
---|
735 | if(jet(f, 0) != 0) |
---|
736 | { |
---|
737 | ERROR("The input polynomial must be contained in maxideal(1)."); |
---|
738 | } |
---|
739 | |
---|
740 | /* get Milnor number before continuing error check */ |
---|
741 | int mu; |
---|
742 | if(size(#) > 0) // read optional parameter |
---|
743 | { |
---|
744 | if(size(#) > 1 || typeof(#[1]) != "int") |
---|
745 | { |
---|
746 | ERROR("Wrong optional parameters."); |
---|
747 | } |
---|
748 | else |
---|
749 | { |
---|
750 | mu = #[1]; |
---|
751 | } |
---|
752 | } |
---|
753 | else // compute Milnor number |
---|
754 | { |
---|
755 | mu = milnornumber(f); |
---|
756 | } |
---|
757 | |
---|
758 | /* continue error check */ |
---|
759 | if(mu < 0) |
---|
760 | { |
---|
761 | ERROR("The Milnor number of the input polynomial must be"+newline |
---|
762 | +"non-negative and finite."); |
---|
763 | } |
---|
764 | |
---|
765 | int k; // an upper bound for the determinacy, |
---|
766 | // we use several methods: |
---|
767 | |
---|
768 | /* Milnor number */ |
---|
769 | k = mu+1; |
---|
770 | f = jet(f, k); |
---|
771 | |
---|
772 | /* highest corner */ |
---|
773 | int hc; |
---|
774 | for(i = 0; i < 3; i++) |
---|
775 | { |
---|
776 | f = jet(f, k); |
---|
777 | hc = deg(highcorner(std(maxideal(i)*jacob(f)))); |
---|
778 | hc = hc+2-i; |
---|
779 | if(hc < k) |
---|
780 | { |
---|
781 | k = hc; |
---|
782 | } |
---|
783 | } |
---|
784 | |
---|
785 | return(k); |
---|
786 | } |
---|
787 | example |
---|
788 | { |
---|
789 | "EXAMPLE:"; |
---|
790 | echo = 2; |
---|
791 | ring r = 0, (x,y), ds; |
---|
792 | poly f = x3+xy3; |
---|
793 | determinacy(f); |
---|
794 | } |
---|
795 | |
---|