1 | ////////////////////////////////////////////////////////////////////////////// |
---|
2 | version="version classify2.lib 4.1.2.0 Feb_2019 "; // $Id$ |
---|
3 | category="Commutative Algebra"; |
---|
4 | info=" |
---|
5 | LIBRARY: classify2.lib Classification of isolated singularities |
---|
6 | AUTHORS: Janko Boehm, email: boehm@mathematik.uni-kl.de |
---|
7 | Magdaleen Marais, email: magdaleen.marais@up.ac.za |
---|
8 | Gerhard Pfister, email: pfister@mathematik.uni-kl.de |
---|
9 | |
---|
10 | |
---|
11 | OVERVIEW: |
---|
12 | We classify isolated singularities of corank <=2 and modality <=2 with respect |
---|
13 | to right-equivalence over the complex numbers according to Arnold's list. We |
---|
14 | determine the type and, for positive modality, the parameter. |
---|
15 | |
---|
16 | V.I. Arnold has described normal forms and has developed a classifier for, in |
---|
17 | particular, all isolated hypersurface singularities over the complex numbers up |
---|
18 | to modality 2. Building on a series of 105 theorems, this classifier determines |
---|
19 | the type of the given singularity. However, for positive modality, this does not |
---|
20 | fix the right equivalence class of the singularity, since the values of the |
---|
21 | moduli parameters are not specified. |
---|
22 | |
---|
23 | This library implements an alternative classification algorithm for isolated |
---|
24 | hypersurface singularities of corank and modality up to two. For a |
---|
25 | singularity given by a polynomial over the rationals, the algorithm determines |
---|
26 | its right equivalence class by specifying a polynomial representative in Arnold's |
---|
27 | list of normal forms. In particular, the algorithm also determines values for the |
---|
28 | moduli parameters. |
---|
29 | |
---|
30 | The implementation is based on the paper |
---|
31 | |
---|
32 | Janko Boehm, Magdaleen Marais, Gerhard Pfister: A Classification Algorithm for |
---|
33 | Complex Singularities of Corank and Modality up to Two, Singularities and |
---|
34 | Computer Algebra - Festschrift for Gert-Martin Greuel on the Occasion of his |
---|
35 | 70th Birthday, Springer 2017, http://arxiv.org/abs/1604.04774, |
---|
36 | https://doi.org/10.1007/978-3-319-28829-1_2 |
---|
37 | |
---|
38 | There are functions for determining a normal form equation and for determining |
---|
39 | the complex type of the singularity. |
---|
40 | |
---|
41 | Acknowledgements: This research was supported by |
---|
42 | the Staff Exchange Bursary Programme of the University of Pretoria, DFG SPP 1489, |
---|
43 | DFG TRR 195. The financial assistance of the National Research Foundation (NRF), |
---|
44 | South Africa, towards this research is hereby acknowledged. Opinions expressed |
---|
45 | and conclusions arrived at are those of the author and are not necessarily to be |
---|
46 | attributed to the National Research Foundation, South Africa. |
---|
47 | |
---|
48 | KEYWORDS: |
---|
49 | singularities; classification |
---|
50 | |
---|
51 | SEE ALSO: realclassify_lib, classify_lib |
---|
52 | |
---|
53 | PROCEDURES: |
---|
54 | complexClassify(I); classifier returning a normal form equation |
---|
55 | complexType(I); classifier returning the type and modality |
---|
56 | |
---|
57 | "; |
---|
58 | |
---|
59 | |
---|
60 | LIB "polyclass.lib"; |
---|
61 | //LIB "sing.lib"; |
---|
62 | LIB "absfact.lib"; |
---|
63 | LIB "primdec.lib"; |
---|
64 | |
---|
65 | static proc mod_init() |
---|
66 | { |
---|
67 | LIB "gfanlib.so"; |
---|
68 | } |
---|
69 | |
---|
70 | |
---|
71 | |
---|
72 | static proc idealcontainment(ideal I,ideal J){ |
---|
73 | int sz=size(reduce(I,std(J))); |
---|
74 | return(sz==0);} |
---|
75 | |
---|
76 | static proc exponentvector(poly f){ |
---|
77 | return(leadexp(f));} |
---|
78 | |
---|
79 | static proc weighteddeg(poly f,intvec w){ |
---|
80 | return(deg(f,w));} |
---|
81 | |
---|
82 | static proc piecewisedegree(poly f, list pw){ |
---|
83 | int d1=weighteddeg(f,pw[1]); |
---|
84 | int d2; |
---|
85 | for (int i=2;i<=size(pw);i++){ |
---|
86 | d2=weighteddeg(f,pw[i]); |
---|
87 | if (d2<d1){d1=d2;} |
---|
88 | } |
---|
89 | return(d1);} |
---|
90 | |
---|
91 | //list NP = newtonPolygon(x^2*y^2+x^5+y^7); |
---|
92 | static proc piecewiseWeightOfPolygon(list NP){ |
---|
93 | if (size(NP)==1){return(list(NP[1][1]));} |
---|
94 | if (size(NP)>2){ERROR("Not implemented for more than two faces";)} |
---|
95 | intvec vt = setintersection(NP[1][2],NP[2][2])[1]; |
---|
96 | int d1=weighteddeg(var(1)^vt[1]*var(2)^vt[2],NP[1][1]); |
---|
97 | int d2=weighteddeg(var(1)^vt[1]*var(2)^vt[2],NP[2][1]); |
---|
98 | int c1 = lcm(d1,d2) div d1; |
---|
99 | int c2 = lcm(d1,d2) div d2; |
---|
100 | return(list(c1*NP[1][1],c2*NP[2][1])); |
---|
101 | } |
---|
102 | //piecewiseWeightOfPolygon(NP); |
---|
103 | |
---|
104 | |
---|
105 | static proc piecewiseJet(poly f, list pw, int d){ |
---|
106 | if (f==0){return(f);} |
---|
107 | list mon = terms(f); |
---|
108 | poly ff; |
---|
109 | int wt; |
---|
110 | for (int i=1;i<=size(mon);i++){ |
---|
111 | wt = piecewisedegree(mon[i],pw); |
---|
112 | if (wt==d){ff=ff+mon[i];} |
---|
113 | } |
---|
114 | return(ff);} |
---|
115 | |
---|
116 | |
---|
117 | static proc monomials(poly f){ |
---|
118 | list L; |
---|
119 | poly m; |
---|
120 | poly p = f; |
---|
121 | int i = 1; |
---|
122 | while (p<>0){ |
---|
123 | L[i]=leadmonom(p); |
---|
124 | p=p-lead(p); |
---|
125 | i=i+1; |
---|
126 | } |
---|
127 | return(L);} |
---|
128 | |
---|
129 | static proc terms(poly f){ |
---|
130 | list L; |
---|
131 | poly m; |
---|
132 | poly p = f; |
---|
133 | int i = 1; |
---|
134 | while (p<>0){ |
---|
135 | L[i]=lead(p); |
---|
136 | p=p-lead(p); |
---|
137 | i=i+1; |
---|
138 | } |
---|
139 | return(L);} |
---|
140 | |
---|
141 | |
---|
142 | |
---|
143 | static proc lowestMonomials(poly f, intvec w){ |
---|
144 | if (f==0){return(list(f));} |
---|
145 | list mon = monomials(f); |
---|
146 | int wt; |
---|
147 | int mindg = weighteddeg(mon[1],w); |
---|
148 | for (int i=2;i<=size(mon);i++){ |
---|
149 | wt = weighteddeg(mon[i],w); |
---|
150 | if (wt<mindg){mindg=wt;} |
---|
151 | } |
---|
152 | list L; |
---|
153 | int j=1; |
---|
154 | for (i=1;i<=size(mon);i++){ |
---|
155 | if (weighteddeg(mon[i],w)==mindg){ |
---|
156 | L[j]=mon[i]; |
---|
157 | j=j+1; |
---|
158 | } |
---|
159 | } |
---|
160 | return(L);} |
---|
161 | |
---|
162 | |
---|
163 | |
---|
164 | static proc lowestTerms(poly f, intvec w){ |
---|
165 | if (f==0){return(list(f));} |
---|
166 | list mon = terms(f); |
---|
167 | int wt; |
---|
168 | int mindg = weighteddeg(mon[1],w); |
---|
169 | for (int i=2;i<=size(mon);i++){ |
---|
170 | wt = weighteddeg(mon[i],w); |
---|
171 | if (wt<mindg){mindg=wt;} |
---|
172 | } |
---|
173 | list L; |
---|
174 | int j=1; |
---|
175 | for (i=1;i<=size(mon);i++){ |
---|
176 | if (weighteddeg(mon[i],w)==mindg){ |
---|
177 | L[j]=mon[i]; |
---|
178 | j=j+1; |
---|
179 | } |
---|
180 | } |
---|
181 | return(L);} |
---|
182 | |
---|
183 | static proc lowestJet(poly f, intvec w){ |
---|
184 | if (f==0){return(f);} |
---|
185 | list mon = terms(f); |
---|
186 | int wt; |
---|
187 | int mindg = weighteddeg(mon[1],w); |
---|
188 | for (int i=2;i<=size(mon);i++){ |
---|
189 | wt = weighteddeg(mon[i],w); |
---|
190 | if (wt<mindg){mindg=wt;} |
---|
191 | } |
---|
192 | return(jet(f,mindg,w));} |
---|
193 | |
---|
194 | |
---|
195 | static proc removeMon(poly f,poly f0,intvec u1,intvec u2,poly t) { |
---|
196 | if (t==0){ERROR("term should not be zero");} |
---|
197 | if (printlevel>1){ |
---|
198 | "Algorithm 4"; |
---|
199 | "f0="+string(f0); |
---|
200 | "u1="+string(u1); |
---|
201 | "u2="+string(u2); |
---|
202 | "t="+string(t); |
---|
203 | } |
---|
204 | poly dfx=diff(f0,var(1)); |
---|
205 | poly dfy=diff(f0,var(2)); |
---|
206 | if (printlevel>1){ |
---|
207 | "dfx = "+string(dfx); |
---|
208 | "dfy = "+string(dfy); |
---|
209 | } |
---|
210 | poly mx=lowestJet(dfx,u2); |
---|
211 | mx=lowestJet(mx,u1); |
---|
212 | poly my=lowestJet(dfy,u1); |
---|
213 | my=lowestJet(my,u2); |
---|
214 | if (printlevel>1){ |
---|
215 | "mx = "+string(mx); |
---|
216 | "my = "+string(my); |
---|
217 | } |
---|
218 | if (NF(t,mx)==0) { |
---|
219 | f=subst(f,x,x-t/mx); |
---|
220 | f=truncateAtHighestCorner(f); |
---|
221 | return(f); |
---|
222 | } |
---|
223 | if (NF(t,my)==0) { |
---|
224 | f=subst(f,y,y-t/my); |
---|
225 | f=truncateAtHighestCorner(f); |
---|
226 | return(f); |
---|
227 | } |
---|
228 | if (printlevel>0){"did nothing";} |
---|
229 | return(f); |
---|
230 | } |
---|
231 | |
---|
232 | static proc poshull(list L){ |
---|
233 | intmat A[size(L)][size(L[1])]=L[1..size(L)]; |
---|
234 | return(coneViaPoints(A)); |
---|
235 | } |
---|
236 | |
---|
237 | static proc facetsAsCones(cone c){ |
---|
238 | intmat R= intmat(rays(c)); |
---|
239 | intmat Fc=intmat(facets(c)); |
---|
240 | intmat M = intmat(facets(c)*transpose(R)); |
---|
241 | intmat r[1][ncols(R)]; |
---|
242 | int j; |
---|
243 | int l=1; |
---|
244 | list F; |
---|
245 | for (int i=1;i<=nrows(M);i++){ |
---|
246 | list L; |
---|
247 | l=1; |
---|
248 | //"Face";i; |
---|
249 | for (j=1;j<=ncols(M);j++){ |
---|
250 | if (M[i,j]==0){ |
---|
251 | //"ray";j; |
---|
252 | r=R[j,1..ncols(R)]; |
---|
253 | L[l]=intvec(r); |
---|
254 | l=l+1; |
---|
255 | } |
---|
256 | } |
---|
257 | F[i]=list(intvec(Fc[i,1..ncols(Fc)]),L);kill L; |
---|
258 | } |
---|
259 | return(F);} |
---|
260 | |
---|
261 | static proc newtonPolygon(poly f){ |
---|
262 | //"f="+string(f); |
---|
263 | list M=monomials(f); |
---|
264 | //print(M); |
---|
265 | list E; |
---|
266 | intvec v; |
---|
267 | for (int i=1;i<=size(M);i++){ |
---|
268 | v=exponentvector(M[i]); |
---|
269 | v[size(v)+1]=1; |
---|
270 | E[i]=v; |
---|
271 | } |
---|
272 | E[size(M)+1]=intvec(1,0,0); |
---|
273 | E[size(M)+2]=intvec(0,1,0); |
---|
274 | cone c = poshull(E); |
---|
275 | //print(facets(c)); |
---|
276 | list F=facetsAsCones(c); |
---|
277 | //"Facets as cones";print(F); |
---|
278 | list F1; |
---|
279 | list ry; |
---|
280 | int l=1; |
---|
281 | intvec nv; |
---|
282 | intvec ry1,ry2; |
---|
283 | list ryy; |
---|
284 | list fi; |
---|
285 | for (i=1;i<=size(F);i++){ |
---|
286 | ry=F[i][2]; |
---|
287 | nv=F[i][1]; |
---|
288 | nv=nv[1..2]; |
---|
289 | ry1=ry[1];ry1=ry1[1..2]; |
---|
290 | ry2=ry[2];ry2=ry2[1..2]; |
---|
291 | if (ry1[1]>ry2[1]) {ryy=list(ry2,ry1);} else {ryy=list(ry1,ry2);} |
---|
292 | if ((ry[1][3]<>0) and (ry[2][3]<>0)){fi[l]=ryy[1][1];F1[l]=list(nv,ryy);l=l+1;} |
---|
293 | } |
---|
294 | intvec fii=sort(fi)[2]; |
---|
295 | list F2; |
---|
296 | int j; |
---|
297 | for (i=1;i<=size(fii);i++){ |
---|
298 | F2[i]=F1[fii[i]]; |
---|
299 | } |
---|
300 | return(F2); |
---|
301 | } |
---|
302 | |
---|
303 | static proc maximumFiltration(poly f){ |
---|
304 | return(ord(f)); |
---|
305 | } |
---|
306 | |
---|
307 | static proc factorLowestJet(poly f){ |
---|
308 | poly g= jet(f,maximumFiltration(f)); |
---|
309 | list L = factorize(g); |
---|
310 | intvec fii=sort(L[2])[2]; |
---|
311 | list LL; |
---|
312 | for (int i=1;i<=size(fii);i++){ |
---|
313 | LL[size(fii)-i+1]=list(L[1][fii[i]],L[2][fii[i]]); |
---|
314 | } |
---|
315 | list LLL=LL[1..size(LL)-1]; |
---|
316 | if (size(LL)==2){ |
---|
317 | return(list(LL[1])); |
---|
318 | } |
---|
319 | return(LLL);} |
---|
320 | |
---|
321 | static proc lintrafo(poly f, ideal T){ |
---|
322 | def R=basering; |
---|
323 | list L=ringlist(R); |
---|
324 | L[2]=L[2]+list("xx","yy"); |
---|
325 | L[3][1][2]=intvec(1,1,1,1); |
---|
326 | L[3][1][1]="lp"; |
---|
327 | def RR=ring(L); |
---|
328 | setring RR; |
---|
329 | ideal T=fetch(R,T); |
---|
330 | ideal I = xx-T[1],yy-T[2]; |
---|
331 | option(redSB); |
---|
332 | ideal Is = std(I); |
---|
333 | poly fx = -subst(Is[2],x,0)/(coeff(Is[2],x)); |
---|
334 | poly fy = -subst(Is[1],y,0)/(coeff(Is[1],y)); |
---|
335 | poly ff = subst(fetch(R,f),x,fx,y,fy); |
---|
336 | list L=ringlist(RR); |
---|
337 | L[2]=list("xx","yy"); |
---|
338 | L[3][1][2]=intvec(1,1); |
---|
339 | def R1=ring(L); |
---|
340 | setring R1; |
---|
341 | poly ff = imap(RR,ff); |
---|
342 | setring R; |
---|
343 | return(fetch(R1,ff)); |
---|
344 | } |
---|
345 | |
---|
346 | static proc coeff(poly f, poly m){ |
---|
347 | poly p = f; |
---|
348 | int i = 1; |
---|
349 | while (p<>0){ |
---|
350 | if (m==leadmonom(p)){return(leadcoef(p));} |
---|
351 | p=p-lead(p); |
---|
352 | } |
---|
353 | return(0);} |
---|
354 | |
---|
355 | |
---|
356 | static proc coeff2(poly f, poly m){ |
---|
357 | poly radm=(radical(m))[1]; |
---|
358 | matrix M =coef(f,radm); |
---|
359 | for (int i=1;i<=ncols(M);i++){ |
---|
360 | if (M[1,i]==m){return(M[2,i]);} |
---|
361 | } |
---|
362 | return(0);} |
---|
363 | |
---|
364 | |
---|
365 | |
---|
366 | static proc doLinearTransformation(Poly F){ |
---|
367 | setring(F.in); |
---|
368 | poly f = F.value; |
---|
369 | poly g= jet(f,maximumFiltration(f)); |
---|
370 | def R= basering; |
---|
371 | def S=absFactorize(g); |
---|
372 | setring S; |
---|
373 | list L = absolute_factors; |
---|
374 | if (printlevel>0){ |
---|
375 | "absolute factors"; |
---|
376 | L; |
---|
377 | } |
---|
378 | number mp; |
---|
379 | for (int i=1;i<=size(L[3]);i++){ |
---|
380 | if (L[3][i]<>a){ |
---|
381 | mp = leadcoef(L[3][i]); |
---|
382 | } |
---|
383 | } |
---|
384 | if (mp==0){ |
---|
385 | setring R; |
---|
386 | } else { |
---|
387 | if (printlevel>0){"minimal polynomial "+string(mp);} |
---|
388 | minpoly = mp; |
---|
389 | poly g=imap(R,g); |
---|
390 | poly f=imap(R,f); |
---|
391 | } |
---|
392 | list L = factorize(g); |
---|
393 | int alllinear=1; |
---|
394 | for (i =1;i<=size(L[1]);i++){ |
---|
395 | if (deg(L[1][i])>1){alllinear=0;} |
---|
396 | } |
---|
397 | if (alllinear==0){ |
---|
398 | return(F); |
---|
399 | } |
---|
400 | intvec fii=sort(L[2])[2]; |
---|
401 | list LL; |
---|
402 | for (i=1;i<=size(fii);i++){ |
---|
403 | LL[size(fii)-i+1]=list(L[1][fii[i]],L[2][fii[i]]); |
---|
404 | } |
---|
405 | list LLL=LL[1..size(LL)-1]; |
---|
406 | // bugfix for inconsistent behaviour of Singular for list assignment |
---|
407 | if (size(LL)==2){ |
---|
408 | LLL=list(LL[1]); |
---|
409 | } |
---|
410 | L=LLL; |
---|
411 | //list L=factorLowestJet(f); |
---|
412 | poly ff; |
---|
413 | if (size(L)==1) { |
---|
414 | if (NF(L[1][1],y)==0) { |
---|
415 | if (printlevel>0){ "case1 ";} |
---|
416 | ff=lintrafo(f,ideal(y,x)); |
---|
417 | Poly FF=makePoly(ff); |
---|
418 | return(FF); |
---|
419 | } else { |
---|
420 | if (printlevel>0){"case2 ";} |
---|
421 | ff=lintrafo(f,ideal(L[1][1],y)); |
---|
422 | Poly FF=makePoly(ff); |
---|
423 | return(FF); |
---|
424 | } |
---|
425 | } |
---|
426 | if (size(L)==2) { |
---|
427 | if (printlevel>0){"case3 ";} |
---|
428 | ff=lintrafo(f,ideal(L[1][1],L[2][1])); |
---|
429 | Poly FF=makePoly(ff); |
---|
430 | return(FF); |
---|
431 | } |
---|
432 | if ((size(L)==3) and (jet(f,3)==0)) { |
---|
433 | if (printlevel>0){"case4 ";} |
---|
434 | setring R; |
---|
435 | poly ff; |
---|
436 | list L = factorLowestJet(g); |
---|
437 | if (NF(L[1][1],y)==0) |
---|
438 | { |
---|
439 | ff=lintrafo(f,ideal(y,x)); |
---|
440 | //Poly FF=makePoly(ff); |
---|
441 | //return(FF); |
---|
442 | } |
---|
443 | else |
---|
444 | { |
---|
445 | ff=lintrafo(f,ideal(L[1][1],y)); |
---|
446 | //Poly FF=makePoly(ff); |
---|
447 | // return(FF); |
---|
448 | } |
---|
449 | number a1=coeff(ff,x^3*y); |
---|
450 | number a2=coeff(ff,x^2*y^2); |
---|
451 | ff=lintrafo(ff,ideal(x,y+(a1/(2*a2))*x)); |
---|
452 | Poly FF=makePoly(ff); |
---|
453 | return(FF); |
---|
454 | } |
---|
455 | return(F);} |
---|
456 | |
---|
457 | static proc isOnPolygon(intvec v,list NP){ |
---|
458 | int tst; |
---|
459 | cone c1,c2; |
---|
460 | v[size(v)+1]=1; |
---|
461 | int j; |
---|
462 | intvec w; |
---|
463 | for (int i=1;i<=size(NP);i++){ |
---|
464 | list L; |
---|
465 | for (j=1;j<=size(NP[i][2]);j++){ |
---|
466 | w=NP[i][2][j]; |
---|
467 | w[size(w)+1]=1; |
---|
468 | L[j]=w; |
---|
469 | } |
---|
470 | c1=poshull(L); |
---|
471 | L[size(L)+1]=v; |
---|
472 | c2=poshull(L); |
---|
473 | if (c1==c2){return(1);} |
---|
474 | kill L; |
---|
475 | } |
---|
476 | return(0);} |
---|
477 | |
---|
478 | static proc latticePoints(list NP){ |
---|
479 | int maxx,maxy; |
---|
480 | for (int i=1;i<=size(NP);i++){ |
---|
481 | if (NP[i][2][1][1]>maxx) {maxx=NP[i][2][1][1];} |
---|
482 | if (NP[i][2][2][1]>maxx) {maxx=NP[i][2][2][1];} |
---|
483 | if (NP[i][2][1][2]>maxy) {maxy=NP[i][2][1][2];} |
---|
484 | if (NP[i][2][2][2]>maxy) {maxy=NP[i][2][2][2];} |
---|
485 | } |
---|
486 | int x,y; |
---|
487 | list LP; |
---|
488 | for (x=0;x<=maxx;x++){ |
---|
489 | for (y=0;y<=maxy;y++){ |
---|
490 | if (isOnPolygon(intvec(x,y),NP)){LP[size(LP)+1]=intvec(x,y);} |
---|
491 | } |
---|
492 | } |
---|
493 | return(LP);} |
---|
494 | |
---|
495 | |
---|
496 | static proc allMonomialsInPiecewiseWeightedDegree(list wt,int d){ |
---|
497 | list L; |
---|
498 | poly mon; |
---|
499 | int ix,iy; |
---|
500 | for (ix=0;ix<=d;ix++){ |
---|
501 | for (iy=0;iy<=d;iy++){ |
---|
502 | mon=var(1)^ix*var(2)^iy; |
---|
503 | if (piecewisedegree(mon,wt)==d){L[size(L)+1]=mon;} |
---|
504 | } |
---|
505 | } |
---|
506 | return(L); |
---|
507 | } |
---|
508 | //allMonomialsInPiecewiseWeightedDegree(list(intvec(3,2)),8); |
---|
509 | |
---|
510 | |
---|
511 | |
---|
512 | static proc latticeToMonomials(list L){ |
---|
513 | list M; |
---|
514 | for (int i=1;i<=size(L);i++){ |
---|
515 | M[size(M)+1]=var(1)^(L[i][1])*var(2)^(L[i][2]); |
---|
516 | } |
---|
517 | return(M);} |
---|
518 | |
---|
519 | static proc termsOnPolygon(poly f1,NP){ |
---|
520 | list HB=latticePoints(NP); |
---|
521 | poly g; |
---|
522 | for (int j=1;j<=size(HB);j++){ |
---|
523 | g=g+coeff(f1,x^HB[j][1]*y^HB[j][2])*x^HB[j][1]*y^HB[j][2]; |
---|
524 | } |
---|
525 | return(g); |
---|
526 | } |
---|
527 | |
---|
528 | static proc factorWeightedLowestJet(poly f,list F){ |
---|
529 | poly g= termsOnPolygon(f,list(F)); |
---|
530 | list L=factorize(g); |
---|
531 | intvec fii=sort(L[2])[2]; |
---|
532 | list LL; |
---|
533 | for (int i=1;i<=size(fii);i++){ |
---|
534 | LL[size(fii)-i+1]=list(L[1][fii[i]],L[2][fii[i]]); |
---|
535 | } |
---|
536 | list LLL=LL[1..size(LL)-1]; |
---|
537 | if (size(LL)==2){ |
---|
538 | return(list(LL[1])); |
---|
539 | } |
---|
540 | return(LLL);} |
---|
541 | |
---|
542 | |
---|
543 | |
---|
544 | static proc doWeightedLinearTransformation(poly f1,list NF1){ |
---|
545 | if (printlevel>0){"++++++++++++++++++ "+string(deg(f1));} |
---|
546 | list fac=factorWeightedLowestJet(f1,NF1); |
---|
547 | if (printlevel>0){print("Fac "+string(fac));} |
---|
548 | //subst(fac[1][1],y,0); |
---|
549 | if (deg(subst(fac[1][1],var(2),0))>1) {return(0);} |
---|
550 | poly rs=subst(fac[1][1],var(1),0); |
---|
551 | f1=subst(f1,var(1),var(1)/(coeff(fac[1][1],var(1)))-rs/(coeff(fac[1][1],var(1)))); |
---|
552 | f1=truncateAtHighestCorner(f1); |
---|
553 | return(f1); |
---|
554 | } |
---|
555 | |
---|
556 | |
---|
557 | |
---|
558 | static proc corner(list NP){ |
---|
559 | return(NP[1][2][2]); |
---|
560 | } |
---|
561 | |
---|
562 | |
---|
563 | |
---|
564 | static proc weightsOfNewtonPolygon(list NP){ |
---|
565 | list wt; |
---|
566 | for (int i =1;i<=size(NP);i++){ |
---|
567 | wt[i]=NP[i][1]; |
---|
568 | } |
---|
569 | return(wt);} |
---|
570 | |
---|
571 | |
---|
572 | static proc twofaceelimination(poly f,intvec co){ |
---|
573 | int mu=milnor(f); |
---|
574 | list NP=newtonPolygon(f); |
---|
575 | list L=latticePoints(NP); |
---|
576 | //intvec co=corner(NP); |
---|
577 | list L1; |
---|
578 | int l=1; |
---|
579 | for (int j=1;j<=size(L);j++){ |
---|
580 | if ((L[j][1]==co[1]-1) or (L[j][2]==co[2]-1)) { |
---|
581 | L1[l]=L[j]; |
---|
582 | l=l+1; |
---|
583 | } |
---|
584 | } |
---|
585 | poly f0=termsOnPolygon(f,NP); |
---|
586 | poly f1=f; |
---|
587 | intvec co1; |
---|
588 | number c1; |
---|
589 | poly mon; |
---|
590 | list wt; |
---|
591 | list NPadj=adjecentFacesOfNewtonPolygon(NP,co); |
---|
592 | for (j=1;j<=size(L1);j++){ |
---|
593 | co1=L1[j]; |
---|
594 | mon=x^co1[1]*y^co1[2]; |
---|
595 | c1=coeff(f0,mon); |
---|
596 | mon=c1*x^co1[1]*y^co1[2]; |
---|
597 | //"removing "+string(mon); |
---|
598 | if (c1<>0){ |
---|
599 | wt=weightsOfNewtonPolygon(NPadj); |
---|
600 | f1=removeMon(f,f0,wt[1],wt[2],mon); |
---|
601 | //f1=jet(f1,mu+1); |
---|
602 | } |
---|
603 | } |
---|
604 | return(f1); |
---|
605 | } |
---|
606 | |
---|
607 | static proc adjecentFacesOfNewtonPolygon(list NP,intvec co){ |
---|
608 | list NP1; |
---|
609 | for (int i=1;i<=size(NP);i++){ |
---|
610 | if (isVertexOfFace(NP[i][2],co)){ |
---|
611 | NP1[size(NP1)+1]=NP[i]; |
---|
612 | } |
---|
613 | } |
---|
614 | return(NP1);} |
---|
615 | |
---|
616 | static proc isVertexOfFace(list F,intvec co){ |
---|
617 | if ((F[1]==co) or (F[2]==co)){return(1);} |
---|
618 | return(0);} |
---|
619 | |
---|
620 | static proc dotwofaceelimination(poly f,intvec si){ |
---|
621 | poly fnew=twofaceelimination(f,si); |
---|
622 | if (fnew==f){ |
---|
623 | return(f); |
---|
624 | } else { |
---|
625 | return(dotwofaceelimination(fnew, si)); |
---|
626 | }; |
---|
627 | }; |
---|
628 | |
---|
629 | |
---|
630 | static proc innerVertices(NP){ |
---|
631 | list L; |
---|
632 | int j,jj; |
---|
633 | intvec co; |
---|
634 | for (j=1;j<=size(NP)-1;j++){ |
---|
635 | L[j]=NP[j][2][2]; |
---|
636 | } |
---|
637 | return(L); |
---|
638 | } |
---|
639 | |
---|
640 | static proc setunion(list A, list B){ |
---|
641 | list L=A; |
---|
642 | int i,j,l; |
---|
643 | l=1; |
---|
644 | int tst; |
---|
645 | for (i=1;i<=size(B);i++){ |
---|
646 | tst=0; |
---|
647 | for (j=1;j<=size(L);j++){ |
---|
648 | if (L[j]==B[i]){tst=1;} |
---|
649 | } |
---|
650 | if (tst==0){L[size(L)+1]=B[i];} |
---|
651 | } |
---|
652 | return(L);} |
---|
653 | |
---|
654 | |
---|
655 | static proc setintersection(list A, list B){ |
---|
656 | list L; |
---|
657 | if ((size(A)==0) or (size(B)==0)){return(L);} |
---|
658 | int i,j,l; |
---|
659 | l=1; |
---|
660 | for (i=1;i<=size(A);i++){ |
---|
661 | for (j=1;j<=size(B);j++){ |
---|
662 | if (B[j]==A[i]){L[l]=A[i];l=l+1;} |
---|
663 | } |
---|
664 | } |
---|
665 | return(L);} |
---|
666 | |
---|
667 | static proc setminus(list A, list B){ |
---|
668 | list L; |
---|
669 | if (size(A)==0){return(A);} |
---|
670 | int i,j,l; |
---|
671 | l=1; |
---|
672 | int tst; |
---|
673 | for (i=1;i<=size(A);i++){ |
---|
674 | tst=0; |
---|
675 | for (j=1;j<=size(B);j++){ |
---|
676 | if (B[j]==A[i]){tst=1;break;} |
---|
677 | } |
---|
678 | if (tst==0){L[size(L)+1]=A[i];} |
---|
679 | } |
---|
680 | return(L);} |
---|
681 | |
---|
682 | |
---|
683 | |
---|
684 | static proc vectorfieldToTransformation(list vf, intvec wt, int n){ |
---|
685 | list L; |
---|
686 | poly sumf=x; |
---|
687 | poly f=x; |
---|
688 | int currentorder=0; |
---|
689 | number fac =1; |
---|
690 | int i=1; |
---|
691 | while (currentorder<n){ |
---|
692 | f=(vf[1]*diff(f,x)+vf[2]*diff(f,y)); |
---|
693 | sumf=sumf+1/fac*f; |
---|
694 | i=i+1; |
---|
695 | fac=fac*i; |
---|
696 | currentorder=deg(f,wt); |
---|
697 | } |
---|
698 | L[1]=sumf; |
---|
699 | sumf=y; |
---|
700 | f=y; |
---|
701 | currentorder=0; |
---|
702 | fac=1; |
---|
703 | i=1; |
---|
704 | while (currentorder<n){ |
---|
705 | f=(vf[1]*diff(f,x)+vf[2]*diff(f,y)); |
---|
706 | sumf=sumf+1/fac*f; |
---|
707 | i=i+1; |
---|
708 | fac=fac*i; |
---|
709 | currentorder=deg(f,wt); |
---|
710 | } |
---|
711 | L[2]=sumf; |
---|
712 | return(L); |
---|
713 | } |
---|
714 | |
---|
715 | //vectorfieldToTransformation(list(3*y^2*c,-2*x*c),intvec(3,2),5); |
---|
716 | |
---|
717 | |
---|
718 | static proc removeMonW(poly f,poly f0,intvec u1,intvec u2,poly t,list T) { |
---|
719 | if (t==0){ERROR("term should not be zero");} |
---|
720 | if (printlevel>0){"Algorithm 7";} |
---|
721 | //"f0="+string(f0); |
---|
722 | //"u1="+string(u1); |
---|
723 | //"u2="+string(u2); |
---|
724 | if (printlevel>0){"t="+string(t);} |
---|
725 | list NP = newtonPolygon(T[2]); |
---|
726 | list wtNP = piecewiseWeightOfPolygon(NP); |
---|
727 | int currentweighteddeg=piecewisedegree(t,wtNP); |
---|
728 | poly oldf = f; |
---|
729 | poly dfx=diff(f0,var(1)); |
---|
730 | poly dfy=diff(f0,var(2)); |
---|
731 | ideal jac=jacob(f0); |
---|
732 | if (idealcontainment(ideal(t),jac)){return(f);} |
---|
733 | matrix co; |
---|
734 | int mu=milnor(f); |
---|
735 | f=removeMon(f,f0,intvec(1,1),intvec(1,1),t); |
---|
736 | //f=jet(f,mu+1); |
---|
737 | if (printlevel>0){">>> after transformation by Algorithm 4 "+string(piecewiseJet(f,wtNP,currentweighteddeg));} |
---|
738 | if (f!=oldf) {return(f);} |
---|
739 | int isparameterdeg=0; |
---|
740 | for (int i=1; i<=size(T[3]);i++){ |
---|
741 | if (currentweighteddeg==piecewisedegree(T[3][i],wtNP)){isparameterdeg=1;break;} |
---|
742 | } |
---|
743 | if ((f==oldf) and isparameterdeg) {return(f);} |
---|
744 | int dprime = piecewisedegree(T[3][1],wtNP); |
---|
745 | if (printlevel>0){"weighted degree of lowest parameter monomial "+string(dprime);} |
---|
746 | poly m = T[3][1]; |
---|
747 | list fac = factorize(f0); |
---|
748 | if (fac[2][2]!=2){ERROR("there should be a quadratic factor");} |
---|
749 | poly g0 = fac[1][2]; |
---|
750 | if (printlevel>0){"g0= "+string(g0);} |
---|
751 | list philist = vectorfieldToTransformation(list(diff(g0,y),-diff(g0,x)),intvec(3,2),5); |
---|
752 | if (printlevel>0){"automorphism from vectorfield "+string(philist);} |
---|
753 | poly s = coeff(f,m)*m; |
---|
754 | phi = basering, philist[1],philist[2]; |
---|
755 | poly tprime = lowestJet(phi(s)-s,intvec(3,2)); |
---|
756 | if (printlevel>0){"tprime = "+string(tprime);} |
---|
757 | tprime = -f0+doEliminationByDegree(tprime+f0,f0,u1,u2,wtNP,currentweighteddeg,mu); |
---|
758 | tprime = lowestJet(tprime,intvec(3,2)); |
---|
759 | if (printlevel>0){"tprime after elimination = "+string(tprime);} |
---|
760 | //matrix cc = lift(t,ideal(tprime)); |
---|
761 | //"lift "+string(cc); |
---|
762 | //poly c = -cc[1,1]; |
---|
763 | poly c = -t/tprime; |
---|
764 | if (printlevel>0){"c= "+string(c);} |
---|
765 | if (c==0){ERROR("c should not be zero");} |
---|
766 | philist = vectorfieldToTransformation(list(c*diff(g0,y),-c*diff(g0,x)),intvec(3,2),5); |
---|
767 | phi = basering, philist[1],philist[2]; |
---|
768 | if (printlevel>0){"f0 = "+string(f0);} |
---|
769 | //"before vector field coefficients of x^2*y^5 and y^8 "+string(coeff(f,x^2*y^5))+", "+string(coeff(f,y^8)); |
---|
770 | poly ff0=f-f0; |
---|
771 | poly h =f0+phi(ff0); |
---|
772 | //h=jet(h,mu+1); |
---|
773 | h=truncateAtHighestCorner(h); |
---|
774 | //"after vector field coefficients of x^2*y^5 and y^8 "+string(coeff(h,x^2*y^5))+", "+string(coeff(h,y^8)); |
---|
775 | h=doEliminationByDegree(h,f0,u1,u2,wtNP,currentweighteddeg,mu); |
---|
776 | //"after elimination coefficients of x^2*y^5 and y^8 "+string(coeff(h,x^2*y^5))+", "+string(coeff(h,y^8)); |
---|
777 | if (printlevel>0){"h= "+string(h);} |
---|
778 | return(h);} |
---|
779 | |
---|
780 | static proc doEliminationByDegree(poly f, poly f0, intvec u1, intvec u2, list wtNP, int d,int mu){ |
---|
781 | //list Lelim=listOfMonomialsInPiecewiseWeightedDegree(f,wtNP,d); |
---|
782 | list Lelim=allMonomialsInPiecewiseWeightedDegree(wtNP,d); |
---|
783 | if (printlevel>0){"monomials to be considered "+string(Lelim);} |
---|
784 | poly me; |
---|
785 | for (int j=1;j<=size(Lelim);j++){ |
---|
786 | me=Lelim[size(Lelim)-j+1]; |
---|
787 | if (printlevel>0){"eliminating "+string(me)+" with coefficient "+string(coeff(f,me));} |
---|
788 | if (coeff(f,me)!=0){ |
---|
789 | if (printlevel>0){"before elimination "+string(listOfMonomialsInPiecewiseWeightedDegree(f,wtNP,d));} |
---|
790 | f=removeMon(f,f0,u1,u2,coeff(f,me)*me); |
---|
791 | //f=jet(f,mu+1); |
---|
792 | if (printlevel>0){"after elimination "+string(listOfMonomialsInPiecewiseWeightedDegree(f,wtNP,d));} |
---|
793 | } |
---|
794 | } |
---|
795 | return(f);} |
---|
796 | |
---|
797 | |
---|
798 | static proc rescaling(Poly G){ |
---|
799 | //"rescaling to be implemented"; |
---|
800 | return(G);} |
---|
801 | |
---|
802 | static proc eliminateonandabove(Poly fl,list T){ |
---|
803 | def R= basering; |
---|
804 | setring(fl.in); |
---|
805 | poly f = fl.value; |
---|
806 | if (printlevel>0){"on and above";} |
---|
807 | if (printlevel>0){"current polynomial = "+string(f);} |
---|
808 | if (printlevel>0){ |
---|
809 | "Type:"; |
---|
810 | T; |
---|
811 | } |
---|
812 | list NP = newtonPolygon(f); |
---|
813 | if (printlevel>0){ |
---|
814 | "Newton polygon = "; |
---|
815 | NP; |
---|
816 | } |
---|
817 | if (printlevel>0){"on Newton polygon = "+string(termsOnPolygon(f,NP));} |
---|
818 | list NPNF = newtonPolygon(T[2]); |
---|
819 | if (printlevel>0){ |
---|
820 | "Newton polygon normal form = "; |
---|
821 | NPNF; |
---|
822 | } |
---|
823 | list monfNPNF=monomials(termsOnPolygon(f,NPNF)); |
---|
824 | list allmonNPNF=latticeToMonomials(latticePoints(NPNF)); |
---|
825 | list nonNFmon=setminus(setminus(allmonNPNF,T[3]),monomials(T[2])); |
---|
826 | list tobeelim=setminus(setminus(monfNPNF,T[3]),monomials(T[2])); |
---|
827 | if (printlevel>0){ |
---|
828 | if (size(nonNFmon)>0) {"Monomials which can be eliminated "+string(nonNFmon);} |
---|
829 | if (size(tobeelim)>0) {"Monomials which have to be eliminated "+string(tobeelim);} |
---|
830 | } |
---|
831 | if (size(tobeelim)>1){ |
---|
832 | if (printlevel>0){"to be implemented";} |
---|
833 | if (T[1]!="X_9"){ERROR("should be X9");} |
---|
834 | Poly F=makePoly(f); |
---|
835 | return(F); |
---|
836 | } |
---|
837 | if (size(tobeelim)==1){ |
---|
838 | tobeelim[1]; |
---|
839 | f=eliminateMonomialOnDiagonal(f,tobeelim[1],NPNF); |
---|
840 | f=eliminateMonomialsAboveNewtonPolygon(f,T); |
---|
841 | Poly F=makePoly(f); |
---|
842 | return(F); |
---|
843 | } |
---|
844 | f=eliminateMonomialsAboveNewtonPolygon(f,T); |
---|
845 | Poly F=makePoly(f); |
---|
846 | return(F); |
---|
847 | } |
---|
848 | |
---|
849 | static proc eliminateMonomialsAboveNewtonPolygon(poly f,list T) |
---|
850 | { |
---|
851 | int Whash; |
---|
852 | if ((T[1][1]=="W") and (T[1][2]=="#")){Whash=1;} |
---|
853 | if (printlevel>0){"eliminate above the newton polygon";} |
---|
854 | list NP = newtonPolygon(T[2]); |
---|
855 | list wtNP = piecewiseWeightOfPolygon(NP); |
---|
856 | if (printlevel>0){"piecewise weight "+string(wtNP);} |
---|
857 | int dprime = piecewisedegree(T[3][size(T[3])],wtNP); |
---|
858 | if (printlevel>0){"weighted degree of highest parameter monomial "+string(dprime);} |
---|
859 | int degdiag = piecewisedegree(terms(T[2])[1],wtNP); |
---|
860 | if (printlevel>0){"weighted degree of diagonal "+string(degdiag);} |
---|
861 | int j; |
---|
862 | poly f0; |
---|
863 | int w1,w2; |
---|
864 | poly me; |
---|
865 | matrix co; |
---|
866 | ideal jac; |
---|
867 | map phi; |
---|
868 | int mu = milnor(f); |
---|
869 | f0=termsOnPolygon(f,NP); |
---|
870 | if (printlevel>0){"On Newton polygon of normal form "+string(f0);} |
---|
871 | int dprimelow = piecewisedegree(T[3][1],wtNP); |
---|
872 | if (printlevel>0){"weighted degree of lowest parameter monomial "+string(dprimelow);} |
---|
873 | for (int d =degdiag+1;d<=dprime;d++){ |
---|
874 | if (printlevel>0){"Considering degree "+string(d);} |
---|
875 | //list Lelim=listOfMonomialsInPiecewiseWeightedDegree(f,wtNP,d); |
---|
876 | list Lelim=allMonomialsInPiecewiseWeightedDegree(wtNP,d); |
---|
877 | if (printlevel>0){"monomials to be eliminated "+string(Lelim);} |
---|
878 | for (j=1;j<=size(Lelim);j++){ |
---|
879 | me=Lelim[size(Lelim)-j+1]; |
---|
880 | if (printlevel>0){"eliminating "+string(me)+" with coefficient "+string(coeff(f,me));} |
---|
881 | if (coeff(f,me)!=0){ |
---|
882 | if (size(NP)==1){ |
---|
883 | if (Whash==1) { |
---|
884 | f=removeMonW(f,f0,intvec(1,1),intvec(1,1),coeff(f,me)*me,T); |
---|
885 | //f=jet(f,mu+1); |
---|
886 | } else { |
---|
887 | f=removeMon(f,f0,intvec(1,1),intvec(1,1),coeff(f,me)*me); |
---|
888 | } |
---|
889 | } |
---|
890 | if (size(NP)==2){ |
---|
891 | jac=jacob(f0); |
---|
892 | if (idealcontainment(ideal(me),jac)){ |
---|
893 | f=removeMon(f,f0,wtNP[2],wtNP[1],coeff(f,me)*me); |
---|
894 | //f=jet(f,mu+1); |
---|
895 | } |
---|
896 | } |
---|
897 | } |
---|
898 | } |
---|
899 | //"afterwards "+string(piecewiseJet(f,wtNP,d)); |
---|
900 | kill Lelim; |
---|
901 | me=piecewiseJet(f,wtNP,d); |
---|
902 | if (printlevel>0){"still to be eliminated or written in Arnold's basis "+string(me);} |
---|
903 | f0=termsOnPolygon(f,NP); |
---|
904 | jac = jacob(f0); |
---|
905 | for (j=1; j<=size(T[3]);j++){ |
---|
906 | if (piecewisedegree(T[3][j],wtNP)==d){jac=jac+ideal(T[3][j]);} |
---|
907 | } |
---|
908 | co=lift(jac,me); |
---|
909 | if (printlevel>0){"coefficients of lift "+string(co);} |
---|
910 | phi = basering, x-co[1,1],y-co[2,1]; |
---|
911 | f=phi(f); |
---|
912 | //f=jet(f,mu+1); |
---|
913 | f=truncateAtHighestCorner(f); |
---|
914 | if (printlevel>0){">>> after handling current degree "+string(piecewiseJet(f,wtNP,d));} |
---|
915 | } |
---|
916 | poly fjet; |
---|
917 | for (d =degdiag;d<=dprime;d++){ |
---|
918 | fjet=fjet+piecewiseJet(f,wtNP,d); |
---|
919 | } |
---|
920 | if (printlevel>0){">>>>>>>>>>>>>>>>>>>>>>>>>>> after elimination above "+string(fjet);} |
---|
921 | return(fjet);} |
---|
922 | |
---|
923 | static proc listOfMonomialsInPiecewiseWeightedDegree(poly f, list wtNP,int d){ |
---|
924 | poly tobeelim=piecewiseJet(f,wtNP,d); |
---|
925 | //list NP = newtonPolygon(f); |
---|
926 | //"On Newton polygon "+string(termsOnPolygon(f,NP)); |
---|
927 | //"Considering jet "+string(tobeelim); |
---|
928 | if (tobeelim==0){return(list());} |
---|
929 | def R=basering; |
---|
930 | list L=ringlist(R); |
---|
931 | L[3][1][1]="dp"; |
---|
932 | def RR=ring(L); |
---|
933 | kill L; |
---|
934 | setring RR; |
---|
935 | poly tobeelim=fetch(R,tobeelim); |
---|
936 | list Lelim=monomials(tobeelim); |
---|
937 | setring R; |
---|
938 | list Lelim=fetch(RR,Lelim); |
---|
939 | return(Lelim); |
---|
940 | } |
---|
941 | |
---|
942 | static proc eliminateMonomialOnDiagonal(poly f,poly mon,list NPNF){ |
---|
943 | if (printlevel>0){"weighted linear transformation on diagonal";} |
---|
944 | def R=basering; |
---|
945 | list L=ringlist(R); |
---|
946 | L[2]=L[2]+list("aa"); |
---|
947 | L[3][1][2]=intvec(1,1,1); |
---|
948 | L[3][1][1]="lp"; |
---|
949 | def RR=ring(L); |
---|
950 | setring RR; |
---|
951 | poly f=fetch(R,f); |
---|
952 | int wx=NPNF[1][1][1]; |
---|
953 | int wy=NPNF[1][1][2]; |
---|
954 | poly ff = subst(f,x,x+aa*y^(wx div wy)); |
---|
955 | poly mon = fetch(R,mon); |
---|
956 | poly c=coeff2(ff,mon); |
---|
957 | ring T =0,(aa),dp; |
---|
958 | poly c=imap(RR,c); |
---|
959 | if (printlevel>0){"must be zero "+string(c);} |
---|
960 | def S=absFactorize(c); |
---|
961 | setring S; |
---|
962 | list L = absolute_factors; |
---|
963 | if (printlevel>0){"possible transformations "+string(L);} |
---|
964 | poly equ; |
---|
965 | for (int i=2;i<=size(L[1]);i++){if (deg(L[1][i])==1){equ=L[1][i];break;}} |
---|
966 | if (equ==0){ERROR("There should be a transformation on the diagonal defined over the current field");} |
---|
967 | number a1=number(subst(equ,aa,0)); |
---|
968 | number a2=number(diff(equ,aa)); |
---|
969 | number valueforaa=-a1/a2; |
---|
970 | setring RR; |
---|
971 | number valueforaa=number(imap(S,valueforaa)); |
---|
972 | if (printlevel>0){"value "+string(valueforaa);} |
---|
973 | ff=subst(ff,aa,valueforaa); |
---|
974 | //ff=truncateAtHighestCorner(ff); |
---|
975 | if (printlevel>0){"after transformation "+string(monomials(termsOnPolygon(ff,NPNF)));} |
---|
976 | setring R; |
---|
977 | return(fetch(RR,ff)); |
---|
978 | } |
---|
979 | |
---|
980 | static proc determinecase(poly f){ |
---|
981 | if (printlevel>0){"Start determinecase";} |
---|
982 | int t=timer; |
---|
983 | int mu=milnor(f); |
---|
984 | if (printlevel>0){ |
---|
985 | "f = "+string(f); |
---|
986 | "Milnor number "+string(mu); |
---|
987 | } |
---|
988 | list case; |
---|
989 | list NP=newtonPolygon(f); |
---|
990 | list NP1; |
---|
991 | //list lowestordinaryjet1; |
---|
992 | //list lowestordinaryjet=lowestJet(f,intvec(1,1)); |
---|
993 | list cases=possibleTypesOfModalityAndCorankUpToTwo(mu); |
---|
994 | int tt=timer; |
---|
995 | //int ttt,tttt,ttt1,tttt1; |
---|
996 | //int lll; |
---|
997 | for (int i=1;i<=size(cases);i++){ |
---|
998 | NP1=newtonPolygon(cases[i][2]); |
---|
999 | //lowestordinaryjet1=lowestJet(cases[i][2],intvec(1,1)); |
---|
1000 | //"-->"+string(compareNormalsOfNewtonPolygons(NP1,NP)); |
---|
1001 | //ttt=timer; |
---|
1002 | //lll=compareNormalsOfNewtonPolygons(NP1,NP); |
---|
1003 | //ttt=timer-ttt; |
---|
1004 | //ttt1=ttt1+ttt; |
---|
1005 | //tttt=timer; |
---|
1006 | //lll=(termsOnPolygon(f,NP1)==termsOnPolygon(f,NP)); |
---|
1007 | //tttt=timer-tttt; |
---|
1008 | //tttt1=tttt1+tttt; |
---|
1009 | //and (termsOnPolygon(f,NP1)==termsOnPolygon(f,NP)) |
---|
1010 | if ((compareNormalsOfNewtonPolygons(NP1,NP))) {case[size(case)+1]=cases[i];} |
---|
1011 | } |
---|
1012 | if (size(case)==0){return(list());} |
---|
1013 | if (size(case)>1){ERROR("Several possible Newton polygons (should not happen) "+string(case));} |
---|
1014 | if (printlevel>0){"Time for determinecase "+string(timer-t);} |
---|
1015 | //+" for compare normals "+string(ttt1)+" for terms "+string(tttt1); |
---|
1016 | return(case[1]); |
---|
1017 | } |
---|
1018 | |
---|
1019 | static proc compareLists(list NP1,list NP2){ |
---|
1020 | if (size(NP1)!=size(NP2)){return(0);} |
---|
1021 | for (int i=1;i<=size(NP1);i++){ |
---|
1022 | if (not (NP1[i]==NP2[i])){return(0);} |
---|
1023 | } |
---|
1024 | return(1);} |
---|
1025 | |
---|
1026 | |
---|
1027 | static proc compareNormalsOfNewtonPolygons(list NP1,list NP2){ |
---|
1028 | if (size(NP1)!=size(NP2)){return(0);} |
---|
1029 | for (int i=1;i<=size(NP1);i++){ |
---|
1030 | if (not (NP1[i][1]==NP2[i][1])){return(0);} |
---|
1031 | } |
---|
1032 | return(1);} |
---|
1033 | |
---|
1034 | static proc possibleTypesOfModalityAndCorankUpToTwo(int mu){ |
---|
1035 | list cases; |
---|
1036 | list individuals = list("E[6]",x^3+y^4,list()), |
---|
1037 | list("E[7]",x^3+x*y^3,list()), |
---|
1038 | list("E[8]",x^3+y^5,list()), |
---|
1039 | list("X[9]",x^4+x^2*y^2+y^4,list(x^2*y^2)), |
---|
1040 | list("J[10]",x^3+x^2*y^2+y^6,list(x^2*y^2)), |
---|
1041 | list("E[12]",x^3+y^7,list(x*y^5)), |
---|
1042 | list("E[13]",x^3+x*y^5,list(y^8)), |
---|
1043 | list("E[14]",x^3+y^8,list(x*y^6)), |
---|
1044 | list("E[18]",x^3+y^10,list(x*y^7,x*y^8)), |
---|
1045 | list("E[19]",x^3+x*y^7,list(y^11,y^12)), |
---|
1046 | list("E[20]",x^3+y^11,list(x*y^8,x*y^9)), |
---|
1047 | list("Z[11]",x^3*y+y^5,list(x*y^4)), |
---|
1048 | list("Z[12]",x^3*y+x*y^4,list(x^2*y^3)), |
---|
1049 | list("Z[13]",x^3*y+y^6,list(x*y^5)), |
---|
1050 | list("Z[17]",x^3*y+y^8,list(x*y^6,x*y^7)), |
---|
1051 | list("Z[18]",x^3*y+x*y^6,list(y^9,y^10)), |
---|
1052 | list("Z[19]",x^3*y+y^9,list(x*y^7,x*y^8)), |
---|
1053 | list("W[12]",x^4+y^5,list(x^2*y^3)), |
---|
1054 | list("W[13]",x^4+x*y^4,list(y^6)), |
---|
1055 | list("W[17]",x^4+x*y^5,list(y^7,y^8)), |
---|
1056 | list("W[18]",x^4+y^7,list(x^2*y^4,x^2*y^5)), |
---|
1057 | list("J[3,0]",x^3+x^2*y^3+y^9,list(x^2*y^3,x*y^7)), |
---|
1058 | list("Z[1,0]",x^3*y+x^2*y^3+y^7,list(x^2*y^3,x*y^6)), |
---|
1059 | list("W[1,0]",x^4+x^2*y^3+y^6,list(x^2*y^3,x^2*y^4)); |
---|
1060 | // |
---|
1061 | for (int i=1;i<=size(individuals);i++){ |
---|
1062 | if (mu==milnor(individuals[i][2])){cases[size(cases)+1]=individuals[i];} |
---|
1063 | } |
---|
1064 | // |
---|
1065 | list families; |
---|
1066 | |
---|
1067 | if (mu>=1){families[size(families)+1]=list("A["+string(mu)+"]", x^2+y^(mu+1),list())}; |
---|
1068 | |
---|
1069 | if (mu>=4){families[size(families)+1]=list("D["+string(mu)+"]", x^2*y+y^(mu-1),list())}; |
---|
1070 | |
---|
1071 | if (mu>=11){families[size(families)+1] = list("J["+string(mu)+"]", x^3+x^2*y^2+y^(mu-4), list(y^(mu-4))); |
---|
1072 | for (i=1;i<mu-9;i++){families[size(families)+1] = list("Y["+string(4+i)+","+string(4+mu-9-i)+"]", x^(4+i)+x^2*y^2+y^(4+(mu-9-i)),list(x^2*y^2));}; |
---|
1073 | } |
---|
1074 | // |
---|
1075 | if (mu>=10){families[size(families)+1]=list("X["+string(mu)+"]", x^4+x^2*y^2+y^(mu-5),list(y^(mu-5)))}; |
---|
1076 | // |
---|
1077 | if (mu>=17){families[size(families)+1]=list("J[3,"+string(mu-16)+"]",x^3+x^2*y^3+y^(mu-7) ,list(y^(mu-7),y^(mu-6)))}; |
---|
1078 | // |
---|
1079 | if(mu>=16){families[size(families)+1]=list("Z[1,"+string(mu-15)+"]",x^3*y+x^2*y^3+y^(mu-8),list(y^(mu-8),y^(mu-7))); |
---|
1080 | // |
---|
1081 | families[size(families)+1]=list("W[1,"+string(mu-15)+"]",x^4+x^2*y^3+y^(mu-9) ,list(y^(mu-9),y^(mu-8))); |
---|
1082 | // |
---|
1083 | if (mu mod 2 == 0) {families[size(families)+1]=list("W#[1,"+string(mu-15)+"]",(x^2+y^3)^2+x*y^((mu div 2)-3)+x*y^((mu div 2)-2) ,list(x*y^((mu div 2)-3),x*y^((mu div 2)-2)));} |
---|
1084 | // |
---|
1085 | if (mu mod 2 == 1) {families[size(families)+1]=list("W#[1,"+string(mu-15)+"]",(x^2+y^3)^2+x^2*y^((mu-9) div 2)+x^2*y^((mu-7) div 2) ,list(x^2*y^((mu-9) div 2),x^2*y^((mu-7) div 2)));} |
---|
1086 | }; |
---|
1087 | // |
---|
1088 | cases=cases+families; |
---|
1089 | return(cases);} |
---|
1090 | |
---|
1091 | |
---|
1092 | |
---|
1093 | static proc findFaces(poly f,list NP,list S0) |
---|
1094 | { |
---|
1095 | if (printlevel>0){"S0 = "+string(S0);} |
---|
1096 | list facescontainingS0; |
---|
1097 | list mNPi; |
---|
1098 | for (int i=1;i<=size(NP);i++) |
---|
1099 | { |
---|
1100 | mNPi=monomials(termsOnPolygon(f,list(NP[i]))); |
---|
1101 | if (printlevel>0){ |
---|
1102 | "monomials on face "+string(mNPi); |
---|
1103 | "intersection "+string(setintersection(mNPi,S0)); |
---|
1104 | } |
---|
1105 | if (size(S0)==size(setintersection(mNPi,S0))) |
---|
1106 | {facescontainingS0[size(facescontainingS0)+1]=NP[i];} |
---|
1107 | } |
---|
1108 | return(facescontainingS0);} |
---|
1109 | |
---|
1110 | |
---|
1111 | |
---|
1112 | static proc highestCornerDeterminacyBound(poly f){ |
---|
1113 | ideal I2 = maxideal(2)*jacob(f); |
---|
1114 | ideal I1 = maxideal(1)*jacob(f); |
---|
1115 | ideal I0 = maxideal(0)*jacob(f); |
---|
1116 | int d2 = deg(highcorner(std(I2))); |
---|
1117 | int d1 = 1+deg(highcorner(std(I1))); |
---|
1118 | int d0 = 2+deg(highcorner(std(I0))); |
---|
1119 | //"-------->"+string(d0)+", "+string(d1)+", "+string(d2); |
---|
1120 | return(min(d0,min(d1,d2))); |
---|
1121 | } |
---|
1122 | |
---|
1123 | static proc truncateAtHighestCorner(def f){ |
---|
1124 | //int t; |
---|
1125 | if (typeof(f)=="Poly") { |
---|
1126 | //t = timer; |
---|
1127 | ideal I2 = maxideal(2)*jacob(f.value); |
---|
1128 | int d2 = deg(highcorner(std(I2))); |
---|
1129 | //int d2=milnor(f.value)+1; |
---|
1130 | //t=timer-t; |
---|
1131 | //"truncation "+string(deg(f.value))+" -> "+string(d2)+" "+string(t); |
---|
1132 | Poly F = jet(f.value,d2); |
---|
1133 | return(F); |
---|
1134 | } else { |
---|
1135 | //t = timer; |
---|
1136 | ideal I2 = maxideal(2)*jacob(f); |
---|
1137 | int d2 = deg(highcorner(std(I2))); |
---|
1138 | //int d2=milnor(f)+1; |
---|
1139 | //t=timer-t; |
---|
1140 | //"truncation "+string(deg(f))+" -> "+string(d2)+" "+string(t); |
---|
1141 | return(jet(f,d2)); |
---|
1142 | } |
---|
1143 | } |
---|
1144 | |
---|
1145 | |
---|
1146 | |
---|
1147 | proc complexClassify(Poly fl) |
---|
1148 | "USAGE: complexClassify(fl); fl Poly |
---|
1149 | ASSUME: fl is a bivariate polynomial defining a curve singularity at (0,0) of modality <=2 and corank <=2. |
---|
1150 | The ordering of the basering must be local.@* |
---|
1151 | |
---|
1152 | RETURN: A normal form equation g for fl represented by an element of type NormalFormEquation.@* |
---|
1153 | |
---|
1154 | Note: So far the final scaling step is not implemented, so to obtain a normal form equation |
---|
1155 | it may be necessary to do a transformation of the form x->lambda1*x, y->lambda2*y@* |
---|
1156 | |
---|
1157 | Note: Case X9 is not implemented yet. In this case an error is returned stating that it is case X9. |
---|
1158 | |
---|
1159 | g.normalFormEquation.in stores the polynomial ring containing the normal form equation.@* |
---|
1160 | |
---|
1161 | g.normalFormEquation.value is the normal form equation. |
---|
1162 | |
---|
1163 | To access the normal form equation as a polynomial do:@* |
---|
1164 | |
---|
1165 | def S = g.normalFormEquation.in;@* |
---|
1166 | setring S;@* |
---|
1167 | poly nf = g.normalFormEquation.value;@* |
---|
1168 | |
---|
1169 | KEYWORDS: Arnold, classification. |
---|
1170 | |
---|
1171 | EXAMPLE: example complexClassify; shows an example |
---|
1172 | |
---|
1173 | " |
---|
1174 | |
---|
1175 | { |
---|
1176 | fl=truncateAtHighestCorner(fl); |
---|
1177 | def curring=basering; |
---|
1178 | int i,j; |
---|
1179 | int mu=milnor(fl.value); |
---|
1180 | Poly F=doLinearTransformation(fl); |
---|
1181 | def S = F.in; |
---|
1182 | setring S; |
---|
1183 | poly f = F.value; |
---|
1184 | poly f0=lowestJet(f,intvec(1,1)); |
---|
1185 | if (printlevel>0){print("f0 = "+string(f0));} |
---|
1186 | list S0 = monomials(f0); |
---|
1187 | list NP,mf0,iv; |
---|
1188 | poly f1,ff; |
---|
1189 | list si,T; |
---|
1190 | list delta; |
---|
1191 | list ef0; |
---|
1192 | list Gamma0; |
---|
1193 | list deltalist; |
---|
1194 | int foundcorner; |
---|
1195 | list si0; |
---|
1196 | while (1==1) { |
---|
1197 | NP=newtonPolygon(f);if (printlevel>0){print("Newton polygon ");print(NP);} |
---|
1198 | if (printlevel>0){print("S0 = "+string(S0));} |
---|
1199 | for (i=1;i<=size(S0);i++){ |
---|
1200 | ef0[i]=exponentvector(S0[i]); |
---|
1201 | } |
---|
1202 | iv=innerVertices(NP);if (printlevel>0){print("inner vertices "+string(iv)+" exponents "+string(ef0));} |
---|
1203 | si0=setintersection(iv,ef0); |
---|
1204 | for (j=1;j<=size(si0);j++){ |
---|
1205 | if ((si0[j][1]!=1) and (si0[j][2]!=1)){si[size(si)+1]=si0[j];} |
---|
1206 | } |
---|
1207 | if (size(si)>0) { |
---|
1208 | if (size(si)>1){ERROR("This set should contain only one element");} |
---|
1209 | if (printlevel>0){print("two face case, corner "+string(si));} |
---|
1210 | if ((si[1]!=intvec(2,2)) and (si[1]!=intvec(2,3))) {ERROR("modality is bigger than 2");} |
---|
1211 | if (printlevel>0){"f="+string(f);} |
---|
1212 | f=dotwofaceelimination(f,si[1]); |
---|
1213 | if (printlevel>0){print("after two face elimination "+string(termsOnPolygon(f,newtonPolygon(f))));} |
---|
1214 | NP=newtonPolygon(f); |
---|
1215 | if (printlevel>0){ |
---|
1216 | "Newton polygon ";NP; |
---|
1217 | "si = "+string(si); |
---|
1218 | } |
---|
1219 | Gamma0 = findFaces(f,NP,list(var(1)^si[1][1]*var(2)^si[1][2])); |
---|
1220 | if (printlevel>0){ |
---|
1221 | "Gamma0 = "+string(Gamma0); |
---|
1222 | } |
---|
1223 | f1=termsOnPolygon(f,Gamma0); |
---|
1224 | if (printlevel>0){ |
---|
1225 | "f1= "+string(f1); |
---|
1226 | } |
---|
1227 | T = determinecase(f1); |
---|
1228 | if (size(T)==0){ERROR("modality is bigger than 2");}; |
---|
1229 | if (printlevel>0){"Type "+T[1];} |
---|
1230 | Poly G = eliminateonandabove(f,T); |
---|
1231 | G=rescaling(G); |
---|
1232 | NormalFormEquation NFE; |
---|
1233 | NFE.vars = ringlist(curring)[2]; |
---|
1234 | NFE.singularityType=T[1]; |
---|
1235 | NFE.milnorNumber =mu; |
---|
1236 | NFE.normalFormEquation=G; |
---|
1237 | NFE.modality=size(T[3]); |
---|
1238 | setring G.in; |
---|
1239 | list para; |
---|
1240 | for (j =1;j<=size(T[3]);j++){para[j]=list(coeff(G.value,T[3][j])*T[3][j]);} |
---|
1241 | NFE.parameters=para; |
---|
1242 | NFE.corank=2; |
---|
1243 | NFE.determinacy=highestCornerDeterminacyBound(G.value); |
---|
1244 | NFE.realCase=0; |
---|
1245 | setring curring; |
---|
1246 | return(NFE); |
---|
1247 | } |
---|
1248 | //delta = face of NP containing S0 |
---|
1249 | deltalist = findFaces(f,NP,S0); |
---|
1250 | //if (size(deltalist)!=1){ERROR("There should be exactly one face containing S0");} |
---|
1251 | delta = deltalist[1]; |
---|
1252 | f1=termsOnPolygon(f,list(delta)); |
---|
1253 | if (printlevel>0){print("f1= ",string(f1)+" milnor "+string(milnor(f1)));} |
---|
1254 | if (milnor(f1)<>-1) { |
---|
1255 | if (printlevel>0){ |
---|
1256 | print("one face case"); |
---|
1257 | "----";f1; |
---|
1258 | newtonPolygon(f1); |
---|
1259 | } |
---|
1260 | T = determinecase(f1); |
---|
1261 | if (size(T)==0){ERROR("modality is bigger than 2");} |
---|
1262 | if (printlevel>0){"Type "+T[1];} |
---|
1263 | if (size(T[3])==0){ |
---|
1264 | Poly G=makePoly(T[2]); |
---|
1265 | NormalFormEquation NFE; |
---|
1266 | NFE.singularityType=T[1]; |
---|
1267 | NFE.milnorNumber =mu; |
---|
1268 | NFE.normalFormEquation=G; |
---|
1269 | NFE.modality=size(T[3]); |
---|
1270 | NFE.parameters=list(); |
---|
1271 | NFE.corank=2; |
---|
1272 | NFE.vars = ringlist(curring)[2]; |
---|
1273 | if (T[1][1]=="A"){NFE.corank=1;} |
---|
1274 | NFE.determinacy=highestCornerDeterminacyBound(G.value); |
---|
1275 | NFE.realCase=0; |
---|
1276 | setring curring; |
---|
1277 | return(NFE); |
---|
1278 | } |
---|
1279 | Poly G=eliminateonandabove(f,T); |
---|
1280 | G=rescaling(G); |
---|
1281 | NormalFormEquation NFE; |
---|
1282 | NFE.singularityType=T[1]; |
---|
1283 | NFE.milnorNumber =mu; |
---|
1284 | NFE.normalFormEquation=G; |
---|
1285 | NFE.modality=size(T[3]); |
---|
1286 | NFE.vars = ringlist(curring)[2]; |
---|
1287 | setring G.in; |
---|
1288 | list para; |
---|
1289 | for (j =1;j<=size(T[3]);j++){para[j]=list(coeff(G.value,T[3][j])*T[3][j]);} |
---|
1290 | NFE.parameters=para; |
---|
1291 | NFE.corank=2; |
---|
1292 | NFE.determinacy=highestCornerDeterminacyBound(G.value); |
---|
1293 | NFE.realCase=0; |
---|
1294 | setring curring; |
---|
1295 | return(NFE); |
---|
1296 | } |
---|
1297 | ff=doWeightedLinearTransformation(f,delta); |
---|
1298 | if (ff==0) { |
---|
1299 | // factor of highest multiplicity is not x-non-linear |
---|
1300 | if (not compareLists(monomials(termsOnPolygon(f,list(delta))),monomials(termsOnPolygon((x^2-y^3)^2,list(delta))))) {ERROR("modality is bigger than 2");} |
---|
1301 | mu=milnor(f); |
---|
1302 | if (mu mod 2 == 0) {T=list("W#[1,"+string(mu-15)+"]",(x^2+y^3)^2+x*y^((mu div 2)-3)+x*y^((mu div 2)-2) ,list(x*y^((mu div 2)-3),x*y^((mu div 2)-2)));} |
---|
1303 | if (mu mod 2 == 1) {T=list("W#[1,"+string(mu-15)+"]",(x^2+y^3)^2+x^2*y^((mu-9) div 2)+x^2*y^((mu-7) div 2) ,list(x^2*y^((mu-9) div 2),x^2*y^((mu-7) div 2)));} |
---|
1304 | if (printlevel>0){"Type "+T[1];} |
---|
1305 | //Whash cases |
---|
1306 | Poly G=eliminateonandabove(f,T); |
---|
1307 | G=rescaling(G); |
---|
1308 | NormalFormEquation NFE; |
---|
1309 | NFE.singularityType=T[1]; |
---|
1310 | NFE.milnorNumber =mu; |
---|
1311 | NFE.normalFormEquation=G; |
---|
1312 | NFE.modality=size(T[3]); |
---|
1313 | NFE.vars = ringlist(curring)[2]; |
---|
1314 | setring G.in; |
---|
1315 | list para; |
---|
1316 | for (j =1;j<=size(T[3]);j++){para[j]=list(coeff(G.value,T[3][j])*T[3][j]);} |
---|
1317 | NFE.parameters=para; |
---|
1318 | NFE.corank=2; |
---|
1319 | NFE.determinacy=highestCornerDeterminacyBound(G.value); |
---|
1320 | NFE.realCase=0; |
---|
1321 | setring curring; |
---|
1322 | return(NFE); |
---|
1323 | } |
---|
1324 | f=ff; |
---|
1325 | //f=jet(f,mu+1); |
---|
1326 | f=truncateAtHighestCorner(f); |
---|
1327 | f0=termsOnPolygon(f,list(delta));if (printlevel>0){print("f0= "+string(f0));} |
---|
1328 | S0=monomials(f0); |
---|
1329 | } |
---|
1330 | } |
---|
1331 | example |
---|
1332 | { "EXAMPLE:"; echo=2; |
---|
1333 | ring R = 0,(x,y),ds; |
---|
1334 | Poly f = -16065*y^9*x-5103*y^9-128*x^7-595*y^9*x^4+6*y^5-21*y^9*x^5-2187*y^7+4*x^2*y^2+12*x*y^3-5102*y^8-4480*x^6*y^2-45359*x^4*y^4-54431*x^2*y^6-19152*x^5*y^3-64258*x^3*y^5-25515*y^7*x-14490*y^9*x^2-4830*y^9*x^3-27208*x*y^8-77490*x^3*y^6-25872*x^5*y^4+20*y^3*x^3+4*x^4*y^2-15116*x^4*y^3+8*x^3*y^2-6384*x^6*y^3+y^6+28*y^3*x^2-6048*x^5*y^2+28*y^4*x+30*y^4*x^2-1344*x^6*y-2835*y^10-51555*x^3*y^7-4235*y^8*x^4-17185*x^4*y^7-22668*x^3*y^4-84*y^5*x^7-5040*y^4*x^6-560*y^3*x^7-280*y^4*x^7-19320*y^5*x^5-2380*y^5*x^6-672*x^6*y^6-40880*x^4*y^6-8610*x^5*y^6-2289*x^5*y^7-38717*y^8*x^2-20384*y^8*x^3+x^8*y^8-x^7*y^7-14*x^7*y^6-105*x^6*y^7-280*x^5*y^8+21*x^6*y^8+8*x^7*y^8+14*y^5*x-20402*y^5*x^2-10204*y^6*x-5670*y^10*x-3234*y^10*x^2-630*y^10*x^3-35*y^10*x^4-189*y^12-945*y^11-y^14-1197*y^11*x-21*y^13-399*y^11*x^2-140*y^12*x-35*y^11*x^3-21*y^12*x^2-7*y^13*x-61803*y^7*x^2-57960*y^5*x^4-448*x^7*y+9*y^4-672*x^7*y^2; |
---|
1335 | NormalFormEquation F = complexClassify(f); |
---|
1336 | F; |
---|
1337 | def S = F.normalFormEquation.in; |
---|
1338 | setring S; |
---|
1339 | poly nf = F.normalFormEquation.value; |
---|
1340 | } |
---|
1341 | |
---|
1342 | |
---|
1343 | |
---|
1344 | |
---|
1345 | |
---|
1346 | proc complexType(Poly fl) |
---|
1347 | "USAGE: complexType(fl); fl Poly |
---|
1348 | ASSUME: fl is a bivariate polynomial defining a curve singularity at (0,0) of modality and corank <=2. |
---|
1349 | The ordering of the basering must be local.@* |
---|
1350 | |
---|
1351 | RETURN: A list with the complex type of fl and the modality.@* |
---|
1352 | |
---|
1353 | |
---|
1354 | KEYWORDS: Arnold, classification. |
---|
1355 | |
---|
1356 | EXAMPLE: example complexType; shows an example |
---|
1357 | |
---|
1358 | " |
---|
1359 | |
---|
1360 | { |
---|
1361 | fl=truncateAtHighestCorner(fl); |
---|
1362 | def curring=basering; |
---|
1363 | int i,j; |
---|
1364 | int mu=milnor(fl.value); |
---|
1365 | Poly F=doLinearTransformation(fl); |
---|
1366 | def S = F.in; |
---|
1367 | setring S; |
---|
1368 | poly f = F.value; |
---|
1369 | poly f0=lowestJet(f,intvec(1,1)); |
---|
1370 | if (printlevel>0){print("f0 = "+string(f0));} |
---|
1371 | list S0 = monomials(f0); |
---|
1372 | list NP,mf0,iv; |
---|
1373 | poly f1,ff; |
---|
1374 | list si,T; |
---|
1375 | list delta; |
---|
1376 | list ef0; |
---|
1377 | list Gamma0; |
---|
1378 | list deltalist; |
---|
1379 | int foundcorner; |
---|
1380 | list si0; |
---|
1381 | while (1==1) { |
---|
1382 | NP=newtonPolygon(f); |
---|
1383 | if (printlevel>0){print("Newton polygon ");print(NP);} |
---|
1384 | if (printlevel>0){print("S0 = "+string(S0));} |
---|
1385 | for (i=1;i<=size(S0);i++){ |
---|
1386 | ef0[i]=exponentvector(S0[i]); |
---|
1387 | } |
---|
1388 | iv=innerVertices(NP); |
---|
1389 | if (printlevel>0){print("inner vertices "+string(iv)+" exponents "+string(ef0));} |
---|
1390 | si0=setintersection(iv,ef0); |
---|
1391 | for (j=1;j<=size(si0);j++){ |
---|
1392 | if ((si0[j][1]!=1) and (si0[j][2]!=1)){si[size(si)+1]=si0[j];} |
---|
1393 | } |
---|
1394 | if (size(si)>0) { |
---|
1395 | if (size(si)>1){ERROR("This set should contain only one element");} |
---|
1396 | if (printlevel>0){print("two face case, corner "+string(si));} |
---|
1397 | if ((si[1]!=intvec(2,2)) and (si[1]!=intvec(2,3))) {ERROR("modality is bigger than 2");} |
---|
1398 | if (printlevel>0){"f="+string(f);} |
---|
1399 | f=dotwofaceelimination(f,si[1]); |
---|
1400 | if (printlevel>0){print("after two face elimination "+string(termsOnPolygon(f,newtonPolygon(f)))+" Milnor "+string(milnor(f)));} |
---|
1401 | NP=newtonPolygon(f); |
---|
1402 | if (printlevel>0){ |
---|
1403 | "Newton polygon ";NP; |
---|
1404 | "si = "+string(si); |
---|
1405 | } |
---|
1406 | Gamma0 = findFaces(f,NP,list(var(1)^si[1][1]*var(2)^si[1][2])); |
---|
1407 | if (printlevel>0){ |
---|
1408 | "Gamma0 = "+string(Gamma0); |
---|
1409 | } |
---|
1410 | f1=termsOnPolygon(f,Gamma0); |
---|
1411 | if (printlevel>0){ |
---|
1412 | "f1= "+string(f1); |
---|
1413 | } |
---|
1414 | T = determinecase(f1); |
---|
1415 | if (size(T)==0){ERROR("modality is bigger than 2");}; |
---|
1416 | if (printlevel>0){"Type "+T[1];} |
---|
1417 | return(list(T[1],size(T[3]))); |
---|
1418 | } |
---|
1419 | //delta = face of NP containing S0 |
---|
1420 | deltalist = findFaces(f,NP,S0); |
---|
1421 | //if (size(deltalist)!=1){ERROR("There should be exactly one face containing S0");} |
---|
1422 | delta = deltalist[1]; |
---|
1423 | f1=termsOnPolygon(f,list(delta)); |
---|
1424 | if (printlevel>0){print("f1= ",string(f1)+" milnor "+string(milnor(f1)));} |
---|
1425 | if (milnor(f1)<>-1) { |
---|
1426 | if (printlevel>0){ |
---|
1427 | print("one face case"); |
---|
1428 | "----";f1; |
---|
1429 | newtonPolygon(f1); |
---|
1430 | } |
---|
1431 | T = determinecase(f1); |
---|
1432 | if (size(T)==0){ERROR("modality is bigger than 2");} |
---|
1433 | if (printlevel>0){"Type "+T[1];} |
---|
1434 | if (size(T[3])==0){ |
---|
1435 | return(list(T[1],size(T[3]))); |
---|
1436 | } |
---|
1437 | return(list(T[1],size(T[3]))); |
---|
1438 | } |
---|
1439 | ff=doWeightedLinearTransformation(f,delta); |
---|
1440 | if (ff==0) { |
---|
1441 | // factor of highest multiplicity is not x-non-linear |
---|
1442 | if (not compareLists(monomials(termsOnPolygon(f,list(delta))),monomials(termsOnPolygon((x^2-y^3)^2,list(delta))))) {ERROR("modality is bigger than 2");} |
---|
1443 | mu=milnor(f); |
---|
1444 | if (mu mod 2 == 0) {T=list("W#[1,"+string(mu-15)+"]",(x^2+y^3)^2+x*y^((mu div 2)-3)+x*y^((mu div 2)-2) ,list(x*y^((mu div 2)-3),x*y^((mu div 2)-2)));} |
---|
1445 | if (mu mod 2 == 1) {T=list("W#[1,"+string(mu-15)+"]",(x^2+y^3)^2+x^2*y^((mu-9) div 2)+x^2*y^((mu-7) div 2) ,list(x^2*y^((mu-9) div 2),x^2*y^((mu-7) div 2)));} |
---|
1446 | if (printlevel>0){"Type "+T[1];} |
---|
1447 | //Whash cases |
---|
1448 | return(list(T[1],size(T[3]))); |
---|
1449 | } |
---|
1450 | f=ff; |
---|
1451 | //f=jet(f,mu+1); |
---|
1452 | f=truncateAtHighestCorner(f); |
---|
1453 | f0=termsOnPolygon(f,list(delta));if (printlevel>0){print("f0= "+string(f0));} |
---|
1454 | S0=monomials(f0); |
---|
1455 | } |
---|
1456 | } |
---|
1457 | example |
---|
1458 | { "EXAMPLE:"; echo=2; |
---|
1459 | ring R = 0,(x,y),ds; |
---|
1460 | Poly f = -16065*y^9*x-5103*y^9-128*x^7-595*y^9*x^4+6*y^5-21*y^9*x^5-2187*y^7+4*x^2*y^2+12*x*y^3-5102*y^8-4480*x^6*y^2-45359*x^4*y^4-54431*x^2*y^6-19152*x^5*y^3-64258*x^3*y^5-25515*y^7*x-14490*y^9*x^2-4830*y^9*x^3-27208*x*y^8-77490*x^3*y^6-25872*x^5*y^4+20*y^3*x^3+4*x^4*y^2-15116*x^4*y^3+8*x^3*y^2-6384*x^6*y^3+y^6+28*y^3*x^2-6048*x^5*y^2+28*y^4*x+30*y^4*x^2-1344*x^6*y-2835*y^10-51555*x^3*y^7-4235*y^8*x^4-17185*x^4*y^7-22668*x^3*y^4-84*y^5*x^7-5040*y^4*x^6-560*y^3*x^7-280*y^4*x^7-19320*y^5*x^5-2380*y^5*x^6-672*x^6*y^6-40880*x^4*y^6-8610*x^5*y^6-2289*x^5*y^7-38717*y^8*x^2-20384*y^8*x^3+x^8*y^8-x^7*y^7-14*x^7*y^6-105*x^6*y^7-280*x^5*y^8+21*x^6*y^8+8*x^7*y^8+14*y^5*x-20402*y^5*x^2-10204*y^6*x-5670*y^10*x-3234*y^10*x^2-630*y^10*x^3-35*y^10*x^4-189*y^12-945*y^11-y^14-1197*y^11*x-21*y^13-399*y^11*x^2-140*y^12*x-35*y^11*x^3-21*y^12*x^2-7*y^13*x-61803*y^7*x^2-57960*y^5*x^4-448*x^7*y+9*y^4-672*x^7*y^2; |
---|
1461 | complexType(f); |
---|
1462 | } |
---|
1463 | |
---|
1464 | |
---|
1465 | |
---|
1466 | |
---|
1467 | |
---|
1468 | |
---|
1469 | |
---|
1470 | |
---|
1471 | |
---|