source: git/Singular/LIB/absfact.lib @ 4173c7

spielwiese
Last change on this file since 4173c7 was 4173c7, checked in by Hans Schoenemann <hannes@…>, 13 years ago
use div instead of /, part 1 git-svn-id: file:///usr/local/Singular/svn/trunk@14191 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 17.8 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.
24SEE ALSO: factorize
25
26PROCEDURES:
27  absFactorize();        absolute factorization of poly
28";
29
30////////////////////////////////////////////////////////////////////
31static proc partialDegree(poly p, int i)
32"USAGE:  partialDegree(p,i);   p poly, i int
33RETURN: int, the degree of p in the i-th variable
34"
35{
36  int n = nvars(basering);
37  intvec tmp;
38  tmp[n] = 0;
39  tmp[i] = 1;
40  return(deg(p,tmp));
41}
42////////////////////////////////////////////////////////////////////
43static proc belongTo(string s, list l)
44"USAGE:  belongTo(s,l);   s string, l list
45RETURN: 1 if s belongs to l, 0 otherwise
46"
47{
48  string tmp;
49  for(int i = 1; i <= size(l); i++) {
50    tmp = l[i];
51    if (tmp == s) {
52      return(1);
53    }
54  }
55  return(0);
56}
57////////////////////////////////////////////////////////////////////
58static proc variableWithSmallestPositiveDegree(poly p)
59"USAGE:  variableWithSmallestPositiveDegree(p);  p poly
60RETURN: int;  0 if p is constant. Otherwise, the index of the
61        variable which has the smallest positive degree in p.
62"
63{
64  int n = nvars(basering);
65  int v = 0;
66  int d = deg(p);
67  int d_loc;
68  for(int i = 1; i <= n; i++) {
69    d_loc = partialDegree(p, i);
70    if (d_loc >= 1 and d_loc <= d) {
71      v = i;
72      d = d_loc;
73    }
74  }
75  return(v);
76}
77////////////////////////////////////////////////////////////////////
78static proc smallestProperSimpleFactor(poly p)
79"USAGE:  smallestProperSimpleFactor(p);   p poly
80RETURN: poly:  a proper irreducible simple factor of p of smallest
81        degree. If no such factor exists, 0 is returned.
82"
83{
84  list p_facts = factorize(p);
85  int s = size(p_facts[1]);
86  int d = deg(p)+1;
87  poly q = 0;
88  poly f;
89  int e;
90  for(int i = 1; i <= s; i++)
91  {
92    f = p_facts[1][i];
93    e = deg(f);
94    if (e >= 1 and e < d and p_facts[2][i] == 1)
95    {
96      q = f / leadcoef(f);
97      d = e;
98    }
99  }
100  return(q);
101}
102////////////////////////////////////////////////////////////////////
103static proc smallestProperFactor(poly p)
104"USAGE:  smallestProperFactor(p);   p poly
105RETURN: poly:  a proper irreducible factor of p of smallest degree.
106        If p is constant, 0 is returned.
107"
108{
109  list p_facts = factorize(p);
110  int s = size(p_facts[1]);
111  int d = deg(p)+1;
112  poly q = 0;
113  poly f;
114  int e;
115  for(int i = 1; i <= s; i++)
116  {
117    f = p_facts[1][i];
118    e = deg(f);
119    if (e >= 1 and e < d)
120    {
121      q = f / leadcoef(f);
122      d = e;
123    }
124  }
125  return(q);
126}
127////////////////////////////////////////////////////////////////////
128static proc extensionContainingSmoothPoint(poly p, int m)
129"USAGE:  extensionContainingSmoothPoint(p,m);  p poly, m int
130RETURN: poly:  an irreducible univariate polynomial that defines an
131        algebraic extension of the current ground field that contains
132        a smooth point of the hypersurface defined by p=0.
133"
134{
135  int n = nvars(basering) - 1;
136  poly q = 0;
137  int i;
138  list a;
139  for(i=1;i<=n+1;i++){a[i] = 0;}
140  a[m] = var(n+1);
141  // The list a is to be taken with random entries in [-e, e].
142  // Every 10 * n trial, e is incremented by 1.
143  int e = 1;
144  int nbtrial = 0;
145  map h;
146  while (q == 0)
147  {
148    h = basering, a[1..n+1];
149    q = smallestProperSimpleFactor(h(p));
150    for(i = 1; i  <= n ; i = i + 1)
151    {
152      if (i != m)
153      {
154        a[i] = random(-e, e);
155      }
156    }
157    nbtrial++;
158    if (nbtrial >= 10 * n)
159    {
160      e = e + 1;
161      nbtrial = 0;
162    }
163  }
164  return(q);
165}
166////////////////////////////////////////////////////////////////////
167static proc RothsteinTragerResultant(poly g, poly f, int m)
168"USAGE:  RothsteinTragerResultant(g,f,m);  g,f poly, m int
169RETURN: poly
170NOTE:   To be called by the RothsteinTrager procedure only.
171"
172{
173  def MPz = basering;
174  int n = nvars(MPz) - 1;
175  int d = partialDegree(f, m);
176  poly df = diff(f, var(m));
177  list a;
178  int i;
179  for(i=1;i<=n+1;i++){ a[i] = 0; }
180  a[m] = var(m);
181  poly q = 0;
182  int e = 1;
183  int nbtrial = 0;
184  map h;
185  while (q == 0)
186  {
187    h = MPz, a[1..n+1];
188    q = resultant(h(f), h(df) * var(n+1) - h(g), var(m));
189    if (deg(q) == d)
190    {
191      return(q/leadcoef(q));
192    }
193    q = 0;
194    for(i = 1; i  <= n ; i++)
195    {
196      if (i != m)
197      {
198        a[i] = random(-e, e);
199      }
200    }
201    nbtrial++;
202    if (nbtrial >= 10 * n)
203    {
204      e++;
205      nbtrial = 0;
206    }
207  }
208}
209////////////////////////////////////////////////////////////////////
210static proc RothsteinTrager(list g, poly p, int m, int expectedDegQ)
211"USAGE:  RothsteinTrager(g,p,m,d);  g list, p poly, m,d int
212RETURN: list L consisting of two entries of type poly
213NOTE:   the return value is the Rothstein-Trager partial fraction
214        decomposition of the quotient s/p, where s is a generic linear
215        combination of the elements of g. The genericity via d
216        (the expected degree of L[1]).
217"
218{
219  def MPz = basering;
220  int n = nvars(MPz) - 1;
221  poly dp = diff(p, var(m));
222  int r = size(g);
223  list a;
224  int i;
225  for(i=1;i<=r;i++){a[i] = 0;}
226  a[r] = 1;
227  int nbtrial = 0;
228  int e = 1;
229  poly s;
230  poly q;
231  while (1)
232  {
233    s = 0;
234    for(i = 1; i <= r; i++){s = s + a[i] * g[i];}
235    q = RothsteinTragerResultant(s, p, m);
236    q = smallestProperFactor(q);
237    if (deg(q) == expectedDegQ)
238    {
239      // Go into the quotient by q(z)=0
240      ring MP_z = (0,var(n+1)), (x(1..n)), dp;
241      list lMP_z = ringlist(MP_z);
242      lMP_z[1][4] = ideal(imap(MPz,q));
243      list tmp = ringlist(MPz)[2];
244      lMP_z[2] = list(tmp[1..n]);
245      def MPq = ring(lMP_z);
246      setring(MPq);
247      poly f = gcd(imap(MPz, p), par(1) * imap(MPz, dp) - imap(MPz, s));
248      f = f / leadcoef(f);
249      setring(MPz);
250      return(list(q, imap(MPq, f)));
251    }
252    for(i = 1; i  <= r ; i++)
253    {
254      a[i] = random(-e, e);
255    }
256    nbtrial++;
257    if (nbtrial >= 10 * r)
258    {
259      e++;
260      nbtrial = 0;
261    }
262  }
263}
264////////////////////////////////////////////////////////////////////
265static proc absFactorizeIrreducible(poly p)
266"USAGE:   absFactorizeIrreducible(p);   p poly
267ASSUME:  p is an irreducible polynomial that does not depend on the last
268         variable @z of the basering.
269RETURN:  list L of two polynomials: q=L[1] is an irreducible polynomial of
270         minimal degree in @z such that p has an absolute factor
271         over K[@z]/<q>, and f represents such an absolute factor.
272"
273{
274  int dblevel = printlevel - voice + 2;
275  dbprint(dblevel,"Entering absfact.lib::absFactorizeIrreducible with ",p);
276  def MPz = basering;
277  int d = deg(p);
278  int n = nvars(MPz) - 1;
279
280  if (d < 1)
281  {
282    return(list(var(n+1), p));
283  }
284  int m = variableWithSmallestPositiveDegree(p);
285  // var(m) is now considered as the main variable.
286
287  poly q = extensionContainingSmoothPoint(p, m);
288  int r = deg(q);
289  if (r == 1)
290  {
291    return(list(var(n+1), p));
292  }
293
294  // Go into the quotient by q(z)=0
295  ring MP_z = (0,var(n+1)), (x(1..n)), dp;
296  list lMP_z = ringlist(MP_z);
297  lMP_z[1][4] = ideal(imap(MPz,q));
298  list tmp = ringlist(MPz)[2];
299  lMP_z[2] = list(tmp[1..n]);
300  def MPq = ring(lMP_z);
301  setring(MPq);
302  dbprint(dblevel-1,"Factoring in algebraic extension");
303  // "Factoring p in the algebraic extension...";
304  poly p_loc = imap(MPz, p);
305  poly f = smallestProperSimpleFactor(p_loc);
306  int degf = deg(f);
307
308  if (degf == d)
309  {
310    setring(MPz);
311    return(list(var(n+1), p));
312  }
313
314  if (degf * r == d)
315  {
316    setring(MPz);
317    return(list(q, imap(MPq, f)));
318  }
319  dbprint(dblevel-1,"Absolutely irreducible factor found");
320  dbprint(dblevel,"Minimizing field extension");
321  // "Need to find a minimal extension";
322  poly co_f = p_loc / f;
323  poly e = diff(f, var(m)) * co_f;
324  setring(MPz);
325  poly e = imap(MPq, e);
326  list g;
327  int i;
328  for(i = 1; i <= r; i++)
329  {
330    g[i] = subst(e, var(n+1), 0);
331    e = diff(e, var(n+1));
332  }
333  return(RothsteinTrager(g, p, m, d div degf));
334}
335
336////////////////////////////////////////////////////////////////////
337proc absFactorize(poly p, list #)
338"USAGE:  absFactorize(p [,s]);   p poly, s string
339ASSUME: coefficient field is the field of rational numbers or a
340        transcendental extension thereof
341RETURN: ring @code{R} which is obtained from the current basering
342        by adding a new parameter (if a string @code{s} is given as a
343        second input, the new parameter gets this string as name). The ring
344        @code{R} comes with a list @code{absolute_factors} with the
345        following entries:
346@format
347    absolute_factors[1]: ideal   (the absolute factors)
348    absolute_factors[2]: intvec  (the multiplicities)
349    absolute_factors[3]: ideal   (the minimal polynomials)
350    absolute_factors[4]: int     (total number of nontriv. absolute factors)
351@end format
352        The entry @code{absolute_factors[1][1]} is a constant, the
353        entry @code{absolute_factors[3][1]} is the parameter added to the
354        current ring.@*
355        Each of the remaining entries @code{absolute_factors[1][j]} stands for
356        a class of conjugated absolute factors. The corresponding entry
357        @code{absolute_factors[3][j]} is the minimal polynomial of the
358        field extension over which the factor is minimally defined (its degree
359        is the number of conjugates in the class). If the entry
360        @code{absolute_factors[3][j]} coincides with @code{absolute_factors[3][1]},
361        no field extension was necessary for the @code{j}th (class of)
362        absolute factor(s).
363NOTE:   All factors are presented denominator- and content-free. The constant
364        factor (first entry) is chosen such that the product of all (!) the
365        (denominator- and content-free) absolute factors of @code{p} equals
366        @code{p / absolute_factors[1][1]}.
367SEE ALSO: factorize, absPrimdecGTZ
368EXAMPLE: example absFactorize; shows an example
369"
370{
371  int dblevel = printlevel - voice + 2;
372  dbprint(dblevel,"Entering absfact.lib::absFactorize with ",p);
373  def MP = basering;
374  int i;
375  if (char(MP) != 0)
376  {
377    ERROR("// absfact.lib::absFactorize is only implemented for "+
378          "characteristic 0");
379  }
380  if(minpoly!=0)
381  {
382    ERROR("// absfact.lib::absFactorize is not implemented for algebraic "
383          +"extensions");
384  }
385
386  int n = nvars(MP);
387  int pa=npars(MP);
388  list lMP= ringlist(MP);
389  intvec vv,vk;
390  for(i=1;i<=n;i++){vv[i]=1;}
391  vk=vv,1;
392
393  //if the basering has parameters, add the parameters to the variables
394  //takes care about coefficients and possible denominators
395  if(pa>0)
396  {
397    poly qh=cleardenom(p);
398    if (p==0)
399    {
400      number cok=0;
401    }
402    else
403    {
404      number cok=leadcoef(p)/leadcoef(qh);
405    }
406    p=qh;
407    string sp;
408    for(i=1;i<=npars(basering);i++)
409    {
410      sp=string(par(i));
411      sp=sp[2..size(sp)-1];
412      lMP[2][n+i]=sp;
413      vv=vv,1;
414    }
415    lMP[1]=0;
416    n=n+npars(MP);
417  }
418
419  // MPz is obtained by adding the new variable @z to MP
420  // ordering is wp(1...1)
421  // All the above subroutines work in MPz
422  string newvar;
423  if(size(#)>0)
424  {
425    if(typeof(#[1])=="string")
426    {
427      newvar=#[1];
428    }
429    else
430    {
431      newvar = "a";
432    }
433  }
434  else
435  {
436    newvar = "a";
437  }
438  if (newvar=="a")
439  {
440    if(belongTo(newvar, lMP[2])||defined(a)){newvar = "b";}
441    if(belongTo(newvar, lMP[2])||defined(b)){newvar = "c";}
442    if(belongTo(newvar, lMP[2])||defined(c)){newvar = "@c";}
443    while(belongTo(newvar, lMP[2]))
444    {
445       newvar = "@" + newvar;
446    }
447  }
448  lMP[2][n+1] = newvar;
449
450  // new ordering
451  vv=vv,1;
452  list orst;
453  orst[1]=list("wp",vv);
454  orst[2]=list("C",0);
455  lMP[3]=orst;
456
457  def MPz = ring(lMP);
458  setring(MPz);
459  poly p=imap(MP,p);
460
461  // special treatment in the homogeneous case, dropping one variable
462  // by substituting the first variable by 1
463  int ho=homog(p);
464  if(ho)
465  {
466    int dh=deg(p);
467    p=subst(p,var(1),1);
468    int di=deg(p);
469  }
470  list rat_facts = factorize(p);
471  int s = size(rat_facts[1]);
472  list tmpf; // absolute factors
473  intvec tmpm; // respective multiplicities
474  tmpf[1] = list(var(n+1), leadcoef(imap(MP,p)));
475  tmpm[1] = 1;
476  poly tmp;
477  for(i = 2; i <= s; i++)
478  {
479    tmp = rat_facts[1][i];
480    tmp = tmp / leadcoef(tmp);
481    tmpf[i] = absFactorizeIrreducible(tmp);
482    tmpm[i] = rat_facts[2][i];
483  }
484  // the homogeneous case, homogenizing the result
485  // the new variable has to have degree 0
486  // need to change the ring
487  if(ho)
488  {
489    list ll=ringlist(MPz);
490    vv[size(vv)]=0;
491    ll[3][1][2]=vv;
492    def MPhelp=ring(ll);
493    setring(MPhelp);
494    list tmpf=imap(MPz,tmpf);
495    for(i=2;i<=size(tmpf);i++)
496    {
497      tmpf[i][2]=homog(tmpf[i][2],var(1));
498    }
499    if(dh>di)
500    {
501      tmpf[size(tmpf)+1]=list(var(n+1),var(1));
502      tmpm[size(tmpm)+1]=dh-di;
503    }
504    setring(MPz);
505    tmpf=imap(MPhelp,tmpf);
506  }
507  // in case of parameters we have to go back to the old ring
508  // taking care about constant factors
509  if(pa)
510  {
511    n=nvars(MP);
512    list lM=ringlist(MP);
513    orst[1]=list("wp",vk);
514    orst[2]=list("C",0);
515    lM[2][n+1] = newvar;
516    lM[3]=orst;
517    def MPout=ring(lM);
518    setring(MPout);
519    list tmpf=imap(MPz,tmpf);
520    number cok=imap(MP,cok);
521    tmpf[1][2]=cok*tmpf[1][2];
522  }
523  else
524  {
525    def MPout=MPz;
526  }
527  // if we want the output as string
528  if(size(#)>0)
529  {
530    if(typeof(#[1])=="int")
531    {
532      if(#[1]==77)
533      {  // undocumented feature for Gerhard's absPrimdecGTZ
534        if (size(tmpf)<2){ list abs_fac = list(var(n+1),poly(1)); }
535        else { list abs_fac=tmpf[2..size(tmpf)]; }
536        abs_fac=abs_fac,newvar;
537        return(string(abs_fac));
538      }
539    }
540  }
541  // preparing the output for SINGULAR standard
542  // a list: factors(ideal),multiplicities(intvec),minpolys(ideal),
543  // number of factors in the absolute factorization
544  // the output(except the coefficient) should have no denominators
545  // and no content
546  ideal facts,minpols;
547  intvec mults;
548  int nfacts;
549  number co=1;
550  minpols[1]=tmpf[1][1];
551  facts[1]=tmpf[1][2];     //the coefficient
552  for(i=2;i<=size(tmpf);i++)
553  {
554    minpols[i]=cleardenom(tmpf[i][1]);
555    facts[i]=cleardenom(tmpf[i][2]);
556    co=co*(leadcoef(tmpf[i][2])/leadcoef(facts[i]))^(deg(minpols[i])*tmpm[i]);
557  }
558  facts[1]=facts[1]*co;
559  for(i=1;i<=size(tmpm);i++)
560  {
561    mults[i]=tmpm[i];
562  }
563  for(i=2;i<=size(mults);i++)
564  {
565    nfacts=nfacts+mults[i]*deg(minpols[i]);
566  }
567  list absolute_factors=facts,mults,minpols,nfacts;
568
569  // create ring with extra parameter `newvar` for output:
570  setring(MP);
571  list Lout=ringlist(MP);
572  if(!pa)
573  {
574    list Lpar=list(char(MP),list(newvar),list(list("lp",intvec(1))),ideal(0));
575  }
576  else
577  {
578    list Lpar=Lout[1];
579    Lpar[2][size(Lpar[2])+1]=newvar;
580    vv=Lpar[3][1][2];
581    vv=vv,1;
582    Lpar[3][1][2]=vv;
583  }
584  Lout[1]=Lpar;
585  def MPo=ring(Lout);
586  setring(MPo);
587  list absolute_factors=imap(MPout,absolute_factors);
588  export absolute_factors;
589  setring(MP);
590
591  dbprint( printlevel-voice+3,"
592// 'absFactorize' created a ring, in which a list absolute_factors (the
593// absolute factors) is stored.
594// To access the list of absolute factors, type (if the name S was assigned
595// to the return value):
596        setring(S); absolute_factors; ");
597  return(MPo);
598}
599example
600{
601  "EXAMPLE:"; echo = 2;
602  ring R = (0), (x,y), lp;
603  poly p = (-7*x^2 + 2*x*y^2 + 6*x + y^4 + 14*y^2 + 47)*(5x2+y2)^3*(x-y)^4;
604  def S = absFactorize(p) ;
605  setring(S);
606  absolute_factors;
607}
608
609/*
610  ring r=0,(x,t),dp;
611  poly p=x^4+(t^3-2t^2-2t)*x^3-(t^5-2t^4-t^2-2t-1)*x^2
612         -(t^6-4t^5+t^4+6t^3+2t^2)*x+(t^6-4t^5+2t^4+4t^3+t^2);
613  def S = absFactorize(p,"s");
614  setring(S);
615  absolute_factors;
616
617  ring r1=(0,a,b),(x,y),dp;
618  poly p=(a3-a2b+27ab3-27b4)/(a+b5)*x2+(a2+27b3)*y;
619  def S = absFactorize(p);
620  setring(S);
621  absolute_factors;
622
623  ring r2=0,(x,y,z,w),dp;
624  poly f=(x2+y2+z2)^2+w4;
625  def S =absFactorize(f);
626  setring(S);
627  absolute_factors;
628
629ring r=0,(x),dp;
630poly p=0;
631def S = absFactorize(p);
632setring(S);
633absolute_factors;
634
635ring r=0,(x),dp;
636poly p=7/11;
637def S = absFactorize(p);
638setring(S);
639absolute_factors;
640
641ring r=(0,a,b),(x,y),dp;
642poly p=0;
643def S = absFactorize(p);
644setring(S);
645absolute_factors;
646
647ring r=(0,a,b),(x,y),dp;
648poly p=(a+1)/b;
649def S = absFactorize(p);
650setring(S);
651absolute_factors;
652
653ring r=(0,a,b),(x,y),dp;
654poly p=(a+1)/b*x;
655def S = absFactorize(p,"s");
656setring(S);
657absolute_factors;
658
659ring r=(0,a,b),(x,y),dp;
660poly p=(a+1)/b*x + 1;
661def S = absFactorize(p,"s");
662setring(S);
663absolute_factors;
664
665ring r=(0,a,b),(x,y),dp;
666poly p=(a+1)/b*x + y;
667def S = absFactorize(p,"s");
668setring(S);
669absolute_factors;
670
671ring r=0,(x,t),dp;
672poly p=x^4+(t^3-2t^2-2t)*x^3-(t^5-2t^4-t^2-2t-1)*x^2
673       -(t^6-4t^5+t^4+6t^3+2t^2)*x+(t^6-4t^5+2t^4+4t^3+t^2);
674def S = absFactorize(p,"s");
675setring(S);
676absolute_factors;
677
678ring r1=(0,a,b),(x,y),dp;
679poly p=(a3-a2b+27ab3-27b4)/(a+b5)*x2+(a2+27b3)*y;
680def S = absFactorize(p);
681setring(S);
682absolute_factors;
683
684ring r2=0,(x,y,z,w),dp;
685poly f=(x2+y2+z2)^2+w4;
686def S =absFactorize(f);
687setring(S);
688absolute_factors;
689
690ring r3=0,(x,y,z,w),dp;
691poly f=(x2+y2+z2)^4+w8;
692def S =absFactorize(f);
693setring(S);
694absolute_factors;
695
696ring r4=0,(x,y),dp;
697poly f=y6-(2x2-2x-14)*y4-(4x3+35x2-6x-47)*y2+14x4-12x3-94x2;
698def S=absFactorize(f);
699setring(S);
700absolute_factors;
701
702ring R1 = 0, x, dp;
703def S1 = absFactorize(x4-2);
704setring(S1);
705absolute_factors;
706
707ring R3 = 0, (x,y), dp;
708poly f = x2y4+y6+2x3y2+2xy4-7x4+7x2y2+14y4+6x3+6xy2+47x2+47y2;
709def S3 = absFactorize(f);
710setring(S3);
711absolute_factors;
712
713ring R4 = 0, (x,y), dp;
714poly f = y4+2*xy2-7*x2+14*y2+6*x+47;
715def S4 = absFactorize(f);
716setring(S4);
717absolute_factors;
718
719*/
Note: See TracBrowser for help on using the repository browser.