source: git/Singular/LIB/puiseuxexpansions.lib @ eb2c4ac

spielwiese
Last change on this file since eb2c4ac was eb2c4ac, checked in by Hans Schoenemann <hannes@…>, 7 months ago
avoid names x,y in puiseuxexpansions.lib
  • Property mode set to 100644
File size: 44.6 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2///////////////////////////////////////////////////////////////////////////////
3///////////////////////////////////////////////////////////////////////////////
4version="version puiseuxexpansion.lib 4.3.1.0 Jul_2022 "; // $Id$
5category="Commutative Algebra";
6info="
7LIBRARY:  puiseuxexpansion.lib  Puiseux expansions over algebraic extensions
8AUTHORS:  J. Boehm, j.boehm at mx.uni-saarland.de
9          W. Decker, decker at mathematik.uni-kl.de
10          S. Laplagne, slaplagn at dm.uba.ar
11          G. Pfister, seelisch at mathematik.uni-kl.de
12
13OVERVIEW:
14This library implements the Newton-Puiseux algorithm to compute Puiseux
15expansions and provides a class and procedures to work with Puiseux
16expansions in K[X][Y], where K is the field Q of rational numbers
17or an algebraic extension of Q.
18
19PROCEDURES:
20 puiseuxList(PP, maxDeg, iVarX, iVarY);     Computes the Puiseux expansion of PP at X = 0
21 makePuiseux(f, denom, fr);       Creates a Puiseux element
22 makePuiseuxFromPoly(poly f);     Converts a polynomail to Puiseux data type
23 printPuiseux(f);                 Prints information for Puiseux elements
24 puiseux(f, maxDeg, atOrigin);    Computes the Puiseux expansions of f
25";
26
27LIB "normal.lib";
28LIB "assprimeszerodim.lib";
29LIB "integralbasis.lib";
30
31
32static proc mod_init()
33{
34  newstruct("Puiseux","ring in, poly numer, poly denom, int fraction");
35  system("install","Puiseux","print",printPuiseux,1);
36  system("install","Puiseux","+",addPuiseux,2);
37  system("install","Puiseux","*",multPuiseux,2);
38  system("install","Puiseux","^",expPuiseux,2);
39  system("install","Puiseux","==",equalPuiseux,2);
40  system("install","Puiseux","=",makePuiseuxFromPoly,1);
41}
42
43// Extended the basering of P to ring S if possible.
44// That is, all the variables used in P must be also present in S.
45static proc extendBaseringPuiseux(Puiseux P, def S)
46{
47  def R = basering;
48  def R1 = P.in;
49  setring R1;
50  poly Pnumer = P.numer;
51  poly Pdenom = P.denom;
52  setring S;
53  Puiseux Q;
54  Q.in = basering;
55  Q.numer = imap(R1, Pnumer);
56  Q.denom = imap(R1, Pdenom);
57  Q.fraction = P.fraction;
58  setring R;
59  return(Q);
60}
61
62// Substitue var(idx) from f by Puiseux expansion g
63static proc substPuiseux(Puiseux f, int idx, Puiseux g)
64{
65  if (idx==1){ERROR("First variable is used for Puiseux series coefficients");}
66
67  def R = basering;
68  def S = f.in;
69  setring S;
70  int lcmfr = lcm(f.fraction,g.fraction);
71  int s1 = lcmfr div (f.fraction);
72  int s2 = lcmfr div (g.fraction);
73  poly ff = subst(f.numer,var(1),var(1)^s1);
74  poly gg = subst(g.numer,var(1),var(1)^s2);
75  poly fd = subst(f.denom,var(1),var(1)^s1);
76  poly gd = subst(g.denom,var(1),var(1)^s2);
77  def S1 =addvarsTo(S,"@h",0);
78  setring S1;
79  intvec wt=intvec(0:(nvars(S1)-1),1);
80  wt[idx]=1;
81  def S2=changeord(list(list("wp",wt)));
82  setring S2;
83  poly ff = fetch(S,ff);
84  poly gg = fetch(S,gg);
85  poly gd = fetch(S,gd);
86  poly ffh = homog(ff,@h);
87  int dg = deg(ffh);
88  ffh = subst(ffh,@h,gd);
89  ffh = subst(ffh,var(idx),gg);
90  kill ff;
91  kill gg;
92
93  setring S;
94  kill ff;
95  kill gg;
96
97  poly ffh = fetch(S2,ffh);
98  Puiseux F;
99  F.in = S;
100
101  F.numer = ffh;
102  kill ffh;
103
104  F.denom = fd*gd^dg;
105  kill fd;
106  kill gd;
107  kill dg;
108
109  F.fraction = lcmfr;
110  kill lcmfr;
111
112  setring R;
113
114  return(F);
115}
116
117proc makePuiseuxFromPoly(poly f)
118"USAGE: makePuiseuxFromPoly(f); f polynomial in two variables
119RETURN: make a ring independent polynomial over Puiseux series
120EXAMPLE: example makePuiseuxFromPoly, shows an example"
121{
122  return(makePuiseux(f,1,1));
123}
124example
125{
126 "EXAMPLE:"; echo=2;
127 ring R=0,(x,y),dp;
128 poly f=3*x^2+x+1;
129 makePuiseuxFromPoly(f);
130}
131
132proc makePuiseux(poly f, poly denom, int fr)
133"USAGE: makePuiseux(f, denom, fr); f polynomial in two variables, denom polynomial in
134        the first variable of the ring, int fr
135RETURN: make a ring independent polynomial over Puiseux series
136EXAMPLE: example makePuiseux, shows an example"
137{
138  Puiseux F;
139  F.in=basering;
140  F.numer=f;
141  F.denom =denom;
142  F.fraction = fr;
143  return(F);
144}
145example
146{
147 "EXAMPLE:"; echo=2;
148 ring R=0,(x,y),dp;
149 poly f=3*x^2+x*y+1;
150 makePuiseux(f,x^2,3);
151}
152
153// Print information for Puiseux elements
154proc printPuiseux(Puiseux f)
155"USAGE: printPuiseux(f); f Puiseux expansion
156RETURN: prints information for Puiseux elements
157EXAMPLE: example printPuiseux, shows an example"
158{
159  def R = basering;
160  def S = f.in;
161  setring S;
162  "(", f.numer, ") / ", f.denom;
163  "Denominator of exponent : ",  f.fraction;
164  if(minpoly <> 0) {"Minimal polynomial: ", minpoly;}
165  setring R;
166}
167example
168{
169 "EXAMPLE:"; echo=2;
170 ring R=0,(x,y),dp;
171 poly f=3*x^2+x*y+1;
172 Puiseux F = makePuiseux(f,x^2,3);
173 printPuiseux(F);
174
175}
176
177// Sum of two Puiseux elements
178static proc addPuiseux(Puiseux f,Puiseux g)
179{
180  def R=basering;
181  def S1=f.in;
182  setring S1;
183  int lcmfr = lcm(f.fraction,g.fraction);
184  int s1 = lcmfr div (f.fraction);
185  int s2 = lcmfr div (g.fraction);
186  poly ff = subst(f.numer,var(1),var(1)^s1);
187  poly gg = subst(g.numer,var(1),var(1)^s2);
188  poly fd = subst(f.denom,var(1),var(1)^s1);
189  poly gd = subst(g.denom,var(1),var(1)^s2);
190  Puiseux F;
191  poly comden = lcm(fd,gd);
192  F.numer  = (comden/fd)*ff+ (comden/gd)*gg ;
193  F.denom = comden;
194  F.fraction = lcmfr;
195  F.in = S1;
196  setring R;
197  return(F);
198}
199
200// Product of two Puiseux elements
201static proc multPuiseux(Puiseux f,Puiseux g)
202{
203  def R=basering;
204  def S1=f.in;
205  setring S1;
206  int lcmfr = lcm(f.fraction,g.fraction);
207  int s1 = lcmfr div (f.fraction);
208  int s2 = lcmfr div (g.fraction);
209  poly ff = subst(f.numer,var(1),var(1)^s1);
210  poly gg = subst(g.numer,var(1),var(1)^s2);
211  poly fd = subst(f.denom,var(1),var(1)^s1);
212  poly gd = subst(g.denom,var(1),var(1)^s2);
213  poly comden = lcm(fd,gd);
214  Puiseux F;
215  F.numer  = ff*gg ;
216  F.denom = fd*gd;
217  F.fraction = lcmfr;
218  F.in = S1;
219  setring R;
220  return(F);
221}
222
223// Compares if two Puiseux elements are equal.
224// It takes commond denominator for the rational functions and
225// also for the denominator of the exponents before comparing.
226static proc equalPuiseux(Puiseux f,Puiseux g)
227{
228  def R=basering;
229  def S1=f.in;
230  setring S1;
231  int lcmfr = lcm(f.fraction,g.fraction);
232  int s1 = lcmfr div (f.fraction);
233  int s2 = lcmfr div (g.fraction);
234  poly ff = subst(f.numer,var(1),var(1)^s1);
235  poly gg = subst(g.numer,var(1),var(1)^s2);
236  poly fd = subst(f.denom,var(1),var(1)^s1);
237  poly gd = subst(g.denom,var(1),var(1)^s2);
238  int result = (ff*gd)==(gg*fd);
239  kill ff,gg,fd,gd; // these stay in f.in for some reason, so delete them
240  setring R;
241  return(result);
242}
243
244// Raise Puiseux expansion f to the power n.
245static proc expPuiseux(Puiseux f, int n)
246{
247  def R=basering;
248  def S1=f.in;
249  setring S1;
250  Puiseux F;
251  F.numer=(f.numer)^n;
252  F.denom =(f.denom)^n;
253  F.fraction = f.fraction;
254  F.in = S1;
255  setring R;
256  return(F);
257}
258
259
260
261//-----------------------------------------------------------------------------
262// Puiseux expansions
263//-----------------------------------------------------------------------------
264
265// Main call por computing Puiseux expansions
266// It returns a list containing information about the expansions
267proc puiseux(poly f, int maxDeg, int atOrigin)
268"USAGE: puiseux(f, maxDeg, atOrigin); f polynomial in two variables, int maxDeg, int atOrigin
269RETURN: the Puiseux expansions of f developed up to degree maxDeg. If atOrigin = 1, only the expansions passing through the origin will be returned.
270EXAMPLE: example puiseux, shows an example"
271{
272  if(atOrigin == 1)
273  {
274    list p = puiseuxMainOneExpansion(f, maxDeg, 0, 1, 1);
275    //list p = puiseuxMain(f, maxDeg, 0, 1, 1);
276  } else
277  {
278    list p = puiseuxMainOneExpansion(f, maxDeg, 0, 0, 1);
279    //list p = puiseuxMain(f, maxDeg, 0, 0, 1);
280  }
281
282  return(p);
283}
284example
285{
286 "EXAMPLE:"; echo=2;
287 ring R=0,(x,y),dp;
288 poly f=y^3 + x^2 + x^8;
289 puiseux(f,3,0);
290}
291
292
293// Finds the coefficient of y^slope in the Puiseux expansion
294// slN is the numerator of the slope
295// slD is the denominator
296static proc puiseuxStep(poly f, int slN, int slD){
297  def R = basering;
298
299  // We treat the negative exponents separately
300  poly f2 = subst(f, var(1), var(1)^slD);
301  if(slN >= 0)
302  {
303    poly f3 = subst(f2, var(2), var(2)*(var(1)^slN));
304    matrix M = coef(f3,var(1));
305    poly c = M[2, ncols(M)];
306  } else
307  {
308    // All this is a workaround to work with negative exponents,
309    // which are not allowed in Singular.
310    int dy = deg(f2, intvec(0,1));
311    int dx3, dy3;
312
313    // We multiply by a large enough power of x so that the next
314    // substitution is possible.
315    poly f3 = f2*var(1)^(-slN*dy);
316
317    // Now we replace x^(|slN|)*y by  x^(|slN|)*x^slN*y = y
318    poly f3Temp = f3;
319    poly f4;
320    poly leadP;
321    int mExp;
322    for(int i = 1; i <= size(f3); i++)
323    {
324      leadP = lead(f3Temp);
325      dy3 = deg(leadP, intvec(0,1));
326      f4 = f4 + leadP / var(1)^(dy3*(-slN));
327      f3Temp = f3Temp - leadP;
328    }
329    // We can now take the equation to vanish the term with
330    // smalles degree in x
331    matrix M = coef(f4,var(1));
332    poly c = M[2, ncols(M)];
333  }
334
335  return(c);
336}
337
338
339// Computes the Puiseux expansions starting with slope >= sN / sD;
340
341// Output:
342// cs[1] = Puiseux expansion
343// cs[2] = denominator of all the exponents
344// cs[6] = code for identifying different classes
345// cs[7] = exponents where new branches appear (this is used for computing
346//         polynomials over the ground field).
347static proc puiseuxMain(poly f, int maxDeg, int sN, int sD, int firstTime)
348{
349//  if(firstTime == 1)
350//  {
351//    "Computing Puiseux expansions for f = ", f;
352//  }
353//  else
354//  {
355//    "PuiseuxMain - not first time - f", f;
356//  }
357
358  int dbg = 0;
359
360  def R = basering;
361  intvec vy = (0,1);
362  int d = deg(f, vy);
363  int h,i,j, ii;
364  int k = 1;
365  int slN, slD;
366  int g;
367  int sizeErg;
368  int stop;
369  int newMaxDeg;
370
371  poly @c;
372  int newExt;
373  list fc;
374  poly minP;
375  list cs;    // coefficients of the p. exp
376  list csT;
377  list out;
378  poly je;    // Minimal equation
379  list fJe;   // Factorization of the initial equation
380  matrix cofs;
381  number c1;
382  poly fNew1, fNew2;
383  poly fTemp;
384  int cod = 0;
385  int mi;
386
387  // Case of Puiseux expansions with finite number of terms
388  list dd = Integralbasis::divideBy(f, var(2));
389  if(dd[2]>0)
390  {
391    for(i = 1; i <= dd[2]; i++)
392    {
393      cs[size(cs) + 1] = list(0, 1);
394      cs[size(cs)][6] = list(cod);
395      cod++;
396    }
397    f = dd[1];
398  }
399
400
401  list l = newtonpoly2(f);
402
403  //list l = newtonpoly(f);
404
405  //if((l[1][1] == 0) && (l[1][2] == 0))
406  //{
407  //  "ERROR: The polynomial must pass through the origin.";
408  //  list cs = list();
409  //} else
410  //{
411
412
413  //for(i = 1; i<size(l); i++)
414  for(i = size(l)-1; i>=1; i--)
415  {
416    slN = l[i+1][1] - l[i][1];
417    slD = l[i][2] - l[i+1][2];
418
419    // We always use positive denominator
420    if(slD < 0){
421      slN = (-1)*slN;
422      slD = (-1)*slD;
423    }
424
425    //if(firstTime == 1)
426    //{
427    //  "Computing for segment ", i, " with slope ", slN, " / ", slD;
428    //} else
429    //{
430    //  "Further developing at segment ", i, " with slope ", slN, " / ", slD;
431    //}
432
433    if(slN != 0)
434    {
435      g = gcd(slD, slN);
436    } else {
437      if(slD != 0)
438      {
439        g = slD;
440      } else
441      {
442        g = 1;
443      }
444    }
445
446    slN = slN div g;
447    slD = slD div g;
448
449
450    // sD = sN = 0 indicates all slopes must be used.
451
452    //if(bigint(sD) * slN > bigint(sN) * slD)
453    //if((bigint(sD) * slN > bigint(sN) * slD) or ((sD == 0) and (sN == 0)))
454    if((slN > 0) or ((sD == 0) and (sN == 0)))
455    {
456      je = minEqNewton(f, slN, slD);
457
458      fJe = factorize(je);
459      fJe = Integralbasis::sortFactors(fJe);
460
461      fNew1 = subst(f, var(1), var(1)^slD);
462      for(ii = 2; ii <= size(fJe[1]); ii++)
463      {
464        if((ii == 2) || (nonRatCoeff(fJe[1][ii]) == 0))  // If the factor is over the ground field then it is a new conjugacy class.
465        {
466          cod++;
467        }
468        fTemp = subst(fJe[1][ii]^(fJe[2][ii]), var(1), var(1)^slD);
469        // Checks if a new extension is needed.
470        // If so, we have a new cutting point for building the factors.
471        if((deg(fJe[1][ii], vy) > 1) or (size(fJe[1]) > 2))
472        {
473          newExt = 1;
474        } else
475        {
476          newExt = 0;
477        }
478
479        @c = puiseuxStep(fTemp, slN, 1);
480        fc = factorize(@c);
481       
482        for(j = 2; j <= size(fc[1]); j++)
483        {
484
485          // Stopping criterium
486          stop = 0;
487          if(maxDeg <= 0){
488            if(fc[2][j] == 1){
489              stop = 1;
490            }
491          } else {
492            if(slN >= slD * maxDeg){
493              stop = 1;
494            }
495          }
496
497          if(fc[1][j] != var(2))
498          {
499            minP = fc[1][j];
500            if(deg(minP)==1)
501            {
502              cofs = coeffs(minP, var(2));
503              c1 = number(-cofs[1,1]/cofs[2,1]);
504              if(defined(erg))
505              {
506                sizeErg = size(erg);
507              } else
508              {
509                sizeErg = 1;
510              }
511              if(stop == 0)
512              {
513                if(slN >= 0)
514                {
515                  fNew2 = subst(fNew1, var(2), (c1+var(2))*(var(1)^slN));
516                } else
517                {
518                  // Negative exponent case
519                  fNew2 = negExp(fNew1, c1, slN);
520                }
521                fNew2= sat(ideal(fNew2), var(1))[1];
522
523                if(maxDeg > 0)
524                {
525                  newMaxDeg = maxDeg * slD - slN;
526                } else {
527                  newMaxDeg = maxDeg;
528                }
529
530                csT = puiseuxMain(fNew2, newMaxDeg, slN, 1, 0);
531
532                for(k = 1; k<=size(csT); k++)
533                {
534                  if(typeof(csT[k]) != "ring")
535                  {
536                    // Case of a polynomial with the corresponding denominator;
537                    if(slN >= 0)
538                    {
539                      cs[size(cs) + 1] = list((csT[k][1] + c1)*(var(1)^(csT[k][2]*slN)), csT[k][2]*slD);
540                    } else
541                    {
542                      cs[size(cs) + 1] = list(csT[k][1] + c1, csT[k][2]*slD);
543                      cs[size(cs)][8] = list("Denominator", var(1) ^ (-slN*csT[k][2]));
544                    }
545
546
547                    // NEW CODING FOR CONJUGATED EXPANSIONS
548                    cs[size(cs)][6] = insert(csT[k][6], cod);
549                   
550                    if(newExt == 1)
551                    {
552                      if(typeof(csT[k][7])=="none")
553                      {
554                        cs[size(cs)][7] = list(0);
555                      } else
556                      {
557                        cs[size(cs)][7] = insert(csT[k][7], 0);
558                      }
559                      cs[size(cs)][7] = sum2All(cs[size(cs)][7], slN * csT[k][2]);
560                    } else
561                    {
562                      if(typeof(csT[k][7]) != "none")
563                      {
564                        cs[size(cs)][7] = sum2All(csT[k][7], slN * csT[k][2]);
565                      }
566                    }
567                  } else
568                  {
569                    def RT = csT[k];
570                    setring RT;
571                    number c1 = fetch(R, c1);
572                    // Change @a by par(1)?
573                    c1 = number(subst(c1, @a, number(erg[sizeErg])));
574
575                    for(h = 1; h <= size(PE); h++)
576                    {
577                      if(slN >= 0)
578                      {
579                        PE[h][1] = (PE[h][1] + c1) * (var(1)^(slN*PE[h][2]));
580                        if(newExt == 1)
581                        {
582                          PE[h][7] = insert(PE[h][7], 0);
583                        }
584                        PE[h][7] = sum2All(PE[h][7], slN * PE[h][2]);
585                        PE[h][2] = PE[h][2] * slD;
586
587                        // NEW CODING FOR CONJUGATED EXPANSIONS
588                        PE[h][6] = insert(PE[h][6], cod);
589                      } else {
590                        PE[h][1] = PE[h][1] + c1;
591                        if(newExt == 1)
592                        {
593                          PE[h][7] = insert(PE[h][7], 0);
594                        }
595
596                        // Check if this is correct
597                        // (before PE[h][8] was defined after PE[h][2] was modified)
598                        PE[h][8] = list("Denominator", var(1) ^ (-slN*PE[h][2]));
599                        PE[h][2] = PE[h][2] * slD;
600
601                        // NEW CODING FOR CONJUGATED EXPANSIONS
602                        PE[h][6] = insert(PE[h][6], cod);
603                      }
604                    }
605                    setring R;
606                    cs[size(cs) + 1] = RT;
607                    kill RT;
608                  }
609                }
610              } else {
611                cs[size(cs) + 1] = list(c1 * var(1)^slN, slD);
612                cs[size(cs)][6] = list(cod);
613                if(newExt == 1)
614                {
615                  cs[size(cs)][7] = list(slN);
616                }
617              }
618            } else {
619              if(npars(R) == 1)
620              {
621                if(!defined(erg))
622                {
623                  list erg;
624                  erg[1] = par(1);
625                }
626                if(!defined(minPolys))
627                {
628                  list minPolys;
629                  minPolys[1] = 1;
630                }
631
632                def S = Integralbasis::splitRingAt(minP, erg);
633                setring S;
634                sizeErg = size(erg);
635                poly newA = erg[size(erg)];  // The root founded. That is, newA is the root of minP in S.
636                poly fNew1 = fetch(R, fNew1);
637
638                if(!defined(minPolys))
639                {
640                  list minPolys = fetch(R, minPolys);
641                }
642
643                // We replace the old a by its image in the new ring.
644                minPolys[size(minPolys)+1] = fetch(R, minP);
645                for(mi = 1; mi <= size(minPolys); mi++)
646                {
647                  minPolys[mi] = subst(minPolys[mi], par(1), erg[size(erg)-1]);
648                }
649                fNew1 = subst(fNew1, par(1), erg[size(erg)-1]);
650
651              } else {
652                def S = Integralbasis::splitRingAt(minP);
653                setring S;
654                list minPolys;
655                minPolys[1] = 1;
656                minPolys[size(minPolys)+1] = fetch(R, minP);
657                poly fNew1 = fetch(R, fNew1);
658                poly newA = par(1);      // In this case, it is the new variable.
659                list erg = par(1);
660                export erg;
661                sizeErg = 1;
662              }
663              export(minPolys);
664
665              if(defined(PES))
666              {
667                PES = list();
668              } else
669              {
670                list PES;
671              }
672              if(stop == 0){
673                // We use the image of the root of minP in the new ring.
674                if(slN >= 0)
675                {
676                  poly fNew2 = subst(fNew1, var(2), (newA + var(2))*var(1)^slN);
677                } else
678                {
679                  poly fNew2 = negExp(fNew1, newA, slN);
680                }
681                fNew2= sat(ideal(fNew2), var(1))[1];
682                if(maxDeg > 0)
683                {
684                  newMaxDeg = maxDeg * slD - slN;
685                } else {
686                  newMaxDeg = maxDeg;
687                }
688                list csTS = puiseuxMain(fNew2, newMaxDeg, slN, 1, 0);
689                for(k = 1; k<=size(csTS); k++)
690                {
691                  if(typeof(csTS[k]) == "ring")
692                  {
693                    def RT = csTS[k];
694                    setring RT;
695                    for(h = 1; h<=size(PE); h++)
696                    {
697                      // The new coefficient is the image of a in the new ring.
698                      if(slN >= 0)
699                      {
700                        PE[h][1] = (PE[h][1] + erg[sizeErg]) * var(1)^(PE[h][2]*slN);
701                        if(newExt == 1)
702                        {
703                          PE[h][7] = insert(PE[h][7], 0);
704                        }
705                        PE[h][7] = sum2All(PE[h][7], slN * PE[h][2]);
706                        PE[h][2] = PE[h][2] * slD;
707
708                        // NEW CODING FOR CONJUGATED EXPANSIONS
709                        PE[h][6] = insert(PE[h][6], cod);
710                      } else
711                      {
712                        PE[h][1] = (PE[h][1] + erg[sizeErg]);
713                        if(newExt == 1)
714                        {
715                          PE[h][7] = insert(PE[h][7], 0);
716                        }
717                        PE[h][7] = sum2All(PE[h][7], slN * PE[h][2]);
718                        PE[h][2] = PE[h][2] * slD;
719
720                        // NEW CODING FOR CONJUGATED EXPANSIONS
721                        PE[h][6] = insert(PE[h][6], cod);
722
723                        PE[h][8] = list("Denominator", var(1) ^ (-slN*PE[h][2]));
724                      }
725                    }
726                    setring R;
727                    cs[size(cs) + 1] = RT;
728                    kill RT;
729                    setring S;
730                  } else {
731                    if(slN >= 0)
732                    {
733                      PES[size(PES) + 1] = list((csTS[k][1] + newA)*var(1)^(csTS[k][2]*slN), csTS[k][2]*slD);
734                    } else {
735                      PES[size(PES) + 1] = list(csTS[k][1] + newA, csTS[k][2]*slD);
736                      PES[size(PES)][8] = list("Denominator", var(1) ^ (-slN*csTS[k][2]));
737                    }
738
739                    // NEW CODING FOR CONJUGATED EXPANSIONS
740                    PES[size(PES)][6] = insert(csTS[k][6], cod);
741
742                    if(newExt == 1)
743                    {
744                      if(typeof(csTS[k][7]) == "none")
745                      {
746                        PES[size(PES)][7] = list(0);
747
748                      } else
749                      {
750                        PES[size(PES)][7] = insert(csTS[k][7], 0);
751                      }
752                      PES[size(PES)][7] = sum2All(PES[size(PES)][7], slN * csTS[k][2]);
753                    } else
754                    {
755                      PES[size(PES)][7] = csTS[k][7];
756                    }
757
758                  }
759                }
760                if(size(PES) > 0){
761                  list PE = PES;
762                  export PE;
763                  setring R;
764                  cs[size(cs) + 1] = S;
765                }
766              } else {
767                poly PE1 = newA*var(1)^slN;
768                list PE = list(list(PE1, slD));
769                PE[1][6] = list(cod);
770                if(newExt == 1)
771                {
772                  PE[1][7] = list(slN);
773                }
774                export PE;
775
776                setring R;
777                cs[size(cs) + 1] = S;
778              }
779              kill S;
780            }
781          }
782          setring R;
783        }
784      }
785    }
786  }
787
788  if(size(cs) == 0){
789    cs = list(list(poly(0),1));
790    cs[1][6] = list(1);
791  }
792
793  return(cs);
794}
795
796
797
798///////////////////////////////////////////////////////////////////////
799
800// Computes the Puiseux expansions starting with slope >= sN / sD;
801
802// Output:
803// cs[1] = Puiseux expansion
804// cs[2] = denominator of all the exponents
805// cs[6] = code for identifying different classes
806// cs[7] = exponents where new branches appear (this is used for computing
807//         polynomials over the ground field).
808
809// It computes only one expansion in each class
810static proc puiseuxMainOneExpansion(poly f, int maxDeg, int sN, int sD, int firstTime)
811{
812  //if(firstTime == 1)
813  //{
814  //  "Computing Puiseux expansions for f = ", f;
815  //}
816  //else
817  //{
818  //  "PuiseuxMain - not first time - f", f;
819  //}
820  //
821  //~;
822
823  int dbg = 0;
824
825  def R = basering;
826  intvec vy = (0,1);
827  int d = deg(f, vy);
828  int h,i,j, ii;
829  int k = 1;
830  int slN, slD;
831  int g;
832  int sizeErg;
833  int stop;
834  int newMaxDeg;
835
836  poly @c;
837  int newExt;
838  list fc;
839  poly minP;
840  list cs;    // coefficients of the p. exp
841  list csT;
842  list out;
843  poly je;    // Minimal equation
844  list fJe;   // Factorization of the initial equation
845  matrix cofs;
846  number c1;
847  poly fNew1, fNew2;
848  poly fTemp;
849  int cod = 0;
850  int mi;
851
852  poly minPBefore;
853  list minPFacs;
854
855  // Case of Puiseux expansions with finite number of terms
856  list dd = Integralbasis::divideBy(f, var(2));
857  if(dd[2]>0)
858  {
859    for(i = 1; i <= dd[2]; i++)
860    {
861      cs[size(cs) + 1] = list(0, 1);
862      cs[size(cs)][6] = list(cod);
863      cod++;
864    }
865    f = dd[1];
866  }
867
868
869  list l = newtonpoly2(f);
870
871  //list l = newtonpoly(f);
872
873  //if((l[1][1] == 0) && (l[1][2] == 0))
874  //{
875  //  "ERROR: The polynomial must pass through the origin.";
876  //  list cs = list();
877  //} else
878  //{
879
880
881  //for(i = 1; i<size(l); i++)
882  for(i = size(l)-1; i>=1; i--)
883  {
884    slN = l[i+1][1] - l[i][1];
885    slD = l[i][2] - l[i+1][2];
886
887    // We always use positive denominator
888    if(slD < 0){
889      slN = (-1)*slN;
890      slD = (-1)*slD;
891    }
892
893    //if(firstTime == 1)
894    //{
895    //  "Computing for segment ", i, " with slope ", slN, " / ", slD;
896    //} else
897    //{
898    //  "Further developing at segment ", i, " with slope ", slN, " / ", slD;
899    //}
900
901    if(slN != 0)
902    {
903      g = gcd(slD, slN);
904    } else {
905      if(slD != 0)
906      {
907        g = slD;
908      } else
909      {
910        g = 1;
911      }
912    }
913
914    slN = slN div g;
915    slD = slD div g;
916
917
918    // sD = sN = 0 indicates all slopes must be used.
919
920    //if(bigint(sD) * slN > bigint(sN) * slD)
921    //if((bigint(sD) * slN > bigint(sN) * slD) or ((sD == 0) and (sN == 0)))
922    if((slN > 0) or ((sD == 0) and (sN == 0)))
923    {
924      je = minEqNewton(f, slN, slD);
925
926      fJe = factorize(je);
927      fJe = Integralbasis::sortFactors(fJe);
928
929      fNew1 = subst(f, var(1), var(1)^slD);
930     
931      for(ii = 2; ii <= size(fJe[1]); ii++)
932      {
933        if((ii == 2) || (nonRatCoeff(fJe[1][ii]) == 0))  // If the factor is over the ground field then it is a new conjugacy class.
934        {
935          cod++;
936        }
937        fTemp = subst(fJe[1][ii]^(fJe[2][ii]), var(1), var(1)^slD);
938        // Checks if a new extension is needed.
939        // If so, we have a new cutting point for building the factors.
940        if((deg(fJe[1][ii], vy) > 1) or (size(fJe[1]) > 2))
941        {
942          newExt = 1;
943        } else
944        {
945          newExt = 0;
946        }
947
948        @c = puiseuxStep(fTemp, slN, 1);
949        fc = factorize(@c);
950       
951        // Stopping criterium
952        stop = 0;
953        if(maxDeg <= 0)
954        {
955          if(fJe[2][ii] == 1)
956          {
957            stop = 1;
958          }
959        } else
960        {
961          if(slN >= slD * maxDeg)
962          {
963            stop = 1;
964          }
965        }
966
967        if(fJe[1][ii] != var(2))
968        {
969          //minP = fc[1][j];
970         
971          // The minpoly of Y^5-1 is Y-1 and the minpoly of Y^5-2 is Y^5-2.
972          minPBefore = subst(fJe[1][ii],var(1),1);
973          minPFacs = factorize(minPBefore);
974          minP = minPFacs[1][2];   // Check if this always the smallest degree
975
976          if(deg(minP)==1)
977          {
978            cofs = coeffs(minP, var(2));
979            c1 = number(-cofs[1,1]/cofs[2,1]);
980            if(defined(erg))
981            {
982              sizeErg = size(erg);
983            } else
984            {
985              sizeErg = 1;
986            }
987            if(stop == 0)
988            {
989              if(slN >= 0)
990              {
991                fNew2 = subst(fNew1, var(2), (c1+var(2))*(var(1)^slN));
992              } else
993              {
994                // Negative exponent case
995                fNew2 = negExp(fNew1, c1, slN);
996              }
997              fNew2= sat(ideal(fNew2), var(1))[1];
998
999              if(maxDeg > 0)
1000              {
1001                newMaxDeg = maxDeg * slD - slN;
1002              } else {
1003                newMaxDeg = maxDeg;
1004              }
1005
1006              csT = puiseuxMainOneExpansion(fNew2, newMaxDeg, slN, 1, 0);
1007
1008              for(k = 1; k<=size(csT); k++)
1009              {
1010                if(typeof(csT[k]) != "ring")
1011                {
1012                  // Case of a polynomial with the corresponding denominator;
1013                  if(slN >= 0)
1014                  {
1015                    cs[size(cs) + 1] = list((csT[k][1] + c1)*(var(1)^(csT[k][2]*slN)), csT[k][2]*slD);
1016                  } else
1017                  {
1018                    cs[size(cs) + 1] = list(csT[k][1] + c1, csT[k][2]*slD);
1019                    cs[size(cs)][8] = list("Denominator", var(1) ^ (-slN*csT[k][2]));
1020                  }
1021
1022
1023                  // NEW CODING FOR CONJUGATED EXPANSIONS
1024                  cs[size(cs)][6] = insert(csT[k][6], cod);
1025                 
1026                  if(newExt == 1)
1027                  {
1028                    if(typeof(csT[k][7])=="none"){
1029                      cs[size(cs)][7] = list(0);
1030                    } else
1031                    {
1032                      cs[size(cs)][7] = insert(csT[k][7], 0);
1033                    }
1034                    cs[size(cs)][7] = sum2All(cs[size(cs)][7], slN * csT[k][2]);
1035                  } else
1036                  {
1037                    if(typeof(csT[k][7]) != "none"){
1038                      cs[size(cs)][7] = sum2All(csT[k][7], slN * csT[k][2]);
1039                    }
1040                  }
1041                } else
1042                {
1043                  def RT = csT[k];
1044                  setring RT;
1045                  number c1 = fetch(R, c1);
1046                  // Change @a by par(1)?
1047                  c1 = number(subst(c1, @a, number(erg[sizeErg])));
1048
1049                  for(h = 1; h <= size(PE); h++)
1050                  {
1051                    if(slN >= 0)
1052                    {
1053                      PE[h][1] = (PE[h][1] + c1) * (var(1)^(slN*PE[h][2]));
1054                      if(newExt == 1)
1055                      {
1056                        PE[h][7] = insert(PE[h][7], 0);
1057                      }
1058                      PE[h][7] = sum2All(PE[h][7], slN * PE[h][2]);
1059                      PE[h][2] = PE[h][2] * slD;
1060
1061                      // NEW CODING FOR CONJUGATED EXPANSIONS
1062                      PE[h][6] = insert(PE[h][6], cod);
1063                    } else {
1064                      PE[h][1] = PE[h][1] + c1;
1065                      if(newExt == 1)
1066                      {
1067                        PE[h][7] = insert(PE[h][7], 0);
1068                      }
1069
1070                      // Check if this is correct
1071                      // (before PE[h][8] was defined after PE[h][2] was modified)
1072                      PE[h][8] = list("Denominator", var(1) ^ (-slN*PE[h][2]));
1073                      PE[h][2] = PE[h][2] * slD;
1074
1075                      // NEW CODING FOR CONJUGATED EXPANSIONS
1076                      PE[h][6] = insert(PE[h][6], cod);
1077                    }
1078                  }
1079                  setring R;
1080                  cs[size(cs) + 1] = RT;
1081                  kill RT;
1082                }
1083              }
1084            } else {
1085              cs[size(cs) + 1] = list(c1 * var(1)^slN, slD);
1086              cs[size(cs)][6] = list(cod);
1087              if(newExt == 1)
1088              {
1089                cs[size(cs)][7] = list(slN);
1090              }
1091            }
1092          } else {
1093            if(npars(R) == 1)
1094            {
1095              if(!defined(erg))
1096              {
1097                list erg;
1098                erg[1] = par(1);
1099              }
1100              if(!defined(minPolys))
1101              {
1102                list minPolys;
1103                minPolys[1] = 1;
1104              }
1105
1106              def S = Integralbasis::splitRingAt(minP, erg);
1107              setring S;
1108              sizeErg = size(erg);
1109              poly newA = erg[size(erg)];  // The root founded. That is, newA is the root of minP in S.
1110              poly fNew1 = fetch(R, fNew1);
1111
1112              if(!defined(minPolys))
1113              {
1114                list minPolys = fetch(R, minPolys);
1115              }
1116
1117              // We replace the old a by its image in the new ring.
1118              minPolys[size(minPolys)+1] = fetch(R, minP);
1119              for(mi = 1; mi <= size(minPolys); mi++)
1120              {
1121                minPolys[mi] = subst(minPolys[mi], par(1), erg[size(erg)-1]);
1122              }
1123              fNew1 = subst(fNew1, par(1), erg[size(erg)-1]);
1124
1125            } else {
1126              def S = Integralbasis::splitRingAt(minP);
1127              setring S;
1128              list minPolys;
1129              minPolys[1] = 1;
1130              minPolys[size(minPolys)+1] = fetch(R, minP);
1131              poly fNew1 = fetch(R, fNew1);
1132              poly newA = par(1);      // In this case, it is the new variable.
1133              list erg = par(1);
1134              export erg;
1135              sizeErg = 1;
1136            }
1137            export(minPolys);
1138
1139            if(defined(PES))
1140            {
1141              PES = list();
1142            } else
1143            {
1144              list PES;
1145            }
1146            if(stop == 0){
1147              // We use the image of the root of minP in the new ring.
1148              if(slN >= 0)
1149              {
1150                poly fNew2 = subst(fNew1, var(2), (newA + var(2))*var(1)^slN);
1151              } else
1152              {
1153                poly fNew2 = negExp(fNew1, newA, slN);
1154              }
1155              fNew2= sat(ideal(fNew2), var(1))[1];
1156              if(maxDeg > 0)
1157              {
1158                newMaxDeg = maxDeg * slD - slN;
1159              } else {
1160                newMaxDeg = maxDeg;
1161              }
1162              list csTS = puiseuxMainOneExpansion(fNew2, newMaxDeg, slN, 1, 0);
1163              for(k = 1; k<=size(csTS); k++)
1164              {
1165                if(typeof(csTS[k]) == "ring")
1166                {
1167                  def RT = csTS[k];
1168                  setring RT;
1169                  for(h = 1; h<=size(PE); h++)
1170                  {
1171                    // The new coefficient is the image of a in the new ring.
1172                    if(slN >= 0)
1173                    {
1174                      PE[h][1] = (PE[h][1] + erg[sizeErg]) * var(1)^(PE[h][2]*slN);
1175                      if(newExt == 1)
1176                      {
1177                        PE[h][7] = insert(PE[h][7], 0);
1178                      }
1179                      PE[h][7] = sum2All(PE[h][7], slN * PE[h][2]);
1180                      PE[h][2] = PE[h][2] * slD;
1181
1182                      // NEW CODING FOR CONJUGATED EXPANSIONS
1183                      PE[h][6] = insert(PE[h][6], cod);
1184                    } else
1185                    {
1186                      PE[h][1] = (PE[h][1] + erg[sizeErg]);
1187                      if(newExt == 1)
1188                      {
1189                        PE[h][7] = insert(PE[h][7], 0);
1190                      }
1191                      PE[h][7] = sum2All(PE[h][7], slN * PE[h][2]);
1192                      PE[h][2] = PE[h][2] * slD;
1193
1194                      // NEW CODING FOR CONJUGATED EXPANSIONS
1195                      PE[h][6] = insert(PE[h][6], cod);
1196
1197                      PE[h][8] = list("Denominator", var(1) ^ (-slN*PE[h][2]));
1198                    }
1199                  }
1200                  setring R;
1201                  cs[size(cs) + 1] = RT;
1202                  kill RT;
1203                  setring S;
1204                } else {
1205                  if(slN >= 0)
1206                  {
1207                    PES[size(PES) + 1] = list((csTS[k][1] + newA)*var(1)^(csTS[k][2]*slN), csTS[k][2]*slD);
1208                  } else {
1209                    PES[size(PES) + 1] = list(csTS[k][1] + newA, csTS[k][2]*slD);
1210                    PES[size(PES)][8] = list("Denominator", var(1) ^ (-slN*csTS[k][2]));
1211                  }
1212
1213                  // NEW CODING FOR CONJUGATED EXPANSIONS
1214                  PES[size(PES)][6] = insert(csTS[k][6], cod);
1215
1216                  if(newExt == 1)
1217                  {
1218                    if(typeof(csTS[k][7]) == "none")
1219                    {
1220                      PES[size(PES)][7] = list(0);
1221
1222                    } else
1223                    {
1224                      PES[size(PES)][7] = insert(csTS[k][7], 0);
1225                    }
1226                    PES[size(PES)][7] = sum2All(PES[size(PES)][7], slN * csTS[k][2]);
1227                  } else
1228                  {
1229                    PES[size(PES)][7] = csTS[k][7];
1230                  }
1231
1232                }
1233              }
1234              if(size(PES) > 0){
1235                list PE = PES;
1236                export PE;
1237                setring R;
1238                cs[size(cs) + 1] = S;
1239              }
1240            } else {
1241              poly PE1 = newA*var(1)^slN;
1242              list PE = list(list(PE1, slD));
1243              PE[1][6] = list(cod);
1244              if(newExt == 1)
1245              {
1246                PE[1][7] = list(slN);
1247              }
1248              export PE;
1249
1250              setring R;
1251              cs[size(cs) + 1] = S;
1252            }
1253            kill S;
1254          }
1255        }
1256        setring R;
1257      }
1258    }
1259  }
1260
1261  if(size(cs) == 0){
1262    cs = list(list(poly(0),1));
1263    cs[1][6] = list(1);
1264  }
1265
1266  return(cs);
1267}
1268
1269
1270///////////////////////////////////////////////////////////////////////////////
1271
1272static proc newtonpoly2(poly f, int #)
1273"USAGE:   newtonpoly2(f);   f poly
1274ASSUME:  basering has exactly two variables; @*
1275         f is convenient, that is, f(x,0) != 0 != f(0,y).
1276RETURN:  list of intvecs (= coordinates x,y of the Newton polygon of f).
1277NOTE:    Procedure uses @code{execute}; this can be avoided by calling
1278         @code{newtonpoly2(f,1)} if the ordering of the basering is @code{ls}.
1279KEYWORDS: Newton polygon
1280EXAMPLE: example newtonpoly2;  shows an example
1281"
1282{
1283  matrix m = coeffs(f, var(2));
1284  matrix my;
1285  list A;
1286  int i;
1287  int j = 1;
1288  for (i=1; i<=nrows(m); i++)
1289  {
1290    if(m[i,1]!=0)
1291    {
1292      my = coef(m[i,1],var(1));
1293      A[j] = intvec(deg(my[1, ncols(my)]), i-1);
1294      j++;
1295    }
1296  }
1297
1298  list l2;
1299  int rep = 1;
1300  while(rep == 1)
1301  {
1302    l2 = list();
1303    l2[1] = A[1];
1304    rep = 0;
1305    j = 2;
1306    for (i=2; i<=size(A)-1; i++)
1307    {
1308      if((A[i][2]-A[i-1][2])*(A[i+1][1]-A[i-1][1]) > (A[i+1][2]-A[i-1][2])*(A[i][1]-A[i-1][1]))
1309      {
1310        l2[j] = A[i];
1311        j++;
1312      } else {
1313        rep = 1;
1314      }
1315    }
1316    l2[j] = A[size(A)];
1317    A = l2;
1318  }
1319//  "newton polygon: ";
1320//  f;
1321//  A;
1322  return(A);
1323}
1324example
1325{
1326 "EXAMPLE:"; echo = 2;
1327  ring r=0,(x,y),ls;
1328  poly f=x5+2x3y-x2y2+3xy5+y6-y7;
1329  newtonpoly2(f);
1330}
1331
1332///////////////////////////////////////////////////////////////////////////////
1333
1334static proc negExp(poly fNew1, poly c1, int slN)
1335{
1336  // All this is a workaround to work with negative exponents,
1337  // which are not allowed in Singular.
1338
1339  //" DIVIDE EXPANSION BY X^", slN;
1340  int dy = deg(fNew1, intvec(0,1));
1341  int dx3;
1342  int dy3;
1343  int dxExpT;
1344
1345  // We multiply by a large enough power of x so that the next
1346  // substitution is possible.
1347  poly f3 = fNew1*(var(1)^(-slN*dy));
1348
1349  //We replace var(2) by var(1)^slN*(c1 + var(2));
1350  poly f3Temp = f3;
1351  poly f4;
1352  poly leadP;
1353
1354  poly expT;
1355  poly expTTemp;
1356  poly leadExpT;
1357
1358  int i;
1359  int j;
1360  poly f4T;
1361  for(i = 1; i <= size(f3); i++)
1362  {
1363    leadP = lead(f3Temp);
1364    dx3 = deg(leadP, intvec(1,0));
1365    dy3 = deg(leadP, intvec(0,1));
1366    expT = (var(1)*(c1 + var(2)))^dy3;
1367    expTTemp = expT;
1368    f4T = 0;
1369    for(j = 1; j <= size(expT); j++)
1370    {
1371      leadExpT = lead(expTTemp);
1372      expTTemp = expTTemp - leadExpT;
1373      dxExpT = deg(leadExpT, intvec(1,0));
1374
1375      // IN k * x^a * leadExpT, substitute the x in leadExpT by x^(-slN)
1376      f4T = f4T + leadP/var(1)^(dxExpT*(-slN))/(var(2)^dy3) * subst(leadExpT, var(1), 1);
1377    }
1378    f4 = f4 + f4T;
1379    f3Temp = f3Temp - leadP;
1380  }
1381  f4 = sat(ideal(f4), var(1))[1];
1382  return(f4);
1383}
1384
1385
1386// Main procedure for computing Puiseux expansions returning a list of
1387// Input: - polynomial PP for which the expansions will be computed
1388//        - integer maxDeg, the maximum degree up to which compute the
1389//          expansions. If maxDeg = -1 computes the singular part
1390//        - int iVarX, the index of the X variable for the ring C{X}[Y].
1391//        - int iVarY, the index of the Y variable for the ring C{X}[Y].
1392proc puiseuxList(Puiseux PP, int maxDeg, int iVarX, int iVarY)
1393"USAGE:  puiseuxList(PP, maxDeg, iVarX, iVarY);   Puiseux expansion PP, int maxDeg
1394         is the integer up to which the Puiseux expansions will be computed
1395         (if maxDeg = -1 computes the singular part),
1396         int iVarX is the index of the X variable for the ring C{X}[Y],
1397         int iVarY, the index of the Y variable for the ring C{X}[Y].
1398ASSUME:  basering has exactly two variables; @*
1399         f is convenient, that is, f(x,0) != 0 != f(0,y).
1400RETURN:  a list with the Puiseux expansions of PP.
1401KEYWORDS: Puiseux expansions
1402EXAMPLE: example puiseuxList;  shows an example
1403"
1404
1405{
1406  def R = basering;
1407  def R_PP = PP.in;
1408  setring R_PP;
1409
1410  poly PP_numer = PP.numer;
1411
1412  // We move to a ring where the only variables are varX and varY
1413  list l = ringlist(R_PP);
1414  list lNew = l;
1415  lNew[2] = list(var(iVarX), var(iVarY));
1416  lNew[3] = list(list("dp", intvec(1,1)), list("C", 0), list("L", 100000));
1417  def S = ring(lNew);
1418  setring S;
1419
1420  int id = 1;
1421  list out;
1422  int i;
1423  int j;
1424
1425  // We compute the Puiseux expansions in this new ring
1426  poly f = imap(R_PP, PP_numer);
1427  list p;
1428
1429  // NOT WORKING - WE NEED TO CONSIDER ALSO THE MULTIPLICITIES
1430  // If we develop to a fix degree we can split the computations.
1431  // If we want to compute the singular parts then we cannot split.
1432  //if(maxDeg != -1)
1433  //{
1434  //  // We compute the Puiseux expansions in this new ring
1435  //  list fFacs = factorize(f);
1436  //  fFacs = Integralbasis::sortFactors(fFacs);
1437  //
1438  //  for(i = 1; i <= size(fFacs[1]); i++)
1439  //  {
1440  //    if((fFacs[1][i] != 1) and (deg(fFacs[1][i], intvec(0,1)) > 0))
1441  //    {
1442  //      p = p + puiseuxMain(fFacs[1][i], maxDeg, 0, 0, 1);
1443  //    }
1444  //  }
1445  //} else
1446  //{
1447  //  p = puiseuxMain(f, maxDeg, 0, 0, 1);
1448  //}
1449
1450  p = puiseuxMain(f, maxDeg, 0, 0, 1);
1451
1452  // We create a list of type Puiseux with one element for each Puiseux
1453  // expansion
1454  for(i = 1; i <= size(p); i++)
1455  {
1456    // Case of Puiseux expansions in an extended ring
1457    // (there can be different expansions in this ring)
1458    if(typeof(p[i])=="ring")
1459    {
1460      def R1 =  p[i];
1461      setring R1;
1462
1463      // We add back the removed variables
1464      list lR1 = ringlist(R1);
1465      list l = imap(R_PP, l);
1466      lR1[2] = l[2];
1467      lR1[3] = l[3];
1468
1469      def S1 = ring(lR1);
1470      setring S1;
1471      list PE = imap(R1, PE);
1472
1473      // We create the new elements of the list, one element
1474      // for each Puiseux expansion in PE.
1475      for(j = 1; j <= size(PE); j++)
1476      {
1477        poly pui = PE[j][1];
1478        Puiseux P = pui;
1479        P.in = basering;
1480        P.numer = pui;
1481
1482        // The denominator of the exponent gets multiplied by the
1483        // the new denominator of the exponent
1484        P.fraction = PE[j][2] * PP.fraction;
1485
1486        // The old denominator is discarded
1487        if(size(PE[j]) >= 8)
1488        {
1489          if(typeof(PE[j][8]) != "none")
1490          {
1491            P.denom = PE[j][8][2];
1492          } else
1493          {
1494            P.denom = poly(1);
1495          }
1496        } else
1497        {
1498          P.denom = poly(1);
1499        }
1500        out[id] = P;
1501        kill P;
1502        kill pui;
1503        id++;
1504      }
1505      setring S;
1506      kill R1;
1507      kill S1;
1508    } else
1509    {
1510      // Case of Puiseux expansions over the basering
1511
1512      // We move back to the original ring with all the variables
1513      setring R_PP;
1514      list p2 = imap(S, p);
1515
1516      // We create the new element of the list
1517      poly pui = p2[i][1];
1518      Puiseux P = pui;
1519      kill pui;
1520
1521      P.in = basering;
1522
1523      // The denominator of the exponent gets multiplied by the
1524      // the new denominator of the exponent
1525      P.fraction = p2[i][2] * PP.fraction;
1526
1527      // The old denominator is discarded
1528      if(size(p2[i]) >= 8)
1529      {
1530        P.denom = p2[i][8][2];
1531      } else
1532      {
1533        P.denom = poly(1);
1534      }
1535      out[id] = P;
1536      kill p2;
1537      kill P;
1538      id++;
1539      setring S;
1540    }
1541  }
1542  setring R_PP;
1543  kill PP_numer;
1544  kill l;
1545  kill lNew;
1546
1547  setring R;
1548  return(out);
1549}
1550example
1551{
1552 "EXAMPLE:"; echo = 2;
1553  ring r=0,(x,y),dp;
1554  Puiseux PP = x5+2x3y-x2y2+3xy5+y6+y7;
1555  puiseuxList(PP, 5, 1, 2);
1556}
1557
1558// The equation corresponding to a segment of the Newton polygon.
1559// Factors of this equation correspong to different groups of
1560// conjugate expansions
1561static proc minEqNewton(poly f, int slN, int slD)
1562{
1563  intvec v = slD, slN;
1564  int d = 0;
1565  poly je = jet(f, d, v);
1566  while(je == 0){
1567    d++;
1568    je = jet(f, d, v);
1569  }
1570  je = Integralbasis::divideBy(je, var(1))[1];
1571  je = Integralbasis::divideBy(je, var(2))[1];
1572  return(je);
1573}
1574
1575///////////////////////////////////////////////////////////////////////////////
1576
1577// Same as splitring, but it allows for a name of the parameter to add
1578// (it will call splitring and replace the parameter, so if a,b,c,o
1579// are in use, it will give an error as in splitring)
1580static proc splitringParname(poly f, string varName, list #)
1581{
1582  def R = basering;
1583  def S = Integralbasis::splitRingAt(f, #);
1584  //def S = splitring(f, #);
1585  setring S;
1586  int isErg = defined(erg);
1587  list l = ringlist(S);
1588  l[1][2][1] = varName;
1589  l[1][4][1] = 0;
1590  poly p =subst(minpoly, par(1), var(1));
1591  def T = ring(l);
1592  setring T;
1593  poly q = fetch(S, p);
1594  number mp = number(subst(q, var(1), par(1)));
1595  minpoly = mp;
1596  if(isErg) {
1597    list erg = fetch(S, erg);
1598    export(erg);
1599  }
1600  setring R;
1601  return(T);
1602}
1603///////////////////////////////////////////////////////////////////////////////
1604
1605// Sum a same number to all elements in the list of integers
1606static proc sum2All(list l, int k)
1607{
1608  for(int i = 1; i <= size(l); i++)
1609  {
1610    l[i] = l[i] + k;
1611  }
1612  return(l);
1613}
1614///////////////////////////////////////////////////////////////////////////////
1615
1616
1617
1618// Checks if the coefficents of a polynomial are rational or contain some parameter
1619static proc nonRatCoeff(poly p)
1620{
1621  matrix cc = coef(p, var(1)*var(2));
1622  int nonRat = 0;
1623  for(int i = 1; i <= ncols(cc); i++)
1624  {
1625    if(pardeg(number(cc[2,i])) > 0)
1626    {
1627      nonRat = 1;
1628    }
1629  }
1630  return(nonRat);
1631}
Note: See TracBrowser for help on using the repository browser.