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

spielwiese
Last change on this file since a96a75 was a96a75, checked in by slap <slaplagne@…>, 2 years ago
Bugs fixed 1. Puiseux expansions were erroneously marked as belonging to different classes. 2. The maximum degree of the expansions is incremented when in one computing Hensel factors to avoid erroneous truncations.
  • Property mode set to 100644
File size: 30.9 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2///////////////////////////////////////////////////////////////////////////////
3///////////////////////////////////////////////////////////////////////////////
4version="version puiseuxexpansion.lib 4.2.1.3 Dec_2021 "; // $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 = puiseuxMain(f, maxDeg, 0, 1, 1);
275  } else
276  {
277    list p = puiseuxMain(f, maxDeg, 0, 0, 1);
278  }
279
280  return(p);
281}
282example
283{
284 "EXAMPLE:"; echo=2;
285 ring R=0,(x,y),dp;
286 poly f=y^3 + x^2 + x^8;
287 puiseux(f,3,0);
288}
289
290
291// Finds the coefficient of y^slope in the Puiseux expansion
292// slN is the numerator of the slope
293// slD is the denominator
294static proc puiseuxStep(poly f, int slN, int slD){
295  def R = basering;
296  poly x = var(1);
297  poly y = var(2);
298
299  // We treat the negative exponents separately
300  poly f2 = subst(f, x, x^slD);
301  if(slN >= 0)
302  {
303    poly f3 = subst(f2, y, y*(x^slN));
304    matrix M = coef(f3,x);
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*x^(-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 / x^(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,x);
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            minP = fc[1][j];
499            if(deg(minP)==1)
500            {
501              cofs = coeffs(minP, var(2));
502              c1 = number(-cofs[1,1]/cofs[2,1]);
503              if(defined(erg)){
504                sizeErg = size(erg);
505              } else {
506                sizeErg = 1;
507              }
508              if(stop == 0)
509              {
510                if(slN >= 0)
511                {
512                  fNew2 = subst(fNew1, var(2), (c1+var(2))*(var(1)^slN));
513                } else
514                {
515                  // Negative exponent case
516                  fNew2 = negExp(fNew1, c1, slN);
517                }
518                fNew2= sat(ideal(fNew2), var(1))[1][1];
519
520                if(maxDeg > 0)
521                {
522                  newMaxDeg = maxDeg * slD - slN;
523                } else {
524                  newMaxDeg = maxDeg;
525                }
526
527                csT = puiseuxMain(fNew2, newMaxDeg, slN, 1, 0);
528
529                for(k = 1; k<=size(csT); k++)
530                {
531                  if(typeof(csT[k]) != "ring")
532                  {
533                    // Case of a polynomial with the corresponding denominator;
534                    if(slN >= 0)
535                    {
536                      cs[size(cs) + 1] = list((csT[k][1] + c1)*(var(1)^(csT[k][2]*slN)), csT[k][2]*slD);
537                    } else {
538                      cs[size(cs) + 1] = list(csT[k][1] + c1, csT[k][2]*slD);
539                      cs[size(cs)][8] = list("Denominator", var(1) ^ (-slN*csT[k][2]));
540                    }
541
542
543                    // NEW CODING FOR CONJUGATED EXPANSIONS
544                    cs[size(cs)][6] = insert(csT[k][6], cod);
545
546                    if(newExt == 1)
547                    {
548                      if(typeof(csT[k][7])=="none"){
549                        cs[size(cs)][7] = list(0);
550                      } else
551                      {
552                        cs[size(cs)][7] = insert(csT[k][7], 0);
553                      }
554                      cs[size(cs)][7] = sum2All(cs[size(cs)][7], slN * csT[k][2]);
555                    } else
556                    {
557                      if(typeof(csT[k][7]) != "none"){
558                        cs[size(cs)][7] = sum2All(csT[k][7], slN * csT[k][2]);
559                      }
560                    }
561                  } else {
562                    def RT = csT[k];
563                    setring RT;
564                    number c1 = fetch(R, c1);
565                    // Change @a by par(1)?
566                    c1 = number(subst(c1, @a, number(erg[sizeErg])));
567
568                    for(h = 1; h <= size(PE); h++)
569                    {
570                      if(slN >= 0)
571                      {
572                        PE[h][1] = (PE[h][1] + c1) * (var(1)^(slN*PE[h][2]));
573                        if(newExt == 1)
574                        {
575                          PE[h][7] = insert(PE[h][7], 0);
576                        }
577                        PE[h][7] = sum2All(PE[h][7], slN * PE[h][2]);
578                        PE[h][2] = PE[h][2] * slD;
579
580                        // NEW CODING FOR CONJUGATED EXPANSIONS
581                        PE[h][6] = insert(PE[h][6], cod);
582                      } else {
583                        PE[h][1] = PE[h][1] + c1;
584                        if(newExt == 1)
585                        {
586                          PE[h][7] = insert(PE[h][7], 0);
587                        }
588
589                        // Check if this is correct
590                        // (before PE[h][8] was defined after PE[h][2] was modified)
591                        PE[h][8] = list("Denominator", var(1) ^ (-slN*PE[h][2]));
592                        PE[h][2] = PE[h][2] * slD;
593
594                        // NEW CODING FOR CONJUGATED EXPANSIONS
595                        PE[h][6] = insert(PE[h][6], cod);
596                      }
597                    }
598                    setring R;
599                    cs[size(cs) + 1] = RT;
600                    kill RT;
601                  }
602                }
603              } else {
604                cs[size(cs) + 1] = list(c1 * var(1)^slN, slD);
605                cs[size(cs)][6] = list(cod);
606                if(newExt == 1)
607                {
608                  cs[size(cs)][7] = list(slN);
609                }
610              }
611            } else {
612              if(npars(R) == 1)
613              {
614                if(!defined(erg))
615                {
616                  list erg;
617                  erg[1] = par(1);
618                }
619                if(!defined(minPolys))
620                {
621                  list minPolys;
622                  minPolys[1] = 1;
623                }
624
625                def S = Integralbasis::splitRingAt(minP, erg);
626                setring S;
627                sizeErg = size(erg);
628                poly newA = erg[size(erg)];  // The root founded. That is, newA is the root of minP in S.
629                poly fNew1 = fetch(R, fNew1);
630
631                if(!defined(minPolys))
632                {
633                  list minPolys = fetch(R, minPolys);
634                }
635
636                // We replace the old a by its image in the new ring.
637                minPolys[size(minPolys)+1] = fetch(R, minP);
638                for(mi = 1; mi <= size(minPolys); mi++)
639                {
640                  minPolys[mi] = subst(minPolys[mi], par(1), erg[size(erg)-1]);
641                }
642                fNew1 = subst(fNew1, par(1), erg[size(erg)-1]);
643
644              } else {
645                def S = Integralbasis::splitRingAt(minP);
646                setring S;
647                list minPolys;
648                minPolys[1] = 1;
649                minPolys[size(minPolys)+1] = fetch(R, minP);
650                poly fNew1 = fetch(R, fNew1);
651                poly newA = par(1);      // In this case, it is the new variable.
652                list erg = par(1);
653                export erg;
654                sizeErg = 1;
655              }
656              export(minPolys);
657
658              if(defined(PES))
659              {
660                PES = list();
661              } else
662              {
663                list PES;
664              }
665              if(stop == 0){
666                // We use the image of the root of minP in the new ring.
667                if(slN >= 0)
668                {
669                  poly fNew2 = subst(fNew1, var(2), (newA + var(2))*var(1)^slN);
670                } else
671                {
672                  poly fNew2 = negExp(fNew1, newA, slN);
673                }
674                fNew2= sat(ideal(fNew2), var(1))[1][1];
675                if(maxDeg > 0)
676                {
677                  newMaxDeg = maxDeg * slD - slN;
678                } else {
679                  newMaxDeg = maxDeg;
680                }
681                list csTS = puiseuxMain(fNew2, newMaxDeg, slN, 1, 0);
682                for(k = 1; k<=size(csTS); k++)
683                {
684                  if(typeof(csTS[k]) == "ring")
685                  {
686                    def RT = csTS[k];
687                    setring RT;
688                    for(h = 1; h<=size(PE); h++)
689                    {
690                      // The new coefficient is the image of a in the new ring.
691                      if(slN >= 0)
692                      {
693                        PE[h][1] = (PE[h][1] + erg[sizeErg]) * var(1)^(PE[h][2]*slN);
694                        if(newExt == 1)
695                        {
696                          PE[h][7] = insert(PE[h][7], 0);
697                        }
698                        PE[h][7] = sum2All(PE[h][7], slN * PE[h][2]);
699                        PE[h][2] = PE[h][2] * slD;
700
701                        // NEW CODING FOR CONJUGATED EXPANSIONS
702                        PE[h][6] = insert(PE[h][6], cod);
703                      } else
704                      {
705                        PE[h][1] = (PE[h][1] + erg[sizeErg]);
706                        if(newExt == 1)
707                        {
708                          PE[h][7] = insert(PE[h][7], 0);
709                        }
710                        PE[h][7] = sum2All(PE[h][7], slN * PE[h][2]);
711                        PE[h][2] = PE[h][2] * slD;
712
713                        // NEW CODING FOR CONJUGATED EXPANSIONS
714                        PE[h][6] = insert(PE[h][6], cod);
715
716                        PE[h][8] = list("Denominator", var(1) ^ (-slN*PE[h][2]));
717                      }
718                    }
719                    setring R;
720                    cs[size(cs) + 1] = RT;
721                    kill RT;
722                    setring S;
723                  } else {
724                    if(slN >= 0)
725                    {
726                      PES[size(PES) + 1] = list((csTS[k][1] + newA)*var(1)^(csTS[k][2]*slN), csTS[k][2]*slD);
727                    } else {
728                      PES[size(PES) + 1] = list(csTS[k][1] + newA, csTS[k][2]*slD);
729                      PES[size(PES)][8] = list("Denominator", var(1) ^ (-slN*csTS[k][2]));
730                    }
731
732                    // NEW CODING FOR CONJUGATED EXPANSIONS
733                    PES[size(PES)][6] = insert(csTS[k][6], cod);
734
735                    if(newExt == 1)
736                    {
737                      if(typeof(csTS[k][7]) == "none")
738                      {
739                        PES[size(PES)][7] = list(0);
740
741                      } else
742                      {
743                        PES[size(PES)][7] = insert(csTS[k][7], 0);
744                      }
745                      PES[size(PES)][7] = sum2All(PES[size(PES)][7], slN * csTS[k][2]);
746                    } else
747                    {
748                      PES[size(PES)][7] = csTS[k][7];
749                    }
750
751                  }
752                }
753                if(size(PES) > 0){
754                  list PE = PES;
755                  export PE;
756                  setring R;
757                  cs[size(cs) + 1] = S;
758                }
759              } else {
760                poly PE1 = newA*var(1)^slN;
761                list PE = list(list(PE1, slD));
762                PE[1][6] = list(cod);
763                if(newExt == 1)
764                {
765                  PE[1][7] = list(slN);
766                }
767                export PE;
768
769                setring R;
770                cs[size(cs) + 1] = S;
771              }
772              kill S;
773            }
774          }
775          setring R;
776        }
777      }
778    }
779  }
780
781  if(size(cs) == 0){
782    cs = list(list(poly(0),1));
783    cs[1][6] = list(1);
784  }
785
786  return(cs);
787}
788
789
790///////////////////////////////////////////////////////////////////////////////
791
792static proc newtonpoly2(poly f, int #)
793"USAGE:   newtonpoly2(f);   f poly
794ASSUME:  basering has exactly two variables; @*
795         f is convenient, that is, f(x,0) != 0 != f(0,y).
796RETURN:  list of intvecs (= coordinates x,y of the Newton polygon of f).
797NOTE:    Procedure uses @code{execute}; this can be avoided by calling
798         @code{newtonpoly2(f,1)} if the ordering of the basering is @code{ls}.
799KEYWORDS: Newton polygon
800EXAMPLE: example newtonpoly2;  shows an example
801"
802{
803  matrix m = coeffs(f, var(2));
804  matrix my;
805  list A;
806  int i;
807  int j = 1;
808  for (i=1; i<=nrows(m); i++)
809  {
810    if(m[i,1]!=0)
811    {
812      my = coef(m[i,1],var(1));
813      A[j] = intvec(deg(my[1, ncols(my)]), i-1);
814      j++;
815    }
816  }
817
818  list l2;
819  int rep = 1;
820  while(rep == 1)
821  {
822    l2 = list();
823    l2[1] = A[1];
824    rep = 0;
825    j = 2;
826    for (i=2; i<=size(A)-1; i++)
827    {
828      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]))
829      {
830        l2[j] = A[i];
831        j++;
832      } else {
833        rep = 1;
834      }
835    }
836    l2[j] = A[size(A)];
837    A = l2;
838  }
839//  "newton polygon: ";
840//  f;
841//  A;
842  return(A);
843}
844example
845{
846 "EXAMPLE:"; echo = 2;
847  ring r=0,(x,y),ls;
848  poly f=x5+2x3y-x2y2+3xy5+y6-y7;
849  newtonpoly2(f);
850}
851
852///////////////////////////////////////////////////////////////////////////////
853
854static proc negExp(poly fNew1, poly c1, int slN)
855{
856  // All this is a workaround to work with negative exponents,
857  // which are not allowed in Singular.
858
859  //" DIVIDE EXPANSION BY X^", slN;
860  int dy = deg(fNew1, intvec(0,1));
861  int dx3;
862  int dy3;
863  int dxExpT;
864
865  // We multiply by a large enough power of x so that the next
866  // substitution is possible.
867  poly f3 = fNew1*(var(1)^(-slN*dy));
868
869  //We replace var(2) by var(1)^slN*(c1 + var(2));
870  poly f3Temp = f3;
871  poly f4;
872  poly leadP;
873
874  poly expT;
875  poly expTTemp;
876  poly leadExpT;
877
878  int i;
879  int j;
880  poly f4T;
881  for(i = 1; i <= size(f3); i++)
882  {
883    leadP = lead(f3Temp);
884    dx3 = deg(leadP, intvec(1,0));
885    dy3 = deg(leadP, intvec(0,1));
886    expT = (var(1)*(c1 + var(2)))^dy3;
887    expTTemp = expT;
888    f4T = 0;
889    for(j = 1; j <= size(expT); j++)
890    {
891      leadExpT = lead(expTTemp);
892      expTTemp = expTTemp - leadExpT;
893      dxExpT = deg(leadExpT, intvec(1,0));
894
895      // IN k * x^a * leadExpT, substitute the x in leadExpT by x^(-slN)
896      f4T = f4T + leadP/var(1)^(dxExpT*(-slN))/(var(2)^dy3) * subst(leadExpT, var(1), 1);
897    }
898    f4 = f4 + f4T;
899    f3Temp = f3Temp - leadP;
900  }
901  f4 = sat(ideal(f4), var(1))[1][1];
902  return(f4);
903}
904
905
906// Main procedure for computing Puiseux expansions returning a list of
907// Input: - polynomial PP for which the expansions will be computed
908//        - integer maxDeg, the maximum degree up to which compute the
909//          expansions. If maxDeg = -1 computes the singular part
910//        - int iVarX, the index of the X variable for the ring C{X}[Y].
911//        - int iVarY, the index of the Y variable for the ring C{X}[Y].
912proc puiseuxList(Puiseux PP, int maxDeg, int iVarX, int iVarY)
913"USAGE:  puiseuxList(PP, maxDeg, iVarX, iVarY);   Puiseux expansion PP, int maxDeg
914         is the integer up to which the Puiseux expansions will be computed
915         (if maxDeg = -1 computes the singular part),
916         int iVarX is the index of the X variable for the ring C{X}[Y],
917         int iVarY, the index of the Y variable for the ring C{X}[Y].
918ASSUME:  basering has exactly two variables; @*
919         f is convenient, that is, f(x,0) != 0 != f(0,y).
920RETURN:  a list with the Puiseux expansions of PP.
921KEYWORDS: Puiseux expansions
922EXAMPLE: example puiseuxList;  shows an example
923"
924
925{
926  def R = basering;
927  def R_PP = PP.in;
928  setring R_PP;
929
930  poly PP_numer = PP.numer;
931
932  // We move to a ring where the only variables are varX and varY
933  list l = ringlist(R_PP);
934  list lNew = l;
935  lNew[2] = list(var(iVarX), var(iVarY));
936  lNew[3] = list(list("dp", intvec(1,1)), list("C", 0), list("L", 100000));
937  def S = ring(lNew);
938  setring S;
939
940  int id = 1;
941  list out;
942  int i;
943  int j;
944
945  // We compute the Puiseux expansions in this new ring
946  poly f = imap(R_PP, PP_numer);
947  list p;
948
949  // NOT WORKING - WE NEED TO CONSIDER ALSO THE MULTIPLICITIES
950  // If we develop to a fix degree we can split the computations.
951  // If we want to compute the singular parts then we cannot split.
952  //if(maxDeg != -1)
953  //{
954  //  // We compute the Puiseux expansions in this new ring
955  //  list fFacs = factorize(f);
956  //  fFacs = Integralbasis::sortFactors(fFacs);
957  //
958  //  for(i = 1; i <= size(fFacs[1]); i++)
959  //  {
960  //    if((fFacs[1][i] != 1) and (deg(fFacs[1][i], intvec(0,1)) > 0))
961  //    {
962  //      p = p + puiseuxMain(fFacs[1][i], maxDeg, 0, 0, 1);
963  //    }
964  //  }
965  //} else
966  //{
967  //  p = puiseuxMain(f, maxDeg, 0, 0, 1);
968  //}
969
970  p = puiseuxMain(f, maxDeg, 0, 0, 1);
971
972  // We create a list of type Puiseux with one element for each Puiseux
973  // expansion
974  for(i = 1; i <= size(p); i++)
975  {
976    // Case of Puiseux expansions in an extended ring
977    // (there can be different expansions in this ring)
978    if(typeof(p[i])=="ring")
979    {
980      def R1 =  p[i];
981      setring R1;
982
983      // We add back the removed variables
984      list lR1 = ringlist(R1);
985      list l = imap(R_PP, l);
986      lR1[2] = l[2];
987      lR1[3] = l[3];
988
989      def S1 = ring(lR1);
990      setring S1;
991      list PE = imap(R1, PE);
992
993      // We create the new elements of the list, one element
994      // for each Puiseux expansion in PE.
995      for(j = 1; j <= size(PE); j++)
996      {
997        poly pui = PE[j][1];
998        Puiseux P = pui;
999        P.in = basering;
1000        P.numer = pui;
1001
1002        // The denominator of the exponent gets multiplied by the
1003        // the new denominator of the exponent
1004        P.fraction = PE[j][2] * PP.fraction;
1005
1006        // The old denominator is discarded
1007        if(size(PE[j]) >= 8)
1008        {
1009          if(typeof(PE[j][8]) != "none")
1010          {
1011            P.denom = PE[j][8][2];
1012          } else
1013          {
1014            P.denom = poly(1);
1015          }
1016        } else
1017        {
1018          P.denom = poly(1);
1019        }
1020        out[id] = P;
1021        kill P;
1022        kill pui;
1023        id++;
1024      }
1025      setring S;
1026      kill R1;
1027      kill S1;
1028    } else
1029    {
1030      // Case of Puiseux expansions over the basering
1031
1032      // We move back to the original ring with all the variables
1033      setring R_PP;
1034      list p2 = imap(S, p);
1035
1036      // We create the new element of the list
1037      poly pui = p2[i][1];
1038      Puiseux P = pui;
1039      kill pui;
1040
1041      P.in = basering;
1042
1043      // The denominator of the exponent gets multiplied by the
1044      // the new denominator of the exponent
1045      P.fraction = p2[i][2] * PP.fraction;
1046
1047      // The old denominator is discarded
1048      if(size(p2[i]) >= 8)
1049      {
1050        P.denom = p2[i][8][2];
1051      } else
1052      {
1053        P.denom = poly(1);
1054      }
1055      out[id] = P;
1056      kill p2;
1057      kill P;
1058      id++;
1059      setring S;
1060    }
1061  }
1062  setring R_PP;
1063  kill PP_numer;
1064  kill l;
1065  kill lNew;
1066
1067  setring R;
1068  return(out);
1069}
1070example
1071{
1072 "EXAMPLE:"; echo = 2;
1073  ring r=0,(x,y),dp;
1074  Puiseux PP = x5+2x3y-x2y2+3xy5+y6+y7;
1075  puiseuxList(PP, 5, 1, 2);
1076}
1077
1078// The equation corresponding to a segment of the Newton polygon.
1079// Factors of this equation correspong to different groups of
1080// conjugate expansions
1081static proc minEqNewton(poly f, int slN, int slD)
1082{
1083  intvec v = slD, slN;
1084  int d = 0;
1085  poly je = jet(f, d, v);
1086  while(je == 0){
1087    d++;
1088    je = jet(f, d, v);
1089  }
1090  je = Integralbasis::divideBy(je, var(1))[1];
1091  je = Integralbasis::divideBy(je, var(2))[1];
1092  return(je);
1093}
1094
1095///////////////////////////////////////////////////////////////////////////////
1096
1097// Same as splitring, but it allows for a name of the parameter to add
1098// (it will call splitring and replace the parameter, so if a,b,c,o
1099// are in use, it will give an error as in splitring)
1100static proc splitringParname(poly f, string varName, list #)
1101{
1102  def R = basering;
1103  def S = Integralbasis::splitRingAt(f, #);
1104  //def S = splitring(f, #);
1105  setring S;
1106  int isErg = defined(erg);
1107  list l = ringlist(S);
1108  l[1][2][1] = varName;
1109  l[1][4][1] = 0;
1110  poly p =subst(minpoly, par(1), var(1));
1111  def T = ring(l);
1112  setring T;
1113  poly q = fetch(S, p);
1114  number mp = number(subst(q, var(1), par(1)));
1115  minpoly = mp;
1116  if(isErg) {
1117    list erg = fetch(S, erg);
1118    export(erg);
1119  }
1120  setring R;
1121  return(T);
1122}
1123///////////////////////////////////////////////////////////////////////////////
1124
1125// Sum a same number to all elements in the list of integers
1126static proc sum2All(list l, int k)
1127{
1128  for(int i = 1; i <= size(l); i++)
1129  {
1130    l[i] = l[i] + k;
1131  }
1132  return(l);
1133}
1134///////////////////////////////////////////////////////////////////////////////
1135
1136
1137
1138// Checks if the coefficents of a polynomial are rational or contain some parameter
1139static proc nonRatCoeff(poly p)
1140{
1141  matrix cc = coef(p, var(1)*var(2));
1142  int nonRat = 0;
1143  for(int i = 1; i <= ncols(cc); i++)
1144  {
1145    if(pardeg(number(cc[2,i])) > 0)
1146    {
1147      nonRat = 1;
1148    }
1149  }
1150  return(nonRat);
1151}
Note: See TracBrowser for help on using the repository browser.