source: git/Singular/LIB/arnoldclassify.lib @ a1f059

spielwiese
Last change on this file since a1f059 was a1f059, checked in by Hans Schoenemann <hannes@…>, 17 months ago
arnoldClassify_to_string
  • Property mode set to 100644
File size: 59.7 KB
RevLine 
[b41ce4]1///////////////////////////////////////////////////////////////////////////////
[a1f059]2version="version arnoldclassify.lib 4.3.1.2 Nov_2022 "; // $Id$
[b41ce4]3category="Singularities";
4info="
5LIBRARY:  arnoldClassify.lib     Arnol'd Classifier of Singularities
6AUTHOR:   Eva Maria Hemmerling,  ehemmerl@rhrk.uni-kl.de
7
8OVERVIEW:
9   A library for classifying isolated hypersurface singularities
10   from the list of V.I. Arnol'd w.r.t. right equivalence up to corank 2.
11   The method relies on Baciu's list of Milnor codes and Newton's rotating
12   ruler method to distinguish the Y- and Z- singularities.
13
14REFERENCES:
15[AVG85] Arnold, Varchenko, Gusein-Zade: Singularities of Differentiable Maps.
16Vol. 1: The classification of critical points caustics and wave fronts.
17Birkh\"auser, Boston 1985
18
19[Bac01] Corina Baciu: The classification of Hypersurface Singularities
20using the Milnor Code, Diplomarbeit, Universit\"at Kaiserslautern, 2001.
21
22[GP12] Greuel, Pfister: A SINGULAR Introduction to Commutative Algebra,
23Springer Science and Business Media, 2012
24
25[Hem17] Eva Maria Hemmerling: Algorithmic Arnol'd Classification in SINGULAR,
26Master Thesis, TU Kaiserslautern, 2017.
27
28SEE ALSO:  classify_lib, realclassify_lib
29
30KEYWORDS: singularities; Arnol'd; classification; Milnor code
31
32PROCEDURES:
33arnoldListAllSeries();   list of all singularity series up to corank 2
34arnoldShowSeries(S);     data defining a singularity series S
35arnoldNormalForm(S,#);   normalform for a singularity series S
36arnoldClassify(f);       singularity class of a power series f
37arnoldCorank(f);         corank of singularity defined by f
38arnoldDeterminacy(f,#);  upper bound for the determinacy of f
39arnoldMilnorCode(f,#);   Milnor Code of a singularity f
40arnoldMorseSplit(f,#);   result of Splitting Lemma applied to f
41";
42
43LIB "inout.lib";
44LIB "elim.lib";
45LIB "sing.lib";
46LIB "findifs.lib";
47
48///////////////////////////////////////////////////////////////////////////////
49static proc mod_init()
50{
51//* define new datastruct singclass
52newstruct("singclass", "string Restrictions, string NormalForm,
53intvec MilnorCode, int Tjurina, int Determinacy, int Milnor, int Corank,
54int Modality, int s, int r, int k, string Class, string Series ");
55
56//* define new datastruct singclass
57newstruct("singseries", " string Restrictions, string MilnorCode,
58 string MilnorNumber, string Corank, string Modality, string SpecialForm,
59 string NormalForm, string Series ");
60
61LIB "SingularityDBM.lib";
62arnold_classify_init();
63}
64
65///////////////////////////////////////////////////////////////////////////////
66proc arnoldClassify( poly fPoly )
67"USAGE:   arnoldClassify (f); f poly
[078a4a]68ASSUME:  The basering is local of characteristic 0 and f defines an
[b41ce4]69         isolated singularity from Arnol'd's list of corank at most 2.
[078a4a]70COMPUTE: singularity class with respect to right equivalence and
[b41ce4]71         invariants used in the process of classification
72RETURN:  Singularity class of f of type singclass containing
73         @* - name of singularity series as listed by arnoldListAllSeries(),
74         @* - name of singularity class,
75         @* - parameters k,r,s defining the singularity class, -1 if not used,
76         @* - modality, corank, Milnor number, determinacy,
77         @* - Tjurina number, -2 if not computed, -1 if infinite,
78         @* - Milnor code, -1 if not computed,
79         @* - normal form of the singularity series from Arnol'd's list,
80         @* - restrictions on parameters as string in SINGULAR syntax.
81EXAMPLE: example arnoldClassify; shows an example
82"
83{
84
85//* initialisation
86//* new datastructure singclass to save intrinsic information
87 singclass f;
88 f = init_newsingclass( f );
89
90//* Check the basering
91 if( check_basering() )
92 {
93  f.Class = "NoClass";
94 }
95
96//* Detect units
97 if(jet(fPoly,0) != 0 )
98   {
99   return( ERROR("The polynomial is a unit." ));
100 }
101
102//* Compute Basic Invariants: Corank, Determinacy and Milnornumber
103 ideal Jf = std(jacob( fPoly ));
104 f.Milnor = vdim( Jf );
105 if( f.Milnor < 0 )
106 { ERROR("Milnornumber of the singularity must be finite.")}
107 f.Determinacy = arnoldDeterminacy( Jf , f.Milnor);
108 f.Corank = arnoldCorank( fPoly );
109
110//* Check if Milnornumber is finite
111 if( f.Milnor < 0 ){
112  f.Class = "NoClass";
113  ERROR("Milnornumber of the singularity must be finite.")
114  return( f );
115 }
116
117//* Singularities with Milnornumber = 0 belong to A[0];
118 if( f.Milnor == 0 ){
119  f.Class = "A[0]";
120  f.Series = "A[k]";
121  f.k = 0;
122  f.r = -1;
123  f.s = -1;
124  f.Modality = 0;
125  f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
126  return( f );
127 }
128
129//* Singularities with Corank <=1 belong to A[k];
130 if( f.Corank <= 1){
131  f.Class = "A["+string(f.Milnor)+"]";
132  f.Series = "A[k]";
133  f.k = f.Milnor;
134  f.r = -1;
135  f.s = -1;
136  f.Modality = 0;
137  f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
138  return( f );
139 }
140
141//* Reduce f to f.Determinacy-jet
142  if(f.Determinacy >0 ){
143   fPoly = jet(fPoly, f.Determinacy);
144  }
145
146//* Classification of singularities of Corank 2
147 if( f.Corank == 2 ){
148   return( Classify_Corank2(f, fPoly));
149 }
150
151//* Classification of singularities of Corank 3
152 if( f.Corank == 3 ){
153  return( Classify_Corank3(f, fPoly));
154 }
155
156//* No classification for singularities of corank > 3
157 f.Class = "NoClass";
158 return( f );
159
160}
161example
162{ "EXAMPLE:"; echo=2;
163  ring r = 0,(x,y),ds;
164  int k = random(3,10);
165  poly g = x4 + x2*y^(2*k+1)+x*y^(3*k+1)+ y^(4*k +1);
166  arnoldClassify(g);
167  map phi=r,y-x^2,x+y;
168  phi(g);
169  arnoldClassify(phi(g));
170  ring C = (0,i), (x,y), ds;
171  minpoly = i2 + 1;
172  poly f =(x2+y2)^2+x5;
173  arnoldClassify(f);
174}
175
176///////////////////////////////////////////////////////////////////////////////
177static proc Classify_Corank2( singclass f, poly fPoly)
178{
179//* Classifies Singularities of Corank 2; Series D,J,E,W,X,Y,Z
180
181//* Zero 4-jet is no listed singularity
182 if( jet( fPoly, 4) == 0 ){
183  f.Class = "NoClass";
184  return(f);
185 }
186
187//* Compute MilnorCode
188 f.MilnorCode = arnoldMilnorCode( fPoly );
189
190//* Classification of the series D,J and E using MilnorCode
191//* (have Nonzero 3-jet and milnorcode of length 3)
192 if (size( f.MilnorCode ) == 3 ){
193  def g = Classify_Corank2_MC_3( f );
194  return( g );
195 }
196
197//* Classification of series W, X and class Z[k,r,s] using MilnorCode
198//* (have zero 3-jet and milnorcode of length 5)
199 if ( size( f.MilnorCode ) == 5 ){
200  def g = Classify_Corank2_MC_5( f, fPoly );
201  if( g.Class != ""){
202   return( g );
203  }
204
205//* Classification of series Y and Z except Z[k,r,s]
206//* Splitting Lemma needed
207  f.k = (f.MilnorCode[3]+1) div 2;
208   return( Classify_Z_Y( f, fPoly ));
209 }
210 ERROR("MilnorCode must have length 3 or 5.");
211}
212
213///////////////////////////////////////////////////////////////////////////////
214static proc Classify_Z_Y( singclass f, poly fPoly )
215{
216
217//* classifies singularity f defined by polynomial fPoly which is of type
218//* Z or Y with zero 3-jet and non-zero 4-jet.
219//* The idea is to transform the Newton Polygon of fPoly until we can decide
[417a505]220//* whether f is of type Z or Y. For further information abouth the algorithm
[b41ce4]221//* refer to [Hem17]
222
223//* Apply Splitting lemma if needed such that fPoly is given in two variables
224 if( nvars(basering) > f.Corank ){
225  fPoly = Morse_arnold( fPoly, f.Determinacy );
226
227  //* Reduce number of variables
[078a4a]228  def @br = basering;
229  list L = ringlist(@br);
[b41ce4]230  if( size(L[1]) == 0 )
[078a4a]231  { ring red@br = char(basering), (x(1..f.Corank)), (c,ds); }
[b41ce4]232  else{
233    number m = leadcoef(L[1][4][1]);
[078a4a]234    ring red@br = (char(basering),t), (x(1..f.Corank)), (c,ds);
235    number m = imap(@br, m);
[b41ce4]236    minpoly = m;
237  }
[078a4a]238   map MapReduce = @br, maxideal(1);
[b41ce4]239  poly fPoly = MapReduce( fPoly );
240 }
241
242//* Newton Polygon transformation, not needed if f.k = 1
243if( f.k > 1 ){
244 fPoly = Transform_Newton_Polygon( f, fPoly );
245}
246
247//* Check if quasihomogeneous part of fPoly ( w.r.t. to (k,1)) has
248//* linear factor of multiplicity 3
249 def b = Check_linear_factor_mult_3( f, fPoly );
250
251if( b==1){return ( Classify_Y( f, fPoly ));}
252else{ return ( Classify_Z( f ) );}
253}
254
255///////////////////////////////////////////////////////////////////////////////
256static proc Check_linear_factor_mult_3( singclass f, poly fPoly)
257{
258 //* set weights;
259 intvec weights = f.k,1;
260
261 //* Compute quasihomogeneous part of fPoly and set y=1
262 poly h1 = jet( fPoly, 4*weights[1], weights);
263 poly h2 = subst(h1, var(2), 1);
264 //* Compute gcd of h and derivative twice
265 h2 = gcd( h2, diff(h2,var(1)));
266 h2 = gcd( h2, diff(h2, var(1)));
267
268 //* Switch coordinates and repeat the last step
269 if( deg(h2) == 0 ){
270  h2 = subst(h1, var(1), 1);
271  //* Compute gcd of h and derivative
272  h2 = gcd( h2, diff(h2,var(2)));
273  h2 = gcd( h2, diff(h2, var(2)));
274 }
275
276 if (deg(h2) > 0){ return( 0 );}
277 else{ return( 1 ); }
278}
279
280///////////////////////////////////////////////////////////////////////////////
281static proc Transform_Newton_Polygon( singclass f, poly fPoly)
282{
283//* Check that 3-jet is zero and 4-jet is not zero
284 if( jet( fPoly, 3 ) != 0 ){ERROR( "Three-jet must be zero.");}
285 if( jet( fPoly, 4 ) == 0 ){ERROR("Four-jet must not be zero.");}
286
287 //* Make sure var(1) is contained in the 4-jet
288 if( diff( jet(fPoly,4) , var(1)) == 0 ){
289  map CoordSwitch = basering, var(2), var(1);
290  fPoly = CoordSwitch(fPoly);
291 }
292
293 //* Eliminate all monomials with weighted degree smaller than 4k (= 4*f.k)
294 return( Classify_Z_Y_EliminateMonomials( f, fPoly ));
295}
296
297///////////////////////////////////////////////////////////////////////////////
298static proc Classify_Z_Y_EliminateMonomials( singclass f, poly fPoly)
299{
300//* Eliminates all monomials with weighted degree smaller than 4k (= 4*f.k)
301
302//* Computes the weights for the first monomial to eliminate
303//* Use Newtons Rulers Method
304intvec weights = Classify_Z_Y_FindWeights( fPoly, f.k);
305
306//* if the weights are (k,1), no monomials to eliminate, return fPoly
307if (weights[1] == f.k && weights[2] ==1 ){ return (fPoly);}
308
309//* Compute quasihomogeneous part w.r.t to weights
310poly fquasi = jet( fPoly, 4*weights[1], weights);
311
312//* Check if fquasi contains monomial of the form y^l
313if ( Coeff (fquasi, var(2), var(2)^(4*weights[1] div weights[2])) == 0 )
314{
315 int counter = 1;
316
317 while ( counter < 10){
318  //* Find genereous transformation
319  int ran = random(-50,50);
320  map Phi = basering, x + ran*t*var(2)^( (d div weights[2] - 2) div 2 ),y;
321  if (Coeff (Phi(fquasi), var(2), var(2)^(d div weights[2])) != 0 )
322  { counter = 11;}
323  counter = counter+1;
324 }
325 fPoly = Phi(fPoly);
326}
327
328//* Compute the fourthzero if fquasi has one
329poly factor = Classify_Z_Y_FourthZero( fquasi );
330
331//* If fquasi has a fourth zero
332if( factor == 0 ){ ERROR("Field extension needed.")}
333
334poly inversefactor =
335number ( 1 /leadcoef(factor)) * (2*lead(factor) - factor);
336map CoordChange = basering, inversefactor, var(2);
337fPoly = CoordChange( fPoly );
338fPoly = jet( fPoly, f.Determinacy);
339
340return (Classify_Z_Y_EliminateMonomials( f, fPoly ) );
341}
342
343/////////////////////////////////////////////////////////////////////////////
344static proc Classify_Z_Y_FindWeights( poly fPoly, int k )
345{
346//* Input poly f with zero 3-jet and non-zero 4-jet with var(1) contained
347//* in 4-jet. Computes the weights for all monomials lying in the Newton
[417a505]348//* Diagram under the line defined by x^4 and x^3y^k and returns the
[b41ce4]349//* smallest one. In case there is no other monomial than x^4, return
350//* weights (k,1)
351
352//* Cut off monomials above the line x^4 and x^3y^k
353 intvec weights = k, 1;
354 fPoly = jet( fPoly, 4*k, weights );
355
356//* Check if leading term of f is x^4;
357 poly firstlead = lead(fPoly);
358 if( leadmonom(firstlead) != var(1)^4 ){ ERROR("Something went wrong");}
359
360//* Compute the the minimal weights
361 fPoly = fPoly - firstlead;
362 intvec weights2;
363 while( fPoly != 0 ){
364  weights2 = qhweight( firstlead + lead( fPoly ));
365  if (weights2[1] * weights[2] < weights2[2] * weights[1]){ weights = weights2;}
366  fPoly = fPoly - lead(fPoly);
367 }
368 return(weights);
369}
370
371///////////////////////////////////////////////////////////////////////////////
372static proc Classify_Z_Y_FourthZero( poly fPoly )
373{
374//* Assume f = (ax + by^l)^4
375//* Differentiate three times w.r.t x ( f''' = 24a^3(ax +by^l) )
376//* and take gcd of f and f'''
377
378poly f3 = diff(diff(diff( fPoly, var(1)), var(1)), var(1));
379poly factor = gcd( f3, fPoly);
380
381//* Check if the factor^4 is equal input f
382if( factor^4 == fPoly){ return(factor); }
383else{ return(0);}
384}
385
386///////////////////////////////////////////////////////////////////////////////
387static proc Classify_Y( singclass f, poly fPoly)
388{
389//* Determines the parameters r and s of the singularity f of type Y[k,r,s]
390//* Case k = 1
391intvec weights;
392  if ( f.k ==1 ){
393   //* Set an upper bound
394   int bound = f.Milnor + 1;
395
396   //* Transform Newton Polygon and compute intersection point of y axis
397   weights, fPoly = Transform_Newton_Polygon_Y( f, fPoly, bound );
398    int b1 = int(( 2*weights[1]+ 2*f.k*weights[2] )div weights[2] - 4) ;
399
400   //* Switch variables and repeat the last step
401   map CoordChange = basering, var(2), var(1);
402   weights, fPoly =
403   Transform_Newton_Polygon_Y( f, CoordChange(fPoly), bound );
404   int b2 = int(( 2*weights[1]+ 2*f.k*weights[2] )div weights[2] - 4) ;
405
406   //* Y-class
407   f.Series = "Y[1,r,s]";
408   f.s = min(b1,b2);
409   f.r =  int(f.Milnor - 9 - f.s);
410   f.Class = "Y["+string(f.k)+","+string(f.r)+","+string(f.s)+"]";
411   f.Modality = 3*f.k - 2;
[078a4a]412   f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
[b41ce4]413   return( f );
414  }
415  else{
416  //* Set an upper bound = r+s;
417   int bound = 4*f.k + (f.Milnor - 12*f.k +3);
418   weights, fPoly = Transform_Newton_Polygon_Y( f, fPoly, bound);
419   int b1 =
420   int((2*weights[1] + 2* f.k*weights[2]) div weights[2]  - 4*f.k);
421   int b2 = f.Milnor - 12*f.k +3 - b1;
422
423   //* Y-class
424   f.Series = "Y[k,r,s]";
425   f.r = max(b1,b2);
426   f.s =  min(b1,b2);
427   f.Class = "Y["+string(f.k)+","+string(f.r)+","+string(f.s)+"]";
428   f.Modality = 3*f.k - 2;
[078a4a]429   f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
[b41ce4]430   return( f );
431  }
432}
433
434///////////////////////////////////////////////////////////////////////////////
435static proc Transform_Newton_Polygon_Y( singclass f, poly fPoly, int bound )
436{
437//* Eliminates Monomials until the Newton Polygon of fPoly coincides
438//* with the Newton Polygon of the normalform
439
440//* Reduce fPoly to monomials which might get eliminated
441  def weights1 = qhweight(var(1)^2*var(2)^(2*f.k) + var(2)^bound);
442  intvec weights2 = 1,0;
443  poly fquasi = jet( jet(fPoly, 2*weights1[1]+2*f.k*weights1[2], weights1),2, weights2) ;
444
445//* Find minimal ("steepest") weights
446  def weights3 = Classify_Y_FindWeights( fquasi, f.k, weights1);
447
448//* if the weights are (k,1), no monomials to eliminate, return fPoly
449  if (weights3 == weights1 )
450  { return (ERROR("Parameters r and s must be greater 0."));}
451
452//* Compute quasihomogeneous jet w.r.t to weight
453  def d =  2*weights3[1]+2*f.k*weights3[2];
454  fquasi = jet(fPoly, d, weights3);
455
456//* Check if monomial of the form y^l is contained in fquasi
457  if ( Coeff (fquasi, var(2), var(2)^(d div weights3[2])) == 0 )
458  {
459   int counter = 1;
460
461   while ( counter < 10){
462    //* Find genereous transformation
463    int ran = random(-50,50);
464    map Phi = basering,x+ran*t*var(2)^((d div weights3[2] - 2)div 2 ),y;
465    if (Coeff (Phi(fquasi), var(2), var(2)^(d div weights3[2])) != 0 )
466    { counter = 11;}
467    counter = counter+1;
468   }
469  fPoly = Phi(fPoly);
470  }
471
472//* Compute the zero of multiplicity 2 if fquasi has one
473
474 poly factor = jet(fPoly,2*weights3[1]+2*f.k*weights3[2],weights3);
475 factor =
476 Classify_Y_TwoZero(factor);
477
478//* If fquasi has a zero with mult 2
479 if( factor == 0 ){ return( weights3, fPoly );}
480
481 poly inversefactor =number(1/leadcoef(factor))*(2*lead(factor)-factor);
482 map CoordChange = basering, inversefactor, var(2);
483 fPoly = CoordChange( fPoly );
484 fPoly = jet( fPoly, f.Determinacy);
485
486 fquasi =
487 jet( jet(fPoly, 2*weights3[1]+2*f.k*weights3[2], weights3), 2, weights2);
488 if( leadmonom(fquasi) != var(1)^2*var(2)^(2*f.k))
489 { return( weights3, fPoly ); }
490
491return (Transform_Newton_Polygon_Y( f, fPoly, bound ) );
492}
493
494/////////////////////////////////////////////////////////////////////////////
495static proc Classify_Y_FindWeights( poly f, int k, intvec weights1 )
496{
497//* Assume f is a polynomial with no polynomials lying in the Newton
[417a505]498//* Diagram above the line defined by x^2y^2 and x^2y^(bound)
[b41ce4]499//* Computes the weight of all monomials lying under the line and
500//* returns the smallest one
501
502//* Check if leading term of f is x^2y^2k;
503 poly firstlead = lead(f);
504
505 if( leadmonom(firstlead) != var(1)^2*var(2)^(2*k) )
506 { ERROR("Something went wrong");}
507
508//* Compute the the minimal weights
509 f = f - firstlead;
510 intvec weights2;
511 while( f != 0 ){
512  weights2 = qhweight( firstlead + lead( f ));
513  if (weights2[1] * weights1[2] < weights2[2] * weights1[1])
514  { weights1 = weights2;}
515  f = f - lead(f);
516 }
517 return( weights1 );
518}
519
520///////////////////////////////////////////////////////////////////////////////
521static proc Classify_Y_TwoZero( poly f )
522{
523//* Assume f = (ax + by^l)^2* x^2
524
525def factors = factorize ( f );
526
527for( int i = 1; i<= size(factors[2]); i++)
528{
529 if( factors[2][i] == 2 )
530 {
531   if ( factors[1][i] != 1 && factors[1][i] - lead(factors[1][i]) != 0)
532  { return( factors[1][i]);}
533 }
534}
535return( 0 );
536}
537
538///////////////////////////////////////////////////////////////////////////////
539static proc Classify_Corank3( singclass f, poly fPoly)
540{
541return( "Not implemented yet.")
542//* todo
543
544}
545
546
547
548///////////////////////////////////////////////////////////////////////////////
549static proc Classify_Corank2_MC_3 ( singclass f )
550{
551//* Classifies the singularities of corank 2 with non-zero 3-jet using
552//* the Milnor code.
553
554 intvec mc = f.MilnorCode;
555
556 if( mc[1] != 1 ){
557  f.Class = "NoClass";
558   return( f ); }
559
560//* Check type "D[k]"
561 if( mc[2] == 1 ){
562    f.Class = "D["+string(mc[3]+3)+"]";
563  f.Series = "D[k]";
564  f.k = mc[3]+3;
565  f.r = -1;
566  f.s = -1;
567  f.Modality = 0;
568  f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
569  return( f );
570 }
571
572//* Check type "J[k,r]"
573 if(mc[3] > mc[2]&& mc[2]>1){
574  f.k = (mc[2]+1) div 2;
575  f.r = f.Milnor - 3*mc[2]-1;
576  f.s = -1;
577  f.Class = "J["+string(f.k)+","+string(f.r)+"]";
578  f.Series = "J[k,r]";
579  f.Modality = f.k-1;
580  f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
581  return( f );
582 }
583
584 else{
585
586//* Check type "J[k,0]"
587  if( f.Milnor mod 6 == 4 ){
588           f.k = (f.Milnor + 2 ) div 6;
589   f.Class ="J[" + string(f.k) + "," + string(0) +"]";
590   f.Series = "J[k,0]";
591   f.r = -1;
592   f.s = -1;
593   f.Modality = f.k-1;
594   f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
595   return( f );
596        }
597
598//* Check type "E[6k]"
599  if(f.Milnor mod 6 == 0 ){
600   f.k = f.Milnor div 6;
601   f.r = -1;
602   f.s = -1;
603   f.Class = "E[" + string(6*f.k)+ "]";
604   f.Series = "E[6k]";
605   f.Modality = f.k - 1;
606   f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
607   return( f );
608  }
609
610//* Check type "E[6k+1]"
611  if( f.Milnor mod 6 == 1 ){
612   f.k = (f.Milnor - 1) div 6;
613   f.r = -1;
614   f.s = -1;
615   f.Class = "E[" + string(6*f.k+1) + "]";
616   f.Series = "E[6k+1]";
617   f.Modality = f.k -1;
618   f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
619   return( f );
620  }
621
622//* Check type "E[6k+2]"
623  if( f.Milnor mod 6 == 2 ){
624   f.k = (f.Milnor - 2) div 6;
625   f.r = -1;
626   f.s = -1;
627   f.Class = "E[" + string(6*f.k+2) + "]";
628   f.Series = "E[6k+2]";
629   f.Modality = f.k - 1;
630   f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
631   return( f );
632  }
633
634 }
635 f.Class = "NoClass";
636 return( f );
637}
638
639///////////////////////////////////////////////////////////////////////////////
640static proc Classify_Corank2_MC_5 ( singclass f , poly fPoly)
641{
642
643//* Classifies the singularity classes of the series W, Z and the singularity
644//* class Z[k,r,s] using the Milnor code. They are of corank 2, have
645//* zero 3-jet, non-zero 4-jet and the Milnor code is of length 5.
646
647 intvec mc = f.MilnorCode;
648 if( mc[1] != 1 || mc[2]!= 1){  f.Class = "NoClass"; return(f);}
649 if( mc[3] mod 2 == 0 )
650 {
651
652//* Check type "W#[k,2r]"
653  if( mc[3] == mc[5] && mc[3] < mc[4] ){
654   f.k = mc[3] div 2;
655   f.r = mc[4]-mc[3];
656   f.s = -1;
657   f.Class = "W#[" + string(f.k) + "," + string (2*f.r) +"]";
658   f.Series = "W#[k,2r]";
659   f.Modality = 3*f.k -1;
660   f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
661   return(f);
662  }
663
664//* Check type "W[k,r]"
665  if( mc[3]==mc[4] && mc[3]<mc[5] ){
666   f.k = mc[3] div 2;
667   f.r = mc[5] - mc[3];
668   f.s = -1;
669   f.Class = "W[" +string(f.k) +"," + string(f.r) + "]";
670   f.Series = "W[k,r]";
671   f.Modality = 3*f.k - 1;
672   f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
673   return(f);
674  }
675
676//* Check type "W#[k,2r-1]"
677  if(mc[5] < mc[3] && mc[3] < mc[4]
678  && f.Milnor mod 2 == 0
679  && 2*mc[4]+mc[5]+2> f.Milnor div 2){
680   f.k = mc[3] div 2;
681   f.r = (f.Milnor - 12*f.k-2) div 2;
682   f.s = -1;
683   f.Class = "W#[" +string(f.k) +"," +string(2*f.r-1) +"]";
684   f.Series = "W#[k,2r-1]";
685   f.Modality = 3*f.k -1;
686   f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
687   return(f);
688  }
689 }
690
691//* Check type "X[k,0]"
692 if( f.Milnor mod 12 == 9 && mc[3]>=mc[5] && mc[5]<=3*mc[3]-2*mc[4] ){
693  f.k = (f.Milnor + 3) div 12;
694  f.r = 0;
695  f.s = -1;
696  f.Class = "X[" + string(f.k) +",0]";
697  if(f.k == 1){f.Series = "X[1,0]";}
698  else{f.Series = "X[k,0]";}
699  f.Modality = 3*f.k -2;
700  f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
701  return(f);
702 }
703
704//* Check type "W[12k+1]"
705 if( f.Milnor mod 12 == 1 && mc[3]==mc[5] && mc[4]<mc[3] ){
706   f.k = (f.Milnor - 1) div 12;
707  f.r = -1;
708  f.s = -1;
709   f.Class = "W[" + string(12*f.k+1) +"]";
710  f.Series = "W[12k+1]";
711  f.Modality = 3*f.k -2;
712  f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
713   return(f);
714 }
715
716//* Check type "W[12k+5]"
717 if( f.Milnor mod 12 == 5 && mc[3]==mc[5] && mc[4]<mc[3]){
718   f.k = (f.Milnor - 5) div 12;
719  f.r = -1;
720  f.s = -1;
721   f.Class = "W[" + string(12*f.k+5) +"]";
722  f.Series = "W[12k+5]";
723  f.Modality = 3*f.k -1;
724  f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
725   return(f);
726 }
727
728//* Check type "W[12k+6]"
729 if( mc[3] >= mc[5] && f.Milnor mod 12 == 6  &&
730  (2*mc[4]+mc[5]+2) < (f.Milnor div 2) && mc[3]> mc[5]){
731   f.k = (f.Milnor-6) div 12;
732   f.r = -1;
733   f.s = -1;
734   f.Class = "W[" +string(12*f.k+6) +"]";
735   f.Series = "W[12k+6]";
736   f.Modality = 3*f.k -1;
737   f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
738   return(f);
739 }
740
741//* Check type "W[k,0]"
742 if( f.Milnor mod 12 == 3  && mc[5]<= (3*mc[3] -2*mc[4])){
743  f.k = (f.Milnor - 3) div 12;
744  f.r = -1;
745  f.s = -1;
746  f.Class = "W["+ string(f.k) +",0]";
747  f.Series = "W[k,0]";
748  f.Modality = 3*f.k -1;
749  f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
750  return(f);
751 }
752
753//* Check type "W[12k]"
754 if(f.Milnor mod 12 == 0  && mc[3]> mc[5] &&
755  2*mc[4]+ mc[5]+ 2 < (f.Milnor div 2)){
756  f.k = f.Milnor div 12;
757  f.r = -1;
758  f.s = -1;
759  f.Class = "W[" + string(12*f.k) +"]";
760  f.Series = "W[12k]";
761  f.Modality = 3*f.k -2;
762  f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
763  return(f);
764 }
765
766//* Check type "X[k,r]"
767 if( (mc[4] < mc[3] && mc[3] < mc[5] ) ||
768 (mc[3]==mc[4] && mc[3] < mc[5] ) ){
769  f.k = (2*mc[3] + mc[4] + 3) div 6;
770  f.r = mc[5] - mc[3];
771  f.s = -1;
772  f.Class = "X[" + string(f.k) +"," + string(f.r)+ "]";
773  if(f.k == 1){f.Series = "X[1,r]";}
774  else{f.Series = "X[k,r]";}
775  f.Modality = 3*f.k -2;
776  f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
777  return(f);
778 }
779
780 //* Need Tjurina number to classify Z[k,r,s]
781 //* Check type "Z[k,r,s]"
782 f.Tjurina = tjurina(fPoly);
783 if( mc[3]<mc[4] && mc[4] < mc[5] &&
784   f.Milnor > ( 3* (mc[3] + mc[4] +1)) &&
785   mc[4] mod 2 == 1 &&
786   f.Tjurina == f.Milnor - (mc[4] + mc[3]) div 2 ){
787
788   f.k = (mc[3]+1) div 2;
789   f.r = (mc[4] - mc[3]) div 2;
790   f.s = mc[5] - mc[4];
791   f.Class = "Z[" + string(f.k) +"," + string(f.r)+ ","+string(f.s)+"]";
792   if(f.k == 1){f.Series = "Z[1,r,s]";}
793   else{f.Series = "Z[k,r,s]";}
794   f.Modality = 3*f.k +f.r -2;
795   f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
796   return(f);
797 }
798
799 if( mc[3] mod 2 == 1){
800  return(f);
801 }
802 f.Class = "NoClass";
803 return(f);
804}
805
806///////////////////////////////////////////////////////////////////////////////
807static proc Classify_Z( singclass f)
808{
809//* If the input is a singularity of the series Z, this procedure determines
810//* the class and the parameters.
811
812 intvec mc = f.MilnorCode;
813 f.k = (mc[3]+1) div 2;
814
815//* Check type "Z[k,r]"
816 if (f.Milnor mod 6 == 3 && mc[5]<=mc[4]){
817  f.r = (mc[4] - mc[3]) div 2;
818  f.s = -1;
819  f.Class = "Z["+string(f.k)+","+string(f.r)+"]";
820  if(f.k == 1){ f.Series = "Z[1,r]";}
821  else {f.Series = "Z[k,r]";}
822  f.Modality = 3*f.k + f.r-2;
823  f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
824  return(f);
825 }
826
827//* Check type "Z[k,12k+6r]"
828 if (f.Milnor mod 6 ==0 && mc[4]==mc[5] && mc[3]<mc[5]){
829  f.r = (f.Milnor - 12*f.k) div 6;
830  f.s = -1;
831     f.Class = "Z["+string(f.k)+","+string(12*f.k+6*f.r)+"]";
832  if(f.k == 1){ f.Series = "Z[1,6r+12]";}
833  else {f.Series = "Z[k,12k+6r]";}
834  f.Modality = 3*f.k + f.r - 2;
835  f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
836  return(f);
837 }
838
839//* Check type "Z[k,12k+6r+1]"
840 if (f.Milnor mod 6 == 1 && mc[5] mod 2 ==1 ){
841  f.r = (f.Milnor - 12*f.k -1) div 6;
842  f.s = -1;
843     f.Class = "Z["+string(f.k)+","+string(12*f.k+6*f.r+1)+"]";
844  if(f.k == 1){ f.Series = "Z[1,6r+13]";}
845  else {f.Series = "Z[k,12k+6r+1]";}
846  f.Modality = 3*f.k +f.r - 2;
847  f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
848  return(f);
849 }
850
851//* Check type "Z[k,12k+6r-1]"
852 if (f.Milnor mod 6 ==5 && mc[4] > mc[5]  ){
853  f.r = (f.Milnor - 12*f.k +1) div 6 ;
854  f.s = -1;
855     f.Class = "Z["+string(f.k)+","+string(12*f.k+6*f.r-1)+"]";
856  if(f.k == 1){ f.Series = "Z[1,6r+11]";}
857  else {f.Series = "Z[k,12k+6r-1]";}
858  f.Modality = 3*f.k + f.r -2;
859  f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
860  return(f);
861 }
862
863 //* Check type "Z[k,r,s]"
864 if( mc[3]<mc[4] && mc[4]<mc[5] &&
865  f.Milnor>3*(mc[4]+mc[3]+1) && mc[4] mod 2 ==1){
866        f.r = (mc[4] - mc[3]) div 2;
867        f.s = mc[5] - mc[4];
868  f.Class = "Z["+string(f.k)+","+string(f.r)+","+string(f.s)+"]";
869  if(f.k == 1){ f.Series = "Z[1,r,s]";}
870  else {f.Series = "Z[k,r,s]";}
871  f.Modality = 3*f.k +f.r -2;
872  f.NormalForm, f.Restrictions = NormalFormAndRestrictionsDB(f.Series);
873  return(f);
874 }
875}
876
877///////////////////////////////////////////////////////////////////////////////
878static proc arnold_classify_init()
879{
880 // Check or create Singularity Database
881   string s;
882   link l="DBM:r Singularitylist";
883   s = read(l,"VERSION");
884   if (s == "" ) {
885     create_singularity_dbm();
886   }
887   close(l);
888   l="DBM:r Singularitylist";
889   s = read(l,"VERSION");
890}
891
892
893///////////////////////////////////////////////////////////////////////////////
894proc arnoldMorseSplit ( poly f, list # )
895"USAGE:   arnoldMorseSplit(f); f poly
896ASSUME:  base ring is local, f in maxideal(2) has isolated critical point at 0
897COMPUTE: result of Splitting Lemma applied to f
[078a4a]898RETURN:  polynomial g in maxideal(3) right equivalent to f
[b41ce4]899EXAMPLE: example arnoldMorseSplit; shows an example
900"
901{
902//* save basering
[078a4a]903 def @br = basering;
[b41ce4]904 int n = nvars(basering);
905
[417a505]906//* Set Determinacy (Needed for Splitting Lemma) whether given in # or compute
[b41ce4]907 int Determinacy;
908 if (size(#) == 0 ){
909  Determinacy = arnoldDeterminacy(std(jacob(f)));
910 }
911 else{
912  Determinacy = #[1];
913 }
914//* define ring_ext to change variable names in x(1), ..., x(n)
915//* requirement for procedure morse_arnold
916 ring ring_ext = char(basering), (x(1..n)), (c, ds);
[078a4a]917 map Phi = @br, maxideal(1);
[b41ce4]918
919
920//* Apply Splitting Lemma
921 poly f = Morse_arnold( Phi( f ), Determinacy );
922
923//* Define inverse map to map f back into original ring
[078a4a]924 setring @br;
[b41ce4]925 map Phi_inverse = ring_ext, maxideal(1);
926
927 return( Phi_inverse(f) );
928
929}
930example
931{ "EXAMPLE"; echo=2;
932   ring r=0,(x,y,z),ds;
933   export r;
934   poly f=(x2+3y-2z)^2+xyz-(x-y3+x2*z3)^3;
935   poly g=arnoldMorseSplit(f);
936   g;
937}
938
939
940///////////////////////////////////////////////////////////////////////////////
941static proc Morse_arnold( poly f , int Determinacy)
942{
943//* find g(x_1, ..., x_s) such that f ~ g + x_s+1^2 + ... + x_n^2
944//* s = Corank(f);
945
946//* initialisation
947 int n, i, j, k;
948 poly f2, imagef, Q, P;
949 ideal B1, B2;
950 map Phi, Phi2;
951
952//* save basering;
[078a4a]953 def @br = basering;
[b41ce4]954 n = nvars(basering);
955
956//* Check if Determinacy is finite
957 if (Determinacy < 0 ){
958  ERROR("Determinacy must be finite.");
959 }
960 f = jet(f, Determinacy);
961
962//* Change order of variables sorted by frequency of appearance ->
963//* improves running time
964 f = ReorderVar(f);
965
966//* B1 defines the map which renumerates the variables such that
967//* g depends on the first s variables
968 B1 = maxideal(1);
969 i = 1;  //Variable x_i
970 j = 1;  //running index for B1, coordinate transformation x_i -> x_j
971
972//* While-loop over the variables,
973//* Write f = Q x_i^2 + P x_i + R,
974//* Either Q=P=0 and g depends on x_i or g is independent from x_i
975//*  and eliminate P
976 while( i <= n )
977 {
978  f2 = jet(f,2);
979  k = i+1; //* Variable for coordinate transformation
980
981  if((f2 - subst( f2, var(i), 0 )) == 0 ){
982   //* g depends on x_i -> send x_i via B1 to x_j
983   B1[i] = var(j);
984   j = j+1;
985  }
986  else{
987   //* Write f = Q x_i^2 + P x_i + R;
988      Q = Coeff( f2, var(i), var(i)^2);
989
990   //* Check if
[417a505]991   if( Q == 0 && i == n ){ ERROR("Something went wrong"); }
[b41ce4]992
993   //* Find coordinate change Phi2 such that Q(0) != 0
[078a4a]994   Phi2 = @br, maxideal(1);  //* Identity
[b41ce4]995   while( Q == 0 && i<n && k <= n ){
996    B2 = maxideal(1);
997    B2[k]= var(k)+ var(i);
[078a4a]998    Phi2 = @br, B2;
[b41ce4]999    imagef = Phi2( f );
1000    Q = Coeff( jet(imagef,2), var(i), var(i)^2);
1001    k = k+1;
1002   }
1003   f = Phi2(f);
1004   f = jet(f, Determinacy);
1005
1006   P = Coeff( f , var(i), var(i) );
1007   //* Apply coordinate changes until P = 0
1008   while( P != 0 ){
1009    //* Raise ord(P) by completing the square until > determinacy
1010     P = P / number(2*Q);
1011     B2 = maxideal(1);
1012     B2[i] = var(i) - P;
[078a4a]1013     Phi2 = @br, B2;
[b41ce4]1014     f = Phi2(f);
1015     f = jet(f, Determinacy);
1016     P = Coeff( f, var(i), var(i));
1017   }
1018  B1[i] = 0;
1019  f = subst(f, var(i), 0);
1020  }
1021 i = i+1;
1022 }
[078a4a]1023 Phi = @br, B1;
[b41ce4]1024 f = Phi(f);
1025 return( f );
1026}
1027
1028///////////////////////////////////////////////////////////////////////////////
1029static proc Coeff(poly f, list #)
1030{
1031//* initialisation
1032  poly   a, term;
1033  int    n, i;
1034  matrix K;
1035
1036  n     = nvars(basering);
1037  i     = 1;
1038  term  = #[2];
1039  K     = coef(f, #[1]);
1040
1041  while( (i<=ncols(K)) && (K[1,i] != term) )
1042  { i++;
1043    if(i>ncols(K)) { break; }
1044  }
1045  if(i<=ncols(K)) { a = K[2,i]; }
1046  if(i>ncols(K)) { a = 0; }
1047
1048  return(a);
1049}
1050
1051///////////////////////////////////////////////////////////////////////////////
1052static proc ReorderVar( poly f )
1053"USAGE:    GetRf();"
1054{
1055//* Initialisation
1056 intvec frequ;
1057 int maxim, maxim_var, i, j, n;
1058 ideal B = maxideal(1);
1059
1060//* save basering;
[078a4a]1061 def @br = basering;
[b41ce4]1062 n = nvars(basering);
1063
1064//* frequ[i] determines the frequency of appearance of x_i in f
1065 for ( i=1; i <= n; i=i+1 ){
1066  frequ[i] = ncols( coef( f, var(i) ));
1067  if( Coeff(f, var(i), 0) == 0 ) { frequ[i] = frequ[i]+1; }
1068 }
1069
1070//* Determines the order of the variables, lowest frequency first
1071 for( i=n; i>0; i=i-1 ){
1072  maxim = 0;
1073  maxim_var = 0;
1074  for (j = 1; j<= n; j=j+1 ){
1075   if(frequ[j] > maxim ){ maxim = frequ[j]; maxim_var = j;}
1076  }
1077  B[maxim_var]= var(i);
1078  frequ[maxim_var]=-1;
1079 }
1080
1081//* Reorder variables
[078a4a]1082 map Phi = @br, B;
[b41ce4]1083 return (Phi(f));
1084}
1085
1086///////////////////////////////////////////////////////////////////////////////
1087static proc Hcode (intvec a)
1088"USAGE:   Hcode(a); a intvec
1089RETURN:  intvec consisting of the numbers of successive repetitions of
1090         entries in a
1091EXAMPLE: example Hcode; shows an example."
1092{
1093  int n=size(a)-1;
1094  a[n+1]=0;
1095  int i,j,k=1,0,k;
1096  int a0=a[i];
1097  intvec c;
1098
1099  while( i <= n )
1100  {
1101    j++;
1102    k=0;
1103    while( i <= n && a0 == a[i] )
1104    {
1105      i++;
1106      k++;
1107    }
1108    c[j]=k;
1109    a0=a[i];
1110  }
1111
1112  return(c);
1113}
1114example
1115{ "EXAMPLE:"; echo=2;
1116  intvec v = 1,3,4,4,4,4,4,4,4,3,9,9,9,9,1;
1117  Hcode(v);
1118}
1119
1120///////////////////////////////////////////////////////////////////////////////
1121proc arnoldMilnorCode (poly f, list #)
1122"USAGE:   arnoldMilnorCode(f[,e]); f poly, e int
1123ASSUME:  basering is local, f has isolated critical point at 0
1124COMPUTE: Milnor code of f consisting of the numbers of successive repetitions
1125         of coefficients of the 2nd Hilbert series of basering/(jacob(f)^e),
1126         see [Bac01].
1127RETURN:  Milnor code of f as intvec where e=1 by default
1128EXAMPLE: example arnoldMilnorCode; shows an example"
1129{
1130  int  e=1;
1131  if(size(#)==1) { e=#[1]; }
1132  ideal jf=std(jacob(f)^e);
1133  intvec v=hilb(jf,2);
1134
1135  return(Hcode(v));
1136}
1137example
1138{ "EXAMPLE:"; echo=2;
1139  ring r=0,(x,y,z),ds;
1140  poly f=x2y+y3+z2;
1141  arnoldMilnorCode(f);
1142  arnoldMilnorCode(f,2);
1143  // a big second argument may result in memory overflow
1144}
1145
1146///////////////////////////////////////////////////////////////////////////////
1147proc arnoldDeterminacy( I, list # )
1148"USAGE:   arnoldDeterminacy( I[, m]); I poly or ideal, m int.
[078a4a]1149ASSUME:  the basering is local, I is the Jacobian ideal of a polynomial f
[b41ce4]1150         with isolated critical point at 0, m is the Milnor number of f
1151COMPUTE: determinacy bound k for f w.r.t. right equivalence
[078a4a]1152RETURN:  integer k s.th. f is right-k-determined, -1 for infinity
[b41ce4]1153NOTE:    uses [Cor. A.9.7,GP12]
1154EXAMPLE: example arnoldDeterminacy; shows an example"
1155{
1156 //* Case: input = poly *//
1157 if(typeof(I)=="poly")
1158 {
1159  ideal J=std(jacob(I));
1160  return(arnoldDeterminacy(J, #));
1161  }
1162 //* Case Input = Ideal *//
1163  if (typeof(I)=="ideal")
1164  {
1165  int k;   //* upper bound of determinacy
1166  int time;
1167
1168  //* Computation of Milnor number *//
1169  if( size(#) > 0){  k = #[1] + 1; }
1170  else { k = vdim(std(I)) + 1;}
1171
1172  //* If milnor number infinite > return -1 otherwise apply A.9.7.
1173  if( k == -1 ){ return (-1); }
1174
1175   int m;
1176   I=std(I);
1177  for(int i=0;i<=2;i++)
1178  {
1179   m = deg(highcorner(I))+2-i;
1180   if(m<k){ k = m;}
1181
1182   //* if the computation of the standard bases for takes > 10 sek
1183   //* algo breaks and returns upper bound found so far
1184   time = timer;
1185   if(i<2){ I = std(maxideal(1)* I); }
1186   if( timer - time > 12 ){ break; }
1187  }
1188    return( k );
1189  }
1190 ERROR("The input has to be a Jacobian ideal or a polynomial");
1191}
1192examples
1193{
1194  ring r=0,(x,y),ds;
1195  poly f=x3+xy3;
1196  ideal I=std(jacob(f));
1197  int k=arnoldDeterminacy(I);
1198  print(k);
1199}
1200
1201
1202///////////////////////////////////////////////////////////////////////////////
1203proc arnoldCorank(poly f)
1204"USAGE:   arnoldCorank(f);  f poly
1205ASSUME:  basering is local, f in maxideal(2) has isolated critical point at 0
1206RETURN:  corank of the Hessian matrix of f
1207EXAMPLE: example arnoldCorank; shows an example"
1208{
1209  matrix M = jacob(jacob(jet(f,2)));
1210  list lba = bareiss(M);
1211  int cr = nvars(basering) - size(module(lba[1]));
1212  return(cr);
1213}
1214example
1215{ "EXAMPLE:"; echo=2;
1216  ring r=0,(x,y,z),ds;
1217  poly f=(x2+3y-2z)^2+xyz-(x-y3+x2*z3)^3;
1218  arnoldCorank(f);
1219}
1220
1221///////////////////////////////////////////////////////////////////////////////
1222static proc check_basering()
1223{
1224  if( char(basering) >= 13 )
1225  {
1226    ERROR("The characteristic of the basering must be at most 13.");
1227    return(1);
1228  }
1229  int n=nvars(basering);
1230  int i=1;
1231  while( i<n && 1>var(i) ) { i++; }
1232  if( 1<var(i) )
1233  {
1234    ERROR("The basering must be local.");
1235    return(1);
1236  }
1237  return(0);  // basering is OK, return (0)
1238}
1239
1240///////////////////////////////////////////////////////////////////////////////
1241static proc init_newsingclass( singclass f )
1242"USAGE:  creates a new singclass, sets the iPoly to the f_input and
1243        the integer values to -1 resp. -2 to indicate that these values
1244        have not been computed yet
1245RETURN: singclass f"
1246{
1247 f.Corank = -1;
1248 f.Milnor = -2;
1249 f.Determinacy = -2;
1250 f.Tjurina = -2;
1251 f.Modality = -1;
1252 intvec mc = -1;
1253 f.MilnorCode = mc;
1254 f.k=-1;
1255 f.r=-1;
1256 f.s=-1;
1257
1258 return(f);
1259}
1260
1261///////////////////////////////////////////////////////////////////////////////
1262proc arnoldShowSeries( string typ )
1263"USAGE:   arnoldShowSeries( S );  S string
1264ASSUME:  S is the name of a singularity series listed by arnoldListAllSeries().
1265RETURN:  data of the singularity series S of type singseries including
1266         @* - Milnor number of S,
1267         @* - Corank of S,
1268         @* - Milnor code of S (see [Bac01]),
1269         @* - normal form of S as string with parameters k,r,s and a,b,c,d,
1270         @* - restrictions on parameters in the normal form in SINGULAR syntax,
1271         @* - normal form with special (valid) parameters.
1272EXAMPLE: example arnoldShowSeries; shows an example
1273"
1274{
1275 string DatabasePath, Database, value;
1276
1277//* Defining the link to the database Singularitylist which stores
1278//* the singularities
1279 string DBMPATH=system("getenv","DBMPATH");
1280 if( DBMPATH != "" ) { DatabasePath = DBMPATH+"/Singularitylist"; }
1281 else { DatabasePath = "Singularitylist"; }
1282 Database="DBM: ",DatabasePath;
1283
1284//* link to the database
1285 link dbmLink=Database;
1286//* database Singularitylist entry for key
1287 value = read(dbmLink, typ);
1288if(value =="")
1289   {ERROR("SingularitySeries not known. Call arnoldListAllSeries() to get
1290    a list of all valid Singularity types.");}
1291 close(dbmLink);
1292 execute(value);
1293 return( f );
1294}
1295example
1296{ "EXAMPLE:"; echo=2;
1297  arnoldShowSeries("Z[k,12k+6r]");
1298}
1299
1300
1301///////////////////////////////////////////////////////////////////////////////
1302static proc normalformDB( typ )
1303{
1304 if( typeof(typ) == "string" )
1305  {
1306  def f = arnoldShowSeries(typ);
1307  return( f.NormalForm );
1308  }
1309
1310 if( typeof(typ) =="singclass")
1311  {
1312  def f = arnoldShowSeries( typ.Series );
1313  return( f.NormalForm );
1314  }
1315}
1316
1317///////////////////////////////////////////////////////////////////////////////
1318static proc NormalFormAndRestrictionsDB( typ )
1319{
1320 if( typeof(typ) == "string" )
1321  {
1322  def f = arnoldShowSeries(typ);
1323  return( f.NormalForm, f.Restrictions);
1324  }
1325
1326 if( typeof(typ) =="singclass")
1327  {
1328  def f = arnoldShowSeries( typ.Series );
1329  return( f.NormalForm, f.Restrictions );
1330  }
1331}
1332
1333///////////////////////////////////////////////////////////////////////////////
1334static proc specialformDB( typ )
1335{
1336 if( typeof(typ) == "string" )
1337  {
1338  def f = arnoldShowSeries(typ);
1339  return( f.SpecialForm );
1340  }
1341
1342 if( typeof(typ) =="singclass")
1343  {
1344  def f = arnoldShowSeries( typ.Series );
1345  return( f.SpecialForm );
1346  }
1347}
1348
1349///////////////////////////////////////////////////////////////////////////////
1350static proc restrictionsDB( typ )
1351{
1352 if( typeof(typ) == "string" )
1353  {
1354  def f = arnoldShowSeries(typ);
1355  return( f.Restrictions );
1356  }
1357
1358 if( typeof(typ) =="singclass")
1359  {
1360  def f = arnoldShowSeries( typ.Series );
1361  return( f.Restrictions );
1362  }
1363}
1364
1365///////////////////////////////////////////////////////////////////////////////
1366proc arnoldNormalForm( typ, list #)
1367"USAGE:   arnoldNormalForm( S [, l]), S string or singclass, l list
1368ASSUME:  If S is of type string, then S is the name of a singularity series
1369         as listed by arnoldListAllSeries() and l may contain suitable
1370         integer parameters k,r,s. Otherwise S of type singclass is a
1371         singularity class from Arnol'd's list.
1372         Optional suitable polynomial parameters a,b,c,d can be appended to l.
1373         If a,b,c,d are not given, valid values are chosen automatically.
1374RETURN:  string NF is the normal form of the series S if no parameters given,
[417a505]1375         or poly NF is the normal form of the class S with parameters k,r,s.
[b41ce4]1376EXAMPLE: example arnoldNormalForm; shows an example"
1377{
1378
1379//* Sort list entries by type
1380 int noparas = 1;
1381 list polys;
1382 for( int i = 1; i<=size(#); i++)
1383 {
1384  if(typeof( #[i] ) == "int" )
1385  {
1386   if( noparas == 1 ){ intvec paras = #[1]; noparas =0; }
1387   else{ paras[size(paras)+1] = #[i];}
1388  }
1389  if( typeof(#[i]) == "poly")
1390  {
1391   polys[size(polys)+1] = #[i];
1392  }
1393  if( typeof(#[i]) != "int" && typeof(#[i]) != "poly")
1394  {ERROR("Entries of list # must be of type poly or int.");}
1395 }
1396
1397//* If no paras and typ = "string" -> string
1398  if( (typeof (typ) == "string") && (noparas == 1) )
1399  { return( specialformDB(typ) ); }
1400
1401//* Set parameter k,r,s
1402 if( typeof( typ ) == "singclass" )
1403 { int k = typ.k; int r = typ.r; int s = typ.s; string Series = typ.Series; }
1404 if( typeof( typ ) == "string" )
1405 { int k,r,s = paras; string Series = typ; }
1406
1407//* Check if basering is defined
1408 if(!defined(basering)){
1409  ERROR("No basering defined!")}
1410
1411//* Check case k=1
1412 if( k == 1 ){
1413  if( Series[1] == "Y"){ Series = replace( Series, "k", string(k)); }
1414  if( Series[1] == "X"){ Series = replace( Series, "k", string(k)); }
1415  if( Series[1] == "Z"){ Series = replace( Series, "k", string(k)); }
1416 }
1417
1418//* If no polynomial parameter Specialform
1419  if(size(polys) != 0 && size(polys) <4 )
1420  {ERROR("Not enough polynomial parameter as input.")}
1421  if(size(polys) == 0 ){
1422    //* Define ring for the NF; NF given in C(x,y)
[078a4a]1423    def @br = basering;
[b41ce4]1424      if(defined(RingNF) != 0 ) { kill RingNF; }
1425    ring RingNF=char(basering),(x,y),(c, ds);
[078a4a]1426    map Conv = @br,maxideal(1);
[b41ce4]1427
1428    string nf_str = specialformDB(Series);
1429    }
1430  else{
1431
1432   //* Set polynomial parameter
1433   poly a,b,c,d;
1434   a = polys[1];
1435   b = polys[2];
1436   c = polys[3];
1437   d = polys[4];
1438
1439   //* Define ring for the NF; NF given in C(x,y)
[078a4a]1440    def @br = basering;
[b41ce4]1441     if(defined(RingNF) != 0 ) { kill RingNF; }
1442     ring RingNF=char(basering),(x,y),(c, ds);
[078a4a]1443     map Conv = @br,maxideal(1);
[b41ce4]1444
1445   //* Map Polynomial parameter
1446   poly a = Conv(a);
1447   poly b = Conv(b);
1448   poly c = Conv(c);
1449   poly d = Conv(d);
1450
1451   //* Get NormalForm from Databank
1452   string nf_str = normalformDB(Series);
1453
1454   //* Replace polynomial parameter
1455   nf_str = replace( nf_str, "a(y)", "("+string(a)+")");
1456   nf_str = replace( nf_str, "b(y)", "("+string(b)+")");
1457   nf_str = replace( nf_str, "c(y)", "("+string(c)+")");
1458   nf_str = replace( nf_str, "d(y)", "("+string(d)+")");
1459
1460   //* Check parameter
1461   int okay = 0;
1462   string res = restrictionsDB(Series);
1463   string str = "if(" + res + "){okay = 1;};";
1464   execute(str);
1465   if(okay!=1)
1466   {ERROR("Parameter do not satisfy restrictions of
1467   the singularity series.")}
1468 }
1469
1470 //* Replace integer parameter
1471 nf_str = replace( nf_str, "k", string(k));
1472 nf_str = replace( nf_str, "r", string(r));
1473 nf_str = replace( nf_str, "s", string(s));
1474
1475 //* Generate polynomial nf
1476 nf_str = "poly nf  = " + nf_str;
1477 execute( nf_str );
1478
1479 //* Map poly nf back to basering;
[078a4a]1480 setring @br;
[b41ce4]1481 map ConvBack = RingNF, maxideal(1);
1482  return( ConvBack(nf) );
1483}
1484example
1485{ "EXAMPLE:"; echo=2;
1486  ring R = 0, (x,y), ds;
1487 poly a,b,c,d;
1488 a= 1+y2;
1489 c= 3+y;
1490 int k = 5;
1491 int r = 4;
1492 int s = 3;
1493 arnoldNormalForm ("W[12k+1]", k,r,s,a,b,c,d);
1494 def f = _;
1495 def sf = arnoldClassify( f );
1496 arnoldNormalForm(sf, a,b,c,d);
1497 arnoldNormalForm("W[12k+1]");
1498 arnoldNormalForm(sf);
1499}
1500
1501///////////////////////////////////////////////////////////////////////////////
1502proc arnoldListAllSeries()
1503"USAGE:    arnoldListAllSeries();
1504RETRURN:  list of names of singularity series listed by Arnol'd up to corank 2
1505EXAMPLE:  example arnoldListAllSeries; shows an example"
1506{
1507list l=
1508"A[k]",
1509"D[k]",
1510"E[6k]",
1511"E[6k+1]",
1512"E[6k+2]",
1513"J[k,0]",
1514"J[k,r]",
1515"W[12k]",
1516"W[12k+1]",
1517"W[12k+5]",
1518"W[12k+6]",
1519"W[k,0]",
1520"W[k,r]",
1521"W#[k,2r]",
1522"W#[k,2r-1]",
1523"X[k,0]",
1524"X[k,r]",
1525"Y[k,r,s]",
1526"Z[k,r]",
1527"Z[k,r,s]",
1528"Z[1,6r+11]",
1529"Z[1,6r+12]",
1530"Z[1,6r+13]",
1531"Z[k,12k+6r]",
1532"Z[k,12k+6r+1]",
1533"Z[k,12k+6r-1]";
1534
1535return(l);
1536}
1537example
1538{ arnoldListAllSeries();
1539}
1540
1541///////////////////////////////////////////////////////////////////////////////
1542///////////////////////////////////////////////////////////////////////////////
1543//* The following part of the library generates two files, Singularitylist.dir
1544//* and Singularitylist.pag containing a data base for singularities up to
1545//* corank 2 listed by Arnol'd.
1546
1547///////////////////////////////////////////////////////////////////////////////
1548static proc makedbm_init()
1549{
1550//* Generates file containing a data base for singularities up to corank 2
1551//* listed by Arnol'd. This file is needed for arnoldclassify.lib.
1552
1553  string s;
1554  link l="DBM:r Singularitylist";
1555  s = read(l,"VERSION");
1556  if (s == "" ) {
1557    "Need to create database...";
1558    create_singularity_dbm();
1559  }
1560  close(l);
1561  l="DBM:r Singularitylist";
1562  s = read(l,"VERSION");
1563  "Creation done. Current version:", s;
1564}
1565///////////////////////////////////////////////////////////////////////////////
1566
1567static proc dbm_read (link l)
1568{
1569  string s="";
1570  s=read(l);
1571  while( s != "" )
1572  {
1573    s,"=",read(l,s);
1574    s=read(l);
1575  }
1576}
1577
1578///////////////////////////////////////////////////////////////////////////////
1579static proc dbm_getnext (link l)
1580{
1581  string s="";
1582  s=read(l);
1583  if( s != "" ) { s,"=",read(l,s); }
1584}
1585
1586///////////////////////////////////////////////////////////////////////////////
1587static proc create_singularity_dbm
1588{
1589  link l="DBM:rw Singularitylist";
1590
1591//*Data typ singseries;
1592  string s;
1593
1594//* A[k]
1595    s = "singseries f;
1596    f.Series = \"A[k]\";
1597    f.Modality = \"0\";
1598    f.Corank = \"1\";
1599    f.MilnorNumber = \"k\";
1600    f.MilnorCode = \"k\";
1601    f.NormalForm = \"x^(k+1)\";
1602    f.SpecialForm = \"x^(k+1)\";
1603    f.Restrictions = \"(k>1)\";";
1604  write(l, "A[k]", s);
1605
1606//* D[k]
1607    s = "singseries f;
1608    f.Series = \"D[k]\";
1609    f.Modality = \"0\";
1610    f.Corank = \"2\";
1611    f.MilnorNumber = \"k\";
1612    f.MilnorCode = \"1,1,k-3\";
1613    f.NormalForm = \"x^2*y+y^(k-1)\";
1614    f.SpecialForm = \"x^2*y+y^(k-1)\";
1615    f.Restrictions = \"(k>=4)\";";
1616  write(l, "D[k]", s);
1617
1618//* J[k,0]
1619    s = "singseries f;
1620    f.Series = \"J[k,0]\";
1621    f.Modality = \"0\";
1622    f.Corank = \"2\";
1623    f.MilnorNumber = \"6*k-2\";
1624    f.MilnorCode = \"1,2*k+j,2*k-2*j-3\";
1625    f.NormalForm = \"x^3 + b(y)*x^2*y^k+c(y)*x*y^(2*k+1)+y^(3*k)\";
1626    f.SpecialForm = \"x^3 + x^2*y^k+y^(3*k)\";
1627    f.Restrictions = \"(k>1)&& (4*b^3 + 27 != 0)&&
1628     (deg(b)==0)&&(deg(c)<=(k-3))&&(k>2||c==0)\";";
1629  write(l, "J[k,0]", s);
1630
1631//* J[k,r]
1632    s = "singseries f;
1633    f.Series = \"J[k,r]\";
1634    f.Modality = \"0\";
1635    f.Corank = \"2\";
1636    f.MilnorNumber = \"6*k-2+r\";
1637    f.MilnorCode = \"1,2*k-1,2*k+r-1\";
1638    f.NormalForm = \"x^3 + x^2*y^k+a(y)*y^(3*k+r)\";
1639    f.SpecialForm = \"x^3 + x^2*y^k+y^(3*k+r)\";
1640    f.Restrictions = \"(k>1)&&(r>0)&&(jet(a,0)!= 0)&&(deg(a)<=(k-2)) \";";
1641  write(l, "J[k,r]", s);
1642
1643//* E[6k]
1644    s = "singseries f;
1645    f.Series = \"E[6k]\";
1646    f.Modality = \"0\";
1647    f.Corank = \"2\";
1648    f.MilnorNumber = \"6*k\";
1649    f.MilnorCode = \"1,2*k+j,2*k-2j-1\";
1650    f.NormalForm = \"x^3 + a(y)*x*y^(2*k+1)+y^(3*k+1)\";
1651    f.SpecialForm = \"x^3+y^(3*k+1)\";
1652    f.Restrictions = \"(k>=1)&&(k>1||a==0)&&(deg(a)<=(k-2))\";";
1653  write(l, "E[6k]", s);
1654
1655//* E[6k+1]
1656    s = "singseries f;
1657    f.Series = \"E[6k+1]\";
1658    f.Modality = \"0\";
1659    f.Corank = \"2\";
1660    f.MilnorNumber = \"6*k+1\";
1661    f.MilnorCode = \"1,2*k,2*k\";
1662    f.NormalForm = \"x^3 + x*y^(2*k+1)+a(y)*y^(3*k+2)\";
1663    f.SpecialForm = \"x^3 + x*y^(2*k+1)\";
1664    f.Restrictions = \"(k>=1)&&(k>1||a==0)&&(deg(a)<=(k-2))\";";
1665  write(l, "E[6k+1]", s);
1666
1667//* E[6k+2]
1668    s = "singseries f;
1669    f.Series = \"E[6k+2]\";
1670    f.Modality = \"0\";
1671    f.Corank = \"2\";
1672    f.MilnorNumber = \"6*k+2\";
1673    f.MilnorCode = \"1,2*k+j+1,2*k-2j-1\";
1674    f.NormalForm = \"x^3 + a(y)*x*y^(2*k+2)+y^(3*k+2)\";
1675    f.SpecialForm = \"x^3 +y^(3*k+2)\";
1676    f.Restrictions = \"(k>=1)&&(k>1||a==0)&&(deg(a)<=(k-2))\";";
1677  write(l, "E[6k+2]", s);
1678
1679//* X[k,0]
1680    s = "singseries f;
1681    f.Series = \"X[k,0]\";
1682    f.Modality = \"3*k-2\";
1683    f.Corank = \"2\";
1684    f.MilnorNumber = \"12*k-3\";
1685    f.MilnorCode = \"1,1,2*k-1+j,2k-1-2*j+t,2*k-1+j-2t\";
1686    f.NormalForm = \"x^4 + b(y)*x^3*y^k + a(y)*x^2*y^(2*k) + x*y^(3*k)\";
1687    f.SpecialForm = \"x^4 + x^3*y^k + x*y^(3*k)\";
1688    f.Restrictions = \"(jet(a,0)*jet(b,0)!=9)&&(k>1)&&(4*(jet(a,0)^3+jet(b,0)^3)
1689     - jet(a,0)^2*jet(b,0)^2-18* jet(a,0)*jet(b,0) + 27 !=0)&&(deg(a)<=(k-2))
1690     &&(deg(b)<=(2*k-2))\";";
1691  write(l, "X[k,0]", s);
1692
1693    //* X[1,0]
1694        s = "singseries f;
1695        f.Series = \"X[1,0]\";
1696        f.Modality = \"1\";
1697        f.Corank = \"2\";
1698        f.MilnorNumber = \"9\";
1699        f.MilnorCode = \"1,1,1+j,1-2*j+t,1+j-2t\";
1700        f.NormalForm = \"x^4 + a(y)*x^2*y^2 + y^4\";
1701        f.SpecialForm = \"x^4 + x^2*y^2 + y^4\";
1702        f.Restrictions = \"(deg(a)==0)&&(jet(a,0)^2!=4)\";";
1703      write(l, "X[1,0]", s);
1704
1705    //* X[k,r]
1706        s = "singseries f;
1707        f.Series = \"X[k,r]\";
1708        f.Modality = \"3*k-2\";
1709        f.Corank = \"2\";
1710        f.MilnorNumber = \"12*k-3+r\";
1711        f.MilnorCode = \"1,1,2*k-1+j,2k-1-2*j,2*k-1+j+r\";
1712        f.NormalForm = \"x4+a(y)*x3*y^(k)+x^2*y^(2*k)+b(y)*y^(4*k+r)\";
1713        f.SpecialForm = \"x4+x3*y^(k)+x^2*y^(2*k)+y^(4*k+r)\";
1714        f.Restrictions = \"(k>1)&&(r>0)&&(deg(a)<=(k-2))&&(jet(a,0)^2!=4)&&
1715        (jet(b,0)!=0)&&(deg(b)<=(2*k-2))\";";
1716      write(l, "X[k,r]", s);
1717
1718    //* X[1,r]
1719        s = "singseries f;
1720        f.Series = \"X[1,r]\";
1721        f.Modality = \"1\";
1722        f.Corank = \"2\";
1723        f.MilnorNumber = \"9+r\";
1724        f.MilnorCode = \"1,1,1+j,1-2*j,1+j+r\";
1725        f.NormalForm = \"x4+x^2*y^2+a(y)*y^(4+r)\";
1726        f.SpecialForm = \"x4+x^2*y^2+y^(4+r)\";
1727        f.Restrictions = \"(deg(a)==0)&&(jet(a,0)!=0)\";";
1728        write(l, "X[1,r]", s);
1729
1730    //* Y[k,r,s]
1731        s = "singseries f;
1732        f.Series = \"Y[k,r,s]\";
1733        f.Modality = \"3*k-2\";
1734        f.Corank = \"2\";
1735        f.MilnorNumber = \"12*k-3+r+s\";
1736        f.MilnorCode = \"1,1,2*k-1,2*k-1+j,2*k-1-2*j+r+s\";
1737        f.NormalForm = \"((x + a(y)*y^k)^2 + b(y)*y^(2*k+s))*(x2 + y^(2*k+r))\";
1738        f.SpecialForm = \"((x + y^k)^2 + y^(2*k+s))*(x2 + y^(2*k+r))\";
1739        f.Restrictions = \"(jet(a,0)!=0)&&(deg(a)<=(k-2))&&(k>1)&&(jet(b,0)!=0)
[078a4a]1740        &&(1<=s)&&(s<=r)\";";
[b41ce4]1741        write(l, "Y[k,r,s]", s);
1742
1743    //* Y[1,r,s]
1744        s = "singseries f;
1745        f.Series = \"Y[1,r,s]\";
1746        f.Modality = \"1\";
1747        f.Corank = \"2\";
1748        f.MilnorNumber = \"9+r+s\";
1749        f.MilnorCode = \"1,1,1,1+j,1-2*j+r+s\";
1750        f.NormalForm = \" x^(4+r)+ a(y)*x2*y2 + y^(4+s)\";
1751        f.SpecialForm = \" x^(4+r)+ x2*y2 + y^(4+s)\";
[078a4a]1752        f.Restrictions = \"(deg(a)==0)&&(jet(a,0)!=0)&&(1<=s)&&(s<=r)\";";
[b41ce4]1753        write(l, "Y[1,r,s]", s);
1754
1755    //* Z[k,r]
1756        s = "singseries f;
1757        f.Series = \"Z[k,r]\";
1758        f.Modality = \"3*k+r-2\";
1759        f.Corank = \"2\";
1760        f.MilnorNumber = \"12*k-3+6*r\";
1761        f.MilnorCode = \"1,1,2*k-1,2*k-1+j,2*k-1+6*r-2*j\";
1762        f.NormalForm = \"(x + a(y)*y^k)*(x^3 + d(y)*x2*y^(k+1) +
1763        c(y)*x*y^(2*k+2*r+1) + y^(3*k+3*r))\";
1764        f.SpecialForm = \"(x + y^k)*(x^3 + 2*y^(k+1) + x*y^(2*k+2*r+1) +
1765        y^(3*k+3*r))\";
1766        f.Restrictions = \"(k>1)&&(r>=0)&&(4*d^3+27!=0)&&(deg(d)==0)&&
1767        (deg(c)<=(2*k+r-3))&&(deg(a)<=(k-2))\";";
1768        write(l, "Z[k,r]", s);
1769
1770    //* Z[1,r]
1771        s = "singseries f;
1772        f.Series = \"Z[1,r]\";
1773        f.Modality = \"1+r\";
1774        f.Corank = \"2\";
1775        f.MilnorNumber = \"9+6*r\";
1776        f.MilnorCode = \"1,1,1,1+j,1+6*r-2*j\";
1777        f.NormalForm = \"y*(x^3 + d(y)*x^2*y^(2) + c(y)*x*y^(2+2*r+1) +
1778        y^(3+3*r))\";
1779        f.SpecialForm = \"y*(x^3 + x^2*y^(2) + x*y^(2+2*r+1) +
1780        y^(3+3*r))\";
1781        f.Restrictions = \"(r>=0)&&(4*d^3+27!=0)&&(deg(d)==0)
1782        &&(deg(c)<=(r-1))\";";
1783        write(l, "Z[1,r]", s);
1784
1785    //* Z[k,r,s]
1786        s = "singseries f;
1787        f.Series = \"Z[k,r,s]\";
1788        f.Modality = \"3*k+r-2\";
1789        f.Corank = \"2\";
1790        f.MilnorNumber = \"12*k+6*r+s-3\";
1791        f.MilnorCode = \"1,1,2*k-1,2*k-1+2*r,2*k-1+2*r-s\";
1792        f.NormalForm = \"(x^2 + a(y)*x*y^k + b(y)*y^(2*k+r))*
1793        (x^2 + y^(2*k+2*r+s))\";
1794        f.SpecialForm = \"(x^2 + x*y^k + y^(2*k+r))*(x^2 + y^(2*k+2*r+s))\";
1795        f.Restrictions = \"(k>1)&&(r>=0)&&(deg(a)<=(k-2))&&(jet(a,0)!=0)&&
1796        (jet(b,0)!=0)&&(deg(b)<=(2*k+r-2))\";";
1797        write(l, "Z[k,r,s]", s);
1798
1799    //* Z[1,r,s]
1800        s = "singseries f;
1801        f.Series = \"Z[1,r,s]\";
1802        f.Modality = \"1+r\";
1803        f.Corank = \"2\";
1804        f.MilnorNumber = \"9+6*r+s\";
1805        f.MilnorCode = \"1,1,1,1+2*r,1+2*r-s\";
1806        f.NormalForm = \"y*(x^3 + x^2*y^(r+1) + b(y)*y^(3*r+s+3))\";
1807        f.SpecialForm = \"y*(x^3 + x^2*y^(r+1) + y^(3*r+s+3))\";
1808        f.Restrictions = \"(r>=0)&&(jet(b,0)!=0)&&(deg(b)<=(2*k+r-2))\";";
1809        write(l, "Z[1,r,s]", s);
1810
1811    //* Z[k,12k+6r-1]
1812        s = "singseries f;
1813        f.Series = \"Z[k,12k+6r-1]\";
1814        f.Modality = \"3*k+r-2\";
1815        f.Corank = \"2\";
1816        f.MilnorNumber = \"12*k+6r-1\";
1817        f.MilnorCode = \"1,1,2k-1,2k-1+j,2k+1+6*r-2*j\";
1818        f.NormalForm = \"(x + a(y)*y^k)*(x^3 + b(y)*x*y^(2*k+2*r+1) +
1819         y^(3*k+3*r+1))\";
1820        f.SpecialForm = \"(x + y^k)*(x^3 + x*y^(2*k+2*r+1) + y^(3*k+3*r+1))\";
1821        f.Restrictions = \" (k>1)&&(r>=0)&&(deg(a)<=(k-2))&&(jet(a,0)!=0)&&
1822        (jet(b,0)!=0)&&(deg(b)<=(2*k+r-2))\";";
1823        write(l, "Z[k,12k+6r-1]", s);
1824
1825        //* Z[1,6r+11]
1826            s = "singseries f;
1827            f.Series = \"Z[1,6r+11]\";
1828            f.Modality = \"1+r\";
1829            f.Corank = \"2\";
1830            f.MilnorNumber = \"6r+11\";
1831            f.MilnorCode = \"1,1,1,1+j,3+6*r-2*j\";
1832            f.NormalForm = \"y*(x^3 + b(y)*x*y^(2+2*r+1) + y^(3+3*r+1))\";
1833            f.SpecialForm = \"y*(x^3 + x*y^(2+2*r+1) + y^(3+3*r+1))\";
1834            f.Restrictions = \"(r>=0)&&(deg(b)<=(r))\";";
1835            write(l, "Z[1,6r+11]", s);
1836
1837    //* Z[k,12k+6r+1]
1838        s = "singseries f;
1839        f.Series = \"Z[k,12k+6r+1]\";
1840        f.Modality = \"3*k+r-2\";
1841        f.Corank = \"2\";
1842        f.MilnorNumber = \"12*k+6r+1\";
1843        f.MilnorCode = \"1,1,2k-1,2k-1+j,2k+3+6*r-2*j\";
1844        f.NormalForm = \"(x + a(y)*y^k)*(x^3 + b(y)*x*y^(2*k+2*r+2) +
1845         y^(3*k+3*r+2))\";
1846        f.SpecialForm = \"(x + y^k)*(x^3 + x*y^(2*k+2*r+2) + y^(3*k+3*r+2))\";
1847        f.Restrictions = \" (k>1)&&(r>=0)&&(deg(a)<=(k-2))&&(jet(a,0)!=0)&&
1848        (jet(b,0)!=0)&&(deg(b)<=(2*k+r-2))\";";
1849        write(l, "Z[k,12k+6r+1]", s);
1850
1851    //* Z[1,6r+13]
1852        s = "singseries f;
1853        f.Series = \"Z[1,6r+13]\";
1854        f.Modality = \"1+r\";
1855        f.Corank = \"2\";
1856        f.MilnorNumber = \"6r+13\";
1857        f.MilnorCode = \"1,1,1,1+j,5+6*r-2*j\";
1858        f.NormalForm = \"y*(x^3 + b(y)*x*y^(2*r+4) + y^(3*r+5))\";
1859        f.SpecialForm = \"y*(x^3 + x*y^(2*r+4) + y^(3*r+5))\";
1860        f.Restrictions = \" (r>=0)&&(deg(b)<=(r))\";";
1861        write(l, "Z[1,6r+13]", s);
1862
1863    //* Z[k,12k+6r]
1864        s = "singseries f;
1865        f.Series = \"Z[k,12k+6r]\";
1866        f.Modality = \"3*k+r-2\";
1867        f.Corank = \"2\";
1868        f.MilnorNumber = \"12*k+6r\";
1869        f.MilnorCode = \"1,1,2k-1,2k-1+2*r,2k+2*r\";
1870        f.NormalForm = \"(x + a(y)*y^k)*(x^3 + x*y^(2*k+2*r+1) +
1871        b(y)* y^(3*k+3*r+2))\";
1872        f.SpecialForm = \"(x + y^k)*(x^3 + x*y^(2*k+2*r+1) +y^(3*k+3*r+2))\";
1873        f.Restrictions = \" (k>1)&&(r>=0)&&(deg(a)<=(k-2))&&(jet(a,0)!=0)&&
1874        (jet(b,0)!=0)&&(deg(b)<=(2*k+r-2))\";";
1875        write(l, "Z[k,12k+6r]", s);
1876
1877
1878    //* Z[1,6r+12]
1879        s = "singseries f;
1880        f.Series = \"Z[1,6r+12]\";
1881        f.Modality = \"1+r\";
1882        f.Corank = \"2\";
1883        f.MilnorNumber = \"6*r+12\";
1884        f.MilnorCode = \"1,1,1,1+2*r,2+2*r\";
1885        f.NormalForm = \"y*(x^3 + x*y^(2*r+3) +b(y)* y^(3*r+5))\";
1886        f.SpecialForm = \"y*(x^3 + x*y^(2*r+3) +y^(3*r+5))\";
1887        f.Restrictions = \"(r>=0)&&(deg(b)<=(r))\";";
1888        write(l, "Z[1,6r+12]", s);
1889
1890
1891    //* W[k,r]
1892        s = "singseries f;
1893        f.Series = \"W[k,r]\";
1894        f.Modality = \"3*k-1\";
1895        f.Corank = \"2\";
1896        f.MilnorNumber = \"12*k+3+r\";
1897        f.MilnorCode = \"1,1,2k,2k,2k+r\";
1898        f.NormalForm = \"x4+a(y)*x^3*y^(k+1)+x^2*y^(2*k+1)+b(y)*y^(4*k+2+r) \";
1899        f.SpecialForm = \"x4+x^2*y^(2*k+1)+y^(4*k+2+r) \";
1900        f.Restrictions = \"(k>=1)&&(r>0)&&(k>1||a==0)&&(deg(a)<=(k-2))&&
1901        (jet(b,0)!=0)&&(deg(b)<=(2*k-1))\";";
1902        write(l, "W[k,r]", s);
1903
1904    //* W[k,0]
1905        s = "singseries f;
1906        f.Series = \"W[k,0]\";
1907        f.Modality = \"3*k-1\";
1908        f.Corank = \"2\";
1909        f.MilnorNumber = \"12*k+3\";
1910        f.MilnorCode = \"1,1,2k+j,2k-2-2*j+t,2k+6+j+2*t\";
1911        f.NormalForm = \"x4+b(y)*x2*y^(2*k+1)+a(y)*x*y^(3*k+2)+y^(4*k+2) \";
1912        f.SpecialForm = \"x4+x2*y^(2*k+1)+y^(4*k+2) \";
1913        f.Restrictions = \" (k>=1)&&(k>1||a==0)&&(deg(a)<=(k-2))&&
1914        (jet(b,0)^2!=4)&&(deg(b)<=(2*k-1))\";";
1915        write(l, "W[k,0]", s);
1916
1917    //* W[12k]
1918        s = "singseries f;
1919        f.Series = \"W[12k]\";
1920        f.Modality = \"3*k-2\";
1921        f.Corank = \"2\";
1922        f.MilnorNumber = \"12*k\";
1923        f.MilnorCode = \"1,1,2k+j,2k-3-2*j+t,2k+3+j-2*t\";
1924        f.NormalForm = \"x4+a(y)*x*y^(3*k+1)+c(y)*x^2*y^(2*k+1)+y^(4*k+1)\";
1925        f.SpecialForm = \"x4+x^2*y^(2*k+1)+y^(4*k+1)\";
1926        f.Restrictions = \"(k>=1)&&(k>1||a==0)&&(deg(a)<=(k-2))&&
1927        (deg(c)<=(2*k-2))\";";
1928        write(l, "W[12k]", s);
1929
1930    //* W[12k+1]
1931        s = "singseries f;
1932        f.Series = \"W[12k+1]\";
1933        f.Modality = \"3*k-2\";
1934        f.Corank = \"2\";
1935        f.MilnorNumber = \"12*k+1\";
1936        f.MilnorCode = \"1,1,2k+j,2k-1-2*j,2k+j\";
1937        f.NormalForm = \"x4+x*y^(3*k+1)+a(y)*x^2*y^(2*k+1)+c(y)*y^(4*k+2) \";
1938        f.SpecialForm = \"x4+x*y^(3*k+1)+y^(4*k+2) \";
1939        f.Restrictions = \"(k>=1)&&(k>1||a==0)&&(deg(a)<=(k-2))&&
1940        (deg(c)<=(2*k-2))\";";
1941        write(l, "W[12k+1]", s);
1942
1943    //* W[12k+5]
1944        s = "singseries f;
1945        f.Series = \"W[12k+5]\";
1946        f.Modality = \"3*k-1\";
1947        f.Corank = \"2\";
1948        f.MilnorNumber = \"12*k+5\";
1949        f.MilnorCode = \"1,1,2k+j,2k+1-2*j,2k+j\";
1950        f.NormalForm = \"x4+x*y^(3*k+2)+a(y)*x^2*y^(2*k+2)+b(y)*y^(4*k+3) \";
1951        f.SpecialForm = \"x4+x*y^(3*k+2)+y^(4*k+3) \";
1952        f.Restrictions = \"(k>=1)&&(k>1||a==0)&&(deg(a)<=(k-2))&&
1953        (deg(b)<=(2*k-1))\";";
1954        write(l, "W[12k+5]", s);
1955
1956    //* W[12k+6]
1957        s = "singseries f;
1958        f.Series = \"W[12k+6]\";
1959        f.Modality = \"3*k-1\";
1960        f.Corank = \"2\";
1961        f.MilnorNumber = \"12*k+6\";
1962        f.MilnorCode = \"1,1,2k+j,2k-3-2*j+t,2k+9+j-2*t\";
1963        f.NormalForm = \"x4+a(y)*x*y^(3*k+3)+b(y)*x^2*y^(2*k+2)+y^(4*k+3) \";
1964        f.SpecialForm = \"x4+x^2*y^(2*k+2)+y^(4*k+3) \";
1965        f.Restrictions = \"(k>=1)&&(k>1||a==0)&&(deg(a)<=(k-2))&&
1966        (deg(b)<=(2*k-1))\";";
1967        write(l, "W[12k+6]", s);
1968
1969    //* W#[k,2r]
1970        s = "singseries f;
1971        f.Series = \"W#[k,2r]\";
1972        f.Modality = \"3*k-1\";
1973        f.Corank = \"2\";
1974        f.MilnorNumber = \"12*k+3+2*r\";
1975        f.MilnorCode = \"1,1,2k,2k+r,2k\";
1976        f.NormalForm = \"(x2+y^(2*k+1))^2+b(y)*x^2*y^(2*k+1+r)+
1977        a(y)*x*y^(3*k+2+r) \";
1978        f.SpecialForm = \"(x2+y^(2*k+1))^2+x^2*y^(2*k+1+r) \";
1979        f.Restrictions = \"(k>=1)&&(r>0)&&(k>1||a==0)&&(deg(a)<=(k-2))&&
1980        (jet(b,0)!=0)&&(deg(b)<=(2*k-1))\";";
1981        write(l, "W#[k,2r]", s);
1982
1983    //* W#[k,2r-1]
1984        s = "singseries f;
1985        f.Series = \"W#[k,2r-1]\";
1986        f.Modality = \"3*k-1\";
1987        f.Corank = \"2\";
1988        f.MilnorNumber = \"12*k+2+2*r\";
1989        f.MilnorCode = \"1,1,2k,2k-3+j,2*k+5+2*r-2*j\";
1990        f.NormalForm = \"(x2+y^(2*k+1))^2+b(y)*x*y^(3*k+1+r)+
1991        a(y)*y^(4*k+2+r)\";
1992        f.SpecialForm = \"(x2+y^(2*k+1))^2+x*y^(3*k+1+r)\";
1993        f.Restrictions = \"(k>=1)&&(r>0)&&(k>1||a==0)&&(deg(a)<=(k-2))
1994        &&(jet(b,0)!=0)&&(deg(b)<=(2*k-1))\";";
1995        write(l, "W#[k,2r-1]", s);
1996
1997 write(l,"VERSION", "1.0");
1998  close(l);
1999}
2000
2001///////////////////////////////////////////////////////////////////////////////
2002static proc read_singularity_db( string typ )
2003{
2004  string DBMPATH=system("getenv","DBMPATH");
2005  string DatabasePath, Database, S, Text, Tp;
2006
2007  if( DBMPATH != "" ) { DatabasePath = DBMPATH+"/Singularitylist"; }
2008  else { DatabasePath = "Singularitylist"; }
2009  Database="DBM: ",DatabasePath;
2010
2011  link dbmLink=Database;
2012  Tp = read(dbmLink, typ);
2013    return(Tp);
2014
2015}
2016
[a1f059]2017proc arnoldClassify_to_string(poly  fPoly)
2018"USAGE:   arnoldClassify_to_string (f); f poly
2019ASSUME:  The basering is local of characteristic 0 and f defines an
2020         isolated singularity from Arnol'd's list of corank at most 2.
2021COMPUTE: singularity class with respect to right equivalence and
2022         invariants used in the process of classification
2023RETURN:  string: separated by |:
2024         @* - name of singularity series as listed by arnoldListAllSeries(),
2025         @* - name of singularity class,
2026         @* - parameters k,r,s defining the singularity class, -1 if not used,
2027         @* - modality, corank, Milnor number, determinacy,
2028         @* - Tjurina number, -2 if not computed, -1 if infinite,
2029         @* - Milnor code, -1 if not computed,
2030         @* - normal form of the singularity series from Arnol'd's list,
2031         @* - restrictions on parameters as string in SINGULAR syntax.
2032EXAMPLE: example arnoldClassify_to_string; shows an example
2033"
2034{
2035  singclass f=arnoldClassify(fPoly);
2036  string s=f.Restrictions+"|"+
2037    f.NormalForm+"|"+
2038    string(f.MilnorCode)+"|"+
2039    string(f.Tjurina)+"|"+
2040    string(f.Determinacy)+"|"+
2041    string(f.Milnor)+"|"+
2042    string(f.Corank)+"|"+
2043    string(f.Modality)+"|"+
2044    string(f.r)+"|"+
2045    string(f.s)+"|"+
2046    f.Class+"|"+
2047    f.Series;
2048  return(s);
2049}
2050example
2051{ "EXAMPLE:"; echo=2;
2052  ring r = 0,(x,y),ds;
2053  int k = random(3,10);
2054  poly g = x4 + x2*y^(2*k+1)+x*y^(3*k+1)+ y^(4*k +1);
2055  arnoldClassify_to_string(g);
2056}
2057
Note: See TracBrowser for help on using the repository browser.