source:git/Singular/LIB/absfact.lib@121d2ae

jengelh-datetimespielwiese
Last change on this file since 121d2ae was 121d2ae, checked in by Martin Lee <martinlee84@…>, 10 years ago
chg: added new system command "absBiFact"
• Property mode set to `100644`
File size: 24.1 KB
Line
1//
2version="\$Id\$";
3category="Factorization";
4info="
5LIBRARY: absfact.lib   Absolute factorization for characteristic 0
6AUTHORS: Wolfram Decker,       decker at math.uni-sb.de
7         Gregoire Lecerf,      lecerf at math.uvsq.fr
8         Gerhard Pfister,      pfister at mathematik.uni-kl.de
9
10OVERVIEW:
11A library for computing the absolute factorization of multivariate
12polynomials f with coefficients in a field K of characteristic zero.
13Using Trager's idea, the implemented algorithm computes an absolutely
14irreducible factor by factorizing over some finite extension field L
15(which is chosen such that V(f) has a smooth point with coordinates in L).
16Then a minimal extension field is determined making use of the
17Rothstein-Trager partial fraction decomposition algorithm.
18
19REFERENCES:
20G. Cheze, G. Lecerf: Lifting and recombination techniques for absolute
21                  factorization. Journal of Complexity, 23(3):380-420, 2007.
22
23KEYWORDS: factorization; absolute factorization.
25
26PROCEDURES:
27  absFactorize();        absolute factorization of poly
28  absBiFactorize();      absolute factorization of bivariate poly
29";
30
31////////////////////////////////////////////////////////////////////
32static proc partialDegree(poly p, int i)
33"USAGE:  partialDegree(p,i);   p poly, i int
34RETURN: int, the degree of p in the i-th variable
35"
36{
37  int n = nvars(basering);
38  intvec tmp;
39  tmp[n] = 0;
40  tmp[i] = 1;
41  return(deg(p,tmp));
42}
43////////////////////////////////////////////////////////////////////
44static proc belongTo(string s, list l)
45"USAGE:  belongTo(s,l);   s string, l list
46RETURN: 1 if s belongs to l, 0 otherwise
47"
48{
49  string tmp;
50  for(int i = 1; i <= size(l); i++) {
51    tmp = l[i];
52    if (tmp == s) {
53      return(1);
54    }
55  }
56  return(0);
57}
58////////////////////////////////////////////////////////////////////
59static proc variableWithSmallestPositiveDegree(poly p)
60"USAGE:  variableWithSmallestPositiveDegree(p);  p poly
61RETURN: int;  0 if p is constant. Otherwise, the index of the
62        variable which has the smallest positive degree in p.
63"
64{
65  int n = nvars(basering);
66  int v = 0;
67  int d = deg(p);
68  int d_loc;
69  for(int i = 1; i <= n; i++) {
70    d_loc = partialDegree(p, i);
71    if (d_loc >= 1 and d_loc <= d) {
72      v = i;
73      d = d_loc;
74    }
75  }
76  return(v);
77}
78////////////////////////////////////////////////////////////////////
79static proc smallestProperSimpleFactor(poly p)
80"USAGE:  smallestProperSimpleFactor(p);   p poly
81RETURN: poly:  a proper irreducible simple factor of p of smallest
82        degree. If no such factor exists, 0 is returned.
83"
84{
85  list p_facts = factorize(p);
86  int s = size(p_facts[1]);
87  int d = deg(p)+1;
88  poly q = 0;
89  poly f;
90  int e;
91  for(int i = 1; i <= s; i++)
92  {
93    f = p_facts[1][i];
94    e = deg(f);
95    if (e >= 1 and e < d and p_facts[2][i] == 1)
96    {
97      q = f / leadcoef(f);
98      d = e;
99    }
100  }
101  return(q);
102}
103////////////////////////////////////////////////////////////////////
104static proc smallestProperFactor(poly p)
105"USAGE:  smallestProperFactor(p);   p poly
106RETURN: poly:  a proper irreducible factor of p of smallest degree.
107        If p is constant, 0 is returned.
108"
109{
110  list p_facts = factorize(p);
111  int s = size(p_facts[1]);
112  int d = deg(p)+1;
113  poly q = 0;
114  poly f;
115  int e;
116  for(int i = 1; i <= s; i++)
117  {
118    f = p_facts[1][i];
119    e = deg(f);
120    if (e >= 1 and e < d)
121    {
122      q = f / leadcoef(f);
123      d = e;
124    }
125  }
126  return(q);
127}
128////////////////////////////////////////////////////////////////////
129static proc extensionContainingSmoothPoint(poly p, int m)
130"USAGE:  extensionContainingSmoothPoint(p,m);  p poly, m int
131RETURN: poly:  an irreducible univariate polynomial that defines an
132        algebraic extension of the current ground field that contains
133        a smooth point of the hypersurface defined by p=0.
134"
135{
136  int n = nvars(basering) - 1;
137  poly q = 0;
138  int i;
139  list a;
140  for(i=1;i<=n+1;i++){a[i] = 0;}
141  a[m] = var(n+1);
142  // The list a is to be taken with random entries in [-e, e].
143  // Every 10 * n trial, e is incremented by 1.
144  int e = 1;
145  int nbtrial = 0;
146  map h;
147  while (q == 0)
148  {
149    h = basering, a[1..n+1];
150    q = smallestProperSimpleFactor(h(p));
151    for(i = 1; i  <= n ; i = i + 1)
152    {
153      if (i != m)
154      {
155        a[i] = random(-e, e);
156      }
157    }
158    nbtrial++;
159    if (nbtrial >= 10 * n)
160    {
161      e = e + 1;
162      nbtrial = 0;
163    }
164  }
165  return(q);
166}
167////////////////////////////////////////////////////////////////////
168static proc RothsteinTragerResultant(poly g, poly f, int m)
169"USAGE:  RothsteinTragerResultant(g,f,m);  g,f poly, m int
170RETURN: poly
171NOTE:   To be called by the RothsteinTrager procedure only.
172"
173{
174  def MPz = basering;
175  int n = nvars(MPz) - 1;
176  int d = partialDegree(f, m);
177  poly df = diff(f, var(m));
178  list a;
179  int i;
180  for(i=1;i<=n+1;i++){ a[i] = 0; }
181  a[m] = var(m);
182  poly q = 0;
183  int e = 1;
184  int nbtrial = 0;
185  map h;
186  while (q == 0)
187  {
188    h = MPz, a[1..n+1];
189    q = resultant(h(f), h(df) * var(n+1) - h(g), var(m));
190    if (deg(q) == d)
191    {
193    }
194    q = 0;
195    for(i = 1; i  <= n ; i++)
196    {
197      if (i != m)
198      {
199        a[i] = random(-e, e);
200      }
201    }
202    nbtrial++;
203    if (nbtrial >= 10 * n)
204    {
205      e++;
206      nbtrial = 0;
207    }
208  }
209}
210////////////////////////////////////////////////////////////////////
211static proc RothsteinTrager(list g, poly p, int m, int expectedDegQ)
212"USAGE:  RothsteinTrager(g,p,m,d);  g list, p poly, m,d int
213RETURN: list L consisting of two entries of type poly
214NOTE:   the return value is the Rothstein-Trager partial fraction
215        decomposition of the quotient s/p, where s is a generic linear
216        combination of the elements of g. The genericity via d
217        (the expected degree of L[1]).
218"
219{
220  def MPz = basering;
221  int n = nvars(MPz) - 1;
222  poly dp = diff(p, var(m));
223  int r = size(g);
224  list a;
225  int i;
226  for(i=1;i<=r;i++){a[i] = 0;}
227  a[r] = 1;
228  int nbtrial = 0;
229  int e = 1;
230  poly s;
231  poly q;
232  while (1)
233  {
234    s = 0;
235    for(i = 1; i <= r; i++){s = s + a[i] * g[i];}
236    q = RothsteinTragerResultant(s, p, m);
237    q = smallestProperFactor(q);
238    if (deg(q) == expectedDegQ)
239    {
240      // Go into the quotient by q(z)=0
241      ring MP_z = (0,var(n+1)), (x(1..n)), dp;
242      list lMP_z = ringlist(MP_z);
243      lMP_z[1][4] = ideal(imap(MPz,q));
244      list tmp = ringlist(MPz)[2];
245      lMP_z[2] = list(tmp[1..n]);
246      def MPq = ring(lMP_z);
247      setring(MPq);
248      poly f = gcd(imap(MPz, p), par(1) * imap(MPz, dp) - imap(MPz, s));
249      f = f / leadcoef(f);
250      setring(MPz);
251      return(list(q, imap(MPq, f)));
252    }
253    for(i = 1; i  <= r ; i++)
254    {
255      a[i] = random(-e, e);
256    }
257    nbtrial++;
258    if (nbtrial >= 10 * r)
259    {
260      e++;
261      nbtrial = 0;
262    }
263  }
264}
265////////////////////////////////////////////////////////////////////
266static proc absFactorizeIrreducible(poly p)
267"USAGE:   absFactorizeIrreducible(p);   p poly
268ASSUME:  p is an irreducible polynomial that does not depend on the last
269         variable @z of the basering.
270RETURN:  list L of two polynomials: q=L[1] is an irreducible polynomial of
271         minimal degree in @z such that p has an absolute factor
272         over K[@z]/<q>, and f represents such an absolute factor.
273"
274{
275  int dblevel = printlevel - voice + 2;
276  dbprint(dblevel,"Entering absfact.lib::absFactorizeIrreducible with ",p);
277  def MPz = basering;
278  int d = deg(p);
279  int n = nvars(MPz) - 1;
280
281  if (d < 1)
282  {
283    return(list(var(n+1), p));
284  }
285  int m = variableWithSmallestPositiveDegree(p);
286  // var(m) is now considered as the main variable.
287
288  poly q = extensionContainingSmoothPoint(p, m);
289  int r = deg(q);
290  if (r == 1)
291  {
292    return(list(var(n+1), p));
293  }
294
295  list tmp = ringlist(MPz)[2];
296  // Go into the quotient by q(z)=0
297  ring MP_z = (0,var(n+1)), (x(1..n)), dp;
298  list lMP_z = ringlist(MP_z);
299  lMP_z[1][4] = ideal(imap(MPz,q));
300  lMP_z[2] = list(tmp[1..n]);
301  def MPq = ring(lMP_z);
302  setring(MPq);
303  dbprint(dblevel-1,"Factoring in algebraic extension");
304  // "Factoring p in the algebraic extension...";
305  poly p_loc = imap(MPz, p);
306  poly f = smallestProperSimpleFactor(p_loc);
307  int degf = deg(f);
308
309  if (degf == d)
310  {
311    setring(MPz);
312    return(list(var(n+1), p));
313  }
314
315  if (degf * r == d)
316  {
317    setring(MPz);
318    return(list(q, imap(MPq, f)));
319  }
320  dbprint(dblevel-1,"Absolutely irreducible factor found");
321  dbprint(dblevel,"Minimizing field extension");
322  // "Need to find a minimal extension";
323  poly co_f = p_loc / f;
324  poly e = diff(f, var(m)) * co_f;
325  setring(MPz);
326  poly e = imap(MPq, e);
327  list g;
328  int i;
329  for(i = 1; i <= r; i++)
330  {
331    g[i] = subst(e, var(n+1), 0);
332    e = diff(e, var(n+1));
333  }
334  return(RothsteinTrager(g, p, m, d div degf));
335}
336
337////////////////////////////////////////////////////////////////////
338proc absFactorize(poly p, list #)
339"USAGE:  absFactorize(p [,s]);   p poly, s string
340ASSUME: coefficient field is the field of rational numbers or a
341        transcendental extension thereof
342RETURN: ring @code{R} which is obtained from the current basering
343        by adding a new parameter (if a string @code{s} is given as a
344        second input, the new parameter gets this string as name). The ring
345        @code{R} comes with a list @code{absolute_factors} with the
346        following entries:
347@format
348    absolute_factors[1]: ideal   (the absolute factors)
349    absolute_factors[2]: intvec  (the multiplicities)
350    absolute_factors[3]: ideal   (the minimal polynomials)
351    absolute_factors[4]: int     (total number of nontriv. absolute factors)
352@end format
353        The entry @code{absolute_factors[1][1]} is a constant, the
354        entry @code{absolute_factors[3][1]} is the parameter added to the
355        current ring.@*
356        Each of the remaining entries @code{absolute_factors[1][j]} stands for
357        a class of conjugated absolute factors. The corresponding entry
358        @code{absolute_factors[3][j]} is the minimal polynomial of the
359        field extension over which the factor is minimally defined (its degree
360        is the number of conjugates in the class). If the entry
361        @code{absolute_factors[3][j]} coincides with @code{absolute_factors[3][1]},
362        no field extension was necessary for the @code{j}th (class of)
363        absolute factor(s).
364NOTE:   All factors are presented denominator- and content-free. The constant
365        factor (first entry) is chosen such that the product of all (!) the
366        (denominator- and content-free) absolute factors of @code{p} equals
367        @code{p / absolute_factors[1][1]}.
369EXAMPLE: example absFactorize; shows an example
370"
371{
372  int dblevel = printlevel - voice + 2;
373  dbprint(dblevel,"Entering absfact.lib::absFactorize with ",p);
374  def MP = basering;
375  int i;
376  if (char(MP) != 0)
377  {
378    ERROR("// absfact.lib::absFactorize is only implemented for "+
379          "characteristic 0");
380  }
381  if(minpoly!=0)
382  {
383    ERROR("// absfact.lib::absFactorize is not implemented for algebraic "
384          +"extensions");
385  }
386
387  int n = nvars(MP);
388  int pa=npars(MP);
389  list lMP= ringlist(MP);
390  list buflMP= lMP;
391  intvec vv,vk;
392  for(i=1;i<=n;i++){vv[i]=1;}
393  vk=vv,1;
394
395  //if the basering has parameters, add the parameters to the variables
396  //takes care about coefficients and possible denominators
397  if(pa>0)
398  {
399    poly qh=cleardenom(p);
400    if (p==0)
401    {
402      number cok=0;
403    }
404    else
405    {
407    }
408    p=qh;
409    string sp;
410    for(i=1;i<=npars(basering);i++)
411    {
412      sp=string(par(i));
413      sp=sp[2..size(sp)-1];
414      lMP[2][n+i]=sp;
415      vv=vv,1;
416    }
417    lMP[1]=0;
418    n=n+npars(MP);
419  }
420
421  // MPz is obtained by adding the new variable @z to MP
422  // ordering is wp(1...1)
423  // All the above subroutines work in MPz
424  string newvar;
425  if(size(#)>0)
426  {
427    if(typeof(#[1])=="string")
428    {
429      newvar=#[1];
430    }
431    else
432    {
433      newvar = "a";
434    }
435  }
436  else
437  {
438    newvar = "a";
439  }
440  if (newvar=="a")
441  {
442    if(belongTo(newvar, lMP[2])||defined(a)){newvar = "b";}
443    if(belongTo(newvar, lMP[2])||defined(b)){newvar = "c";}
444    if(belongTo(newvar, lMP[2])||defined(c)){newvar = "@c";}
445    while(belongTo(newvar, lMP[2]))
446    {
447       newvar = "@" + newvar;
448    }
449  }
450  lMP[2][n+1] = newvar;
451
452  // new ordering
453  vv=vv,1;
454  list orst;
455  orst[1]=list("wp",vv);
456  orst[2]=list("C",0);
457  lMP[3]=orst;
458
459  def MPz = ring(lMP);
460  setring(MPz);
461  poly p=imap(MP,p);
462
463  // special treatment in the homogeneous case, dropping one variable
464  // by substituting the first variable by 1
465  int ho=homog(p);
466  if(ho)
467  {
468    int dh=deg(p);
469    p=subst(p,var(1),1);
470    int di=deg(p);
471  }
472  list rat_facts = factorize(p);
473  int s = size(rat_facts[1]);
474  list tmpf; // absolute factors
475  intvec tmpm; // respective multiplicities
477  tmpm[1] = 1;
478  poly tmp;
479  for(i = 2; i <= s; i++)
480  {
481    tmp = rat_facts[1][i];
482    tmp = tmp / leadcoef(tmp);
483    tmpf[i] = absFactorizeIrreducible(tmp);
484    tmpm[i] = rat_facts[2][i];
485  }
486  // the homogeneous case, homogenizing the result
487  // the new variable has to have degree 0
488  // need to change the ring
489  if(ho)
490  {
491    list ll=ringlist(MPz);
492    vv[size(vv)]=0;
493    ll[3][1][2]=vv;
494    def MPhelp=ring(ll);
495    setring(MPhelp);
496    list tmpf=imap(MPz,tmpf);
497    for(i=2;i<=size(tmpf);i++)
498    {
499      tmpf[i][2]=homog(tmpf[i][2],var(1));
500    }
501    if(dh>di)
502    {
503      tmpf[size(tmpf)+1]=list(var(n+1),var(1));
504      tmpm[size(tmpm)+1]=dh-di;
505    }
506    setring(MPz);
507    tmpf=imap(MPhelp,tmpf);
508  }
509  // in case of parameters we have to go back to the old ring
510  // taking care about constant factors
511  if(pa)
512  {
513    setring(MP);
514    n=nvars(MP);
515    list lM=ringlist(MP);
516    orst[1]=list("wp",vk);
517    orst[2]=list("C",0);
518    lM[2][n+1] = newvar;
519    lM[3]=orst;
520    def MPout=ring(lM);
521    setring(MPout);
522    list tmpf=imap(MPz,tmpf);
523    number cok=imap(MP,cok);
524    tmpf[1][2]=cok*tmpf[1][2];
525  }
526  else
527  {
528    def MPout=MPz;
529  }
530  // if we want the output as string
531  if(size(#)>0)
532  {
533    if(typeof(#[1])=="int")
534    {
535      if(#[1]==77)
536      {  // undocumented feature for Gerhard's absPrimdecGTZ
537        if (size(tmpf)<2){ list abs_fac = list(var(n+1),poly(1)); }
538        else { list abs_fac=tmpf[2..size(tmpf)]; }
539        abs_fac=abs_fac,newvar;
540        return(string(abs_fac));
541      }
542    }
543  }
544  // preparing the output for SINGULAR standard
545  // a list: factors(ideal),multiplicities(intvec),minpolys(ideal),
546  // number of factors in the absolute factorization
547  // the output(except the coefficient) should have no denominators
548  // and no content
549  ideal facts,minpols;
550  intvec mults;
551  int nfacts;
552  number co=1;
553  minpols[1]=tmpf[1][1];
554  facts[1]=tmpf[1][2];     //the coefficient
555  for(i=2;i<=size(tmpf);i++)
556  {
557    minpols[i]=cleardenom(tmpf[i][1]);
558    facts[i]=cleardenom(tmpf[i][2]);
560  }
561  facts[1]=facts[1]*co;
562  for(i=1;i<=size(tmpm);i++)
563  {
564    mults[i]=tmpm[i];
565  }
566  for(i=2;i<=size(mults);i++)
567  {
568    nfacts=nfacts+mults[i]*deg(minpols[i]);
569  }
570  list absolute_factors=facts,mults,minpols,nfacts;
571
572  // create ring with extra parameter `newvar` for output:
573  setring(MP);
574  list Lout=ringlist(MP);
575  if(!pa)
576  {
577    list Lpar=list(char(MP),list(newvar),list(list("lp",intvec(1))),ideal(0));
578  }
579  else
580  {
581    list Lpar=Lout[1];
582    Lpar[2][size(Lpar[2])+1]=newvar;
583    vv=Lpar[3][1][2];
584    vv=vv,1;
585    Lpar[3][1][2]=vv;
586  }
587  Lout[1]=Lpar;
588  def MPo=ring(Lout);
589  setring(MPo);
590  list absolute_factors=imap(MPout,absolute_factors);
591  export absolute_factors;
592  setring(MP);
593
594  dbprint( printlevel-voice+3,"
595// 'absFactorize' created a ring, in which a list absolute_factors (the
596// absolute factors) is stored.
597// To access the list of absolute factors, type (if the name S was assigned
598// to the return value):
599        setring(S); absolute_factors; ");
600  return(MPo);
601}
602example
603{
604  "EXAMPLE:"; echo = 2;
605  ring R = (0), (x,y), lp;
606  poly p = (-7*x^2 + 2*x*y^2 + 6*x + y^4 + 14*y^2 + 47)*(5x2+y2)^3*(x-y)^4;
607  def S = absFactorize(p) ;
608  setring(S);
609  absolute_factors;
610}
611
612////////////////////////////////////////////////////////////////////////////////
613proc absBiFactorize(poly p, list #)
614"USAGE:  absBiFactorize(p [,s]);   p poly, s string
615ASSUME: coefficient field is the field of rational numbers or a
616        transcendental extension thereof
617RETURN: ring @code{R} which is obtained from the current basering
618        by adding a new parameter (if a string @code{s} is given as a
619        second input, the new parameter gets this string as name). The ring
620        @code{R} comes with a list @code{absolute_factors} with the
621        following entries:
622@format
623    absolute_factors[1]: ideal   (the absolute factors)
624    absolute_factors[2]: intvec  (the multiplicities)
625    absolute_factors[3]: ideal   (the minimal polynomials)
626    absolute_factors[4]: int     (total number of nontriv. absolute factors)
627@end format
628        The entry @code{absolute_factors[1][1]} is a constant, the
629        entry @code{absolute_factors[3][1]} is the parameter added to the
630        current ring.@*
631        Each of the remaining entries @code{absolute_factors[1][j]} stands for
632        a class of conjugated absolute factors. The corresponding entry
633        @code{absolute_factors[3][j]} is the minimal polynomial of the
634        field extension over which the factor is minimally defined (its degree
635        is the number of conjugates in the class). If the entry
636        @code{absolute_factors[3][j]} coincides with @code{absolute_factors[3][1]},
637        no field extension was necessary for the @code{j}th (class of)
638        absolute factor(s).
639NOTE:   The input polynomial p is assumed to be bivariate or to contain one para-
640        meter and one variable. All factors are presented denominator- and
641        content-free. The constant factor (first entry) is chosen such that the
642        product of all (!) the (denominator- and content-free) absolute factors
643        of @code{p} equals @code{p / absolute_factors[1][1]}.
645EXAMPLE: example absBiFactorize; shows an example
646"
647{
648  int dblevel = printlevel - voice + 2;
649  dbprint(dblevel,"Entering absfact.lib::absBiFactorize with ",p);
650  def MP = basering;
651  int i;
652  if (char(MP) != 0)
653  {
654    ERROR("// absfact.lib::absBiFactorize is only implemented for "+
655          "characteristic 0");
656  }
657  if(minpoly!=0)
658  {
659    ERROR("// absfact.lib::absBiFactorize is not implemented for algebraic "
660          +"extensions");
661  }
662
663  int n = nvars(MP);
664  int pa=npars(MP);
665  list lMP= ringlist(MP);
666  intvec vv,vk;
667  for(i=1;i<=n;i++){vv[i]=1;}
668  vk=vv,1;
669
670  //if the basering has parameters, add the parameters to the variables
671  //takes care about coefficients and possible denominators
672  if(pa>0)
673  {
674    poly qh=cleardenom(p);
675    if (p==0)
676    {
677      number cok=0;
678    }
679    else
680    {
682    }
683    p=qh;
684    string sp;
685    for(i=1;i<=npars(basering);i++)
686    {
687      sp=string(par(i));
688      sp=sp[2..size(sp)-1];
689      lMP[2][n+i]=sp;
690      vv=vv,1;
691    }
692    lMP[1]=0;
693    n=n+npars(MP);
694  }
695
696  // MPz is obtained by adding the new variable @z to MP
697  // ordering is wp(1...1)
698  // All the above subroutines work in MPz
699  string newvar;
700  if(size(#)>0)
701  {
702    if(typeof(#[1])=="string")
703    {
704      newvar=#[1];
705    }
706    else
707    {
708      newvar = "a";
709    }
710  }
711  else
712  {
713    newvar = "a";
714  }
715  if (newvar=="a")
716  {
717    if(belongTo(newvar, lMP[2])||defined(a)){newvar = "b";}
718    if(belongTo(newvar, lMP[2])||defined(b)){newvar = "c";}
719    if(belongTo(newvar, lMP[2])||defined(c)){newvar = "@c";}
720    while(belongTo(newvar, lMP[2]))
721    {
722       newvar = "@" + newvar;
723    }
724  }
725
726  // create ring with extra parameter `newvar` for output:
727  setring(MP);
728  list Lout=ringlist(MP);
729  if(!pa)
730  {
731    list Lpar=list(char(MP),list(newvar),list(list("lp",intvec(1))),ideal(0));
732  }
733  else
734  {
735    list Lpar=Lout[1];
736    Lpar[2][size(Lpar[2])+1]=newvar;
737    vv=Lpar[3][1][2];
738    vv=vv,1;
739    Lpar[3][1][2]=vv;
740  }
741  Lout[1]=Lpar;
742  def MPo=ring(Lout);
743  setring(MPo);
744
745  poly p=imap(MP,p);
746
747  if (size (variables (p)) > 2)
748  {
749    ERROR("// absfact.lib::absBiFactorize is not implemented for polynomials "
750          +"with more than 2 variables");
751  }
752
753  // special treatment in the homogeneous case, dropping one variable
754  // by substituting the first variable by 1
755  int ho=homog(p);
756  if(ho)
757  {
758    int dh=deg(p);
759    p=subst(p,var(1),1);
760    int di=deg(p);
761  }
762
763  list tmpf=system ("absBiFact", p);
764
765  // the homogeneous case, homogenizing the result
766  // the new variable has to have degree 0
767  // need to change the ring
768  if(ho)
769  {
770    list ll=ringlist(MPo);
771    vv[size(vv)]=0;
772    ll[3][1][2]=vv;
773    def MPhelp=ring(ll);
774    setring(MPhelp);
775    list tmpf=imap(MPo,tmpf);
776    for(i=2;i<=size(tmpf[1]);i++)
777    {
778      tmpf[1][i]=homog(tmpf[1][i],var(1));
779    }
780    if(dh>di)
781    {
782      tmpf[1][size(tmpf[1])+1]=var(1);
783      tmpf[2][size(tmpf[2])+1]=dh-di;
784      tmpf[3][size(tmpf[3])+1]=par(npars(MPo));
785      tmpf[4]= tmpf[4]+dh-di;
786    }
787    setring(MPo);
788    tmpf=imap(MPhelp,tmpf);
789  }
790
791  if (pa)
792  {
793    number cok=imap(MP,cok);
794    tmpf[1][1]=cok*tmpf[1][1];
795  }
796
797  // if we want the output as string
798  if(size(#)>0)
799  {
800    if(typeof(#[1])=="int")
801    {
802      if(#[1]==77)
803      {  // undocumented feature for Gerhard's absPrimdecGTZ
804        if (size(tmpf[1])<2){ list abs_fac = list(var(n+1),poly(1)); }
805        else
806        {
807          list abs_fac= tmpf[3][2];
808          abs_fac= abs_fac, tmpf[1][2];
809          for (i= 3; i <= size(tmpf[1]); i++)
810          {
811            abs_fac=abs_fac,tmpf[3][i];
812            abs_fac=abs_fac,tmpf[1][i];
813          }
814        }
815        abs_fac=abs_fac,newvar;
816        return(string(abs_fac));
817      }
818    }
819  }
820
821  list absolute_factors= tmpf;
822  export absolute_factors;
823  setring(MP);
824
825  dbprint( printlevel-voice+3,"
826// 'absBiFactorize' created a ring, in which a list absolute_factors (the
827// absolute factors) is stored.
828// To access the list of absolute factors, type (if the name S was assigned
829// to the return value):
830        setring(S); absolute_factors; ");
831  return(MPo);
832}
833example
834{
835  "EXAMPLE:"; echo = 2;
836  ring R = (0), (x,y), lp;
837  poly p = (-7*x^2 + 2*x*y^2 + 6*x + y^4 + 14*y^2 + 47)*(5x2+y2)^3*(x-y)^4;
838  def S = absBiFactorize(p) ;
839  setring(S);
840  absolute_factors;
841}
842
843
844/*
845  ring r=0,(x,t),dp;
846  poly p=x^4+(t^3-2t^2-2t)*x^3-(t^5-2t^4-t^2-2t-1)*x^2
847         -(t^6-4t^5+t^4+6t^3+2t^2)*x+(t^6-4t^5+2t^4+4t^3+t^2);
848  def S = absFactorize(p,"s");
849  setring(S);
850  absolute_factors;
851
852  ring r1=(0,a,b),(x,y),dp;
853  poly p=(a3-a2b+27ab3-27b4)/(a+b5)*x2+(a2+27b3)*y;
854  def S = absFactorize(p);
855  setring(S);
856  absolute_factors;
857
858  ring r2=0,(x,y,z,w),dp;
859  poly f=(x2+y2+z2)^2+w4;
860  def S =absFactorize(f);
861  setring(S);
862  absolute_factors;
863
864ring r=0,(x),dp;
865poly p=0;
866def S = absFactorize(p);
867setring(S);
868absolute_factors;
869
870ring r=0,(x),dp;
871poly p=7/11;
872def S = absFactorize(p);
873setring(S);
874absolute_factors;
875
876ring r=(0,a,b),(x,y),dp;
877poly p=0;
878def S = absFactorize(p);
879setring(S);
880absolute_factors;
881
882ring r=(0,a,b),(x,y),dp;
883poly p=(a+1)/b;
884def S = absFactorize(p);
885setring(S);
886absolute_factors;
887
888ring r=(0,a,b),(x,y),dp;
889poly p=(a+1)/b*x;
890def S = absFactorize(p,"s");
891setring(S);
892absolute_factors;
893
894ring r=(0,a,b),(x,y),dp;
895poly p=(a+1)/b*x + 1;
896def S = absFactorize(p,"s");
897setring(S);
898absolute_factors;
899
900ring r=(0,a,b),(x,y),dp;
901poly p=(a+1)/b*x + y;
902def S = absFactorize(p,"s");
903setring(S);
904absolute_factors;
905
906ring r=0,(x,t),dp;
907poly p=x^4+(t^3-2t^2-2t)*x^3-(t^5-2t^4-t^2-2t-1)*x^2
908       -(t^6-4t^5+t^4+6t^3+2t^2)*x+(t^6-4t^5+2t^4+4t^3+t^2);
909def S = absFactorize(p,"s");
910setring(S);
911absolute_factors;
912
913ring r1=(0,a,b),(x,y),dp;
914poly p=(a3-a2b+27ab3-27b4)/(a+b5)*x2+(a2+27b3)*y;
915def S = absFactorize(p);
916setring(S);
917absolute_factors;
918
919ring r2=0,(x,y,z,w),dp;
920poly f=(x2+y2+z2)^2+w4;
921def S =absFactorize(f);
922setring(S);
923absolute_factors;
924
925ring r3=0,(x,y,z,w),dp;
926poly f=(x2+y2+z2)^4+w8;
927def S =absFactorize(f);
928setring(S);
929absolute_factors;
930
931ring r4=0,(x,y),dp;
932poly f=y6-(2x2-2x-14)*y4-(4x3+35x2-6x-47)*y2+14x4-12x3-94x2;
933def S=absFactorize(f);
934setring(S);
935absolute_factors;
936
937ring R1 = 0, x, dp;
938def S1 = absFactorize(x4-2);
939setring(S1);
940absolute_factors;
941
942ring R3 = 0, (x,y), dp;
943poly f = x2y4+y6+2x3y2+2xy4-7x4+7x2y2+14y4+6x3+6xy2+47x2+47y2;
944def S3 = absFactorize(f);
945setring(S3);
946absolute_factors;
947
948ring R4 = 0, (x,y), dp;
949poly f = y4+2*xy2-7*x2+14*y2+6*x+47;
950def S4 = absFactorize(f);
951setring(S4);
952absolute_factors;
953
954*/
Note: See TracBrowser for help on using the repository browser.