source: git/Singular/LIB/realclassify.lib @ 2473c9

spielwiese
Last change on this file since 2473c9 was 2473c9, checked in by Hans Schoenemann <hannes@…>, 5 years ago
fix: shadowing variables (k:int/poly, i:int/number, n:int/number) in realclassify.lib
  • Property mode set to 100644
File size: 153.5 KB
Line 
1////////////////////////////////////////////////////////////////////////////
2version="version realclassify.lib 4.1.2.0 Feb_2019 "; // $Id$
3category="Singularities";
4info="
5LIBRARY:  realclassify.lib   Classification of real singularities
6AUTHOR:   Janko Boehm,       boehm@mathematik.uni-kl.de
7          Magdaleen Marais,  magdaleen@aims.ac.za
8          Andreas Steenpass, steenpass@mathematik.uni-kl.de
9
10OVERVIEW:
11   A library for classifying isolated hypersurface singularities over the reals
12   w.r.t. right equivalence, based on the determinator of singularities by
13   V.I. Arnold. This library is based on classify2.lib by the first and second author and G. Pfister, but
14   handles the real case, while classify2.lib does the complex classification.
15
16REFERENCES:
17Arnold, Varchenko, Gusein-Zade: Singularities of Differentiable Maps.
18Vol. 1: The classification of critical points caustics and wave fronts.
19Birkh\"auser, Boston 1985
20
21J. Boehm, M.S. Marais, A. Steenpass: The Classification of Real Singularities Using Singular.
22Part III: Unimodal Singularities of Corank 2, https://arxiv.org/abs/1512.09028
23
24Greuel, Lossen, Shustin: Introduction to singularities and deformations.
25Springer, Berlin 2007
26
27M.S. Marais, A. Steenpass: The Classification of Real Singularities Using SINGULAR. Part I:
28Splitting Lemma and Simple Singularities, J. Symb. Comput. 68 (2015), 61-71
29
30M.S. Marais, A. Steenpass: The Classification of Real Singularities Using SINGULAR. Part II:
31The Structure of the Equivalence Classes of the Unimodal Singularities,
32J. Symb. Comput. 74 (2016), 346-366
33
34
35
36Acknowledgements: This research was supported by
37the Staff Exchange Bursary Programme of the University of Pretoria, DFG SPP 1489, and
38DFG TRR 195. The financial assistance of the National Research Foundation (NRF),
39South Africa, towards this research is hereby acknowledged. Opinions expressed
40and conclusions arrived at are those of the author and are not necessarily to be
41attributed to the National Research Foundation, South Africa.
42
43
44PROCEDURES:
45 realclassify(f);               real classification of singularities of modality 0 and 1 up to stable equivalence
46 realmorsesplit(f);             splitting lemma in the real case
47 milnornumber(f);               Milnor number
48 determinacy(f);                an upper bound for the determinacy
49 addnondegeneratevariables(f);  find a right equivalent normal form by adding the non-degenerate variables
50";
51
52LIB "linalg.lib";
53LIB "elim.lib";
54LIB "classify.lib";
55LIB "rootsur.lib";
56LIB "rootsmr.lib";
57LIB "atkins.lib";
58LIB "solve.lib";
59LIB "general.lib";
60LIB "classify2.lib";
61LIB "rootisolation.lib";
62static proc mod_init()
63{
64  option(noredefine);
65}
66
67///////////////////////////////////////////////////////////////////////////////
68proc realclassify(poly f, list #)
69"
70USAGE:    realclassify(f); f poly
71RETURN:   A list of elements NF of type NormalFormEquation with the following keys:
72          @* NF.singularityType    = the type of the singularity as a string
73          @* NF.normalFormEquation = an element F of type Poly with the normal form equation
74                                     (a representative of the stable equivalence class of f,
75                                      as chosen by Arnold)
76          @* NF.modality           = the modality of f
77          @* NF.parameters         = a list where each entry is a list consisting out of the parameter term in
78                                     the normalform equation and an interval
79          @* NF.corank             = the corank of f (of type int)
80          @* NF.inertiaIndex       = the interia index of f (of type int)
81          @* NF.milnorNumber       = the milnor number of f (of type int)
82          @* NF.determinacy        = an lower bound for the determinacy of f (of type int)
83
84
85NOTE:     The classification is done over the real numbers, so in contrast to
86          classify.lib, the signs of coefficients of monomials where even
87          exponents occur matter.
88          @* The ground field must be Q (the rational numbers). No field
89             extensions of any kind nor floating point numbers are allowed.
90          @* The monomial order must be local.
91          @* The input polynomial must be contained in maxideal(2) and must be
92             an isolated singularity of modality 0 or 1. The Milnor number is
93             checked for being finite.
94          @* In case of NF.modality=1, the parameter values are given as an element
95             of a number field QQ[a]/(minpoly), the coefficient ring of F in terms
96             of the unique root of minpoly of the ring F.in in the interval specified
97             in NF.parameters in conjunction with the parameter term.
98SEE ALSO: classify
99KEYWORDS: Classification of singularities
100EXAMPLE:  example realclassify; shows an example"
101{
102  /* auxiliary variables */
103  int i, j;
104
105  /* name for the basering */
106  def br1 = basering;
107  list br1vars=ringlist(br1)[2];
108  def br = basering;
109  list L = ringlist(br);
110  L[3][1][1] = "ds";
111  br = ring(L);
112  setring br;
113  poly f = fetch(br1,f);
114
115
116  /* error check */
117  if((charstr(br) != "QQ") and (charstr(br) != "0"))
118  {
119    ERROR("The ground field must be QQ (the rational numbers).");
120  }
121  int n = nvars(br);
122  for(i = 1; i <= n; i++)
123  {
124    if(var(i) > 1)
125    {
126      ERROR("The monomial order must be local.");
127    }
128  }
129  if(jet(f, 1) != 0)
130  {
131    ERROR("The input polynomial must be contained in maxideal(2).");
132  }
133
134  /* compute Milnor number before continuing the error check */
135  int mu = milnornumber(f);
136  /* continue error check */
137  if(mu < 1)
138  {
139    ERROR("The Milnor number of the input polynomial must be"+newline
140      +"positive and finite.");
141  }
142
143
144  /* apply splitting lemma */
145  list morse = realmorsesplit(f, mu);
146  int cr = morse[1];
147  int lambda = morse[2];
148  int d = morse[3];
149  poly rf = morse[4];
150
151  /* determine the type */
152  string typeofsing;
153  Poly nf;
154  list listNFE;
155  if(cr > 2)
156  {
157    ERROR("This case is not yet implemented.");
158  }
159  if(cr == 0)   // case A[1]
160  {
161    setring br1;
162    poly rf = fetch(br, rf);
163    typeofsing, nf = caseA1(rf, n);
164    int m =0;
165  }
166  if(cr == 1)   // case A[k], k > 1
167  {
168    setring br1;
169    poly rf = fetch(br, rf);
170    typeofsing, nf = caseAk(rf, n);
171    int m = 0;
172  }
173  if(cr == 2)
174  {
175    list LL = ringlist(br);
176    LL[2]=list(LL[2][1],LL[2][2]);
177    def Rcc = ring(LL);
178    setring Rcc;
179    poly rf = fetch(br,rf);
180    list dataFromClassify = complexType(rf);
181    int m = dataFromClassify[2];                // the modality of f
182    string complextype = dataFromClassify[1];   // the complex type of f
183    if(m > 1)
184    {
185     ERROR("The input polynomial must be a singularity of modality 0 or 1.");
186    }
187    if(complextype[1,2] == "D[")   // case D[k]
188    {
189      if(mu == 4)   // case D[4]
190      {
191        setring br1;
192        poly rf = fetch(br, rf);
193        typeofsing, nf = caseD4(rf);
194      }
195      else   // case D[k], k > 4
196      {
197        setring br1;
198        poly rf = fetch(br, rf);
199        typeofsing, nf = caseDk(rf, mu);
200      }
201    }
202    if(complextype == "E[6]")   // case E[6]
203    {
204      setring br1;
205      poly rf = fetch(br, rf);
206      typeofsing, nf = caseE6(rf);
207    }
208    if(complextype == "E[7]")   // case E[7]
209    {
210      setring br1;
211      poly rf = fetch(br, rf);
212      typeofsing, nf = caseE7(rf);
213    }
214    if(complextype == "E[8]")   // case E[8]
215    {
216      setring br1;
217      poly rf = fetch(br, rf);
218      typeofsing, nf = caseE8(rf);
219    }
220    if(complextype == "Z[11]")  // case Z[11]
221    {
222      listNFE = caseZ11(rf,br1);
223    }
224    if(complextype == "Z[12]")  // case Z[12]
225    {
226      listNFE = caseZ12(rf, br1);
227    }
228    if(complextype == "Z[13]")  // case Z[12]
229    {
230      listNFE = caseZ13(rf,br1);
231    }
232    if(complextype == "W[12]")  // case W[12]
233    {
234      listNFE = caseW12(rf,br1);
235    }
236    if(complextype == "W[13]")  // case W[13]
237    {
238      listNFE = caseW13(rf,br1);
239    }
240    if(complextype == "E[12]")  // case E[12]
241    {
242      listNFE = caseE12(rf, br1);
243    }
244    if(complextype == "E[13]")  // case E[13]
245    {
246      listNFE = caseE13(rf, br1);
247    }
248    if(complextype == "E[14]")  // case E[14]
249    {
250      listNFE = caseE14(rf, br1);
251    }
252    if((complextype[1,2] == "X[") and (complextype <> "X[9]"))  // case X9k
253    {
254      listNFE = caseX9k(rf,mu,br1);
255    }
256    if((complextype[1,2] == "J[") and (complextype<>"J[10]"))  // case J[10+k]
257    {
258      listNFE = caseJ10k(rf,mu, br1);
259    }
260    if(complextype == "J[10]")   // case J[10]
261    {
262      listNFE = caseJ10(rf, br1);
263    }
264    if(complextype == "X[9]")  // case X9
265    {
266      listNFE = caseX9(rf, br1);
267    }
268    if(complextype[1,2] == "Y[")
269    {
270      listNFE = caseYrs(rf,mu, br1);
271    }
272  }
273  if (m==0){
274       NormalFormEquation NFE;
275       NFE.normalFormEquation = nf;
276       NFE.parameters = list();
277       NFE.singularityType = typeofsing;
278       listNFE = list(NFE);
279  }
280
281
282  if (size(listNFE)==0)
283  {
284      ERROR("This case is not yet implemented.");
285  }
286
287
288    for (i=1;i<=size(listNFE);i++){
289    listNFE[i].vars = br1vars;
290    listNFE[i].realCase = 1;
291    listNFE[i].corank = cr;
292    listNFE[i].inertiaIndex = lambda;
293    listNFE[i].determinacy = d;
294    listNFE[i].modality = m;
295    listNFE[i].determinacy = d;
296    listNFE[i].milnorNumber = mu;
297}
298
299return(listNFE);
300
301}
302example
303{
304  "EXAMPLE:";
305  echo = 2;
306  ring r = 0, (x,y,z), ds;
307  poly f = (x2+3y-2z)^2+xyz-(x-y3+x2z3)^3;
308  realclassify(f);
309  map phi = r, x+2y+y^2+x*y,x+y+y^2+x^2,z;
310
311  poly f  = x^2+y^2-z^2;  // A[1]
312  f=phi(f);
313  realclassify(f);
314
315  poly f  = x^3+y^2-z^2;  // A[2]
316  f=phi(f);
317  realclassify(f);
318
319  poly f  = x^2*y-y^3+z^2;  // D[4]-
320  f=phi(f);
321  realclassify(f);
322
323
324  poly f  = x^3-y^4-z^2;  // E[6]-
325  f=phi(f);
326  realclassify(f);
327
328  poly f  = x^3+x*y^3+z^2;  // E[7]
329  f=phi(f);
330  realclassify(f);
331
332  poly f  = x^3+y^5-z^2;  // E[8]
333  f=phi(f);
334  realclassify(f);
335
336  poly f  = x^3+3*x^2*y^2+x*y^4-z^2;  //J[10]+
337  f=phi(f);
338  realclassify(f);
339
340  poly f  = x^3+x^2*y^2+4*y^9+z^2;  //J[13]+
341  f=phi(f);
342  realclassify(f);
343
344  poly f  = -x^4-x^2*y^2+3y^9-z^2;  //X[14]--
345  f=phi(f);
346  realclassify(f);
347
348  poly f  = -x^2*y^2-x^7+4*y^8+z^2;  //Y[7,8]--
349  f=phi(f);
350  realclassify(f);
351
352  poly f  = (x^2+y^2)^2+5*x^9-z^2;  // tilde Y[9]
353  f=phi(f);
354  realclassify(f);
355
356  poly f  = x^3+y^7+3*x*y^5+z^2;  // E[12]
357  f=phi(f);
358  realclassify(f);
359
360  poly f  = x^3+x*y^5+4*y^8-z^2;  // E[13]
361  f=phi(f);
362  realclassify(f);
363
364  poly f  = x^3+y^8+2*x*y^6+z^2;  // E[14]+
365  f=phi(f);
366  realclassify(f);
367
368  poly f  = x^3*y+y^5+5*x*y^4-z^2;  // Z[11]
369  f=phi(f);
370  realclassify(f);
371
372  poly f  = x^3*y+x*y^4+6*x^2*y^3+z^2;  // Z[12]
373  f=phi(f);
374  realclassify(f);
375
376  poly f  = x^3*y-y^6+2*x*y^5-z^2;  // Z[13]-
377  f=phi(f);
378  realclassify(f);
379
380  poly f  = x^4+y^5+x^2*y^3+z^2;  // W[12]+
381  f=phi(f);
382  realclassify(f);
383
384
385  poly f  = -x^4+x*y^4+y^6-z^2;  // W[13]-
386  f=phi(f);
387  realclassify(f);
388
389  poly p = x^4-x^2*y^2+5*y^4+x*y^3+x^3*y+z^2;  //X9++
390  f=phi(f);
391  realclassify(p);
392
393}
394////////////////////////////////////////////////////////////////////////////////
395static proc coeff(poly f, poly m){
396poly p = f;
397int i = 1;
398while (p<>0){
399  if (m==leadmonom(p)){return(leadcoef(p));}
400  p=p-lead(p);
401}
402return(0);}
403///////////////////////////////////////////////////////////////////////////////
404static proc maximumFiltration(poly f)
405{
406  def br = basering;
407  ring R = 0,(x,y),dp;
408  poly f = fetch(br,f);
409  int n = ord(f);
410  setring br;
411  return(ord(f));
412}
413//////////////////////////////////////////////////////////////////////////////////////
414static proc WeightedJet(poly f, intvec v, int k)
415{
416  def br = basering;
417  ring R = 0,(x,y),wp(v);
418  poly f = fetch(br,f);
419  while(deg(f,v)>k)
420  {
421    f=f-lead(f);
422  }
423  setring br;
424  poly jetf = fetch(R,f);
425  return(jetf)
426}
427//////////////////////////////////////////////////////////////////////////////////////
428static proc lowestJet(poly rf)   //Apply to cases of type: exceptional, J[10] and J[10+k].
429{
430  map phi;
431  def br = basering;
432  int n = maximumFiltration(rf);
433  poly g = jet(rf,n);
434  list factorization = factorize(g);
435  list factors = factorization[1][2];
436  poly f = factorization[1][2];
437  if(size(factorization[1])==3)
438  {
439    if(factorization[2][2] > factorization[2][3])
440    {
441      factors = insert(factors, factorization[1][3], 1);
442    }
443    else
444    {
445      factors = insert(factors, factorization[1][3]);
446    }
447  matrix T[2][2] = factors[1]/var(1), factors[1]/var(2),
448                   factors[2]/var(1), factors[2]/var(2);
449  phi = br, luinverse(T)[2]*matrix(ideal(var(1), var(2)), 2, 1);
450  rf = phi(rf);
451  }
452  else
453  {
454    if(coeff(f,var(1))<>0)
455    {
456      matrix T[2][2] = factors[1]/var(1),factors[1]/var(2),
457                       0, 1;
458      phi = br, luinverse(T)[2]*matrix(ideal(var(1), var(2)), 2, 1);
459      rf = phi(rf);
460    }
461    else
462    {
463      phi = br,var(2),var(1);
464      rf = phi(rf);
465    }
466  }
467  return(rf)
468}
469///////////////////////////////////////////////////////////////////////////////
470static proc RealRootFactor(poly MP, number k1, number k2)
471{
472  def br = basering;
473  ring R = 0,x,dp;
474  number k1 = fetch(br,k1);
475  number k2 = fetch(br,k2);
476  poly MP = fetch(br,MP);
477  list factorization =factorize(MP);
478  while(size(factorization[1])>2 || factorization[2][2]>1)
479  {
480    for(int i=2;i<=size(factorization[1]);i=i+1)
481    {
482      if(sturm(factorization[1][i],k1,k2)>0)
483      {
484        MP=factorization[1][i];
485        factorization =factorize(MP);
486      }
487    }
488  }
489  int n = deg(MP);
490  setring br;
491  MP = fetch(R,MP);
492  number c = (coeff(MP,var(1)^n))^(-1);
493  MP = c*MP;
494  return(MP)
495}
496///////////////////////////////////////////////////////////////////////////////
497static proc eliminateTildeY(poly rf, poly q, int n1)
498{
499  def br = basering;
500  int i = 5;
501  poly p = (jet(rf,i)-jet(rf,4))/(coeff(rf,var(1)^4)*4*q);
502  while(i<>n1)
503  {
504    rf = eliminateDegreeTildeY(rf,p,q,i,n1);
505    poly rf1 = jet(rf,i)-jet(rf,4);
506    if(rf1 <> 0)
507    {
508      ERROR("this remainder should be zero");
509    }
510    i = i+1;
511    p = (jet(rf,i)-jet(rf,4))/(coeff(rf,var(1)^4)*4*q);
512  }
513  rf = eliminateDegreeTildeY(rf,p,q,i,n1);
514  return(rf);
515}
516///////////////////////////////////////////////////////////////////////////////
517static proc eliminateDegreeTildeY(poly rf, poly p, poly q, int i, int n1)
518{
519   poly px = subst(p,var(2),0);
520   poly py = p-px;
521   map phi = br, var(1)-px/var(1), var(2)- py/(coeff(q,var(2)^2)*var(2));
522   rf  = phi(rf);
523   return(jet(rf,n1));
524}
525///////////////////////////////////////////////////////////////////////////////
526static proc caseA1(poly rf, int n)
527{
528  string typeofsing = "A[1]";
529  poly nf = 0;
530  return(typeofsing, nf);
531}
532
533///////////////////////////////////////////////////////////////////////////////
534static proc caseAk(poly rf, int n)
535{
536  /* preliminaries */
537  string typeofsing;
538  poly nf;
539
540  int k = deg(lead(rf), 1:n)-1;
541  if(k%2 == 0)
542  {
543    nf = var(1)^(k+1);
544    typeofsing = "A["+string(k)+"]";
545  }
546  else
547  {
548    if(leadcoef(rf) > 0)
549    {
550      nf = var(1)^(k+1);
551      typeofsing = "A["+string(k)+"]+";
552    }
553    else
554    {
555      nf = -var(1)^(k+1);
556      typeofsing = "A["+string(k)+"]-";
557    }
558  }
559  return(typeofsing, nf);
560}
561
562///////////////////////////////////////////////////////////////////////////////
563static proc caseD4(poly rf)
564{
565  /* preliminaries */
566  string typeofsing;
567  poly nf;
568  def br = basering;
569  map phi;
570
571  rf = jet(rf, 3);
572  number s1 = number(rf/(var(1)^3));
573  number s2 = number(rf/(var(2)^3));
574  if(s2 == 0 && s1 != 0)
575  {
576    phi = br, var(2), var(1);
577    rf = phi(rf);
578  }
579  if(s1 == 0 && s2 == 0)
580  {
581    number t1 = number(rf/(var(1)^2*var(2)));
582    number t2 = number(rf/(var(2)^2*var(1)));
583    if(t1+t2 == 0)
584    {
585      phi = br, var(1)+2*var(2), var(2);
586      rf = phi(rf);
587    }
588    else
589    {
590      phi = br, var(1)+var(2), var(2);
591      rf = phi(rf);
592    }
593  }
594  ring R = 0, y, dp;
595  map phi = br, 1, y;
596  poly rf = phi(rf);
597  int k = nrroots(rf);
598  setring(br);
599  if(k == 3)
600  {
601    nf = var(1)^2*var(2)-var(2)^3;
602    typeofsing = "D[4]-";
603  }
604  else
605  {
606    nf = var(1)^2*var(2)+var(2)^3;
607    typeofsing = "D[4]+";
608  }
609  return(typeofsing, nf);
610}
611
612///////////////////////////////////////////////////////////////////////////////
613static proc caseDk(poly rf, int mu)
614{
615  /* preliminaries */
616  string typeofsing;
617  poly nf;
618  def br = basering;
619  map phi;
620  poly MP;
621
622  rf = jet(rf, mu-1);
623  list factorization = factorize(jet(rf, 3));
624  list factors = factorization[1][2];
625  if(factorization[2][2] == 2)
626  {
627    factors = insert(factors, factorization[1][3], 1);
628  }
629  else
630  {
631    factors = insert(factors, factorization[1][3]);
632  }
633  factors[2] = factorization[1][1]*factors[2];
634  matrix T[2][2] = factors[1]/var(1), factors[1]/var(2),
635         factors[2]/var(1), factors[2]/var(2);
636  phi = br, luinverse(T)[2]*matrix(ideal(var(1), var(2)), 2, 1);
637  rf = phi(rf);
638  rf = jet(rf, mu-1);
639  poly g;
640  int i;
641  for(i = 4; i < mu; i++)
642  {
643    g = jet(rf, i) - var(1)^2*var(2);
644    if(g != 0)
645    {
646      phi = br, var(1)-(g/(var(1)*var(2)))/2,
647          var(2)-(g/var(1)^i)*var(1)^(i-2);
648      rf = phi(rf);
649      rf = jet(rf, mu-1);
650    }
651  }
652  number a = number(rf/var(2)^(mu-1));
653  if(a > 0)
654  {
655    typeofsing = "D["+string(mu)+"]+";
656    nf = var(1)^2*var(2)+var(2)^(mu-1);
657  }
658  else
659  {
660    typeofsing = "D["+string(mu)+"]-";
661    nf = var(1)^2*var(2)-var(2)^(mu-1);
662  }
663  return(typeofsing, nf);
664}
665
666///////////////////////////////////////////////////////////////////////////////
667static proc caseE6(poly rf)
668{
669  /* preliminaries */
670  string typeofsing;
671  poly nf;
672  def br = basering;
673  map phi;
674
675  poly g = jet(rf,3);
676  number s = number(g/(var(1)^3));
677  if(s == 0)
678  {
679    phi = br, var(2), var(1);
680    rf = phi(rf);
681    g = jet(rf,3);
682  }
683  list Factors = factorize(g);
684  poly g1 = Factors[1][2];
685  phi = br, (var(1)-(g1/var(2))*var(2))/(g1/var(1)), var(2);
686  rf = phi(rf);
687  rf = jet(rf,4);
688  number w = number(rf/(var(2)^4));
689  if(w > 0)
690  {
691    typeofsing = "E[6]+";
692    nf = var(1)^3+var(2)^4;
693  }
694  else
695  {
696    typeofsing = "E[6]-";
697    nf = var(1)^3-var(2)^4;
698  }
699  return(typeofsing, nf);
700}
701
702///////////////////////////////////////////////////////////////////////////////
703static proc caseE7(poly rf)
704{
705  string typeofsing = "E[7]";
706  poly nf = var(1)^3+var(1)*var(2)^3;
707  return(typeofsing, nf);
708}
709
710///////////////////////////////////////////////////////////////////////////////
711static proc caseE8(poly rf)
712{
713  string typeofsing = "E[8]";
714  poly nf = var(1)^3+var(2)^5;
715  return(typeofsing, nf);
716}
717///////////////////////////////////////////////////////////////////////////////
718static proc caseZ11(poly rf, br1)
719{
720  /*preliminaries*/
721  string typeofsing;
722  def br = basering;
723  poly nf;
724  map phi;
725  poly MP;
726  list intervalP;
727
728  /* II. Eliminate the monomials in supp(rf) underneath Gamma(T) */
729  rf = lowestJet(rf);
730
731  /* IV. Read off the desired information */
732  number b = 1/(coeff(rf,var(1)^3*var(2)));
733  phi = br, var(1),b*var(2);
734  rf = phi(rf);
735  b=coeff(rf,var(1)*var(2)^4)^(15)*coeff(rf,var(2)^5)^(-11);
736  MP = var(1)^(15)-b;
737  number n = maxabs(MP);
738  MP = RealRootFactor(MP,-(n+1),n+1);
739  n = maxabs(MP);
740  if(deg(MP)==1)
741  {
742    number a = -number(subst(MP,var(1),0));
743    setring br1;
744    number a = fetch(br,a);
745    number n = fetch(br,n);
746  }
747  else
748  {
749    string mord = ringlist(br1)[3][1][1];
750    ring R1 = (0,a),(var(1),var(2)),ds;
751    list L1 = ringlist(R1);
752    L1[3][1][1]=mord;
753    def R2 = ring(L1);
754    setring R2;
755    kill R1;
756    poly ff = fetch(br,MP);
757    minpoly = number(subst(ff,var(1),a));
758    number n = fetch(br,n);
759  }
760  Poly nfe = var(1)^3+var(2)^5+a*var(1)*var(2)^4;
761  intervalP = list(-n-1,n+1);
762  typeofsing = "Z[11]";
763  NormalFormEquation F;
764  F.normalFormEquation = nfe;
765  F.parameters = list(list(a*var(1)*var(2)^4,intervalP));
766  F.singularityType = typeofsing;
767  setring br;
768  return(list(F));
769}
770///////////////////////////////////////////////////////////////////////////////
771static proc caseZ12(poly rf, br1)
772{
773  /*preliminaries*/
774  string typeofsing;
775  poly nf;
776  def br = basering;
777  map phi;
778  poly MP;
779  list intervalP;
780
781  /* II. Eliminate the monomials in supp(rf) underneath Gamma(T) */
782  rf = lowestJet(rf);
783
784  /* III. Eliminate the monomials in supp(rf) above \Gamma(T) which are not in supp(T)*/
785  /* Write the sum of the terms of rf above \Gamma(T) in terms of Arnold's system */
786  rf = rf-coeff(rf,var(2)^6)*var(2)^6-coeff(rf,var(2)^6)*3*(coeff(rf,var(1)*var(2)^4)^(-1))*var(1)^2*var(2)^3;
787  /* IV. Read off the desired information */
788  phi = br, var(1),(1/(coeff(rf,var(1)^3*var(2))))*var(2);
789  rf = phi(rf);
790  number b=(coeff(rf,var(1)^2*var(2)^3))^(11)*(coeff(rf,var(1)*var(2)^4))^(-7);
791  MP = var(1)^(11)-b;
792
793  number n = maxabs(MP);
794  MP = RealRootFactor(MP,-(n+1),n+1);
795  n = maxabs(MP);
796  if(deg(MP)==1)
797  {
798    number a = -number(subst(MP,var(1),0));
799    setring br1;
800    number a = fetch(br,a);
801    number n = fetch(br,n);
802  }
803  else
804  {
805    string mord = ringlist(br1)[3][1][1];
806    ring R1 = (0,a),(var(1),var(2)),ds;
807    list L1 = ringlist(R1);
808    L1[3][1][1]=mord;
809    def R2 = ring(L1);
810    setring R2;
811    kill R1;
812    poly ff = fetch(br,MP);
813    minpoly = number(subst(ff,var(1),a));
814    number n = fetch(br,n);
815  }
816  Poly nfe =var(1)^3*var(2)+var(1)*var(2)^4+a*var(1)^2*var(2)^3;
817  intervalP = list(-n-1,n+1);
818  typeofsing = "Z[12]";
819  NormalFormEquation F;
820  F.normalFormEquation = nfe;
821  F.parameters = list(list(a*var(1)^2*var(2)^3,intervalP));
822  F.singularityType = typeofsing;
823  setring br;
824  return(list(F));
825}
826///////////////////////////////////////////////////////////////////////////////
827static proc caseZ13(poly rf, br1)
828{
829  /*preliminaries*/
830  string typeofsing;
831  poly nf;
832  def br = basering;
833  map phi;
834  poly MP;
835  list intervalP;
836
837  /* II. Eliminate the monomials in supp(rf) underneath or on \Gamma(T)*/
838  rf = lowestJet(rf);
839
840  /* III. Eliminate the monomials above \Gamma(T) which are not in supp(T)*/
841  number c = coeff(rf,var(1)^3*var(2));
842  poly g = 1/(coeff(rf,var(1)^3*var(2)))*rf;
843  intvec v = 5,3;
844
845  /* Consider the terms of g of w-degree j=19 and Apply Algorithm 4*/
846  poly h = (WeightedJet(g,v,19)-WeightedJet(g,v,18))/(3*var(1)^2*var(2));
847  phi = br,var(1)-h,var(2);
848  rf = phi(rf);
849
850  /* IV. Read off the desired information*/
851  phi = br, var(1), (1/c)*var(2);
852  rf = phi(rf);
853  number s = coeff(rf,var(2)^6);
854  int sg = sign(s);
855  if(coeff(rf,var(2)^6)>0)
856  {
857    typeofsing = "Z[13]+";
858  }
859  else
860  {
861    typeofsing = "Z[13]-";
862  }
863  number b = (coeff(rf,var(1)*var(2)^5))^9*(absValue(coeff(rf,var(2)^6)))^(-7);
864  MP = var(1)^(9)-b;
865  number n = maxabs(MP);
866  MP = RealRootFactor(MP,-(n+1),n+1);
867  n = maxabs(MP);
868  if(deg(MP)==1)
869  {
870    number a = -number(subst(MP,var(1),0));
871    setring br1;
872    number a = fetch(br,a);
873    number n = fetch(br,n);
874  }
875  else
876  {
877    string mord = ringlist(br1)[3][1][1];
878    ring R1 = (0,a),(var(1),var(2)),ds;
879    list L1 = ringlist(R1);
880    L1[3][1][1]=mord;
881    def R2 = ring(L1);
882    setring R2;
883    kill R1;
884    poly ff = fetch(br,MP);
885    minpoly = number(subst(ff,var(1),a));
886    number n = fetch(br,n);
887  }
888  Poly nfe =var(1)^3*var(2)+number(sg)*var(2)^6+a*var(1)*var(2)^5;
889  intervalP = list(-n-1,n+1);
890  NormalFormEquation F;
891  F.normalFormEquation = nfe;
892  F.parameters = list(list(a*var(1)*var(2)^5,intervalP));
893  F.singularityType = typeofsing;
894  setring br;
895  return(list(F));
896}
897
898///////////////////////////////////////////////////////////////////////////////
899static proc caseW12(poly rf, br1)
900{
901  /*preliminaries*/
902  string typeofsing;
903  poly nf;
904  def br = basering;
905  map phi;
906  poly MP;
907  list intervalP;
908
909  /* II. Eliminate the monomials in supp(rf) underneath or on \Gamma(T)*/
910  rf = lowestJet(rf);
911
912  /* III. Eliminate the monomials above \Gamma(T) which are not in supp(T)*/
913  poly g = 1/(coeff(rf,var(2)^5))*rf;
914  intvec v = 5,4;
915
916  /* Consider the terms of g of w-degree j=21 and Apply Algorithm 4 */
917  poly h = (WeightedJet(g,v,21)-WeightedJet(g,v,20))/(5*var(2)^4);
918  phi = br,var(1),var(2)-h;
919  rf = phi(rf);
920
921  /* IV. Read off the desired information*/
922  if(coeff(rf,var(2)^5)<0)
923  {
924    phi = br, var(1),-var(2);
925    rf = phi(rf);
926  }
927  number s = coeff(rf,var(1)^4);
928  int sg = sign(s);
929  if(coeff(rf,var(1)^4)>0)
930  {
931    typeofsing = "W[12]+";
932  }
933  else
934  {
935    typeofsing = "W[12]-";
936  }
937  MP = var(1)^(10)-coeff(rf,var(1)^2*var(2)^3)^(10)*(absValue(coeff(rf,var(1)^4)))^(-5)*coeff(rf,var(2)^5)^(-6);
938  if(coeff(rf,var(1)^2*var(2)^3)>0)
939  {
940    number n = maxabs(MP);
941    MP = RealRootFactor(MP,0,n+1);
942    number n1 = 0;
943    number n2 = maxabs(MP)+1;
944  }
945  else
946  {
947    number n = maxabs(MP);
948    MP = RealRootFactor(MP,-(n+1),0);
949    number n1 = -maxabs(MP)-1;
950    number n2 = 0;
951  }
952
953 if(deg(MP)==1)
954  {
955    number a = -number(subst(MP,var(1),0));
956    setring br1;
957    number a = fetch(br,a);
958    number n = fetch(br,n);
959  }
960  else
961  {
962    string mord = ringlist(br1)[3][1][1];
963    ring R1 = (0,a),(var(1),var(2)),ds;
964    list L1 = ringlist(R1);
965    L1[3][1][1]=mord;
966    def R2 = ring(L1);
967    setring R2;
968    kill R1;
969    poly ff = fetch(br,MP);
970    minpoly = number(subst(ff,var(1),a));
971  }
972  number n1 = fetch(br,n1);
973  number n2 = fetch(br,n2);
974  Poly nfe =number(sg)*var(1)^4+var(2)^5+a*var(1)^2*var(2)^3;
975  intervalP = list(n1,n2);
976  NormalFormEquation F;
977  F.normalFormEquation = nfe;
978  F.parameters = list(list(a*var(1)^2*var(2)^3,intervalP));
979  F.singularityType = typeofsing;
980  setring br;
981  return(list(F));
982}
983///////////////////////////////////////////////////////////////////////////////
984static proc caseW13(poly rf, br1)
985{
986  /*preliminaries*/
987  string typeofsing;
988  poly nf;
989  def br = basering;
990  map phi;
991  poly MP;
992  list intervalP;
993
994  /* II. Eliminate the monomials in supp(rf) underneath or on \Gamma(T) */
995  rf = lowestJet(rf);
996
997  /* III. Eliminate the monomials above \Gamma(T) which are not in supp(T)*/
998  poly g = 1/(coeff(rf,var(1)*var(2)^4))*rf;
999  intvec v = 4,3;
1000
1001  /* Consider the terms of w-degree j=17 and Apply Algorithm 4*/
1002  poly h = (WeightedJet(g,v,17)-WeightedJet(g,v,16))/(4*var(1)*var(2)^3);
1003  phi = br,var(1),var(2)-h;
1004  rf = phi(rf);
1005
1006  /* Write the sum of the terms of rf above \Gamma(T) in terms of Arnold's system */
1007  number b = 1/(4*coeff(rf,var(1)^4));
1008  rf = rf-coeff(rf,var(1)^3*var(2)^2)*var(1)^3*var(2)^2-b*coeff(rf,var(1)^3*var(2)^2)*var(2)^6;
1009
1010  /* Read off the desired information*/
1011  phi = br,coeff(rf,var(1)*var(2)^4)^(-(1))*var(1),var(2);
1012  rf = phi(rf);
1013  number s = coeff(rf,var(1)^4);
1014  int sg = sign(s);
1015  if(coeff(rf,var(1)^4)>0)
1016  {
1017      typeofsing = "W[13]+";
1018  }
1019  else
1020  {
1021    typeofsing = "W[13]-";
1022  }
1023  MP = var(1)^(8)-coeff(rf,var(2)^6)^8*(absValue(coeff(rf,var(1)^4)))^3;
1024  if(coeff(rf,var(2)^6)>0)
1025  {
1026    number n = maxabs(MP);
1027    MP = RealRootFactor(MP,0,n+1);
1028    number n1 = 0;
1029    number n2 = maxabs(MP);
1030  }
1031  else
1032  {
1033    number n = maxabs(MP);
1034    MP = RealRootFactor(MP,-(n+1),0);
1035    number n1 = -maxabs(MP)-1;
1036    number n2 = 0;
1037  }
1038
1039  if(deg(MP)==1)
1040  {
1041    number a = -number(subst(MP,var(1),0));
1042    setring br1;
1043    number a = fetch(br,a);
1044    number n = fetch(br,n);
1045  }
1046  else
1047  {
1048    string mord = ringlist(br1)[3][1][1];
1049    ring R1 = (0,a),(var(1),var(2)),ds;
1050    list L1 = ringlist(R1);
1051    L1[3][1][1]=mord;
1052    def R2 = ring(L1);
1053    setring R2;
1054    kill R1;
1055    poly ff = fetch(br,MP);
1056    minpoly = number(subst(ff,var(1),a));
1057  }
1058  number n1 = fetch(br,n1);
1059  number n2 = fetch(br,n2);
1060  Poly nfe =sg*var(1)^4+var(1)*var(2)^4+a*var(2)^6;
1061  intervalP = list(n1,n2);
1062  NormalFormEquation F;
1063  F.normalFormEquation = nfe;
1064  F.parameters = list(list(a*var(2)^6,intervalP));
1065  F.singularityType = typeofsing;
1066  setring br;
1067  return(list(F));
1068  }
1069////////////////////////////////////////////////////////////////////////////////
1070static proc caseE12(poly rf, br1)
1071{
1072  /*preliminaries*/
1073  string typeofsing;
1074  poly nf;
1075  def br = basering;
1076  map phi;
1077  poly MP;
1078  list intervalP;
1079
1080  /* II. Eliminate the monomials in supp(rf) underneath or on \Gamma(T)*/
1081  /* Apply Algorithm 1*/
1082  rf = lowestJet(rf);
1083
1084  poly g1 = coeff(rf,var(1)^3)^(-1)*rf;
1085
1086  /* Consider the terms of g1 of standard degree j=4 and Apply Algorithm 4*/
1087  poly g2 = (jet(g1,4)-var(1)^3)/(3*var(1)^2);
1088  phi = br, var(1)-g2,var(2);
1089  g1 = phi(g1);
1090
1091  /* Consider the terms of g1 of standard degree j=5 and Apply Algorithm 4*/
1092  poly g3 =  (jet(g1,5)-var(1)^3)/(3*var(1)^2);
1093  phi = br, var(1)-g3, var(2);
1094  g1 = phi(g1);
1095  rf = coeff(rf,var(1)^3)*g1;
1096  /* IV. Read off the desired information */
1097  typeofsing = "E[12]";
1098  MP = var(1)^(21)-coeff(rf,var(1)^3)^(-7)*coeff(rf,var(2)^7)^(-15)*coeff(rf,var(1)*var(2)^5)^(21);
1099  number n = maxabs(MP);
1100  MP = RealRootFactor(MP,-(n+1),n+1);
1101  n = maxabs(MP);
1102
1103  if(deg(MP)==1)
1104  {
1105    number a = -number(subst(MP,var(1),0));
1106    setring br1;
1107    number a = fetch(br,a);
1108  }
1109  else
1110  {
1111    string mord = ringlist(br1)[3][1][1];
1112    ring R1 = (0,a),(var(1),var(2)),ds;
1113    list L1 = ringlist(R1);
1114    L1[3][1][1]=mord;
1115    def R2 = ring(L1);
1116    setring R2;
1117    kill R1;
1118    poly ff = fetch(br,MP);
1119    minpoly = number(subst(ff,var(1),a));
1120  }
1121  Poly nfe =var(1)^3+var(2)^7+a*var(1)*var(2)^5;
1122  number n = fetch(br,n);
1123  intervalP = list(-(n+1),n+1);
1124  NormalFormEquation F;
1125  F.normalFormEquation = nfe;
1126  F.parameters = list(list(a*var(1)*var(2)^5,intervalP));
1127  F.singularityType = typeofsing;
1128  setring br;
1129  return(list(F));
1130}
1131////////////////////////////////////////////////////////////////////////////////
1132static proc caseE13(poly rf, br1)
1133{
1134  /*preliminaries*/
1135  string typeofsing;
1136  poly nf;
1137  def br = basering;
1138  map phi;
1139  poly MP;
1140  list intervalP;
1141
1142  /* II. Eliminate the monomials in supp(rf) underneath or on \Gamma(T)*/
1143  /* Apply Algorithm 1*/
1144  rf = lowestJet(rf);
1145
1146  poly g1 = coeff(rf,var(1)^3)^(-1)*rf;
1147
1148  /* Consider the terms of g of standard degree j=4 and Apply Algorithm 4*/
1149  poly g2 = (jet(g1,4)-var(1)^3)/(3*var(1)^2);
1150  phi = br, var(1)-g2,var(2);
1151  g1 = phi(g1);
1152
1153  /* Consider the terms of g of standard degree j=5 and Apply Algorithm 4*/
1154  poly g3 =(jet(g1,5)-var(1)^3)/(3*var(1)^2);
1155  phi = br, var(1)-g3, var(2);
1156  g1 = phi(g1);
1157
1158  rf = coeff(rf,var(1)^3)*g1;
1159
1160  /* Read off the desired information */
1161  phi  = br, coeff(rf,var(1)*var(2)^5)^(-1)*var(1), var(2);
1162  rf = phi(rf);
1163  if(coeff(rf,var(1)^3)<0)
1164  {
1165    phi = br, -var(1),-var(2);
1166    rf = phi(rf);
1167  }
1168  rf = rf-coeff(rf,var(1)^2*var(2)^3)*(var(1)^2*var(2)^3+(3*coeff(rf,var(1)^3))^(-1)*var(2)^(8));
1169  typeofsing = "E[13]";
1170  MP = var(1)^(15)-coeff(rf,var(1)^3)^8*coeff(rf,var(2)^8)^(15);
1171  number n = maxabs(MP);
1172  MP = RealRootFactor(MP,-(n+1),n+1);
1173
1174  n = maxabs(MP);
1175  if(deg(MP)==1)
1176  {
1177    number a = -number(subst(MP,var(1),0));
1178    setring br1;
1179    number a = fetch(br,a);
1180  }
1181  else
1182  {
1183    string mord = ringlist(br1)[3][1][1];
1184    ring R1 = (0,a),(var(1),var(2)),ds;
1185    list L1 = ringlist(R1);
1186    L1[3][1][1]=mord;
1187    def R2 = ring(L1);
1188    setring R2;
1189    kill R1;
1190    poly ff = fetch(br,MP);
1191    minpoly = number(subst(ff,var(1),a));
1192  }
1193  number n = fetch(br,n);
1194  Poly nfe =var(1)^3+var(1)*var(2)^5+a*var(2)^8;
1195  intervalP = list(-(n+1),n+1);
1196  NormalFormEquation F;
1197  F.normalFormEquation = nfe;
1198  F.parameters = list(list(a*var(2)^8,intervalP));
1199  F.singularityType = typeofsing;
1200  setring br;
1201  return(list(F));
1202}
1203
1204////////////////////////////////////////////////////////////////////////////////////////////////
1205static proc caseE14(poly rf, br1)
1206{
1207  /*preliminaries*/
1208  string typeofsing = "E[14]";
1209  poly nf;
1210  def br = basering;
1211  map phi;
1212  poly MP;
1213  list intervalP;
1214
1215  /* II. Eliminate the monomials in supp(f) underneath or on \Gamma(T)*/
1216  /* Apply Algorithm 1*/
1217  rf = lowestJet(rf);
1218
1219  poly g1 = coeff(rf,var(1)^3)^(-1)*rf;
1220
1221  /* Consider the terms of g1 of standard degree j=4 and Apply Algorithm 4*/
1222  poly g2 = (jet(g1,4)-var(1)^3)/(3*var(1)^2);
1223  phi = br, var(1)-g2,var(2);
1224  g1 = phi(g1);
1225
1226  /*Consider the terms of g1 of standard degree j=5 and Apply Algorithm 4*/
1227  poly g3 =  (jet(g1,5)-var(1)^3)/(3*var(1)^2);
1228  phi = br, var(1)-g3, var(2);
1229  g1 = phi(g1);
1230
1231  /*Eliminate the monomials above \Gamma(T) which are not in supp(T)*/
1232
1233  intvec v = 8,3;
1234  /*Consider the terms of g1 of w-degree j=25 and Apply Algorithm 4*/
1235  poly h = (WeightedJet(g1,v,25)-WeightedJet(g1,v,24))/(3*var(1)^2);
1236  phi = br,var(1),var(2)-h;
1237  g1 = phi(g1);
1238
1239  rf = coeff(rf,var(1)^3)*g1;
1240
1241  /*Read off the desired information */
1242  if(coeff(rf,var(1)^3)<0)
1243  {
1244    phi = rf, -var(1),var(2);
1245    rf = phi(rf);
1246  }
1247  int s =  sign(coeff(rf,var(1)^8));
1248  if(coeff(rf,var(1)^8)<0)
1249  {
1250    typeofsing = "E[14]-";
1251  }
1252  else
1253  {
1254    typeofsing = "E[14]+";
1255  }
1256  MP = var(1)^(12)-coeff(rf,var(1)^3)^(-4)*absValue(coeff(rf,var(2)^8))^(-9)*coeff(rf,var(1)*var(2)^6)^(12);
1257  if(coeff(rf,var(1)*var(2)^6)>0)
1258  {
1259    number n = maxabs(MP);
1260    MP = RealRootFactor(MP,0,n+1);
1261    n = maxabs(MP);
1262    number n1 = 0;
1263    number n2 = n+1;
1264  }
1265  else
1266  {
1267    number n = maxabs(MP);
1268    MP = RealRootFactor(MP,-(n+1),0);
1269    n = maxabs(MP);
1270    number n1 = -(n+1);
1271    number n2 = 0;
1272  }
1273  if(deg(MP)==1)
1274  {
1275    number a = -number(subst(MP,var(1),0));
1276    setring br1;
1277    number a = fetch(br,a);
1278  }
1279  else
1280  {
1281    string mord = ringlist(br1)[3][1][1];
1282    ring R1 = (0,a),(var(1),var(2)),ds;
1283    list L1 = ringlist(R1);
1284    L1[3][1][1]=mord;
1285    def R2 = ring(L1);
1286    setring R2;
1287    kill R1;
1288    poly ff = fetch(br,MP);
1289    minpoly = number(subst(ff,var(1),a));
1290  }
1291  number n1 = fetch(br,n1);
1292  number n2 = fetch(br,n2);
1293  Poly nfe =var(1)^3+number(s)*var(2)^8+a*var(1)*var(2)^6;
1294  intervalP = list(n1,n2);
1295  NormalFormEquation F;
1296  F.normalFormEquation = nfe;
1297  F.parameters = list(list(a*var(1)*var(2)^6,intervalP));
1298  F.singularityType = typeofsing;
1299  setring br;
1300  return(list(F));
1301}
1302/////////////////////////////////////////////////////////////////////////////////////////
1303static proc caseJ10(poly rf, br1)
1304{
1305  /*preliminaries*/
1306  map phi;
1307  poly MP;
1308  poly MP1;
1309  poly MP2;
1310  poly MP3;
1311  def br = basering;
1312  string typeofsing;
1313  string typeofsing1;
1314  string typeofsing2;
1315  string typeofsing3;
1316  string intervalP;
1317  list interval1;
1318  list interval2;
1319  list interval3;
1320  intvec v = 2,1;
1321  rf = jet(rf,6);
1322  rf = lowestJet(rf);
1323  number n = coeff(rf,var(1)^3);
1324  rf = (1/n)*rf;
1325  n = (coeff(rf,var(1)^2*var(2)^2)*(1/3));
1326  phi = br,var(1)-n*var(2)^2,var(2);
1327  rf = phi(rf);
1328  number e = coeff(rf,var(2)^6);
1329  number d = coeff(rf,var(1)*var(2)^4);
1330  poly pplus = (4*d^3+27*e^2)*var(1)^6+(-36*d^3-243*e^2)*var(1)^4+(81*d^3+729*e^2)*var(1)^2-729*e^2;
1331  poly pminus = -(4*d^3+27*e^2)*var(1)^6+(-36*d^3-243*e^2)*var(1)^4-(81*d^3+729*e^2)*var(1)^2-729*e^2;
1332  poly k = var(1)^3+d*var(1)+e;
1333  ring R = 0,x,dp;
1334  poly k = fetch(br,k);
1335  int nn = nrroots(k);
1336  kill R;
1337  setring br;
1338  if(nn==1)
1339  {
1340    typeofsing1 = "J[10]+";
1341    if(e<0)
1342    {
1343      MP = pplus;
1344      number n00 = maxabs(MP);
1345      MP = RealRootFactor(MP,0,n00+1);
1346      number n0 = maxabs(MP);
1347      number n1 = 0;
1348      number n2 = n0+1;
1349    }
1350    if(e>0)
1351    {
1352      MP = pplus;
1353      number n00 = maxabs(MP);
1354      MP = RealRootFactor(MP,-(n00+1),0);
1355      number n0 = maxabs(MP);
1356      number n1 = -(n0+1);
1357      number n2 = 0;
1358    }
1359    if(e==0)
1360    {
1361      MP=var(1);
1362      number n1 = 0;
1363      number n2 = 0;
1364    }
1365    if(deg(MP)==1)
1366    {
1367      number a = -number(subst(MP,var(1),0));
1368      setring br1;
1369      number a = fetch(br,a);
1370    }
1371    else
1372    {
1373      string mord = ringlist(br1)[3][1][1];
1374      ring R1 = (0,a),(var(1),var(2)),ds;
1375      list L1 = ringlist(R1);
1376      L1[3][1][1]=mord;
1377      def R2 = ring(L1);
1378      setring R2;
1379      kill R1;
1380      poly ff = fetch(br,MP);
1381      minpoly = number(subst(ff,var(1),a));
1382    }
1383    Poly nfe1 = var(1)^3+a*var(1)^2*var(2)^2+var(1)*var(2)^4;
1384    number n1 = fetch(br,n1);
1385    number n2 = fetch(br,n2);
1386    interval1 = list(n1,n2);
1387    NormalFormEquation F1;
1388    F1.normalFormEquation = nfe1;
1389    F1.parameters = list(list(a*var(1)^2*var(2)^2,interval1));
1390    F1.singularityType = typeofsing1;
1391    setring br;
1392    return(F1);
1393  }
1394  else
1395  {
1396    phi = br,0,0;
1397    k = phi(k);
1398    if(k==0)
1399    {
1400      typeofsing1 = "J[10]+";
1401      typeofsing2 = "J[10]-";
1402      typeofsing3 = "J[10]+";
1403      string mord = ringlist(br1)[3][1][1];
1404
1405      ring R1 = 0,(var(1),var(2)),ds;
1406      list L1 = ringlist(R1);
1407      L1[3][1][1]=mord;
1408      def R2 = ring(L1);
1409      setring R2;
1410      kill R1;
1411      Poly nfe2 =var(1)^3-var(1)*var(2)^4;
1412      interval2 = list(0,0);
1413      NormalFormEquation F2;
1414      F2.normalFormEquation = nfe2;
1415      F2.parameters = list(list(0*var(1)^2*var(2)^2,interval2));
1416      F2.singularityType = typeofsing2;
1417
1418      ring R1 = (0,a),(var(1),var(2)),ds;
1419      list L1 = ringlist(R1);
1420      L1[3][1][1]=mord;
1421      def R2 = ring(L1);
1422      setring R2;
1423      kill R1;
1424      minpoly = a^2-(9/2);
1425      Poly nfe1 =var(1)^3+a*var(1)^2*var(2)^2+var(1)*var(2)^4;
1426      Poly nfe3 =var(1)^3+a*var(1)^2*var(2)^2+var(1)*var(2)^4;
1427      interval1 = list(0,3/2);
1428      interval3 = list(-3/2,0);
1429      NormalFormEquation F1;
1430      F1.normalFormEquation = nfe1;
1431      F1.parameters = list(list(a*var(1)^2*var(2)^2,interval1));
1432      F1.singularityType = typeofsing1;
1433      NormalFormEquation F3;
1434      F3.normalFormEquation = nfe3;
1435      F3.parameters = list(list(a*var(1)^2*var(2)^2,interval3));
1436      F3.singularityType = typeofsing3;
1437      setring br;
1438    }
1439    else
1440    {
1441      int epsilon = 1;
1442      ring R1 = 0,x,dp;
1443      poly pplus=fetch(br,pplus);
1444      ring R = (real,30,i),x,dp;
1445      poly pplus = fetch(br,pplus);
1446      list r;
1447      while(size(r)<6)
1448      {
1449         epsilon = epsilon+1;
1450         r = laguerre_solve(pplus,epsilon);
1451      }
1452      list r1;
1453      for(int j = 1;j<=size(r);j++)
1454      {
1455        if(impart(r[j])==0)
1456        {
1457          bigint s = int(r[j]*10^(epsilon));
1458          r1 = insert(r1,s);
1459          kill s;
1460        }
1461      }
1462      kill r;
1463      kill j;
1464      int kk = size(r1);
1465      setring R1;
1466      list r2;
1467      for(int ii=1; ii<=kk;ii++)
1468      {
1469         number m;
1470         m=number(r1[ii])/(10^(epsilon));
1471         r2 = insert(r2,m);
1472         kill m;
1473      }
1474      kill ii;
1475      list L1 = ringlist(br1);
1476      number t = (number(epsilon))^(-1);
1477      number n1 = sturm(pplus,r2[1]-t,r2[1]+t);
1478      number n2 = sturm(pplus,r2[2]-t,r2[2]+t);
1479      number n3 = sturm(pplus,r2[3]-t,r2[3]+t);
1480      number n4 = sturm(pplus,r2[4]-t,r2[4]+t);
1481      while(n1>1 || n2>1 || n3>1 || n4>1)
1482      {
1483        kill n1;
1484        kill n2;
1485        kill n3;
1486        kill n4;
1487        kill t;
1488        kill r1;
1489        kill r2;
1490        epsilon = epsilon+1;
1491        setring R;
1492        pplus = fetch(br,pplus);
1493        list r = laguerre_solve(pplus,epsilon);
1494        list r1;
1495        for(int j = 1;j<=size(r);j++)
1496        {
1497          if(impart(r[j])==0)
1498          {
1499            bigint s = int(r[j]*10^(epsilon));
1500            r1 = insert(r1,s);
1501            kill s;
1502          }
1503        }
1504        kill r;
1505        kill j;
1506        setring R1;
1507        kk = size(r1);
1508        list r2;
1509        for(int ii=1; ii<=kk;ii++)
1510        {
1511           number m;
1512           m=number(r1[ii])/(10^(epsilon));
1513           r2 = insert(r2,m);
1514           kill m;
1515        }
1516        kill ii;
1517        number t = (number(epsilon))^(-1);
1518        number n1 = sturm(pplus,r2[1]-t,r2[1]+t);
1519        number n2 = sturm(pplus,r2[2]-t,r2[2]+t);
1520        number n3 = sturm(pplus,r2[3]-t,r2[3]+t);
1521        number n4 = sturm(pplus,r2[4]-t,r2[4]+t);
1522      }
1523      setring br;
1524      kill R1;
1525      list r2;
1526      kk = size(r1);
1527      for(int ii=1; ii<=kk;ii++)
1528      {
1529         number m;
1530         m=number(r1[ii])/(10^(epsilon));
1531         r2 = insert(r2,m);
1532         kill m;
1533      }
1534      kill ii;
1535      number t = (number(epsilon))^(-1);
1536      typeofsing1 = "J[10]-";
1537      typeofsing2 = "J[10]+";
1538      typeofsing3 = "J[10]+";
1539
1540      if(e>0)
1541      {
1542        MP3 = RealRootFactor(pplus, r2[4]-t, r2[4]+t);
1543        number n1 = r2[4]-t;
1544        number n2 = r2[4]+t;
1545      }
1546      else
1547      {
1548        MP3 = RealRootFactor(pplus, r2[3]-t, r2[3]+t);
1549        number n1 = r2[3]-t;
1550        number n2 = r2[3]+t;
1551      }
1552      if(deg(MP3)==1)
1553      {
1554        number a = -number(subst(MP3,var(1),0));
1555        def R2 =  br1;
1556        setring R2;
1557        number a = fetch(br,a);
1558      }
1559      else
1560      {
1561        string mord = ringlist(br1)[3][1][1];
1562        ring R1 = (0,a),(var(1),var(2)),ds;
1563        list L1 = ringlist(R1);
1564        L1[3][1][1]=mord;
1565        def R2 = ring(L1);
1566        setring R2;
1567        kill mord;
1568        kill R1;
1569        poly ff = fetch(br,MP3);
1570        minpoly = number(subst(ff,var(1),a));
1571      }
1572      number n1 = fetch(br,n1);
1573      number n2 = fetch(br,n2);
1574      Poly nfe3 =var(1)^3+a*var(1)^2*var(2)^2+var(1)*var(2)^4;
1575      interval3 = list(n1,n2);
1576      NormalFormEquation F3;
1577      F3.normalFormEquation = nfe3;
1578      F3.parameters = list(list(a*var(1)^2*var(2)^2,interval3));
1579      F3.singularityType = typeofsing3;
1580      setring br;
1581      kill R2;
1582
1583      if(e>0)
1584      {
1585        MP2 = RealRootFactor(pplus, r2[2]-t,r2[2]+t);
1586        n1 = r2[2]-t;
1587        n2 = r2[2]+t;
1588      }
1589      else
1590      {
1591        MP2 = RealRootFactor(pplus, r2[1]-t,r2[1]+t);
1592        n1 = r2[1]-t;
1593        n2 = r2[1]+t;
1594      }
1595      if(deg(MP2)==1)
1596      {
1597        number a = -number(subst(MP2,var(1),0));
1598        def R2 =  br1;
1599        setring R2;
1600        number a = fetch(br,a);
1601      }
1602      else
1603      {
1604        string mord = ringlist(br1)[3][1][1];
1605        ring R1 = (0,a),(var(1),var(2)),ds;
1606        list L1 = ringlist(R1);
1607        L1[3][1][1]=mord;
1608        def R2 = ring(L1);
1609        setring R2;
1610        kill mord;
1611        kill R1;
1612        poly ff = fetch(br,MP2);
1613        minpoly = number(subst(ff,var(1),a));
1614      }
1615      number n1 = fetch(br,n1);
1616      number n2 = fetch(br,n2);
1617      Poly nfe2 =var(1)^3+a*var(1)^2*var(2)^2+var(1)*var(2)^4;
1618      interval2 = list(n1,n2);
1619      NormalFormEquation F2;
1620      F2.normalFormEquation = nfe2;
1621      F2.parameters = list(list(a*var(1)^2*var(2)^2,interval2));
1622      F2.singularityType = typeofsing2;
1623      setring br;
1624      kill R2;
1625      number n0 = maxabs(pminus);
1626      if(e>0)
1627      {
1628        MP1 = RealRootFactor(pminus, 0,n0+1);
1629        n0  = maxabs(pminus);
1630        n1 = 0;
1631        n2 = n0+1;
1632      }
1633      else
1634      {
1635        MP1 = RealRootFactor(pminus, -(n0+1),0);
1636        n0 = maxabs(pminus);
1637        n1 = -n0-1;
1638        n2 = 0;
1639      }
1640      if(deg(MP1)==1)
1641      {
1642        number a = -number(subst(MP1,var(1),0));
1643        setring br1;
1644        number a = fetch(br,a);
1645      }
1646      else
1647      {
1648        string mord = ringlist(br1)[3][1][1];
1649        ring R1 = (0,a),(var(1),var(2)),ds;
1650        list L1 = ringlist(R1);
1651        L1[3][1][1]=mord;
1652        def R2 = ring(L1);
1653        setring R2;
1654        kill mord;
1655        kill R1;
1656        poly ff = fetch(br,MP1);
1657        minpoly = number(subst(ff,var(1),a));
1658      }
1659      number n1 = fetch(br,n1);
1660      number n2 = fetch(br,n2);
1661      Poly nfe1 =var(1)^3+a*var(1)^2*var(2)^2-var(1)*var(2)^4;
1662      interval1 = list(n1,n2);
1663      NormalFormEquation F1;
1664      F1.normalFormEquation = nfe1;
1665      F1.parameters = list(list(a*var(1)^2*var(2)^2,interval1));
1666      F1.singularityType = typeofsing1;
1667     }
1668     return(list(F1,F2,F3));
1669  }
1670}
1671
1672///////////////////////////////////////////////////////////////////////////////////////////////
1673static proc caseYrs(poly rf, int mu ,br1)
1674{
1675  //preliminaries
1676  def R = basering;
1677  def br = basering;
1678  map phi;
1679  poly MP1;
1680  poly MP2;
1681  poly MP3;
1682  poly MP4;
1683  poly MP5;
1684  poly MP6;
1685  poly MP7;
1686  poly MP8;
1687  poly nf1;
1688  poly nf2;
1689  poly nf3;
1690  poly nf4;
1691  poly nf5;
1692  poly nf6;
1693  poly nf7;
1694  poly nf8;
1695  list interval1;
1696  list interval2;
1697  list interval3;
1698  list interval4;
1699  list interval5;
1700  list interval6;
1701  list interval7;
1702  list interval8;
1703  string typeofsing1;
1704  string typeofsing2;
1705  string typeofsing3;
1706  string typeofsing4;
1707  string typeofsing5;
1708  string typeofsing6;
1709  string typeofsing7;
1710  string typeofsing8;
1711  int d = determinacy(rf,mu);
1712  rf = jet(rf,d);
1713  poly f0 = rf;
1714  int n1 = (mu-1) div 2;
1715  poly rts = subst(jet(rf,4),var(2),1);
1716  ring S = 0,(x),dp;
1717  poly rts = fetch(br,rts);
1718  int rt = nrroots(rts);
1719  setring br;
1720  kill S;
1721  if((n1 mod 2 == 0) && (rt==0))
1722  {
1723     list factorization = factorize(jet(rf,4));
1724     poly rf1 = factorization[1][2];
1725     map phi = R, var(1)-coeff(rf1,var(1)*var(2))/(2*coeff(rf1,var(1)^2))*var(2),var(2);
1726     rf = phi(rf);
1727     poly q = phi(rf1);
1728     number b = (coeff(rf,var(1)^4));
1729     int sigmab = sign(b);
1730     rf = eliminateTildeY(rf,q,n1);
1731     b = b^(-1);
1732     poly P = var(1)^4-(coeff(rf,var(1)^(n1)))^4*b^(n1);
1733     number c = coeff(rf,var(1)^(n1));
1734     int sigmac = sign(c);
1735     if(sign(b)>0)
1736     {
1737       typeofsing1 = "tilde(Y)["+string(n1)+"]+";
1738     }
1739     else
1740     {
1741        typeofsing1 = "tilde(Y)["+string(n1)+"]-";
1742     }
1743     number m = maxabs(P);
1744     if(sign(c)>0)
1745     {
1746        MP1 = RealRootFactor(P,0,m+1);
1747     }
1748     else
1749     {
1750       MP1 = RealRootFactor(P,-m-1,0);
1751     }
1752     number m1 = maxabs(MP1);
1753     if(deg(MP1)==1)
1754     {
1755       number a = -number(subst(MP1,var(1),0));
1756       def R2 = br1;
1757       setring R2;
1758       number a = fetch(R,a);
1759     }
1760     else
1761     {
1762       string mord = ringlist(br1)[3][1][1];
1763       ring R1 = (0,a),(var(1),var(2)),ds;
1764       list L1 = ringlist(R1);
1765       L1[3][1][1]=mord;
1766       def R2 = ring(L1);
1767       setring R2;
1768       kill mord;
1769       kill R1;
1770       poly ff = fetch(R,MP1);
1771       minpoly = number(subst(ff,var(1),a));
1772     }
1773     Poly nfe1 =sigmab*(var(1)^2+var(2)^2)^2+a*var(1)^(n1);
1774     number m1 = fetch(R,m1);
1775     if(sigmac<0)
1776     {
1777       interval1 = list(-(m1+1),0);
1778     }
1779     else
1780     {
1781       interval1 = list(0,(m1+1));
1782     }
1783     NormalFormEquation F1;
1784     F1.normalFormEquation = nfe1;
1785     F1.parameters = list(list(a*var(1)^(n1),interval1));
1786     F1.singularityType = typeofsing1;
1787
1788     setring br;
1789     return(F1);
1790  }
1791  else
1792  {
1793  if(coeff(rf,var(1)^4)>0)
1794  {
1795    int sigma =1;
1796  }
1797  else
1798  {
1799    int sigma =-1;
1800  }
1801
1802  // Eliminate the monomials in supp(rf) underneath Gamma(T)
1803
1804  poly h = jet(rf,4);
1805  //Implement line 3 in Algorithm 5
1806  int sum;
1807  for(int i=2; i<=size(factorize(h)[2]);i++)
1808  {
1809    sum = sum+factorize(h)[2][i];
1810  }
1811  if(sum==4)
1812  {
1813    list factorization = factorize(h);
1814     list factors = factorization[1][2];
1815         factors = insert(factors,factorization[1][3]);
1816    matrix T[2][2] = factors[1]/var(1), factors[1]/var(2),
1817                     factors[2]/var(1), factors[2]/var(2);
1818    phi = br, luinverse(T)[2]*matrix(ideal(var(1), var(2)), 2, 1);
1819    rf = phi(rf);
1820    rf = jet(rf, d);
1821  }
1822  else
1823  {
1824    def S = absFactorize(h);
1825    setring S;
1826    poly h = fetch(br,h);
1827    poly rf = fetch(br,rf);
1828    poly f =absolute_factors[3][2];
1829    poly g = absolute_factors[1][2];
1830    kill br;
1831    ring br = (0,a),(var(1),var(2)),ds;
1832    poly f = fetch(S,f);
1833    minpoly=number(f);
1834    poly h = fetch(S,h);
1835    //poly rf = fetch(br,rf);
1836    list factors = factorize(h);
1837    matrix T[2][2] = factors[1][2]/var(1), factors[1][2]/var(2),
1838                     factors[1][3]/var(1),factors[1][3]/var(2);
1839    map phi = br, luinverse(T)[2]*matrix(ideal(var(1), var(2)), 2, 1);
1840    poly rf = fetch(S,rf);
1841    rf = phi(rf);
1842    kill S;
1843  }
1844  int n1 = 3;
1845  int n2 = 3;
1846  number c1 = coeff(rf,var(1)^(n1));
1847  number c2 = coeff(rf,var(2)^(n2));
1848  number s1 = 0;
1849  number s2 = 0;
1850  number b;
1851  while(c1==0 || c2==0)
1852  {
1853    if(coeff(rf,var(1)^(n1))==0 && s1==0)
1854    {
1855      n1 = n1+1;
1856      b = coeff(rf,var(1)^(2)*var(2)^(2));
1857      phi = br, var(1), var(2)-(coeff(rf,var(1)^(n1)*var(2))/(2*b))*var(1)^(n1-2);
1858      rf = phi(rf);
1859      c1 = coeff(rf,var(1)^(n1));
1860    }
1861    else
1862    {
1863      s1 = 1;
1864    }
1865    if(coeff(rf,var(2)^(n2))==0 && s2==0)
1866    {
1867      n2 = n2+1;
1868      phi = br, var(1)-(coeff(rf,var(2)^(n2)*var(1))/(2*b))*var(2)^(n2-2), var(2);
1869      rf = phi(rf);
1870      c2 = coeff(rf,var(2)^(n2));
1871    }
1872    else
1873    {
1874      s2 = 1;
1875    }
1876    rf = jet(rf,d);
1877  }
1878
1879  kill c1;
1880  kill c2;
1881  kill s1;
1882  kill s2;
1883  kill b;
1884  //Eliminate the monomials above Gamma(T)
1885  //Read off the desired information
1886  ring S = 0,(x,y),dp;
1887  poly f0 = fetch(R,f0);
1888  poly f0x = subst(jet(f0,4),y,1);
1889  int n = nrroots(f0x);
1890  setring br;
1891  kill S;
1892  if(n>0)
1893  {
1894    if(minpoly == 0) //h has a linear factor over Q
1895    {
1896      if((n1 mod 2 <> 0) && (n2 mod 2 ==0))  //r is odd and s is even
1897      {
1898        phi = br,var(2),var(1);
1899        rf = phi(rf);
1900        int s = n1;
1901        n1 = n2;
1902        n2 = s;
1903      }
1904      poly P = var(1)^(2*n1)-(absValue(coeff(rf,var(1)^2*var(2)^2)))^(-n1*n2)*coeff(rf,var(1)^n1)^(2*n2)*coeff(rf,var(2)^n2)^(2*n1);
1905      poly P1 =  var(1)^(2*n2)-(absValue(coeff(rf,var(1)^2*var(2)^2)))^(-n1*n2)*coeff(rf,var(1)^n1)^(2*n2)*coeff(rf,var(2)^n2)^(2*n1);
1906
1907     if((n1 mod 2 == 0) && (n2 mod 2 ==0))  // r is even and s is even
1908      {
1909        if(coeff(rf,var(1)^n1)>0) //d>0
1910        {
1911          number m = maxabs(P1);
1912          MP2 = RealRootFactor(P1,0,m+1);
1913          number m2 = maxabs(MP2);
1914          if(coeff(rf,var(2)^n2)>0) //e>0
1915          {
1916            m = maxabs(P);
1917            MP1 = RealRootFactor(P,0,m+1);
1918            number m1 = maxabs(MP1);
1919            if(coeff(rf,var(1)^2*var(2)^2)>0)   //b>0
1920            {
1921              typeofsing1  = "Y["+string(n1)+","+string(n2)+"]++";
1922              typeofsing2  = "Y["+string(n2)+","+string(n1)+"]++";
1923              if(deg(MP1)==1)
1924              {
1925                 number a = -number(subst(MP1,var(1),0));
1926                 def R2 =  br1;
1927                 setring R2;
1928                 number a = fetch(br,a);
1929              }
1930              else
1931              {
1932                string mord = ringlist(br1)[3][1][1];
1933                ring R1 = (0,a),(var(1),var(2)),ds;
1934                list L1 = ringlist(R1);
1935                L1[3][1][1]=mord;
1936                def R2 = ring(L1);
1937                kill mord;
1938                setring R2;
1939                kill R1;
1940                poly ff = fetch(br,MP1);
1941                minpoly = number(subst(ff,var(1),a));
1942              }
1943              Poly nfe1 =var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n2);
1944              number m1 = fetch(br,m1);
1945              interval1 = list(0,m1+1);
1946              NormalFormEquation F1;
1947              F1.normalFormEquation = nfe1;
1948              F1.parameters = list(list(a*var(2)^(n2),interval1));
1949              F1.singularityType = typeofsing1;
1950              setring br;
1951              kill R2;
1952
1953              if(deg(MP2)==1)
1954              {
1955                 number a = -number(subst(MP2,var(1),0));
1956                 def R2 = br1;
1957                 setring R2;
1958                 number a = fetch(br,a);
1959              }
1960              else
1961              {
1962                string mord = ringlist(br1)[3][1][1];
1963                ring R1 = (0,a1),(var(1),var(2)),ds;
1964                list L1 = ringlist(R1);
1965                L1[3][1][1]=mord;
1966                def R2 = ring(L1);
1967                setring R2;
1968                kill mord;
1969                kill R1;
1970                poly ff = fetch(br,MP2);
1971                minpoly = number(subst(ff,var(1),a));
1972              }
1973              Poly nfe2 =var(1)^2*var(2)^2+var(1)^(n2)+a*var(2)^(n1);
1974              number m2 = fetch(br,m2);
1975              interval2 = list(0,m2+1);
1976              NormalFormEquation F2;
1977              F2.normalFormEquation = nfe2;
1978              F2.parameters = list(list(a*var(2)^(n1),interval2));
1979              F2.singularityType = typeofsing2;
1980              return(F1,F2);
1981            }
1982            else //b<0
1983            {
1984              typeofsing1  = "Y["+string(n1)+","+string(n2)+"]-+";
1985              typeofsing2  = "Y["+string(n2)+","+string(n1)+"]-+";
1986              if(deg(MP1)==1)
1987              {
1988                 number a = -number(subst(MP1,var(1),0));
1989                 def R2 = br1;
1990                 setring R2;
1991                 number a = fetch(br,a);
1992              }
1993              else
1994              {
1995                string mord = ringlist(br1)[3][1][1];
1996                ring R1 = (0,a),(var(1),var(2)),ds;
1997                list L1 = ringlist(R1);
1998                L1[3][1][1]=mord;
1999                def R2 = ring(L1);
2000                setring R2;
2001                kill R1;
2002                kill mord;
2003                poly ff = fetch(br,MP1);
2004                minpoly = number(subst(ff,var(1),a));
2005              }
2006              Poly nfe1 =-var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n2);
2007              number m1 = fetch(br,m1);
2008              interval1 = list(0,m1+1);
2009              NormalFormEquation F1;
2010              F1.normalFormEquation = nfe1;
2011              F1.parameters = list(list(a*var(2)^(n2),interval1));
2012              F1.singularityType = typeofsing1;
2013              setring br;
2014              kill R2;
2015              if(deg(MP2)==1)
2016              {
2017                 number a = -number(subst(MP2,var(1),0));
2018                 def R2 = br1;
2019                 setring R2;
2020                 number a = fetch(br,a);
2021              }
2022              else
2023              {
2024                string mord = ringlist(br1)[3][1][1];
2025                ring R1 = (0,a),(var(1),var(2)),ds;
2026                list L1 = ringlist(R1);
2027                L1[3][1][1]=mord;
2028                def R2 = ring(L1);
2029                kill mord;
2030                setring R2;
2031                kill R1;
2032                poly ff = fetch(br,MP2);
2033                minpoly = number(subst(ff,var(1),a));
2034              }
2035              Poly nfe2 =-var(1)^2*var(2)^2+var(1)^(n2)+a*var(2)^(n1);
2036              number m2 = fetch(br,m2);
2037              interval2 = list(0,m2+1);
2038              NormalFormEquation F2;
2039              F2.normalFormEquation = nfe2;
2040              F2.parameters = list(list(a*var(2)^(n1),interval2));
2041              F2.singularityType = typeofsing2;
2042              setring br;
2043            }
2044          }
2045          else //e<0
2046          {
2047           number m = maxabs(P);
2048           MP1 = RealRootFactor(P,-m-1,0);
2049           number m1 = maxabs(MP1);
2050           if(coeff(rf,var(1)^2*var(2)^2)>0)  //b>0
2051           {
2052              typeofsing1  = "Y["+string(n1)+","+string(n2)+"]++";
2053              typeofsing2  = "Y["+string(n2)+","+string(n1)+"]+-";
2054              if(deg(MP1)==1)
2055              {
2056                 number a = -number(subst(MP1,var(1),0));
2057                 def R2 = br1;
2058                 setring R2;
2059                 number a = fetch(br,a);
2060              }
2061              else
2062              {
2063                string mord = ringlist(br1)[3][1][1];
2064                ring R1 = (0,a),(var(1),var(2)),ds;
2065                list L1 = ringlist(R1);
2066                L1[3][1][1]=mord;
2067                def R2 = ring(L1);
2068                kill mord;
2069                setring R2;
2070                kill R1;
2071                poly ff = fetch(br,MP1);
2072                minpoly = number(subst(ff,var(1),a));
2073              }
2074              Poly nfe1 =var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n2);
2075              number m1 = fetch(br,m1);
2076              interval1 = list(-(m1+1),0);
2077              NormalFormEquation F1;
2078              F1.normalFormEquation = nfe1;
2079              F1.parameters = list(list(a*var(2)^(n2),interval1));
2080              F1.singularityType = typeofsing1;
2081              setring br;
2082              kill R2;
2083
2084              if(deg(MP2)==1)
2085              {
2086                 number a = -number(subst(MP2,var(1),0));
2087                 def R2 = br1;
2088                 setring R2;
2089                 number a = fetch(br,a);
2090              }
2091              else
2092              {
2093                string mord = ringlist(br1)[3][1][1];
2094                ring R1 = (0,a),(var(1),var(2)),ds;
2095                list L1 = ringlist(R1);
2096                L1[3][1][1]=mord;
2097                def R2 = ring(L1);
2098                kill mord;
2099                setring R2;
2100                kill R1;
2101                poly ff = fetch(br,MP2);
2102                minpoly = number(subst(ff,var(1),a));
2103              }
2104              Poly nfe2 = var(1)^2*var(2)^2-var(1)^(n2)+a*var(2)^(n1);
2105              number m2 = fetch(br,m2);
2106              interval2 = list(0,m2+1);
2107              NormalFormEquation F2;
2108              F2.normalFormEquation = nfe2;
2109              F2.parameters = list(list(a*var(2)^(n1),interval2));
2110              F2.singularityType = typeofsing2;
2111              setring br;
2112              kill R2;
2113           }
2114           else //b<0
2115           {
2116             typeofsing1  = "Y["+string(n1)+","+string(n2)+"]-+";
2117             typeofsing2  = "Y["+string(n2)+","+string(n1)+"]--";
2118             if(deg(MP1)==1)
2119             {
2120                number a = -number(subst(MP1,var(1),0));
2121                def R2 = br1;
2122                setring R2;
2123                number a = fetch(br,a);
2124             }
2125             else
2126             {
2127               string mord = ringlist(br1)[3][1][1];
2128               ring R1 = (0,a),(var(1),var(2)),ds;
2129               list L1 = ringlist(R1);
2130               L1[3][1][1]=mord;
2131               def R2 = ring(L1);
2132               kill mord;
2133               setring R2;
2134               kill R1;
2135               poly ff = fetch(br,MP1);
2136               minpoly = number(subst(ff,var(1),a));
2137             }
2138             Poly nfe1 =-var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n2);
2139             number m1 = fetch(br,m1);
2140             interval1 = list(-(m1+1),0);
2141             NormalFormEquation F1;
2142             F1.normalFormEquation = nfe1;
2143             F1.parameters = list(list(a*var(2)^(n2),interval1));
2144             F1.singularityType = typeofsing1;
2145             setring br;
2146             kill R2;
2147
2148             if(deg(MP2)==1)
2149             {
2150                number a = -number(subst(MP2,var(1),0));
2151                def R2 = br1;
2152                setring R2;
2153                number a = fetch(br,a);
2154             }
2155             else
2156             {
2157               string mord = ringlist(br1)[3][1][1];
2158               ring R1 = (0,a),(var(1),var(2)),ds;
2159               list L1 = ringlist(R1);
2160               L1[3][1][1]=mord;
2161               def R2 = ring(L1);
2162               setring R2;
2163               kill mord;
2164               kill R1;
2165               poly ff = fetch(br,MP2);
2166               minpoly = number(subst(ff,var(1),a));
2167             }
2168             Poly nfe2 = -var(1)^2*var(2)^2-var(1)^(n2)+a*var(2)^(n1);
2169             number m2 = fetch(br,m2);
2170             interval2 = list(0,m2+1);
2171             NormalFormEquation F2;
2172             F2.normalFormEquation = nfe2;
2173             F2.parameters = list(list(a*var(2)^(n1),interval2));
2174             F2.singularityType = typeofsing2;
2175             setring br;
2176           }
2177         }
2178       }
2179       else //d<0
2180       {
2181         number m = maxabs(P1);
2182         MP2 = RealRootFactor(P1,-m-1,0);
2183         number m2 = maxabs(MP2);
2184         if(coeff(rf,var(2)^n2)>0) //e>0
2185         {
2186           number m = maxabs(P);
2187           MP1 = RealRootFactor(P,0,m+1);
2188           number m1 = maxabs(MP1);
2189           if(coeff(rf,var(1)^2*var(2)^2)>0)   //b>0
2190           {
2191             typeofsing1  = "Y["+string(n1)+","+string(n2)+"]+-";
2192             typeofsing2  = "Y["+string(n2)+","+string(n1)+"]++";
2193
2194             if(deg(MP1)==1)
2195             {
2196                number a = -number(subst(MP1,var(1),0));
2197                def R2 = br1;
2198                setring R2;
2199                number a = fetch(br,a);
2200             }
2201             else
2202             {
2203               string mord = ringlist(br1)[3][1][1];
2204               ring R1 = (0,a),(var(1),var(2)),ds;
2205               list L1 = ringlist(R1);
2206               L1[3][1][1]=mord;
2207               def R2 = ring(L1);
2208               setring R2;
2209               kill mord;
2210               kill R1;
2211               poly ff = fetch(br,MP1);
2212               minpoly = number(subst(ff,var(1),a));
2213             }
2214             Poly nfe1 =var(1)^2*var(2)^2-var(1)^(n1)+a*var(2)^(n2);
2215             number m1 = fetch(br,m1);
2216             interval1 = list(0,m1+1);
2217             NormalFormEquation F1;
2218             F1.normalFormEquation = nfe1;
2219             F1.parameters = list(list(a*var(2)^(n2),interval1));
2220             F1.singularityType = typeofsing1;
2221
2222             setring br;
2223             kill R2;
2224
2225             if(deg(MP2)==1)
2226             {
2227                 number a = -number(subst(MP2,var(1),0));
2228                 def R2 = br1;
2229                 setring R2;
2230                 number a = fetch(br,a);
2231
2232             }
2233             else
2234             {
2235               string mord = ringlist(br1)[3][1][1];
2236               ring R1 = (0,a),(var(1),var(2)),ds;
2237               list L1 = ringlist(R1);
2238               L1[3][1][1]=mord;
2239               def R2 = ring(L1);
2240               setring R2;
2241               kill R1;
2242               kill mord;
2243               poly ff = fetch(br,MP2);
2244               minpoly = number(subst(ff,var(1),a));
2245             }
2246             Poly nfe2 = var(1)^2*var(2)^2+var(1)^(n2)+a*var(2)^(n1);
2247             number m2 = fetch(br,m2);
2248             interval2 = list(-(m2+1),0);
2249             NormalFormEquation F2;
2250             F2.normalFormEquation = nfe2;
2251             F2.parameters = list(list(a*var(2)^(n1),interval2));
2252             F2.singularityType = typeofsing2;
2253             setring br;
2254           }
2255           else //b<0
2256           {
2257             typeofsing1  = "Y["+string(n1)+","+string(n2)+"]--";
2258             typeofsing2  = "Y["+string(n2)+","+string(n1)+"]-+";
2259
2260
2261             if(deg(MP1)==1)
2262             {
2263                number a = -number(subst(MP1,var(1),0));
2264                def R2 = br1;
2265                setring R2;
2266                number a = fetch(br,a);
2267             }
2268             else
2269             {
2270               string mord = ringlist(br1)[3][1][1];
2271               ring R1 = (0,a),(var(1),var(2)),ds;
2272               list L1 = ringlist(R1);
2273               L1[3][1][1]=mord;
2274               def R2 = ring(L1);
2275               setring R2;
2276               kill R1;
2277               kill mord;
2278               poly ff = fetch(br,MP1);
2279               minpoly = number(subst(ff,var(1),a));
2280             }
2281             Poly nfe1 =-var(1)^2*var(2)^2-var(1)^(n1)+a*var(2)^(n2);
2282             number m1 = fetch(br,m1);
2283             interval1 = list(0,m1+1);
2284             NormalFormEquation F1;
2285             F1.normalFormEquation = nfe1;
2286             F1.parameters = list(list(a*var(2)^(n2),interval1));
2287             F1.singularityType = typeofsing1;
2288
2289             setring br;
2290             kill R2;
2291
2292             if(deg(MP2)==1)
2293             {
2294                 number a = -number(subst(MP2,var(1),0));
2295                 def R2 = br1;
2296                 setring R2;
2297                 number a = fetch(br,a);
2298             }
2299             else
2300             {
2301               string mord = ringlist(br1)[3][1][1];
2302               ring R1 = (0,a),(var(1),var(2)),ds;
2303               list L1 = ringlist(R1);
2304               L1[3][1][1]=mord;
2305               def R2 = ring(L1);
2306               setring R2;
2307               kill R1;
2308               kill mord;
2309               poly ff = fetch(br,MP2);
2310               minpoly = number(subst(ff,var(1),a));
2311             }
2312             Poly nfe2 = -var(1)^2*var(2)^2+var(1)^(n2)+a*var(2)^(n1);
2313             number m2 = fetch(br,m2);
2314             interval2 = list(-(m2+1),0);
2315             NormalFormEquation F2;
2316             F2.normalFormEquation = nfe2;
2317             F2.parameters = list(list(a*var(2)^(n1),interval2));
2318             F2.singularityType = typeofsing2;
2319             setring br;
2320           }
2321         }
2322         else  //e<0
2323         {
2324           number m = maxabs(P);
2325           MP1 = RealRootFactor(P,-m-1,0);
2326           number m1 = maxabs(MP1);
2327           if(coeff(rf,var(1)^2*var(2)^2)>0)       //b>0
2328           {
2329             typeofsing1  = "Y["+string(n1)+","+string(n2)+"]+-";
2330             typeofsing2  = "Y["+string(n2)+","+string(n1)+"]+-";
2331             if(deg(MP1)==1)
2332             {
2333                number a = -number(subst(MP1,var(1),0));
2334                def R2 = br1;
2335                setring R2;
2336                number a = fetch(br,a);
2337             }
2338             else
2339             {
2340               string mord = ringlist(br1)[3][1][1];
2341               ring R1 = (0,a),(var(1),var(2)),ds;
2342               list L1 = ringlist(R1);
2343               L1[3][1][1]=mord;
2344               def R2 = ring(L1);
2345               setring R2;
2346               kill mord;
2347               kill R1;
2348               poly ff = fetch(br,MP1);
2349               minpoly = number(subst(ff,var(1),a));
2350             }
2351             Poly nfe1 =var(1)^2*var(2)^2-var(1)^(n1)+a*var(2)^(n2);
2352             number m1 = fetch(br,m1);
2353             interval1 = list(-(m1+1),0);
2354             NormalFormEquation F1;
2355             F1.normalFormEquation = nfe1;
2356             F1.parameters = list(list(a*var(2)^(n2),interval1));
2357             F1.singularityType = typeofsing1;
2358
2359             setring br;
2360             kill R2;
2361
2362             if(deg(MP2)==1)
2363             {
2364                 number a = -number(subst(MP2,var(1),0));
2365                 def R2 = br1;
2366                 setring R2;
2367                 number a = fetch(br,a);
2368             }
2369             else
2370             {
2371               string mord = ringlist(br1)[3][1][1];
2372               ring R1 = (0,a),(var(1),var(2)),ds;
2373               list L1 = ringlist(R1);
2374               L1[3][1][1]=mord;
2375               def R2 = ring(L1);
2376               setring R2;
2377               kill mord;
2378               kill R1;
2379               poly ff = fetch(br,MP2);
2380               minpoly = number(subst(ff,var(1),a));
2381             }
2382             Poly nfe2 = var(1)^2*var(2)^2-var(1)^(n2)+a*var(2)^(n1);
2383             number m2 = fetch(br,m2);
2384             interval2 = list(-(m2+1),0);
2385             NormalFormEquation F2;
2386             F2.normalFormEquation = nfe2;
2387             F2.parameters = list(list(a*var(2)^(n1),interval2));
2388             F2.singularityType = typeofsing2;
2389             setring br;
2390           }
2391           else //b<0
2392           {
2393             typeofsing1  = "Y["+string(n1)+","+string(n2)+"]--";
2394             typeofsing2  = "Y["+string(n2)+","+string(n1)+"]--";
2395
2396             if(deg(MP1)==1)
2397             {
2398                number a = -number(subst(MP1,var(1),0));
2399                def R2 = br1;
2400                setring R2;
2401                number a = fetch(br,a);
2402             }
2403             else
2404             {
2405               string mord = ringlist(br1)[3][1][1];
2406               ring R1 = (0,a),(var(1),var(2)),ds;
2407               list L1 = ringlist(R1);
2408               L1[3][1][1]=mord;
2409               def R2 = ring(L1);
2410               setring R2;
2411               kill mord;
2412               kill R1;
2413               poly ff = fetch(br,MP1);
2414               minpoly = number(subst(ff,var(1),a));
2415             }
2416             Poly nfe1 =-var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n2);
2417             number m1 = fetch(br,m1);
2418             interval1 = list(-(m1+1),0);
2419             NormalFormEquation F1;
2420             F1.normalFormEquation = nfe1;
2421             F1.parameters = list(list(a*var(2)^(n2),interval1));
2422             F1.singularityType = typeofsing1;
2423
2424             setring br;
2425             kill R2;
2426
2427             if(deg(MP2)==1)
2428             {
2429                 number a = -number(subst(MP2,var(1),0));
2430                 def R2 = br1;
2431                 setring R2;
2432                 number a = fetch(br,a);
2433             }
2434             else
2435             {
2436               string mord = ringlist(br1)[3][1][1];
2437               ring R1 = (0,a),(var(1),var(2)),ds;
2438               list L1 = ringlist(R1);
2439               L1[3][1][1]=mord;
2440               def R2 = ring(L1);
2441               setring R2;
2442               kill mord;
2443               kill R1;
2444               poly ff = fetch(br,MP2);
2445               minpoly = number(subst(ff,var(1),a));
2446             }
2447             Poly nfe2 = -var(1)^2*var(2)^2-var(1)^(n2)+a*var(2)^(n1);
2448             number m2 = fetch(br,m2);
2449             interval2 = list(-(m2+1),0);
2450             NormalFormEquation F2;
2451             F2.normalFormEquation = nfe2;
2452             F2.parameters = list(list(a*var(2)^(n1),interval2));
2453             F2.singularityType = typeofsing2;
2454             setring br;
2455           }
2456         }
2457       }
2458       if(n1==n2 && sign(coeff(rf,var(1)^n1))==sign(coeff(rf,var(2)^n2)))
2459       {
2460         return(F1);
2461       }
2462       else
2463       {
2464         return(F1,F2);
2465       }
2466      }
2467      if((n1 mod 2 == 0) && (n2 mod 2<>0))  //r is even and s is odd
2468      {
2469        if(coeff(rf,var(1)^2*var(2)^2)>0)  //b>0
2470        {
2471          if(coeff(rf,var(1)^(n1))>0)    //d>0
2472          {
2473            number m = maxabs(P);
2474            MP1 = RealRootFactor(P,-m-1,0);
2475            number m1 = maxabs(MP1);
2476            MP2 = RealRootFactor(P,0,m+1);
2477            number m2 = maxabs(MP2);
2478            MP3 = RealRootFactor(P1,0,m+1);
2479            number m3 = maxabs(MP3);
2480            MP4 = RealRootFactor(P1,0,m+1);
2481            number m4 = maxabs(MP4);
2482
2483            typeofsing1 = "Y["+string(n1)+","+string(n2)+"]++";
2484            typeofsing2 = "Y["+string(n1)+","+string(n2)+"]++";
2485            typeofsing3 = "Y["+string(n2)+","+string(n1)+"]+-";
2486            typeofsing4 = "Y["+string(n2)+","+string(n1)+"]++";
2487
2488
2489             if(deg(MP1)==1)
2490             {
2491                number a = -number(subst(MP1,var(1),0));
2492                def R2 = br1;
2493                setring R2;
2494                number a = fetch(br,a);
2495             }
2496             else
2497             {
2498               string mord = ringlist(br1)[3][1][1];
2499               ring R1 = (0,a),(var(1),var(2)),ds;
2500               list L1 = ringlist(R1);
2501               L1[3][1][1]=mord;
2502               def R2 = ring(L1);
2503               setring R2;
2504               kill mord;
2505               kill R1;
2506               poly ff = fetch(br,MP1);
2507               minpoly = number(subst(ff,var(1),a));
2508             }
2509             Poly nfe1 =var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n2);
2510             number m1 = fetch(br,m1);
2511             interval1 = list(-(m1+1),0);
2512             NormalFormEquation F1;
2513             F1.normalFormEquation = nfe1;
2514             F1.parameters = list(list(a*var(2)^(n2),interval1));
2515             F1.singularityType = typeofsing1;
2516
2517             setring br;
2518             kill R2;
2519
2520             if(deg(MP2)==1)
2521             {
2522                 number a = -number(subst(MP2,var(1),0));
2523                 def R2 = br1;
2524                 setring R2;
2525                 number a = fetch(br,a);
2526             }
2527             else
2528             {
2529               string mord = ringlist(br1)[3][1][1];
2530               ring R1 = (0,a),(var(1),var(2)),ds;
2531               list L1 = ringlist(R1);
2532               L1[3][1][1]=mord;
2533               def R2 = ring(L1);
2534               setring R2;
2535               kill mord;
2536               kill R1;
2537               poly ff = fetch(br,MP2);
2538               minpoly = number(subst(ff,var(1),a));
2539             }
2540             Poly nfe2 = var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n2);
2541             number m2 = fetch(br,m2);
2542             interval2 = list(0,(m2+1));
2543             NormalFormEquation F2;
2544             F2.normalFormEquation = nfe2;
2545             F2.parameters = list(list(a*var(2)^(n2),interval2));
2546             F2.singularityType = typeofsing2;
2547             setring br;
2548             kill R2;
2549
2550             if(deg(MP3)==1)
2551             {
2552                number a = -number(subst(MP3,var(1),0));
2553                def R2 = br1;
2554                setring R2;
2555                number a = fetch(br,a);
2556
2557             }
2558             else
2559             {
2560               string mord = ringlist(br1)[3][1][1];
2561               ring R1 = (0,a),(var(1),var(2)),ds;
2562               list L1 = ringlist(R1);
2563               L1[3][1][1]=mord;
2564               def R2 = ring(L1);
2565               setring R2;
2566               kill mord;
2567               kill R1;
2568               poly ff = fetch(br,MP3);
2569               minpoly = number(subst(ff,var(1),a));
2570             }
2571             Poly nfe3 =var(1)^2*var(2)^2-var(1)^(n2)+a*var(2)^(n1);
2572             number m3 = fetch(br,m3);
2573             interval3 = list(0,(m3+1));
2574             NormalFormEquation F3;
2575             F3.normalFormEquation = nfe3;
2576             F3.parameters = list(list(a*var(2)^(n1),interval3));
2577             F3.singularityType = typeofsing3;
2578
2579             setring br;
2580             kill R2;
2581
2582             if(deg(MP4)==1)
2583             {
2584                 number a = -number(subst(MP4,var(1),0));
2585                 def R2 = br1;
2586                 setring R2;
2587                 number a = fetch(br,a);
2588             }
2589             else
2590             {
2591               string mord = ringlist(br1)[3][1][1];
2592               ring R1 = (0,a),(var(1),var(2)),ds;
2593               list L1 = ringlist(R1);
2594               L1[3][1][1]=mord;
2595               def R2 = ring(L1);
2596               setring R2;
2597               kill mord;
2598               kill R1;
2599               poly ff = fetch(br,MP4);
2600               minpoly = number(subst(ff,var(1),a));
2601             }
2602             Poly nfe4 = var(1)^2*var(2)^2+var(1)^(n2)+a*var(2)^(n1);
2603             number m4 = fetch(br,m4);
2604             interval4 = list(0,(m4+1));
2605             NormalFormEquation F4;
2606             F4.normalFormEquation = nfe4;
2607             F4.parameters = list(list(a*var(2)^(n1),interval4));
2608             F4.singularityType = typeofsing4;
2609             setring br;
2610
2611            return(F1,F2,F3,F4);
2612          }
2613          else  //d<0
2614          {
2615            number m = maxabs(P);
2616            MP1 = RealRootFactor(P,-m-1,0);
2617            number m1 = maxabs(MP1);
2618            MP2 = RealRootFactor(P,0,m+1);
2619            number m2 = maxabs(MP2);
2620            MP3 = RealRootFactor(P1,-m-1,0);
2621            number m3 = maxabs(MP3);
2622            MP4 = RealRootFactor(P1,-m-1,0);
2623            number m4 = maxabs(MP4);
2624
2625            typeofsing1 = "Y["+string(n1)+","+string(n2)+"]+-";
2626            typeofsing2 = "Y["+string(n1)+","+string(n2)+"]+-";
2627            typeofsing3 = "Y["+string(n2)+","+string(n1)+"]+-";
2628            typeofsing4 = "Y["+string(n2)+","+string(n1)+"]++";
2629
2630
2631
2632             if(deg(MP1)==1)
2633             {
2634                number a = -number(subst(MP1,var(1),0));
2635                def R2 = br1;
2636                setring R2;
2637                number a = fetch(br,a);
2638
2639             }
2640             else
2641             {
2642               string mord = ringlist(br1)[3][1][1];
2643               ring R1 = (0,a),(var(1),var(2)),ds;
2644               list L1 = ringlist(R1);
2645               L1[3][1][1]=mord;
2646               def R2 = ring(L1);
2647               setring R2;
2648               kill mord;
2649               kill R1;
2650               poly ff = fetch(br,MP1);
2651               minpoly = number(subst(ff,var(1),a));
2652             }
2653             Poly nfe1 =var(1)^2*var(2)^2-var(1)^(n1)+a*var(2)^(n2);
2654             number m1 = fetch(br,m1);
2655             interval1 = list(-(m1+1),0);
2656             NormalFormEquation F1;
2657             F1.normalFormEquation = nfe1;
2658             F1.parameters = list(list(a*var(2)^(n2),interval1));
2659             F1.singularityType = typeofsing1;
2660
2661             setring br;
2662             kill R2;
2663
2664             if(deg(MP2)==1)
2665             {
2666                 number a = -number(subst(MP2,var(1),0));
2667                 def R2 = br1;
2668                 setring R2;
2669                 number a = fetch(br,a);
2670
2671             }
2672             else
2673             {
2674               string mord = ringlist(br1)[3][1][1];
2675               ring R1 = (0,a),(var(1),var(2)),ds;
2676               list L1 = ringlist(R1);
2677               L1[3][1][1]=mord;
2678               def R2 = ring(L1);
2679               setring R2;
2680               kill mord;
2681               kill R1;
2682               poly ff = fetch(br,MP2);
2683               minpoly = number(subst(ff,var(1),a));
2684             }
2685             Poly nfe2 = var(1)^2*var(2)^2-var(1)^(n1)+a*var(2)^(n2);
2686             number m2 = fetch(br,m2);
2687             interval2 = list(0,(m2+1));
2688             NormalFormEquation F2;
2689             F2.normalFormEquation = nfe2;
2690             F2.parameters = list(list(a*var(2)^(n2),interval2));
2691             F2.singularityType = typeofsing2;
2692             setring br;
2693             kill R2;
2694
2695
2696             if(deg(MP3)==1)
2697             {
2698                number a = -number(subst(MP3,var(1),0));
2699                def R2 = br1;
2700                setring R2;
2701                number a = fetch(br,a);
2702
2703             }
2704             else
2705             {
2706               string mord = ringlist(br1)[3][1][1];
2707               ring R1 = (0,a),(var(1),var(2)),ds;
2708               list L1 = ringlist(R1);
2709               L1[3][1][1]=mord;
2710               def R2 = ring(L1);
2711               setring R2;
2712               kill mord;
2713               kill R1;
2714               poly ff = fetch(br,MP3);
2715               minpoly = number(subst(ff,var(1),a));
2716             }
2717             Poly nfe3 =var(1)^2*var(2)^2-var(1)^(n2)+a*var(2)^(n1);
2718             number m3 = fetch(br,m3);
2719             interval3 = list(-(m3+1),0);
2720             NormalFormEquation F3;
2721             F3.normalFormEquation = nfe3;
2722             F3.parameters = list(list(a*var(2)^(n2),interval3));
2723             F3.singularityType = typeofsing3;
2724
2725             setring br;
2726             kill R2;
2727
2728             if(deg(MP4)==1)
2729             {
2730                 number a = -number(subst(MP4,var(1),0));
2731                 def R2 = br1;
2732                 setring R2;
2733                 number a = fetch(br,a);
2734
2735             }
2736             else
2737             {
2738               string mord = ringlist(br1)[3][1][1];
2739               ring R1 = (0,a),(var(1),var(2)),ds;
2740               list L1 = ringlist(R1);
2741               L1[3][1][1]=mord;
2742               def R2 = ring(L1);
2743               setring R2;
2744               kill mord;
2745               kill R1;
2746               poly ff = fetch(br,MP4);
2747               minpoly = number(subst(ff,var(1),a));
2748             }
2749             Poly nfe4 = var(1)^2*var(2)^2+var(1)^(n2)+a*var(2)^(n1);
2750             number m4 = fetch(br,m4);
2751             interval4 = list(-(m4+1),0);
2752             NormalFormEquation F4;
2753             F4.normalFormEquation = nfe4;
2754             F4.parameters = list(list(a*var(2)^(n1),interval4));
2755             F4.singularityType = typeofsing4;
2756             setring br;
2757             kill R2;
2758
2759            return(F1,F2,F3,F4);
2760          }
2761        }
2762        else
2763        {
2764         if(coeff(rf,var(1)^(n1))<0)    //d<0
2765         {
2766            number m = maxabs(P);
2767            MP1 = RealRootFactor(P,-m-1,0);
2768            number m1 = maxabs(MP1);
2769            MP2 = RealRootFactor(P,0,m+1);
2770            number m2 = maxabs(MP2);
2771            MP3 = RealRootFactor(P1,-(m+1),0);
2772            number m3 = maxabs(MP3);
2773            MP4 = RealRootFactor(P1,-(m+1),0);
2774            number m3 = maxabs(MP3);
2775            typeofsing1 = "Y["+string(n1)+","+string(n2)+"]-+";
2776            typeofsing2 = "Y["+string(n1)+","+string(n2)+"]-+";
2777            typeofsing3 = "Y["+string(n2)+","+string(n1)+"]--";
2778            typeofsing4 = "Y["+string(n2)+","+string(n1)+"]-+";
2779
2780
2781
2782             if(deg(MP1)==1)
2783             {
2784                number a = -number(subst(MP1,var(1),0));
2785                def R2 = br1;
2786                setring R2;
2787                number a = fetch(br,a);
2788
2789             }
2790             else
2791             {
2792               string mord = ringlist(br1)[3][1][1];
2793               ring R1 = (0,a),(var(1),var(2)),ds;
2794               list L1 = ringlist(R1);
2795               L1[3][1][1]=mord;
2796               def R2 = ring(L1);
2797               setring R2;
2798               kill mord;
2799               kill R1;
2800               poly ff = fetch(br,MP1);
2801               minpoly = number(subst(ff,var(1),a));
2802             }
2803             Poly nfe1 =-var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n2);
2804             number m1 = fetch(br,m1);
2805             interval1 = list(-(m1+1),0);
2806             NormalFormEquation F1;
2807             F1.normalFormEquation = nfe1;
2808             F1.parameters = list(list(a*var(2)^(n2),interval1));
2809             F1.singularityType = typeofsing1;
2810
2811             setring br;
2812             kill R2;
2813
2814             if(deg(MP2)==1)
2815             {
2816                 number a = -number(subst(MP2,var(1),0));
2817                 def R2 = br1;
2818                 setring R2;
2819                 number a = fetch(br,a);
2820
2821             }
2822             else
2823             {
2824               string mord = ringlist(br1)[3][1][1];
2825               ring R1 = (0,a),(var(1),var(2)),ds;
2826               list L1 = ringlist(R1);
2827               L1[3][1][1]=mord;
2828               def R2 = ring(L1);
2829               setring R2;
2830               kill mord;
2831               kill R1;
2832               poly ff = fetch(br,MP2);
2833               minpoly = number(subst(ff,var(1),a));
2834             }
2835             Poly nfe2 = -var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n2);
2836             number m2 = fetch(br,m2);
2837             interval2 = list(0,(m2+1));
2838             NormalFormEquation F2;
2839             F2.normalFormEquation = nfe2;
2840             F2.parameters = list(list(a*var(2)^(n2),interval2));
2841             F2.singularityType = typeofsing2;
2842             setring br;
2843             kill R2;
2844
2845             if(deg(MP3)==1)
2846             {
2847                number a = -number(subst(MP3,var(1),0));
2848                def R2 = br1;
2849                setring R2;
2850                number a = fetch(br,a);
2851
2852             }
2853             else
2854             {
2855               string mord = ringlist(br1)[3][1][1];
2856               ring R1 = (0,a),(var(1),var(2)),ds;
2857               list L1 = ringlist(R1);
2858               L1[3][1][1]=mord;
2859               def R2 = ring(L1);
2860               setring R2;
2861               kill mord;
2862               kill R1;
2863               poly ff = fetch(br,MP3);
2864               minpoly = number(subst(ff,var(1),a));
2865             }
2866             Poly nfe3 =-var(1)^2*var(2)^2-var(1)^(n2)+a*var(2)^(n1);
2867             number m3 = fetch(br,m3);
2868             interval3 = list(0,(m3+1));
2869             NormalFormEquation F3;
2870             F3.normalFormEquation = nfe3;
2871             F3.parameters = list(list(a*var(2)^(n1),interval3));
2872             F3.singularityType = typeofsing3;
2873
2874             setring br;
2875             kill R2;
2876
2877             if(deg(MP4)==1)
2878             {
2879                 number a = -number(subst(MP4,var(1),0));
2880                 def R2 = br1;
2881                 setring R2;
2882                 number a = fetch(br,a);
2883
2884             }
2885             else
2886             {
2887               string mord = ringlist(br1)[3][1][1];
2888               ring R1 = (0,a),(var(1),var(2)),ds;
2889               list L1 = ringlist(R1);
2890               L1[3][1][1]=mord;
2891               def R2 = ring(L1);
2892               setring R2;
2893               kill mord;
2894               kill R1;
2895               poly ff = fetch(br,MP4);
2896               minpoly = number(subst(ff,var(1),a));
2897             }
2898             Poly nfe4 = -var(1)^2*var(2)^2+var(1)^(n2)+a*var(2)^(n1);
2899             number m4 = fetch(br,m4);
2900             interval4 = list(0,(m4+1));
2901             NormalFormEquation F4;
2902             F4.normalFormEquation = nfe4;
2903             F4.parameters = list(list(a*var(2)^(n1),interval4));
2904             F4.singularityType = typeofsing4;
2905             setring br;
2906            return(F1,F2,F3,F4);
2907         }
2908         else  //d>0
2909         {
2910            number m = maxabs(P);
2911            MP1 = RealRootFactor(P,-m-1,0);
2912            number m1 = maxabs(MP1);
2913            MP2 = RealRootFactor(P,0,m+1);
2914            number m2 = maxabs(MP2);
2915            MP3 = RealRootFactor(P1,0,m+1);
2916            number m3 = maxabs(MP3);
2917            MP4 = RealRootFactor(P1,0,m+1);
2918            number m4 = maxabs(MP4);
2919            interval3 = "(-infinity,0)";
2920            interval4 = "(-infinity,0)";
2921
2922            typeofsing1 = "Y["+string(n1)+","+string(n2)+"]-+";
2923            typeofsing2 = "Y["+string(n1)+","+string(n2)+"]-+";
2924            typeofsing3 = "Y["+string(n2)+","+string(n1)+"]--";
2925            typeofsing4 = "Y["+string(n2)+","+string(n1)+"]-+";
2926
2927             if(deg(MP1)==1)
2928             {
2929                number a = -number(subst(MP1,var(1),0));
2930                def R2 = br1;
2931                setring R2;
2932                number a = fetch(br,a);
2933
2934             }
2935             else
2936             {
2937               string mord = ringlist(br1)[3][1][1];
2938               ring R1 = (0,a),(var(1),var(2)),ds;
2939               list L1 = ringlist(R1);
2940               L1[3][1][1]=mord;
2941               def R2 = ring(L1);
2942               setring R2;
2943               kill mord;
2944               kill R1;
2945               poly ff = fetch(br,MP1);
2946               minpoly = number(subst(ff,var(1),a));
2947             }
2948             Poly nfe1 =-var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n2);
2949             number m1 = fetch(br,m1);
2950             interval1 = list(-(m1+1),0);
2951             NormalFormEquation F1;
2952             F1.normalFormEquation = nfe1;
2953             F1.parameters = list(list(a*var(2)^(n2),interval1));
2954             F1.singularityType = typeofsing1;
2955
2956             setring br;
2957             kill R2;
2958
2959             if(deg(MP2)==1)
2960             {
2961                 number a = -number(subst(MP2,var(1),0));
2962                 def R2 = br1;
2963                 setring R2;
2964                 number a = fetch(br,a);
2965
2966             }
2967             else
2968             {
2969               string mord = ringlist(br1)[3][1][1];
2970               ring R1 = (0,a),(var(1),var(2)),ds;
2971               list L1 = ringlist(R1);
2972               L1[3][1][1]=mord;
2973               def R2 = ring(L1);
2974               setring R2;
2975               kill mord;
2976               kill R1;
2977               poly ff = fetch(br,MP2);
2978               minpoly = number(subst(ff,var(1),a));
2979             }
2980             Poly nfe2 =-var(1)^2*var(2)^2+var(1)^(n1)+var(2)^(n2);
2981             number m2 = fetch(br,m2);
2982             interval2 = list(0,(m2+1));
2983             NormalFormEquation F2;
2984             F2.normalFormEquation = nfe2;
2985             F2.parameters = list(list(a*var(2)^(n2),interval2));
2986             F2.singularityType = typeofsing2;
2987             setring br;
2988             kill R2;
2989
2990
2991             if(deg(MP3)==1)
2992             {
2993                number a = -number(subst(MP3,var(1),0));
2994                def R2 = br1;
2995                setring R2;
2996                number a = fetch(br,a);
2997
2998             }
2999             else
3000             {
3001               string mord = ringlist(br1)[3][1][1];
3002               ring R1 = (0,a),(var(1),var(2)),ds;
3003               list L1 = ringlist(R1);
3004               L1[3][1][1]=mord;
3005               def R2 = ring(L1);
3006               setring R2;
3007               kill mord;
3008               kill R1;
3009               poly ff = fetch(br,MP3);
3010               minpoly = number(subst(ff,var(1),a));
3011             }
3012             Poly nfe3 =-var(1)^2*var(2)^2-var(1)^(n2)+a*var(2)^(n1);
3013             number m3 = fetch(br,m3);
3014             interval3 = list(-(m3+1),0);
3015             NormalFormEquation F3;
3016             F3.normalFormEquation = nfe3;
3017             F3.parameters = list(list(a*var(2)^(n1),interval3));
3018             F3.singularityType = typeofsing3;
3019
3020             setring br;
3021             kill R2;
3022
3023             if(deg(MP4)==1)
3024             {
3025                 number a = -number(subst(MP4,var(1),0));
3026                 def R2 = br1;
3027                 setring R2;
3028                 number a = fetch(br,a);
3029
3030             }
3031             else
3032             {
3033               string mord = ringlist(br1)[3][1][1];
3034               ring R1 = (0,a),(var(1),var(2)),ds;
3035               list L1 = ringlist(R1);
3036               L1[3][1][1]=mord;
3037               def R2 = ring(L1);
3038               setring R2;
3039               kill mord;
3040               kill R1;
3041               poly ff = fetch(br,MP4);
3042               minpoly = number(subst(ff,var(1),a));
3043             }
3044             Poly nfe4 = -var(1)^2*var(2)^2+var(1)^(n2)+a*var(2)^(n1);
3045             number m4 = fetch(br,m4);
3046             interval4 = list(-(m4+1),0);
3047             NormalFormEquation F4;
3048             F4.normalFormEquation = nfe4;
3049             F4.parameters = list(list(a*var(2)^(n1),interval4));
3050             F4.singularityType = typeofsing4;
3051             setring br;
3052
3053            return(F1,F2,F3,F4);
3054         }
3055       }
3056     }
3057     if((n1 mod 2 <> 0) && (n2 mod 2 <> 0))  // r and s are odd
3058     {
3059       if(coeff(rf,var(1)^2*var(2)^2)>0)
3060       {
3061            number m = maxabs(P);
3062            MP1 = RealRootFactor(P,-m-1,0);
3063            number m1 = maxabs(MP1);
3064            MP3 = RealRootFactor(P,-m-1,0);
3065            number m3 = maxabs(MP3);
3066            MP2 = RealRootFactor(P,0,m+1);
3067            number m2 = maxabs(MP2);
3068            MP4 = RealRootFactor(P,0,m+1);
3069            number m4 = maxabs(MP4);
3070            m = maxabs(P1);
3071            MP5 = RealRootFactor(P1,-m-1,0);
3072            number m5 = maxabs(MP5);
3073            MP7 = RealRootFactor(P1,-m-1,0);
3074            number m7 = maxabs(MP7);
3075            MP6 = RealRootFactor(P1,0,m+1);
3076            number m6 = maxabs(MP6);
3077            MP8 = RealRootFactor(P1,0,m+1);
3078            number m8 = maxabs(MP8);
3079
3080            typeofsing1 = "Y["+string(n1)+","+string(n2)+"]+-";
3081            typeofsing2 = "Y["+string(n1)+","+string(n2)+"]+-";
3082            typeofsing3 = "Y["+string(n1)+","+string(n2)+"]++";
3083            typeofsing4 = "Y["+string(n1)+","+string(n2)+"]++";
3084            typeofsing5 = "Y["+string(n2)+","+string(n1)+"]+-";
3085            typeofsing6 = "Y["+string(n2)+","+string(n1)+"]+-";
3086            typeofsing7 = "Y["+string(n2)+","+string(n1)+"]++";
3087            typeofsing8 = "Y["+string(n2)+","+string(n1)+"]++";
3088
3089             if(deg(MP1)==1)
3090             {
3091                number a = -number(subst(MP1,var(1),0));
3092                def R2 = br1;
3093                setring R2;
3094                number a = fetch(br,a);
3095
3096             }
3097             else
3098             {
3099               string mord = ringlist(br1)[3][1][1];
3100               ring R1 = (0,a),(var(1),var(2)),ds;
3101               list L1 = ringlist(R1);
3102               L1[3][1][1]=mord;
3103               def R2 = ring(L1);
3104               setring R2;
3105               kill mord;
3106               kill R1;
3107               poly ff = fetch(br,MP1);
3108               minpoly = number(subst(ff,var(1),a));
3109             }
3110             Poly nfe1 =var(1)^2*var(2)^2-var(1)^(n1)+a*var(2)^(n2);
3111             number m1 = fetch(br,m1);
3112             interval1 = list(-(m1+1),0);
3113             NormalFormEquation F1;
3114             F1.normalFormEquation = nfe1;
3115             F1.parameters = list(list(a*var(2)^(n2),interval1));
3116             F1.singularityType = typeofsing1;
3117
3118             setring br;
3119             kill R2;
3120
3121             if(deg(MP2)==1)
3122             {
3123                 number a = -number(subst(MP2,var(1),0));
3124                 def R2 = br1;
3125                 setring R2;
3126                 number a = fetch(br,a);
3127
3128             }
3129             else
3130             {
3131               string mord = ringlist(br1)[3][1][1];
3132               ring R1 = (0,a),(var(1),var(2)),ds;
3133               list L1 = ringlist(R1);
3134               L1[3][1][1]=mord;
3135               def R2 = ring(L1);
3136               setring R2;
3137               kill mord;
3138               kill R1;
3139               poly ff = fetch(br,MP2);
3140               minpoly = number(subst(ff,var(1),a));
3141             }
3142             Poly nfe2 =var(1)^2*var(2)^2-var(1)^(n1)+a*var(2)^(n2);
3143             number m2 = fetch(br,m2);
3144             interval2 = list(0,(m2+1));
3145             NormalFormEquation F2;
3146             F2.normalFormEquation = nfe2;
3147             F2.parameters = list(list(a*var(2)^(n2),interval2));
3148             F2.singularityType = typeofsing2;
3149             setring br;
3150             kill R2;
3151
3152
3153             if(deg(MP3)==1)
3154             {
3155                number a = -number(subst(MP3,var(1),0));
3156                def R2 = br1;
3157                setring R2;
3158                number a = fetch(br,a);
3159
3160             }
3161             else
3162             {
3163               string mord = ringlist(br1)[3][1][1];
3164               ring R1 = (0,a),(var(1),var(2)),ds;
3165               list L1 = ringlist(R1);
3166               L1[3][1][1]=mord;
3167               def R2 = ring(L1);
3168               setring R2;
3169               kill mord;
3170               kill R1;
3171               poly ff = fetch(br,MP3);
3172               minpoly = number(subst(ff,var(1),a));
3173             }
3174             Poly nfe3 =var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n2);
3175             number m3 = fetch(br,m3);
3176             interval3 = list(-(m3+1),0);
3177             NormalFormEquation F3;
3178             F3.normalFormEquation = nfe3;
3179             F3.parameters = list(list(a*var(2)^(n2),interval3));
3180             F3.singularityType = typeofsing3;
3181
3182             setring br;
3183             kill R2;
3184
3185             if(deg(MP4)==1)
3186             {
3187                 number a = -number(subst(MP2,var(1),0));
3188                 def R2 = br1;
3189                setring R2;
3190                number a = fetch(br,a);
3191
3192             }
3193             else
3194             {
3195               string mord = ringlist(br1)[3][1][1];
3196               ring R1 = (0,a),(var(1),var(2)),ds;
3197               list L1 = ringlist(R1);
3198               L1[3][1][1]=mord;
3199               def R2 = ring(L1);
3200               setring R2;
3201               kill mord;
3202               kill R1;
3203               poly ff = fetch(br,MP4);
3204               minpoly = number(subst(ff,var(1),a));
3205             }
3206             Poly nfe4 = var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n2);
3207             number m4 = fetch(br,m4);
3208             interval4 = list(0,(m4+1));
3209             NormalFormEquation F4;
3210             F4.normalFormEquation = nfe4;
3211             F4.parameters = list(list(a*var(2)^(n2),interval4));
3212             F4.singularityType = typeofsing4;
3213             setring br;
3214             kill R2;
3215
3216             if(deg(MP5)==1)
3217             {
3218                number a = -number(subst(MP5,var(1),0));
3219                def R2 = br1;
3220                setring R2;
3221                number a = fetch(br,a);
3222
3223             }
3224             else
3225             {
3226               string mord = ringlist(br1)[3][1][1];
3227               ring R1 = (0,a),(var(1),var(2)),ds;
3228               list L1 = ringlist(R1);
3229               L1[3][1][1]=mord;
3230               def R2 = ring(L1);
3231               setring R2;
3232               kill mord;
3233               kill R1;
3234               poly ff = fetch(br,MP5);
3235               minpoly = number(subst(ff,var(1),a));
3236             }
3237             Poly nfe5 =var(1)^2*var(2)^2-var(1)^(n2)+a*var(2)^(n1);
3238             number m5 = fetch(br,m5);
3239             interval5 = list(-(m5+1),0);
3240             NormalFormEquation F5;
3241             F5.normalFormEquation = nfe5;
3242             F5.parameters = list(list(a*var(2)^(n1),interval5));
3243             F5.singularityType = typeofsing5;
3244
3245             setring br;
3246             kill R2;
3247
3248             if(deg(MP6)==1)
3249             {
3250                 number a = -number(subst(MP6,var(1),0));
3251                 def R2 = br1;
3252                setring R2;
3253                number a = fetch(br,a);
3254
3255             }
3256             else
3257             {
3258               string mord = ringlist(br1)[3][1][1];
3259               ring R1 = (0,a),(var(1),var(2)),ds;
3260               list L1 = ringlist(R1);
3261               L1[3][1][1]=mord;
3262               def R2 = ring(L1);
3263               setring R2;
3264               kill mord;
3265               kill R1;
3266               poly ff = fetch(br,MP6);
3267               minpoly = number(subst(ff,var(1),a));
3268             }
3269             Poly nfe6 =var(1)^2*var(2)^2-var(1)^(n2)+a*var(2)^(n1);
3270             number m6 = fetch(br,m6);
3271             interval6 = list(0,(m6+1));
3272             NormalFormEquation F6;
3273             F6.normalFormEquation = nfe6;
3274             F6.parameters = list(list(a*var(2)^(n1),interval6));
3275             F6.singularityType = typeofsing6;
3276             setring br;
3277             kill R2;
3278
3279
3280             if(deg(MP7)==1)
3281             {
3282                number a = -number(subst(MP7,var(1),0));
3283                def R2 = br1;
3284                setring R2;
3285                number a = fetch(br,a);
3286
3287             }
3288             else
3289             {
3290               string mord = ringlist(br1)[3][1][1];
3291               ring R1 = (0,a),(var(1),var(2)),ds;
3292               list L1 = ringlist(R1);
3293               L1[3][1][1]=mord;
3294               def R2 = ring(L1);
3295               setring R2;
3296               kill mord;
3297               kill R1;
3298               poly ff = fetch(br,MP7);
3299               minpoly = number(subst(ff,var(1),a));
3300             }
3301             Poly nfe7 =var(1)^2*var(2)^2+var(1)^(n2)+a*var(2)^(n1);
3302             number m7 = fetch(br,m7);
3303             interval7 = list(-(m7+1),0);
3304             NormalFormEquation F7;
3305             F7.normalFormEquation = nfe7;
3306             F7.parameters = list(list(a*var(2)^(n1),interval7));
3307             F7.singularityType = typeofsing7;
3308
3309             setring br;
3310             kill R2;
3311
3312             if(deg(MP8)==1)
3313             {
3314                 number a = -number(subst(MP8,var(1),0));
3315                 def R2 = br1;
3316                setring R2;
3317                number a = fetch(br,a);
3318
3319             }
3320             else
3321             {
3322               string mord = ringlist(br1)[3][1][1];
3323               ring R1 = (0,a),(var(1),var(2)),ds;
3324               list L1 = ringlist(R1);
3325               L1[3][1][1]=mord;
3326               def R2 = ring(L1);
3327               setring R2;
3328               kill mord;
3329               kill R1;
3330               poly ff = fetch(br,MP8);
3331               minpoly = number(subst(ff,var(1),a));
3332             }
3333             Poly nfe8 = var(1)^2*var(2)^2+var(1)^(n2)+a*var(2)^(n1);
3334             number m8 = fetch(br,m8);
3335             interval8 = list(0,(m8+1));
3336             NormalFormEquation F8;
3337             F8.normalFormEquation = nfe8;
3338             F8.parameters = list(list(a*var(2)^(n1),interval8));
3339             F8.singularityType = typeofsing8;
3340             setring br;
3341       }
3342       else
3343       {
3344            number m = maxabs(P);
3345            MP1 = RealRootFactor(P,-m-1,0);
3346            number m1 = maxabs(MP1);
3347            MP3 = RealRootFactor(P,-m-1,0);
3348            number m2 = maxabs(MP2);
3349            MP2 = RealRootFactor(P,0,m+1);
3350            number m3 = maxabs(MP3);
3351            MP4 = RealRootFactor(P,0,m+1);
3352            number m4 = maxabs(MP4);
3353            MP5 = RealRootFactor(P1,-m-1,0);
3354            number m5 = maxabs(MP5);
3355            MP7 = RealRootFactor(P1,-m-1,0);
3356            number m7 = maxabs(MP7);
3357            MP6 = RealRootFactor(P1,0,m+1);
3358            number m6 = maxabs(MP6);
3359            MP8 = RealRootFactor(P1,0,m+1);
3360            number m8 = maxabs(MP8);
3361            typeofsing1 = "Y["+string(n1)+","+string(n2)+"]--";
3362            typeofsing2 = "Y["+string(n1)+","+string(n2)+"]--";
3363            typeofsing3 = "Y["+string(n1)+","+string(n2)+"]-+";
3364            typeofsing4 = "Y["+string(n1)+","+string(n2)+"]-+";
3365            typeofsing5 = "Y["+string(n2)+","+string(n1)+"]--";
3366            typeofsing6 = "Y["+string(n2)+","+string(n1)+"]--";
3367            typeofsing7 = "Y["+string(n2)+","+string(n1)+"]-+";
3368            typeofsing8 = "Y["+string(n2)+","+string(n1)+"]-+";
3369             if(deg(MP1)==1)
3370             {
3371                number a = -number(subst(MP1,var(1),0));
3372                def R2 = br1;
3373                setring R2;
3374                number a = fetch(br,a);
3375
3376             }
3377             else
3378             {
3379               string mord = ringlist(br1)[3][1][1];
3380               ring R1 = (0,a),(var(1),var(2)),ds;
3381               list L1 = ringlist(R1);
3382               L1[3][1][1]=mord;
3383               def R2 = ring(L1);
3384               setring R2;
3385               kill mord;
3386               kill R1;
3387               poly ff = fetch(br,MP1);
3388               minpoly = number(subst(ff,var(1),a));
3389             }
3390             Poly nfe1 =-var(1)^2*var(2)^2-var(1)^(n1)+a*var(2)^(n2);
3391             number m1 = fetch(br,m1);
3392             interval1 = list(-(m1+1),0);
3393             NormalFormEquation F1;
3394             F1.normalFormEquation = nfe1;
3395             F1.parameters = list(list(a*var(2)^(n2),interval1));
3396             F1.singularityType = typeofsing1;
3397
3398             setring br;
3399             kill R2;
3400
3401             if(deg(MP2)==1)
3402             {
3403                 number a = -number(subst(MP2,var(1),0));
3404                 def R2 = br1;
3405                setring R2;
3406                number a = fetch(br,a);
3407
3408             }
3409             else
3410             {
3411               string mord = ringlist(br1)[3][1][1];
3412               ring R1 = (0,a),(var(1),var(2)),ds;
3413               list L1 = ringlist(R1);
3414               L1[3][1][1]=mord;
3415               def R2 = ring(L1);
3416               setring R2;
3417               kill mord;
3418               kill R1;
3419               poly ff = fetch(br,MP2);
3420               minpoly = number(subst(ff,var(1),a));
3421             }
3422             Poly nfe2 =-var(1)^2*var(2)^2-var(1)^(n1)+a*var(2)^(n2);
3423             number m2 = fetch(br,m2);
3424             interval2 = list(0,(m2+1));
3425             NormalFormEquation F2;
3426             F2.normalFormEquation = nfe2;
3427             F2.parameters = list(list(a*var(2)^(n2),interval2));
3428             F2.singularityType = typeofsing2;
3429             setring br;
3430             kill R2;
3431
3432
3433             if(deg(MP3)==1)
3434             {
3435                number a = -number(subst(MP3,var(1),0));
3436                def R2 = br1;
3437                setring R2;
3438                number a = fetch(br,a);
3439
3440             }
3441             else
3442             {
3443               string mord = ringlist(br1)[3][1][1];
3444               ring R1 = (0,a),(var(1),var(2)),ds;
3445               list L1 = ringlist(R1);
3446               L1[3][1][1]=mord;
3447               def R2 = ring(L1);
3448               setring R2;
3449               kill mord;
3450               kill R1;
3451               poly ff = fetch(br,MP3);
3452               minpoly = number(subst(ff,var(1),a));
3453             }
3454             Poly nfe3 =-var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n2);
3455             number m3 = fetch(br,m3);
3456             interval3 = list(-(m3+1),0);
3457             NormalFormEquation F3;
3458             F3.normalFormEquation = nfe3;
3459             F3.parameters = list(list(a*var(2)^(n2),interval3));
3460             F3.singularityType = typeofsing3;
3461
3462             setring br;
3463             kill R2;
3464
3465             if(deg(MP4)==1)
3466             {
3467                 number a = -number(subst(MP2,var(1),0));
3468                 def R2 = br1;
3469                setring R2;
3470                number a = fetch(br,a);
3471
3472             }
3473             else
3474             {
3475               string mord = ringlist(br1)[3][1][1];
3476               ring R1 = (0,a),(var(1),var(2)),ds;
3477               list L1 = ringlist(R1);
3478               L1[3][1][1]=mord;
3479               def R2 = ring(L1);
3480               setring R2;
3481               kill mord;
3482               kill R1;
3483               poly ff = fetch(br,MP4);
3484               minpoly = number(subst(ff,var(1),a));
3485             }
3486             Poly nfe4 = -var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n2);
3487             number m4 = fetch(br,m4);
3488             interval4 = list(0,(m4+1));
3489             NormalFormEquation F4;
3490             F4.normalFormEquation = nfe4;
3491             F4.parameters = list(list(a*var(2)^(n2),interval4));
3492             F4.singularityType = typeofsing4;
3493             setring br;
3494             kill R2;
3495
3496             if(deg(MP5)==1)
3497             {
3498                number a = -number(subst(MP5,var(1),0));
3499                def R2 = br1;
3500                setring R2;
3501                number a = fetch(br,a);
3502
3503             }
3504             else
3505             {
3506               string mord = ringlist(br1)[3][1][1];
3507               ring R1 = (0,a),(var(1),var(2)),ds;
3508               list L1 = ringlist(R1);
3509               L1[3][1][1]=mord;
3510               def R2 = ring(L1);
3511               setring R2;
3512               kill mord;
3513               kill R1;
3514               poly ff = fetch(br,MP5);
3515               minpoly = number(subst(ff,var(1),a));
3516             }
3517             Poly nfe5 =-var(1)^2*var(2)^2-var(1)^(n2)+a*var(2)^(n1);
3518             number m5 = fetch(br,m5);
3519             interval5 = list(-(m5+1),0);
3520             NormalFormEquation F5;
3521             F5.normalFormEquation = nfe5;
3522             F5.parameters = list(list(a*var(2)^(n1),interval5));
3523             F5.singularityType = typeofsing5;
3524
3525             setring br;
3526             kill R2;
3527
3528             if(deg(MP6)==1)
3529             {
3530                 number a = -number(subst(MP6,var(1),0));
3531                 def R2 = br1;
3532                setring R2;
3533                number a = fetch(br,a);
3534
3535             }
3536             else
3537             {
3538               string mord = ringlist(br1)[3][1][1];
3539               ring R1 = (0,a),(var(1),var(2)),ds;
3540               list L1 = ringlist(R1);
3541               L1[3][1][1]=mord;
3542               def R2 = ring(L1);
3543               setring R2;
3544               kill mord;
3545               kill R1;
3546               poly ff = fetch(br,MP6);
3547               minpoly = number(subst(ff,var(1),a));
3548             }
3549             Poly nfe6 =-var(1)^2*var(2)^2-var(1)^(n2)+a*var(2)^(n1);
3550             number m6 = fetch(br,m6);
3551             interval6 = list(0,(m6+1));
3552             NormalFormEquation F6;
3553             F6.normalFormEquation = nfe6;
3554             F6.parameters = list(list(a*var(2)^(n1),interval6));
3555             F6.singularityType = typeofsing6;
3556             setring br;
3557             kill R2;
3558
3559
3560             if(deg(MP7)==1)
3561             {
3562                number a = -number(subst(MP7,var(1),0));
3563                def R2 = br1;
3564                setring R2;
3565                number a = fetch(br,a);
3566             }
3567             else
3568             {
3569               string mord = ringlist(br1)[3][1][1];
3570               ring R1 = (0,a),(var(1),var(2)),ds;
3571               list L1 = ringlist(R1);
3572               L1[3][1][1]=mord;
3573               def R2 = ring(L1);
3574               setring R2;
3575               kill mord;
3576               kill R1;
3577               poly ff = fetch(br,MP7);
3578               minpoly = number(subst(ff,var(1),a));
3579             }
3580             Poly nfe7 =-var(1)^2*var(2)^2+var(1)^(n2)+a*var(2)^(n1);
3581             number m7 = fetch(br,m7);
3582             interval7 = list(-(m7+1),0);
3583             NormalFormEquation F7;
3584             F7.normalFormEquation = nfe7;
3585             F7.parameters = list(list(a*var(2)^(n1),interval7));
3586             F7.singularityType = typeofsing7;
3587
3588             setring br;
3589             kill R2;
3590
3591             if(deg(MP8)==1)
3592             {
3593                 number a = -number(subst(MP8,var(1),0));
3594                 def R2 = br1;
3595                setring R2;
3596                number a = fetch(br,a);
3597
3598             }
3599             else
3600             {
3601               string mord = ringlist(br1)[3][1][1];
3602               ring R1 = (0,a),(var(1),var(2)),ds;
3603               list L1 = ringlist(R1);
3604               L1[3][1][1]=mord;
3605               def R2 = ring(L1);
3606               setring R2;
3607               kill mord;
3608               kill R1;
3609               poly ff = fetch(br,MP8);
3610               minpoly = number(subst(ff,var(1),a));
3611             }
3612             Poly nfe8 = -var(1)^2*var(2)^2+var(1)^(n2)+a*var(2)^(n1);
3613             number m8 = fetch(br,m8);
3614             interval8 = list(0,(m8+1));
3615             NormalFormEquation F8;
3616             F8.normalFormEquation = nfe8;
3617             F8.parameters = list(list(a*var(2)^(n1),interval8));
3618             F8.singularityType = typeofsing8;
3619             setring br;
3620             kill R2;
3621       }
3622       if(n1==n2)  // r=s
3623       {
3624            return(F1,F2,F3,F4);
3625       }
3626       else  // r not equal to s
3627       {
3628            return(F1,F2,F3,F4,F5,F6,F7,F8);
3629       }
3630     }
3631    }
3632    else // h has no linear factors over Q
3633    {
3634      if(sigma == 1)
3635      {
3636        poly de = (coeff(rf,var(1)^n1)*coeff(rf,var(2)^n1));
3637        ring S = 0,(x,y),ds;
3638        poly de = fetch(br,de);
3639        number t0 = number(de);
3640        if(n1 mod 2 <>0)
3641        {
3642          int r = 1;
3643        }
3644        else
3645        {
3646          int r = 0;
3647        }
3648        if(t0>0)
3649        {
3650          int t = 1;
3651        }
3652        else
3653        {
3654          int t = -1;
3655        }
3656        setring br;
3657        if(r==1)  //r odd
3658        {
3659         "r odd";
3660         number s = (coeff(rf,var(1)^2*var(2)^2));
3661         s = s^(-n1);
3662         poly P = var(1)^2-s*(coeff(rf,var(1)^(n1))*coeff(rf,var(2)^(n1)))^2;
3663         setring R;
3664         poly P = fetch(br,P);
3665         typeofsing1 = "Y["+string(n1)+","+string(n1)+"]++";
3666         typeofsing2 = "Y["+string(n1)+","+string(n1)+"]++";
3667         typeofsing3 = "Y["+string(n1)+","+string(n1)+"]+-";
3668         typeofsing4 = "Y["+string(n1)+","+string(n1)+"]+-";
3669         number m = maxabs(P);
3670         MP1 = RealRootFactor(P,-m-1,0);
3671         number m1 = maxabs(MP1);
3672         MP2 = RealRootFactor(P,0,m+1);
3673         number m2 = maxabs(MP2);
3674         MP3 = RealRootFactor(P,-m-1,0);
3675         number m3 = maxabs(MP3);
3676         MP4 = RealRootFactor(P,0,m+1);
3677         number m4 = maxabs(MP4);
3678         if(deg(MP1)==1)
3679         {
3680              number a = -number(subst(MP1,var(1),0));
3681              def R2 = br1;
3682              setring R2;
3683              number a = fetch(R,a);
3684         }
3685         else
3686         {
3687             string mord = ringlist(br1)[3][1][1];
3688             ring R1 = (0,a),(var(1),var(2)),ds;
3689             list L1 = ringlist(R1);
3690             L1[3][1][1]=mord;
3691             def R2 = ring(L1);
3692               setring R2;
3693               kill mord;
3694               kill R1;
3695             poly ff = fetch(R,MP1);
3696             minpoly = number(subst(ff,var(1),a));
3697         }
3698         Poly nfe1 =var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n1);
3699         number m1 = fetch(R,m1);
3700         interval1 = list(-(m1+1),0);
3701         NormalFormEquation F1;
3702         F1.normalFormEquation = nfe1;
3703         F1.parameters = list(list(a*var(2)^(n1),interval1));
3704         F1.singularityType = typeofsing1;
3705
3706         setring R;
3707         kill R2;
3708
3709         if(deg(MP2)==1)
3710         {
3711             number a = -number(subst(MP2,var(1),0));
3712             def R2 = br1;
3713              setring R2;
3714              number a = fetch(R,a);
3715         }
3716         else
3717         {
3718            string mord = ringlist(br1)[3][1][1];
3719            ring R1 = (0,a),(var(1),var(2)),ds;
3720            list L1 = ringlist(R1);
3721            L1[3][1][1]=mord;
3722            def R2 = ring(L1);
3723               setring R2;
3724               kill mord;
3725               kill R1;
3726            poly ff = fetch(R,MP2);
3727            minpoly = number(subst(ff,var(1),a));
3728          }
3729          Poly nfe2 =var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n1);
3730          number m2 = fetch(R,m2);
3731          interval2 = list(0,(m2+1));
3732          NormalFormEquation F2;
3733          F2.normalFormEquation = nfe2;
3734          F2.parameters = list(list(a*var(2)^(n1),interval2));
3735          F2.singularityType = typeofsing2;
3736          setring R;
3737          kill R2;
3738
3739          if(deg(MP3)==1)
3740          {
3741             number a = -number(subst(MP3,var(1),0));
3742              def R2 = br1;
3743              setring R2;
3744              number a = fetch(R,a);
3745          }
3746          else
3747          {
3748            string mord = ringlist(br1)[3][1][1];
3749            ring R1 = (0,a),(var(1),var(2)),ds;
3750            list L1 = ringlist(R1);
3751            L1[3][1][1]=mord;
3752            def R2 = ring(L1);
3753               setring R2;
3754               kill mord;
3755               kill R1;
3756            poly ff = fetch(R,MP3);
3757            minpoly = number(subst(ff,var(1),a));
3758          }
3759          Poly nfe3 =var(1)^2*var(2)^2-var(1)^(n1)+a*var(2)^(n1);
3760          number m3 = fetch(R,m3);
3761          interval3 = list(-(m3+1),0);
3762          NormalFormEquation F3;
3763          F3.normalFormEquation = nfe3;
3764          F3.parameters = list(list(a*var(2)^(n1),interval3));
3765          F3.singularityType = typeofsing3;
3766
3767          setring R;
3768          kill R2;
3769
3770          if(deg(MP4)==1)
3771          {
3772              number a = -number(subst(MP2,var(1),0));
3773              def R2 = br1;
3774              setring R2;
3775              number a = fetch(R,a);
3776          }
3777          else
3778          {
3779            string mord = ringlist(br1)[3][1][1];
3780            ring R1 = (0,a),(var(1),var(2)),ds;
3781            list L1 = ringlist(R1);
3782            L1[3][1][1]=mord;
3783            def R2 = ring(L1);
3784               setring R2;
3785               kill mord;
3786               kill R1;
3787            poly ff = fetch(R,MP4);
3788            minpoly = number(subst(ff,var(1),a));
3789          }
3790          Poly nfe4 = var(1)^2*var(2)^2-var(1)^(n1)+a*var(2)^(n1);
3791          number m4 = fetch(R,m4);
3792          interval4 = list(0,(m4+1));
3793          NormalFormEquation F4;
3794          F4.normalFormEquation = nfe4;
3795          F4.parameters = list(list(a*var(2)^(n1),interval4));
3796          F4.singularityType = typeofsing4;
3797          setring br;
3798          return(F1,F2,F3,F4);
3799        }
3800        else  //r even
3801        {
3802          //"r even";
3803          if(t==-1) //de<0
3804          {
3805            //"de<0";
3806            number s = number(n1);
3807            s = -s/2;
3808            int ns = int(s);
3809            s = coeff(rf,var(1)^2*var(2)^2)^ns;
3810            kill ns;
3811            poly P = var(1)+s*(coeff(rf,var(1)^n1)*coeff(rf,var(2)^n1));
3812            setring R;
3813            MP1 = fetch(br,P);
3814            MP2 = MP1;
3815            number m = maxabs(MP1);
3816            typeofsing1 = "Y["+string(n1)+","+string(n1)+"]+-";
3817            typeofsing2 = "Y["+string(n1)+","+string(n1)+"]++";
3818
3819             if(deg(MP1)==1)
3820             {
3821                number a = -number(subst(MP1,var(1),0));
3822                def R2 = br1;
3823                setring R2;
3824                number a = fetch(R,a);
3825             }
3826             else
3827             {
3828               string mord = ringlist(br1)[3][1][1];
3829               ring R1 = (0,a),(var(1),var(2)),ds;
3830               list L1 = ringlist(R1);
3831               L1[3][1][1]=mord;
3832               def R2 = ring(L1);
3833               setring R2;
3834               kill mord;
3835               kill R1;
3836               poly ff = fetch(R,MP1);
3837               minpoly = number(subst(ff,var(1),a));
3838             }
3839             Poly nfe1 =var(1)^2*var(2)^2-var(1)^(n1)+a*var(2)^(n1);
3840             number m = fetch(R,m);
3841             interval1 = list(-(m+1),m+1);
3842             NormalFormEquation F1;
3843             F1.normalFormEquation = nfe1;
3844             F1.parameters = list(list(a*var(2)^(n1),interval1));
3845             F1.singularityType = typeofsing1;
3846
3847             setring R;
3848             kill R2;
3849
3850             if(deg(MP2)==1)
3851             {
3852                 number a = -number(subst(MP2,var(1),0));
3853                 def R2 = br1;
3854                 setring R2;
3855                 number a = fetch(R,a);
3856             }
3857             else
3858             {
3859               string mord = ringlist(br1)[3][1][1];
3860               ring R1 = (0,a),(var(1),var(2)),ds;
3861               list L1 = ringlist(R1);
3862               L1[3][1][1]=mord;
3863               def R2 = ring(L1);
3864               setring R2;
3865               kill mord;
3866               kill R1;
3867               poly ff = fetch(R,MP2);
3868               minpoly = number(subst(ff,var(1),a));
3869             }
3870             Poly nfe2 =var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n1);
3871             number m = fetch(R,m);
3872             interval2 = list(-(m+1),(m+1));
3873             NormalFormEquation F2;
3874             F2.normalFormEquation = nfe2;
3875             F2.parameters = list(list(a*var(2)^(n1),interval2));
3876             F2.singularityType = typeofsing2;
3877             setring R;
3878
3879            return(F1,F2);
3880           }
3881           else  //de>0
3882           {
3883             //"de>0";
3884             number s = number(n1);
3885             s = -s/2;
3886             int ns = int(s);
3887             s = coeff(rf,var(1)^2*var(2)^2)^ns;
3888             poly P1 = var(1)+s*(coeff(rf,var(1)^n1)*coeff(rf,var(2)^n1));
3889             poly P2 = var(1)-s*(coeff(rf,var(1)^n1)*coeff(rf,var(2)^n1));
3890             int signE = determineSignExpression(coeff(rf,var(1)^n1));
3891             //"e="+string(coeff(rf,var(1)^n1));
3892             //setring R;
3893             //poly f0x = subst(f0,var(2),0)+1;
3894             //poly f0y = subst(f0,var(1),0)+1;
3895             //ring S0 = 0,(var(1),var(2)),dp;
3896             //poly f0x = fetch(R,f0x);
3897             //poly f0y = fetch(R,f0y);
3898             //int n3 = nrroots(f0x); // change
3899             //int n4 = nrroots(f0y);
3900             setring R;
3901             //kill S0;
3902             if(signE<0)
3903             {
3904               typeofsing1 = "Y["+string(n1)+","+string(n1)+"]+-";
3905               nf1 = var(1)^2*var(2)^2-var(1)^(n1)+var(2)^(n1);
3906               MP1 = fetch(br,P1);
3907               number m = maxabs(MP1);
3908               if(deg(MP1)==1)
3909               {
3910                  number a = -number(subst(MP1,var(1),0));
3911                  def R2 = br1;
3912                  setring R2;
3913                  number a = fetch(R,a);
3914               }
3915               else
3916               {
3917                 string mord = ringlist(br1)[3][1][1];
3918                 ring R1 = (0,a),(var(1),var(2)),ds;
3919                 list L1 = ringlist(R1);
3920                 L1[3][1][1]=mord;
3921                 def R2 = ring(L1);
3922               setring R2;
3923               kill mord;
3924               kill R1;
3925                 poly ff = fetch(R,MP1);
3926                 minpoly = number(subst(ff,var(1),a));
3927               }
3928               Poly nfe1 = var(1)^2*var(2)^2-var(1)^(n1)+a*var(2)^(n1);
3929               number m = fetch(R,m);
3930               interval1 = list(-(m+1),m+1);
3931               NormalFormEquation F1;
3932               F1.normalFormEquation = nfe1;
3933               F1.parameters = list(list(a*var(2)^(n1),interval1));
3934               F1.singularityType = typeofsing1;
3935
3936               setring R;
3937               return(F1);
3938             }
3939             else
3940             {
3941               typeofsing1 = "Y["+string(n1)+","+string(n1)+"]++";
3942               MP1 = fetch(br,P2);
3943               number m = maxabs(MP1);
3944               if(deg(MP1)==1)
3945               {
3946                  number a = -number(subst(MP1,var(1),0));
3947                  def R2 = br1;
3948                  setring R2;
3949                  number a = fetch(R,a);
3950               }
3951               else
3952               {
3953                 string mord = ringlist(br1)[3][1][1];
3954                 ring R1 = (0,a),(var(1),var(2)),ds;
3955                 list L1 = ringlist(R1);
3956                 L1[3][1][1]=mord;
3957                 def R2 = ring(L1);
3958               setring R2;
3959               kill mord;
3960               kill R1;
3961                 poly ff = fetch(R,MP1);
3962                 minpoly = number(subst(ff,var(1),a));
3963               }
3964               Poly nfe1 = var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n1);
3965               number m = fetch(R,m);
3966               interval1 = list(-(m+1),m+1);
3967               NormalFormEquation F1;
3968               F1.normalFormEquation = nfe1;
3969               F1.parameters = list(list(a*var(2)^(n1),interval1));
3970               F1.singularityType = typeofsing1;
3971
3972               setring R;
3973               return(F1);
3974             }
3975           }
3976         }}
3977         else   //sigma<1
3978         {
3979            poly de = (coeff(rf,var(1)^n1)*coeff(rf,var(2)^n1));
3980            poly w0 = coeff(rf,var(1)^2*var(2)^2);
3981            ring S = 0,(x,y),ds;
3982            poly w0 = fetch(br,w0);
3983            poly de = fetch(br,de);
3984            number t0 = number(de);
3985            int w = sign(number(w0));
3986            if(n1 mod 2 <>0)
3987           {
3988             int r = 1;
3989           }
3990           else
3991           {
3992             int r = 0;
3993           }
3994           if(t0>0)
3995           {
3996             int t = 1;
3997           }
3998           else
3999           {
4000             int t = -1;
4001           }
4002           setring br;
4003           if(r==1)   //r odd
4004           {
4005             number c = w*(coeff(rf,var(1)^2*var(2)^2));
4006             number s = c^(-n1);
4007             poly P = var(1)^2-s*(coeff(rf,var(1)^(n1))*coeff(rf,var(2)^(n1)))^2;
4008             setring R;
4009             poly P = fetch(br,P);
4010             typeofsing1 = "Y["+string(n1)+","+string(n1)+"]-+";
4011             typeofsing2 = "Y["+string(n1)+","+string(n1)+"]-+";
4012             typeofsing3 = "Y["+string(n1)+","+string(n1)+"]--";
4013             typeofsing4 = "Y["+string(n1)+","+string(n1)+"]--";
4014             number m = maxabs(P);
4015             poly MP1 = RealRootFactor(P,-m-1,0);
4016             number m1 = maxabs(MP1);
4017             poly MP2 = RealRootFactor(P,0,m+1);
4018             number m2 = maxabs(MP2);
4019             poly MP3 = RealRootFactor(P,-m-1,0);
4020             number m3 = maxabs(MP3);
4021             poly MP4 = RealRootFactor(P,0,m+1);
4022             number m4 = maxabs(MP4);
4023
4024            if(deg(MP1)==1)
4025             {
4026                number a = -number(subst(MP1,var(1),0));
4027                def R2 = br1;
4028                setring R2;
4029                number a = fetch(R,a);
4030             }
4031             else
4032             {
4033               string mord = ringlist(br1)[3][1][1];
4034               ring R1 = (0,a),(var(1),var(2)),ds;
4035               list L1 = ringlist(R1);
4036               L1[3][1][1]=mord;
4037               def R2 = ring(L1);
4038               setring R2;
4039               kill mord;
4040               kill R1;
4041               poly ff = fetch(R,MP1);
4042               minpoly = number(subst(ff,var(1),a));
4043             }
4044             Poly nfe1 =-var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n1);
4045             number m1 = fetch(R,m1);
4046             interval1 = list(-(m1+1),0);
4047             NormalFormEquation F1;
4048             F1.normalFormEquation = nfe1;
4049             F1.parameters = list(list(a*var(2)^(n1),interval1));
4050             F1.singularityType = typeofsing1;
4051
4052             setring R;
4053             kill R2;
4054
4055             if(deg(MP2)==1)
4056             {
4057                 number a = -number(subst(MP2,var(1),0));
4058                 def R2 = br1;
4059                 setring R2;
4060                 number a = fetch(R,a);
4061             }
4062             else
4063             {
4064               string mord = ringlist(br1)[3][1][1];
4065               ring R1 = (0,a),(var(1),var(2)),ds;
4066               list L1 = ringlist(R1);
4067               L1[3][1][1]=mord;
4068               def R2 = ring(L1);
4069               setring R2;
4070               kill mord;
4071               kill R1;
4072               poly ff = fetch(R,MP2);
4073               minpoly = number(subst(ff,var(1),a));
4074             }
4075             Poly nfe2 =-var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n1);
4076             number m2 = fetch(R,m2);
4077             interval2 = list(0,(m2+1));
4078             NormalFormEquation F2;
4079             F2.normalFormEquation = nfe2;
4080             F2.parameters = list(list(a*var(2)^(n1),interval2));
4081             F2.singularityType = typeofsing2;
4082             setring R;
4083             kill R2;
4084
4085
4086             if(deg(MP3)==1)
4087             {
4088                number a = -number(subst(MP3,var(1),0));
4089                def R2 = br1;
4090                setring R2;
4091                number a = fetch(R,a);
4092             }
4093             else
4094             {
4095               string mord = ringlist(br1)[3][1][1];
4096               ring R1 = (0,a),(var(1),var(2)),ds;
4097               list L1 = ringlist(R1);
4098               L1[3][1][1]=mord;
4099               def R2 = ring(L1);
4100               setring R2;
4101               kill mord;
4102               kill R1;
4103               poly ff = fetch(R,MP3);
4104               minpoly = number(subst(ff,var(1),a));
4105             }
4106             Poly nfe3 =-var(1)^2*var(2)^2-var(1)^(n1)+a*var(2)^(n1);
4107             number m3 = fetch(R,m3);
4108             interval3 = list(-(m3+1),0);
4109             NormalFormEquation F3;
4110             F3.normalFormEquation = nfe3;
4111             F3.parameters = list(list(a*var(2)^(n1),interval3));
4112             F3.singularityType = typeofsing3;
4113
4114             setring R;
4115             kill R2;
4116
4117             if(deg(MP4)==1)
4118             {
4119                 number a = -number(subst(MP2,var(1),0));
4120                 def R2 = br1;
4121                setring R2;
4122                number a = fetch(R,a);
4123             }
4124             else
4125             {
4126               string mord = ringlist(br1)[3][1][1];
4127               ring R1 = (0,a),(var(1),var(2)),ds;
4128               list L1 = ringlist(R1);
4129               L1[3][1][1]=mord;
4130               def R2 = ring(L1);
4131               setring R2;
4132               kill mord;
4133               kill R1;
4134               poly ff = fetch(R,MP4);
4135               minpoly = number(subst(ff,var(1),a));
4136             }
4137             Poly nfe4 = -var(1)^2*var(2)^2-var(1)^(n1)+a*var(2)^(n1);
4138             number m4 = fetch(R,m4);
4139             interval4 = list(0,(m4+1));
4140             NormalFormEquation F4;
4141             F4.normalFormEquation = nfe4;
4142             F4.parameters = list(list(a*var(2)^(n1),interval4));
4143             F4.singularityType = typeofsing4;
4144             setring R;
4145            return(F1,F2,F3,F4);
4146
4147         }
4148         else    //r even
4149         {
4150           if(t==-1)  //de<0
4151           {
4152             number s = number(n1);
4153             s = -s/2;
4154             int n = int(s);
4155             number s = coeff(rf,var(1)^2*var(2)^2)^n;
4156             poly P = var(1)+s*(coeff(rf,var(1)^n1)*coeff(rf,var(2)^n1));
4157             setring R;
4158             MP1 = fetch(br,P);
4159             MP2 = MP1;
4160             number m1 = maxabs(MP1);
4161             number m2 = maxabs(MP2);
4162             typeofsing1 = "Y["+string(n1)+","+string(n1)+"]--";
4163             typeofsing2 = "Y["+string(n1)+","+string(n1)+"]-+";
4164
4165             if(deg(MP1)==1)
4166             {
4167                number a = -number(subst(MP1,var(1),0));
4168                def R2 = br1;
4169                setring R2;
4170                number a = fetch(R,a);
4171             }
4172             else
4173             {
4174               string mord = ringlist(br1)[3][1][1];
4175               ring R1 = (0,a),(var(1),var(2)),ds;
4176               list L1 = ringlist(R1);
4177               L1[3][1][1]=mord;
4178               def R2 = ring(L1);
4179               setring R2;
4180               kill mord;
4181               kill R1;
4182               poly ff = fetch(R,MP1);
4183               minpoly = number(subst(ff,var(1),a));
4184             }
4185             Poly nfe1 =-var(1)^2*var(2)^2-var(1)^(n1)+a*var(2)^(n1);
4186             number m1 = fetch(R,m1);
4187             interval1 = list(-(m1+1),m1+1);
4188             NormalFormEquation F1;
4189             F1.normalFormEquation = nfe1;
4190             F1.parameters = list(list(a*var(2)^(n1),interval1));
4191             F1.singularityType = typeofsing1;
4192
4193             setring R;
4194             kill R2;
4195
4196             if(deg(MP2)==1)
4197             {
4198                 number a = -number(subst(MP2,var(1),0));
4199                 def R2 = br1;
4200                setring R2;
4201                number a = fetch(R,a);
4202             }
4203             else
4204             {
4205               string mord = ringlist(br1)[3][1][1];
4206               ring R1 = (0,a),(var(1),var(2)),ds;
4207               list L1 = ringlist(R1);
4208               L1[3][1][1]=mord;
4209               def R2 = ring(L1);
4210               setring R2;
4211               kill mord;
4212               kill R1;
4213               poly ff = fetch(R,MP2);
4214               minpoly = number(subst(ff,var(1),a));
4215             }
4216             Poly nfe2 =-var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n1);
4217             number m2 = fetch(R,m2);
4218             interval2 = list(-(m2+1),(m2+1));
4219             NormalFormEquation F2;
4220             F2.normalFormEquation = nfe2;
4221             F2.parameters = list(list(a*var(2)^(n1),interval2));
4222             F2.singularityType = typeofsing2;
4223             setring R;
4224
4225             return(F1,F2);
4226           }
4227           else      //de>0
4228           {
4229             number s = number(n1);
4230             s = -s/2;
4231             int n = int(s);
4232             number s = coeff(rf,var(1)^2*var(2)^2)^n;
4233             setring R;
4234             poly s0 = fetch(br,s);
4235             if(s0>0)
4236             {
4237               int k = 1;
4238             }
4239             else
4240             {
4241               int k = -1;
4242             }
4243             setring br;
4244             poly P1 = var(1)-k*s*(coeff(rf,var(1)^n1)*coeff(rf,var(2)^n2));
4245             poly P2 = var(1)+k*s*(coeff(rf,var(1)^n1)*coeff(rf,var(2)^n2));
4246             //"e="+string(coeff(rf,var(1)^n1));
4247             int signE = determineSignExpression(coeff(rf,var(1)^n1));
4248             //setring R;
4249             //poly f0x = subst(f0,var(2),0)-1;
4250             //poly f0y = subst(f0,var(1),0)-1;
4251             //ring S = 0,(var(1),var(2)),dp;
4252             //poly f0x = fetch(R,f0x);
4253             //poly f0y = fetch(R,f0y);
4254             //int n3 = nrroots(f0x); //change
4255             //int n4 = nrroots(f0y);
4256             setring R;
4257             if(signE>0)
4258             {
4259               typeofsing1 = "Y["+string(n1)+","+string(n1)+"]-+";
4260               nf1 = -var(1)^2*var(2)^2+var(1)^(n1)+var(2)^(n1);
4261               MP1 = fetch(br,P1);
4262               number m1 = maxabs(MP1);
4263
4264               if(deg(MP1)==1)
4265               {
4266                  number a = -number(subst(MP1,var(1),0));
4267                  def R2 = br1;
4268                  setring R2;
4269                  number a = fetch(R,a);
4270               }
4271               else
4272               {
4273                 string mord = ringlist(br1)[3][1][1];
4274                 ring R1 = (0,a),(var(1),var(2)),ds;
4275                 list L1 = ringlist(R1);
4276                 L1[3][1][1]=mord;
4277                 def R2 = ring(L1);
4278               setring R2;
4279               kill mord;
4280               kill R1;
4281                 poly ff = fetch(R,MP1);
4282                 minpoly = number(subst(ff,var(1),a));
4283               }
4284               Poly nfe1 =-var(1)^2*var(2)^2+var(1)^(n1)+a*var(2)^(n1);
4285               number m1 = fetch(R,m1);
4286               interval1 = list(-(m1+1),m1+1);
4287               NormalFormEquation F1;
4288               F1.normalFormEquation = nfe1;
4289               F1.parameters = list(list(a*var(2)^(n1),interval1));
4290               F1.singularityType = typeofsing1;
4291
4292               setring R;
4293               return(F1);
4294             }
4295             else
4296             {
4297               typeofsing1 = "Y["+string(n1)+","+string(n1)+"]--";
4298               MP1 = fetch(br,P2);
4299               maxabs(MP1);
4300               if(deg(MP1)==1)
4301               {
4302                  number a = -number(subst(MP1,var(1),0));
4303                  def R2 = br1;
4304                  setring R2;
4305                  number a = fetch(R,a);
4306               }
4307               else
4308               {
4309                 string mord = ringlist(br1)[3][1][1];
4310                 ring R1 = (0,a),(var(1),var(2)),ds;
4311                 list L1 = ringlist(R1);
4312                 L1[3][1][1]=mord;
4313                 def R2 = ring(L1);
4314               setring R2;
4315               kill mord;
4316               kill R1;
4317                 poly ff = fetch(R,MP1);
4318                 minpoly = number(subst(ff,var(1),a));
4319               }
4320               Poly nfe1 =-var(1)^2*var(2)^2-var(1)^(n1)+a*var(2)^(n1);
4321               number m1 = fetch(R,m1);
4322               interval1 = list(-(m1+1),m1+1);
4323               NormalFormEquation F1;
4324               F1.normalFormEquation = nfe1;
4325               F1.parameters = list(list(a*var(2)^(n1),interval1));
4326               F1.singularityType = typeofsing1;
4327
4328               setring R;
4329               return(F1);
4330           }
4331         }
4332       }
4333     }
4334    }
4335  }
4336  else         //case \tilde Y[r]--h has no linear factors over R
4337  {
4338    setring br;
4339      number b = (coeff(rf,var(1)^2*var(2)^2));
4340      b = b^(-1);
4341      number c1 = coeff(rf,var(1)^(n1));
4342      number c2 = coeff(rf,var(2)^(n1));
4343      poly P = var(1)^8-(c1*c2)^4*(4^2*b)^(2*n1);
4344      kill b;
4345      setring R;
4346      poly P = fetch(br,P);
4347      if(sigma>0)
4348      {
4349        typeofsing1 = "tilde(Y)["+string(n1)+"]+";
4350        typeofsing2 = "tilde(Y)["+string(n1)+"]+";
4351      }
4352      else
4353      {
4354        typeofsing1 = "tilde(Y)["+string(n1)+"]-";
4355        typeofsing2 = "tilde(Y)["+string(n1)+"]-";
4356      }
4357      number m = maxabs(P);
4358      MP1 = RealRootFactor(P,0,m+1);
4359      number m1 = maxabs(MP1);
4360      MP2 = RealRootFactor(P,-m-1,0);
4361      number m2 = maxabs(MP2);
4362      if(deg(MP1)==1)
4363      {
4364         number a = -number(subst(MP1,var(1),0));
4365         def R2 = br1;
4366         setring R2;
4367         number a = fetch(R,a);
4368      }
4369      else
4370      {
4371        string mord = ringlist(br1)[3][1][1];
4372        ring R1 = (0,a),(var(1),var(2)),ds;
4373        list L1 = ringlist(R1);
4374        L1[3][1][1]=mord;
4375        def R2 = ring(L1);
4376               setring R2;
4377               kill mord;
4378               kill R1;
4379        poly ff = fetch(R,MP1);
4380        minpoly = number(subst(ff,var(1),a));
4381      }
4382      Poly nfe1 =sigma*(var(1)^2+var(2)^2)^2+a*var(1)^(n1);
4383      number m1 = fetch(R,m1);
4384      interval1 = list(0,m1+1);
4385      NormalFormEquation F1;
4386      F1.normalFormEquation = nfe1;
4387      F1.parameters = list(list(a*var(1)^(n1),interval1));
4388      F1.singularityType = typeofsing1;
4389
4390      setring R;
4391      kill R2;
4392      if(deg(MP2)==1)
4393      {
4394         number a = -number(subst(MP2,var(1),0));
4395         def R2 = br1;
4396         setring R2;
4397         number a = fetch(R,a);
4398      }
4399      else
4400      {
4401        string mord = ringlist(br1)[3][1][1];
4402        ring R1 = (0,a),(var(1),var(2)),ds;
4403        list L1 = ringlist(R1);
4404        L1[3][1][1]=mord;
4405        def R2 = ring(L1);
4406               setring R2;
4407               kill mord;
4408               kill R1;
4409        poly ff = fetch(R,MP2);
4410        minpoly = number(subst(ff,var(1),a));
4411      }
4412      Poly nfe2 =sigma*(var(1)^2+var(2)^2)^2+a*var(1)^(n1);
4413      number m2 = fetch(R,m2);
4414      interval2 = list(-(m2+1),0);
4415      NormalFormEquation F2;
4416      F2.normalFormEquation = nfe2;
4417      F2.parameters = list(list(a*var(1)^(n1),interval2));
4418      F2.singularityType = typeofsing2;
4419      setring R;
4420      return(F1,F2);
4421    }
4422  }
4423}
4424
4425static proc determineSignExpression(number e){
4426def R = basering;
4427poly MP = subst(minpoly,a,var(1));
4428poly E = subst(e,a,var(1));
4429ring RS =0,var(1),dp;
4430poly MP = fetch(R,MP);
4431poly E = fetch(R,E);
4432number ze = -coeff(MP,var(1))/(2*coeff(MP,var(1)^2));
4433return(sign(number(subst(E,var(1),ze))));
4434//if (coeff(E,var(1))==0){return(sign(coeff(E,1)));}
4435//number zeroOfE=-coeff(E,1)/coeff(E,var(1));
4436//int signE=sign(zeroOfE*sign(coeff(E,var(1))));
4437//return(signE);
4438}
4439
4440
4441//////////////////////////////////////////////////////////////////////////////////////////////
4442static proc caseJ10k(poly rf, int mu, br1)
4443{
4444/*preliminaries*/
4445  string typeofsing;
4446  def br = basering;
4447  map phi;
4448  poly MP;
4449  list intervalP;
4450  list interval1;
4451  list interval2;
4452  string typeofsing1;
4453  string typeofsing2;
4454
4455  /* II. Eliminate the monomials in supp(rf) underneath or on \Gamma(T)*/
4456  /* Apply Algorithm 1 to rf */
4457  poly g = jet(rf,3);
4458  list factorization = factorize(g);
4459  poly m = factorization[1][2];
4460  if(leadmonom(m)==y)
4461  {
4462    phi = br,var(2),var(1);
4463    rf = phi(rf);
4464  }
4465  else
4466  {
4467    matrix T[2][2] = factorization[1][2]/var(1),factorization[1][2]/var(2) ,
4468                     0, 1;
4469    phi = br, luinverse(T)[2]*matrix(ideal(var(1), var(2)), 2, 1);
4470    rf = phi(rf);
4471  }
4472
4473  /* Apply Algorithm 3 */
4474  intvec v = 2,1;
4475  poly p = WeightedJet(rf,v,6);
4476  p = subst(p,var(2),1);
4477  kill factorization;
4478  list factorization = factorize(p);
4479
4480  for(int i =2; i<=size(factorization[2]);i++)
4481  {
4482    if(factorization[2][i] == 2)
4483    {
4484      poly s = factorization[1][i];
4485    }
4486  }
4487
4488  number r = coeff(s,var(1));
4489  number q = number(-subst(s,var(1),0))/r;
4490  phi = br, var(1)+q*var(2)^2, var(2);
4491  rf = phi(rf);
4492
4493  /*Use Algorithm 4 iteratively*
4494  number c = coeff(rf,var(1)^2*var(2)^2);
4495  number n = (mu-2)*(1/2);
4496  if(round(n)>n)
4497  {
4498    n = round(n)-1;
4499  }
4500  else
4501  {
4502    n = round(n);
4503  }
4504  for(int i=4; i<=n;i++)
4505  {
4506    rf = jet(rf,mu+1);
4507    number t = coeff(rf,var(1)*var(2)^(i));
4508    phi = br, var(1)-(t/(2*c))*var(2)^(i-2),var(2);
4509    rf = phi(rf);
4510  }
4511
4512  /* III. and IV. Discard the monomials above \Gamma(T) and read off the desired information*/
4513
4514  number b = coeff(rf,var(1)^3);
4515  number c = coeff(rf,var(1)^2*var(2)^2);
4516  number t = coeff(rf,var(2)^(mu-4));
4517  MP = var(1)^6-(b^2/(absValue(c))^3)^(mu-4)*t^6;
4518  number n = maxabs(MP);
4519  int sc = sign(c);
4520  if(t>0)
4521  {
4522     MP = RealRootFactor(MP,0,n+1);
4523     n = maxabs(MP);
4524     number n1 = 0;
4525     number n2 = n+1;
4526     number n3 = -(n+1);
4527     number n4 = 0;
4528  }
4529  else
4530  {
4531    MP = RealRootFactor(MP,0,n+1);
4532    n = maxabs(MP);
4533    number n1 = -(n+1);
4534    number n2 = 0;
4535    number n3 = 0;
4536    number n4 = (n+1);
4537
4538  }
4539  if(c>0)
4540  {
4541    typeofsing1 = "J["+string(mu)+"]+";
4542  }
4543  else
4544  {
4545    typeofsing1 = "J["+string(mu)+"]-";
4546  }
4547
4548  if(deg(MP)==1)
4549  {
4550    number a = -number(subst(MP,var(1),0));
4551    def R2 =  br1;
4552    setring R2;
4553    number a = fetch(br,a);
4554
4555  }
4556  else
4557  {
4558    string mord = ringlist(br1)[3][1][1];
4559    ring R1 = (0,a),(var(1),var(2)),ds;
4560    list L1 = ringlist(R1);
4561    L1[3][1][1]=mord;
4562    def R2 = ring(L1);
4563    setring R2;
4564    poly ff = fetch(br,MP);
4565    minpoly = number(subst(ff,var(1),a));
4566  }
4567  number n1 = fetch(br,n1);
4568  number n2 = fetch(br,n2);
4569  Poly nfe1 =var(1)^3+number(sc)*var(1)^2*var(2)^2+a*var(2)^(mu-4);
4570  interval1 = list(n1,n2);
4571  NormalFormEquation F1;
4572  F1.normalFormEquation = nfe1;
4573  F1.parameters = list(list(a*var(2)^(mu-4),interval1));
4574  F1.singularityType = typeofsing1;
4575  number b=-a;
4576  number n3 = fetch(br,n3);
4577  number n4 = fetch(br,n4);
4578  Poly nfe2 =var(1)^3+number(sc)*var(1)^2*var(2)^2+b*var(2)^(mu-4);
4579  interval2 = list(n3,n4);
4580  NormalFormEquation F2;
4581  F2.normalFormEquation = nfe2;
4582  F2.parameters = list(list(b*var(2)^(mu-4),interval2));
4583  F2.singularityType = typeofsing1;
4584  setring br;
4585  if((mu-4) mod 2 == 0)
4586  {
4587    return(list(F1));
4588  }
4589  else
4590  {
4591    return(list(F1,F2));
4592  }
4593}
4594
4595static proc findX9solutions(poly f, number s1, number s2, int expnum){
4596matrix P[4][4] = 1,2,0,0, 0,1,3,0, 0,0,1,5, 1,0,0,-3;
4597list sol;
4598int cfs=0;
4599while (size(sol)<>expnum){
4600  if (rank(P)==4){
4601    sol=findX9solutionsP(f, s1, s2,P);
4602    if (size(sol)<expnum){ERROR("case X9, wrong number of equivalent normal form equations");}
4603    cfs++;
4604  }
4605  P=randommat(4,4,maxideal(0),cfs);
4606}
4607return(sol);}
4608
4609// make static
4610static proc findX9solutionsP(poly f, number s1, number s2,matrix pj){
4611//"entering findX9solutions input "+string(f)+", "+string(s1)+", "+string(s2);
4612def R0 = basering;
4613ring R=0,(@a,@b,@c,@d,var(1),var(2)),lp;
4614poly f = imap(R0,f);
4615map phi = R, @a,@b,@c,@d,@a*var(5)+@b*var(6), @c*var(5)+@d*var(6);
4616poly ft=phi(f);
4617//ft;
4618matrix C=coef(ft,var(5)*var(6));
4619//C;
4620matrix Cc[1][ncols(C)]=C[2,1..ncols(C)];
4621ring R1=0,(@a,@b,@c,@d,t),dp;
4622matrix Cc = imap(R,Cc);
4623number x4coeff = fetch(R0,s1);
4624number y4coeff = fetch(R0,s2);
4625ideal I = Cc[1,1]-x4coeff, Cc[1,2]-0, Cc[1,3]-t, Cc[1,4]-0, Cc[1,5]-y4coeff;
4626short =0;
4627//I=std(radical(I));
4628I=std(I);
4629//"ideal to be solved "+string(I);
4630//"elimination ideal";
4631matrix pj = imap(R0,pj);
4632matrix pjim = pj*matrix([@a,@b,@c,@d]);
4633map phi =basering,pjim[1,1],pjim[2,1],pjim[3,1],pjim[4,1],t;
4634I=phi(I);
4635I =eliminate(I,@a*@b*@c);
4636//elimI;
4637ring R2=0,(@d,t),dp;
4638ideal I = imap(R1,I);
4639//I;
4640//"starting root isolation";
4641//list result = rootIsolation(elimI);
4642//result;
4643//setring R1;
4644list result = rootIsolation(I);
4645//"finished";
4646//list result=0,result2;
4647setring R1;
4648poly mp = eliminate(I,@a*@b*@c*@d)[1];
4649mp = subst(mp,t,@a);
4650setring R;
4651poly mp = imap(R1,mp);
4652mp = subst(mp,@a,var(5));
4653setring R0;
4654poly mp = imap(R,mp);
4655poly mpfac;
4656number i1,i2;
4657//setring R2;
4658setring R1;
4659list intervals;
4660int i,j;
4661int tst;
4662//"result of root isolation ";
4663//result;
4664int ia = 2;
4665//int ia = 5;
4666for (i=1;i<=size(result[2]);i++){
4667   tst=0;
4668   for (j=1;j<=size(intervals);j++){
4669     if (intersects(intervals[j],result[2][i][ia])){
4670        //intersectIntervals(intervals[j],result[2][i][ia]);
4671        intervals[j]=intersectIntervals(intervals[j],result[2][i][ia]);
4672        //"-making interval j smaller ";intervals[j];
4673        tst=1;
4674     }
4675   }
4676   if (tst==0){
4677        //"adding ";string(result[2][i][ia]);
4678        intervals[size(intervals)+1]=result[2][i][ia];
4679   }
4680}
4681//"number of intervals "+string(size(intervals));
4682//intervals;
4683list ret;
4684number i1,i2;
4685for (i=1;i<=size(intervals);i++){
4686   setring R;
4687   //"---";
4688   //intervals[i];
4689   number i1 = intervals[i][1];
4690   number i2 = intervals[i][2];
4691   setring R0;
4692   i1=fetch(R,i1);
4693   i2=fetch(R,i2);
4694   mpfac = RealRootFactor(mp,i1,i2);
4695   ret[i]=list(mpfac,list(i1,i2));
4696}
4697return(ret);
4698}
4699
4700
4701static proc intersects(interval I,interval J){
4702number i1 = I[1];
4703number i2 = I[2];
4704number j1 = J[1];
4705number j2 = J[2];
4706if (((i1<=j1) and (j1<=i2)) or ((i1<=j2) and (j2<=i2)) or ((j1<=i1) and (i1<=j2)) or ((j1<=i2) and (i2<=j2))){
4707    return(1);
4708}
4709return(0);}
4710
4711static proc intersectIntervals(interval I,interval J){
4712number i1 = I[1];
4713number i2 = I[2];
4714number j1 = J[1];
4715number j2 = J[2];
4716interval L;
4717if ((i1<=j1) and (j1<=i2)){
4718  L = j1,min(i2,j2);
4719}
4720if ((i1<=j2) and (j2<=i2)){
4721  L = max(i1,j1),j2;
4722}
4723if ((j1<=i1) and (i1<=j2)){
4724  L = i1,min(j2,i2);
4725}
4726if ((j1<=i2) and (i2<=j2)){
4727  L = max(j1,i1),i2;
4728}
4729return(L);}
4730
4731
4732static proc caseX9(poly rf, def br1){
4733    string typeofsing;
4734    poly nf;
4735    def br = basering;
4736    poly MP2;
4737
4738    poly f4 = jet(rf, 4);
4739    poly f=f4;
4740    map phi;
4741    if (coeff(f,var(1)^4)==0){
4742      phi = br,var(2),var(1);
4743      f=phi(f4);
4744      if (coeff(f,var(1)^4)==0){
4745         phi = br,var(1),var(1)+var(2);
4746         f=phi(f4);
4747         if (coeff(f,var(1)^4)==0){
4748           phi = br,var(1),2*var(1)+var(2);
4749           f=phi(f4);
4750           if (coeff(f,x^4)==0){
4751             phi = br,var(1),3*var(1)+var(2);
4752             f=phi(f4);
4753           }
4754         }
4755      }
4756    }
4757    phi = br, var(1) - coeff(f,var(1)^3*var(2))/coeff(f,var(1)^4)*var(2),var(2);
4758    f=phi(f);
4759    //"no x^3*y term";
4760    //f;
4761    number b = coeff(f,var(1)^4);
4762    number c = coeff(f,var(1)^2*var(2)^2);
4763    number d = coeff(f,var(1)*var(2)^3);
4764    number e = coeff(f,var(2)^4);
4765
4766    list interval1,interval2;
4767    poly nfe1,nfe2;
4768    poly f1 = subst(f,var(2),1);
4769    ring Rgl = 0,var(1),dp;
4770    poly f1 = fetch(br,f1);
4771    int r = nrroots(f1);
4772    setring br;
4773    poly MP1,MP2;
4774    list sol;
4775    string typeofsing1,typeofsing2;
4776    if (r==0){
4777      if (b>0){
4778        sol = findX9solutions(f,1,1,2);
4779        //if (size(sol)<>2){ERROR("case X9, wrong number of normal form equations");}
4780        MP1 = sol[1][1];
4781        MP2 = sol[2][1];
4782        interval1 = sol[1][2];
4783        interval2 = sol[2][2];
4784        nfe1 = var(1)^4+var(2)^4;
4785        nfe2 = var(1)^4+var(2)^4;
4786        typeofsing1="X[9]++";
4787        typeofsing2=typeofsing1;
4788      } else {
4789        sol = findX9solutions(f,-1,-1,2);
4790        //if (size(sol)<>2){ERROR("case X9, wrong number of normal form equations");}
4791        MP1 = sol[1][1];
4792        MP2 = sol[2][1];
4793        interval1 = sol[1][2];
4794        interval2 = sol[2][2];
4795        nfe1 = -var(1)^4-var(2)^4;
4796        nfe2 = -var(1)^4-var(2)^4;
4797        typeofsing1="X[9]--";
4798        typeofsing2=typeofsing1;
4799      }
4800    }
4801    if (r==2){
4802        sol = findX9solutions(f,1,-1,1);
4803        //if (size(sol)<>1){ERROR("case X9, wrong number of normal form equations");}
4804        MP1 = sol[1][1];
4805        MP2 = MP1;
4806        interval1 = sol[1][2];
4807        interval2 = interval1;
4808        nfe1 = var(1)^4-var(2)^4;
4809        nfe2 = -var(1)^4+var(2)^4;
4810        typeofsing1="X[9]+-";
4811        typeofsing2="X[9]-+";
4812    }
4813    if (r==4){
4814        sol = findX9solutions(f,1,1,1);
4815        //if (size(sol)<>1){ERROR("case X9, wrong number of normal form equations");}
4816        sol = sol + findX9solutions(f,-1,-1,1);
4817        //if (size(sol)<>2){ERROR("case X9, wrong number of normal form equations");}
4818        MP1 = sol[1][1];
4819        MP2 = sol[2][1];
4820        interval1 = sol[1][2];
4821        interval2 = sol[2][2];
4822        nfe1 = var(1)^4+var(2)^4;
4823        nfe2 = -var(1)^4-var(2)^4;
4824        typeofsing1="X[9]++";
4825        typeofsing2="X[9]--";
4826    }
4827    if(deg(MP1)==1)
4828    {
4829      number a = -number(subst(MP1,var(1),0));
4830      def R2 =  br1;
4831      setring R2;
4832      number a = fetch(br,a);
4833    }
4834    else
4835    {
4836      string mord = ringlist(br1)[3][1][1];
4837      ring R1 = (0,a),(var(1),var(2)),ds;
4838      list L1 = ringlist(R1);
4839      L1[3][1][1]=mord;
4840      def R2 = ring(L1);
4841      setring R2;
4842      kill R1;
4843      poly ff = fetch(br,MP1);
4844      minpoly = number(subst(ff,var(1),a));
4845    }
4846
4847list interval1 = fetch(br,interval1);
4848poly nfe1 = fetch(br,nfe1);
4849NormalFormEquation F1;
4850Poly nfe1a = nfe1+a*var(1)^2*var(2)^2;
4851F1.normalFormEquation = nfe1a;
4852F1.parameters = list(list(a*var(1)^2*var(2)^2,interval1));
4853F1.singularityType = typeofsing1;
4854
4855setring br;
4856
4857    if(deg(MP2)==1)
4858    {
4859      number a = -number(subst(MP2,var(1),0));
4860      def R2 =  br1;
4861      setring R2;
4862      number a = fetch(br,a);
4863    }
4864    else
4865    {
4866      string mord = ringlist(br1)[3][1][1];
4867      ring R1 = (0,a),(var(1),var(2)),ds;
4868      list L1 = ringlist(R1);
4869      L1[3][1][1]=mord;
4870      def R2 = ring(L1);
4871      setring R2;
4872      kill R1;
4873      poly ff = fetch(br,MP2);
4874      minpoly = number(subst(ff,var(1),a));
4875    }
4876
4877list interval2 = fetch(br,interval2);
4878poly nfe2 = fetch(br,nfe2);
4879NormalFormEquation F2;
4880Poly nfe2a = nfe2+a*var(1)^2*var(2)^2;
4881F2.normalFormEquation = nfe2a;
4882F2.parameters = list(list(a*var(1)^2*var(2)^2,interval2));
4883F2.singularityType = typeofsing2;
4884
4885setring br1;
4886return(list(F1,F2));
4887}
4888
4889
4890
4891//////////////////////////////////////////////////////////////////////////////////////////
4892static proc caseX9k(poly rf,int mu, br1)
4893{
4894  /*preliminaries*/
4895  string typeofsing;
4896  poly nf;
4897  def br = basering;
4898  number a;
4899  map phi;
4900  poly nf1;
4901  poly nf2;
4902  poly MP;
4903  poly MP1;
4904  poly MP2;
4905  list intervalP;
4906  list interval1;
4907  list interval2;
4908  int aconst;
4909
4910  /* II. Eliminate the monomials in supp(rf) underneath or on \Gamma(T)*/
4911  /* Apply Algorithm 2 to rf */
4912
4913  poly g = jet(rf,4);
4914  list factorization = factorize(g);
4915  int i=2;
4916  while(factorization[2][i]<>2)
4917  {
4918    i = i+1;
4919  }
4920  poly m = factorization[1][i];
4921  if(leadmonom(m)==y)
4922  {
4923    phi = br,var(2),var(1);
4924    rf = phi(rf);
4925  }
4926  else
4927  {
4928    matrix T[2][2] = factorization[1][i]/var(1),factorization[1][i]/var(2) ,
4929                     0, 1;
4930    phi = br, luinverse(T)[2]*matrix(ideal(var(1), var(2)), 2, 1);
4931    rf = phi(rf);
4932  }
4933  phi = br,var(1), var(2)-coeff(rf,var(1)^3*var(2))*(1/2)*coeff(rf,var(1)^2*var(2)^2)^(-1)*var(1);
4934  rf = phi(rf);
4935  kill i;
4936  /*Use Algorithm 4 iteratively*/
4937
4938  number c = coeff(rf,var(1)^2*var(2)^2);
4939  number n = number(mu-3)/2;
4940  if(round(n)>n)
4941  {
4942    n = round(n)-1;
4943  }
4944  else
4945  {
4946    n = round(n);
4947  }
4948  for(int i=4; i<=n;i++)
4949  {
4950    rf = jet(rf,mu+1);
4951    number t = coeff(rf,var(1)*var(2)^(i));
4952    phi = br, var(1)-(t/(2*c))*var(2)^(i-2),var(2);
4953    rf = phi(rf);
4954    kill t;
4955  }
4956  kill i;
4957/* III. and IV. Discard the monomials above \Gamma(T) and read off the desired information*/
4958
4959  int b0 = sign(coeff(rf,var(1)^4));
4960  int b1 = sign(coeff(rf,var(1)^2*var(2)^2));
4961  if(b0<0 && b1<0)
4962  {
4963    typeofsing = "X["+string(mu)+"]--";
4964  }
4965  if(b0<0 && b1>0)
4966  {
4967    typeofsing = "X["+string(mu)+"]-+";
4968  }
4969  if(b0>0 && b1<0)
4970  {
4971    typeofsing = "X["+string(mu)+"]+-";
4972  }
4973  if(b0>0 && b1>0)
4974  {
4975    typeofsing = "X["+string(mu)+"]++";
4976  }
4977  MP = var(1)^4-coeff(rf,var(2)^(mu-5))^4*(absValue(coeff(rf,var(1)^4))*coeff(rf,var(1)^2*var(2)^2)^(-2))^(mu-5);
4978  if(mu mod 2 <>0)
4979  {
4980    n = maxabs(MP);
4981    MP = RealRootFactor(MP,-n-1,n+1);
4982    n = maxabs(MP);
4983    if(deg(MP)==1)
4984    {
4985      number a = -number(subst(MP,var(1),0));
4986      def R2 =  br1;
4987      setring R2;
4988      number a = fetch(br,a);
4989    }
4990    else
4991    {
4992      string mord = ringlist(br1)[3][1][1];
4993      ring R1 = (0,a),(var(1),var(2)),ds;
4994      list L1 = ringlist(R1);
4995      L1[3][1][1]=mord;
4996      def R2 = ring(L1);
4997      setring R2;
4998      kill R1;
4999      poly ff = fetch(br,MP);
5000      minpoly = number(subst(ff,var(1),a));
5001    }
5002    number n = fetch(br,n);
5003    Poly nfe = number(b0)*var(1)^4+number(b1)*var(1)^2*var(2)^2+a*var(2)^(mu-5);
5004    intervalP = list(-(n+1),n+1);
5005    NormalFormEquation F;
5006    F.normalFormEquation = nfe;
5007    F.parameters = list(list(a*var(2)^(mu-5),intervalP));
5008    F.singularityType = typeofsing;
5009    return(list(F));
5010  }
5011  else
5012  {
5013    n = maxabs(MP);
5014    MP1 = RealRootFactor(MP,0,n+1);
5015    MP2 = RealRootFactor(MP,-n-1,0);
5016    if(deg(MP1)==1)
5017    {
5018      a = -number(subst(MP1,var(1),0));
5019      def R2 =  br1;
5020      setring R2;
5021      number a = fetch(br,a);
5022      aconst=1;
5023    }
5024    else
5025    {
5026      n = maxabs(MP1);
5027      string mord = ringlist(br1)[3][1][1];
5028      ring R1 = (0,a),(var(1),var(2)),ds;
5029      list L1 = ringlist(R1);
5030      L1[3][1][1]=mord;
5031      def R2 = ring(L1);
5032      setring R2;
5033      kill R1;
5034      poly ff = fetch(br,MP1);
5035      minpoly = number(subst(ff,var(1),a));
5036    }
5037    number n = fetch(br,n);
5038    Poly nfe1 = number(b0)*var(1)^4+number(b1)*var(1)^2*var(2)^2+a*var(2)^(mu-5);
5039    interval1 = list(0,n+1);
5040    NormalFormEquation F1;
5041    F1.normalFormEquation = nfe1;
5042    F1.parameters = list(list(a*var(2)^(mu-5),interval1));
5043    F1.singularityType = typeofsing;
5044    if (aconst){kill a;}
5045    kill n;
5046    setring br;
5047    kill R2;
5048    aconst=0;
5049
5050    if(deg(MP2)==1)
5051    {
5052      number a = -number(subst(MP2,var(1),0));
5053      def R2 =  br1;
5054      setring R2;
5055      number a = fetch(br,a);
5056      aconst = 1;
5057    }
5058    else
5059    {
5060      n = maxabs(MP2);
5061      string mord = ringlist(br1)[3][1][1];
5062      ring R1 = (0,a),(var(1),var(2)),ds;
5063      list L1 = ringlist(R1);
5064      L1[3][1][1]=mord;
5065      def R2 = ring(L1);
5066      setring R2;
5067      kill R1;
5068      poly ff = fetch(br,MP2);
5069      minpoly = number(subst(ff,var(1),a));
5070    }
5071    number n = fetch(br,n);
5072    Poly nfe2 = number(b0)*var(1)^4+number(b1)*var(1)^2*var(2)^2+a*var(2)^(mu-5);
5073    interval2 = list(-n-1,0);
5074    NormalFormEquation F2;
5075    F2.normalFormEquation = nfe2;
5076    F2.parameters = list(list(a*var(2)^(mu-5),interval2));
5077    F2.singularityType = typeofsing;
5078    if(aconst){kill a;}
5079    return(list(F1,F2));
5080  }
5081}
5082//////////////////////////////////////////////////////////////////////////////////////////////
5083
5084
5085proc addnondegeneratevariables(NormalFormEquation NFR)
5086"
5087USAGE:    addnondegeneratevariables(NFR); NFR NormalFormEquation
5088RETURN:   Adds squares of the non-degenerate variables (i.e. var(cr+1), ...,
5089          var(nvars(basering)) for corank cr) to the normalform nf,
5090          with signs according to the inertia index.
5091SEE ALSO: realmorsesplit
5092KEYWORDS: Morse lemma; Splitting lemma
5093EXAMPLE:  example addnondegeneratevariables; shows an example
5094"
5095{
5096  def curring = basering;
5097  Poly nf1 = NFR.normalFormEquation;
5098  def S = nf1.in;
5099  setring S;
5100  poly nf = nf1.value;
5101  list L = ringlist(S);
5102  L[2]=NFR.vars;
5103  def Rnew=ring(L);
5104  setring Rnew;
5105  poly nf = imap(S,nf);
5106  int n = size(NFR.vars);
5107  int lambda = NFR.inertiaIndex;
5108  int cr = NFR.corank;
5109  int i;
5110  for(i = cr+1; i <= n-lambda; i++)
5111  {
5112    nf = nf+var(i)^2;
5113  }
5114  for(i = n-lambda+1; i <= n ; i++)
5115  {
5116    nf = nf-var(i)^2;
5117  }
5118  Poly nfp=nf;
5119  NFR.normalFormEquation = nfp;
5120  setring curring;
5121  return(NFR);
5122}
5123example
5124{
5125  "EXAMPLE:";
5126  echo = 2;
5127  ring r = 0, (x,y,z), ds;
5128  poly f = (x2+3y-2z)^2+xyz-(x-y3+x2z3)^3;
5129  list NFR = realclassify(f);
5130  NFR[1];
5131  addnondegeneratevariables(NFR[1]);
5132}
5133
5134///////////////////////////////////////////////////////////////////////////////
5135proc realmorsesplit(poly f, list #)
5136"
5137USAGE:    realmorsesplit(f[, mu]); f poly, mu int
5138RETURN:   a list consisting of the corank of f, the inertia index, an upper
5139          bound for the determinacy, and the residual form of f
5140NOTE:     The characteristic of the basering must be zero, the monomial order
5141          must be local, f must be contained in maxideal(2) and the Milnor
5142          number of f must be finite.
5143          @* The Milnor number of f can be provided as an optional parameter in
5144             order to avoid that it is computed again.
5145SEE ALSO: morsesplit
5146KEYWORDS: Morse lemma; Splitting lemma
5147EXAMPLE:  example morsesplit; shows an example"
5148{
5149  int i;
5150
5151  /* error check */
5152  if(char(basering) != 0)
5153  {
5154    ERROR("The characteristic must be zero.");
5155  }
5156  int n = nvars(basering);
5157  for(i = 1; i <= n; i++)
5158  {
5159    if(var(i) > 1)
5160    {
5161      ERROR("The monomial order must be local.");
5162    }
5163  }
5164  if(jet(f, 1) != 0)
5165  {
5166    ERROR("The input polynomial must be contained in maxideal(2).");
5167  }
5168
5169  /* get Milnor number before continuing error check */
5170  int mu;
5171  if(size(#) > 0)   // read optional parameter
5172  {
5173    if(size(#) > 1 || typeof(#[1]) != "int")
5174    {
5175      ERROR("Wrong optional parameters.");
5176    }
5177    else
5178    {
5179      mu = #[1];
5180    }
5181  }
5182  else              // compute Milnor number
5183  {
5184    mu = milnornumber(f);
5185  }
5186
5187  /* continue error check */
5188  if(mu < 0)
5189  {
5190    ERROR("The Milnor number of the input polynomial must be"+newline
5191      +"non-negative and finite.");
5192  }
5193
5194  /* compute the determinacy */
5195  int k = determinacy(f, mu);
5196  f = jet(f, k);
5197
5198  /* get jet(f, 2) right */
5199  matrix H = concat(jet(jacob(jacob(f)), 0)/2, unitmat(n));
5200  H = sym_reduce(H);
5201  intvec perm_zero;
5202  intvec perm_neg;
5203  intvec perm_pos;
5204  int c;
5205  int lambda;
5206  for(i = 1; i <= n; i++)
5207  {
5208    if(H[i, i] == 0)
5209    {
5210      perm_zero = perm_zero, i;
5211      c++;
5212    }
5213    if(H[i, i] < 0)
5214    {
5215      perm_neg = perm_neg, i;
5216      lambda++;
5217    }
5218    if(H[i, i] > 0)
5219    {
5220      perm_pos = perm_pos, i;
5221    }
5222  }
5223  intvec perm;
5224  if(size(perm_zero) > 1)
5225  {
5226    perm = perm, perm_zero[2..size(perm_zero)];
5227  }
5228  if(size(perm_neg) > 1)
5229  {
5230    perm = perm, perm_neg[2..size(perm_neg)];
5231  }
5232  if(size(perm_pos) > 1)
5233  {
5234    perm = perm, perm_pos[2..size(perm_pos)];
5235  }
5236  perm = perm[2..size(perm)];
5237  matrix T[n][n];
5238  matrix D[1][n];
5239  for(i = 1; i <= n; i++)
5240  {
5241    T[1..n, i] = H[perm[i], (n+1)..(2*n)];
5242    D[1, i] = H[perm[i], perm[i]];
5243  }
5244  map phi = basering, matrix(maxideal(1))*transpose(T);
5245  f = phi(f);
5246  f = jet(f, k);
5247
5248  /* separate the variables */
5249  phi = basering, maxideal(1);
5250  map corank_part = basering, maxideal(1);
5251  for(i = c+1; i <= n; i++)
5252  {
5253    corank_part[i] = 0;
5254  }
5255  poly h = f-jet(f, 2)-corank_part(f);
5256  poly hi;
5257  while(h != 0)
5258  {
5259    for(i = c+1; i <= n; i++)
5260    {
5261      hi = h/var(i);
5262      phi[i] = var(i)-hi/(2*D[1, i]);
5263      h = h-hi*var(i);
5264    }
5265    f = phi(f);
5266    f = jet(f, k);
5267    h = f-jet(f, 2)-corank_part(f);
5268  }
5269  poly g = f-jet(f, 2);
5270
5271  return(list(c, lambda, k, g));
5272}
5273example
5274{
5275  "EXAMPLE:";
5276  echo = 2;
5277  ring r = 0, (x,y,z), ds;
5278  poly f = (x2+3y-2z)^2+xyz-(x-y3+x2z3)^3;
5279  realmorsesplit(f);
5280}
5281
5282///////////////////////////////////////////////////////////////////////////////
5283/*
5284   symmetric Gauss algorithm
5285
5286   If A is not a square matrix, then the largest upper or left submatrix
5287   is assumed to be symmetric.
5288*/
5289static proc sym_reduce(matrix A)
5290{
5291  int r = nrows(A);
5292  int c = ncols(A);
5293  int n = r;
5294  if(n > c)
5295  {
5296    n = c;
5297  }
5298  poly q;
5299  int i, j;
5300  for(i = 1; i <= n; i++)
5301  {
5302    for(j = i+1; j <= n; j++)
5303    {
5304      if(A[i, j] != 0)
5305      {
5306        while(A[i, i] == 0)
5307        {
5308          A[1..r, i] = A[1..r, i]+A[1..r, j];
5309          A[i, 1..c] = A[i, 1..c]+A[j, 1..c];
5310        }
5311        q = A[i, j]/A[i, i];
5312        A[1..r, j] = A[1..r, j]-q*A[1..r, i];
5313        A[j, 1..c] = A[j, 1..c]-q*A[i, 1..c];
5314      }
5315    }
5316  }
5317  return(A);
5318}
5319
5320///////////////////////////////////////////////////////////////////////////////
5321/*
5322   - apply jet(f, k)
5323   - rewrite f as f = a*var(i)^2+p*var(i)+r with
5324     var(i)-free p and r
5325*/
5326static proc rewriteformorsesplit(poly f, int k, int i)
5327{
5328  f = jet(f, k);
5329  matrix C = coeffs(f, var(i));
5330  poly r = C[1,1];
5331  poly p = C[2,1];
5332  poly a = (f-r-p*var(i))/var(i)^2;
5333  return(f, a, p, r);
5334}
5335
5336///////////////////////////////////////////////////////////////////////////////
5337proc milnornumber(poly f)
5338"
5339USAGE:    milnornumber(f); f poly
5340RETURN:   Milnor number of f, or -1 if the Milnor number is not finite
5341KEYWORDS: Milnor number
5342NOTE:     The monomial order must be local.
5343EXAMPLE:  example milnornumber; shows an example"
5344{
5345  /* error check */
5346  int i;
5347  for(i = nvars(basering); i > 0; i--)
5348  {
5349    if(var(i) > 1)
5350    {
5351      ERROR("The monomial order must be local.");
5352    }
5353  }
5354
5355  return(vdim(std(jacob(f))));
5356}
5357example
5358{
5359  "EXAMPLE:";
5360  echo = 2;
5361  ring r = 0, (x,y), ds;
5362  poly f = x3+y4;
5363  milnornumber(f);
5364}
5365
5366///////////////////////////////////////////////////////////////////////////////
5367proc determinacy(poly f, list #)
5368"
5369USAGE:    determinacy(f[, mu]); f poly, mu int
5370RETURN:   an upper bound for the determinacy of f
5371NOTE:     The characteristic of the basering must be zero, the monomial order
5372          must be local, f must be contained in maxideal(1) and the Milnor
5373          number of f must be finite.
5374          @* The Milnor number of f can be provided as an optional parameter in
5375             order to avoid that it is computed again.
5376SEE ALSO: milnornumber, highcorner
5377KEYWORDS: Determinacy
5378EXAMPLE:  example determinacy; shows an example"
5379{
5380  /* auxiliary variables */
5381  int i;
5382  def br = basering;
5383
5384  /* error check */
5385  if(char(br) != 0)
5386  {
5387    ERROR("The characteristic must be zero.");
5388  }
5389  int n = nvars(br);
5390  for(i = 1; i <= n; i++)
5391  {
5392    if(var(i) > 1)
5393    {
5394      ERROR("The monomial order must be local.");
5395    }
5396  }
5397  if(jet(f, 0) != 0)
5398  {
5399    ERROR("The input polynomial must be contained in maxideal(1).");
5400  }
5401
5402  /* get Milnor number before continuing error check */
5403  int mu;
5404  if(size(#) > 0)   // read optional parameter
5405  {
5406    if(size(#) > 1 || typeof(#[1]) != "int")
5407    {
5408      ERROR("Wrong optional parameters.");
5409    }
5410    else
5411    {
5412      mu = #[1];
5413    }
5414  }
5415  else              // compute Milnor number
5416  {
5417    mu = milnornumber(f);
5418  }
5419
5420  /* continue error check */
5421  if(mu < 0)
5422  {
5423    ERROR("The Milnor number of the input polynomial must be"+newline
5424      +"non-negative and finite.");
5425  }
5426
5427  int k;   // an upper bound for the determinacy,
5428           // we use several methods:
5429
5430  /* Milnor number */
5431  k = mu+1;
5432  f = jet(f, k);
5433
5434  /* highest corner */
5435  int hc;
5436  list lbr = ringlist(br);
5437  if((ordstr(br)[1,2] != "ds")||(size(lbr[3])!=2))
5438  {
5439    lbr[3] = list(list("ds", 1:nvars(br)), list("C", 0));
5440    def br_ds = ring(lbr);
5441    setring(br_ds);
5442    poly f = fetch(br, f);
5443  }
5444  for(i = 0; i < 3; i++)
5445  {
5446    hc = deg(highcorner(std(maxideal(i)*jacob(f))));
5447    hc = hc+2-i;
5448    if(hc < k)
5449    {
5450      k = hc;
5451      f = jet(f, k);
5452    }
5453  }
5454  setring(br);
5455  return(k);
5456}
5457example
5458{
5459  "EXAMPLE:";
5460  echo = 2;
5461  ring r = 0, (x,y), ds;
5462  poly f = x3+xy3;
5463  determinacy(f);
5464}
5465
5466
5467
Note: See TracBrowser for help on using the repository browser.