source: git/Singular/LIB/absfact.lib @ e8eba7

//
version="$Id$";
category="Factorization";
info="
LIBRARY: absfact.lib   Absolute factorization for characteristic 0
AUTHORS: Wolfram Decker,       decker at math.uni-sb.de
         Gregoire Lecerf,      lecerf at math.uvsq.fr
         Gerhard Pfister,      pfister at mathematik.uni-kl.de

OVERVIEW:
A library for computing the absolute factorization of multivariate
polynomials f with coefficients in a field K of characteristic zero.
Using Trager's idea, the implemented algorithm computes an absolutely
irreducible factor by factorizing over some finite extension field L
(which is chosen such that V(f) has a smooth point with coordinates in L).
Then a minimal extension field is determined making use of the
Rothstein-Trager partial fraction decomposition algorithm.

REFERENCES:
G. Cheze, G. Lecerf: Lifting and recombination techniques for absolute
                  factorization. Journal of Complexity, 23(3):380-420, 2007.

KEYWORDS: factorization; absolute factorization.
SEE ALSO: factorize

PROCEDURES:
  absFactorize();        absolute factorization of poly
";

////////////////////////////////////////////////////////////////////
static proc partialDegree(poly p, int i)
"USAGE:  partialDegree(p,i);   p poly, i int
RETURN: int, the degree of p in the i-th variable
"
{
  int n = nvars(basering);
  intvec tmp;
  tmp[n] = 0;
  tmp[i] = 1;
  return(deg(p,tmp));
}
////////////////////////////////////////////////////////////////////
static proc belongTo(string s, list l)
"USAGE:  belongTo(s,l);   s string, l list
RETURN: 1 if s belongs to l, 0 otherwise
"
{
  string tmp;
  for(int i = 1; i <= size(l); i++) {
    tmp = l[i];
    if (tmp == s) {
      return(1);
    }
  }
  return(0);
}
////////////////////////////////////////////////////////////////////
static proc variableWithSmallestPositiveDegree(poly p)
"USAGE:  variableWithSmallestPositiveDegree(p);  p poly
RETURN: int;  0 if p is constant. Otherwise, the index of the
        variable which has the smallest positive degree in p.
"
{
  int n = nvars(basering);
  int v = 0;
  int d = deg(p);
  int d_loc;
  for(int i = 1; i <= n; i++) {
    d_loc = partialDegree(p, i);
    if (d_loc >= 1 and d_loc <= d) {
      v = i;
      d = d_loc;
    }
  }
  return(v);
}
////////////////////////////////////////////////////////////////////
static proc smallestProperSimpleFactor(poly p)
"USAGE:  smallestProperSimpleFactor(p);   p poly
RETURN: poly:  a proper irreducible simple factor of p of smallest
        degree. If no such factor exists, 0 is returned.
"
{
  list p_facts = factorize(p);
  int s = size(p_facts[1]);
  int d = deg(p)+1;
  poly q = 0;
  poly f;
  int e;
  for(int i = 1; i <= s; i++)
  {
    f = p_facts[1][i];
    e = deg(f);
    if (e >= 1 and e < d and p_facts[2][i] == 1)
    {
      q = f / leadcoef(f);
      d = e;
    }
  }
  return(q);
}
////////////////////////////////////////////////////////////////////
static proc smallestProperFactor(poly p)
"USAGE:  smallestProperFactor(p);   p poly
RETURN: poly:  a proper irreducible factor of p of smallest degree.
        If p is constant, 0 is returned.
"
{
  list p_facts = factorize(p);
  int s = size(p_facts[1]);
  int d = deg(p)+1;
  poly q = 0;
  poly f;
  int e;
  for(int i = 1; i <= s; i++)
  {
    f = p_facts[1][i];
    e = deg(f);
    if (e >= 1 and e < d)
    {
      q = f / leadcoef(f);
      d = e;
    }
  }
  return(q);
}
////////////////////////////////////////////////////////////////////
static proc extensionContainingSmoothPoint(poly p, int m)
"USAGE:  extensionContainingSmoothPoint(p,m);  p poly, m int
RETURN: poly:  an irreducible univariate polynomial that defines an
        algebraic extension of the current ground field that contains
        a smooth point of the hypersurface defined by p=0.
"
{
  int n = nvars(basering) - 1;
  poly q = 0;
  int i;
  list a;
  for(i=1;i<=n+1;i++){a[i] = 0;}
  a[m] = var(n+1);
  // The list a is to be taken with random entries in [-e, e].
  // Every 10 * n trial, e is incremented by 1.
  int e = 1;
  int nbtrial = 0;
  map h;
  while (q == 0)
  {
    h = basering, a[1..n+1];
    q = smallestProperSimpleFactor(h(p));
    for(i = 1; i  <= n ; i = i + 1)
    {
      if (i != m)
      {
        a[i] = random(-e, e);
      }
    }
    nbtrial++;
    if (nbtrial >= 10 * n)
    {
      e = e + 1;
      nbtrial = 0;
    }
  }
  return(q);
}
////////////////////////////////////////////////////////////////////
static proc RothsteinTragerResultant(poly g, poly f, int m)
"USAGE:  RothsteinTragerResultant(g,f,m);  g,f poly, m int
RETURN: poly
NOTE:   To be called by the RothsteinTrager procedure only.
"
{
  def MPz = basering;
  int n = nvars(MPz) - 1;
  int d = partialDegree(f, m);
  poly df = diff(f, var(m));
  list a;
  int i;
  for(i=1;i<=n+1;i++){ a[i] = 0; }
  a[m] = var(m);
  poly q = 0;
  int e = 1;
  int nbtrial = 0;
  map h;
  while (q == 0)
  {
    h = MPz, a[1..n+1];
    q = resultant(h(f), h(df) * var(n+1) - h(g), var(m));
    if (deg(q) == d)
    {
      return(q/leadcoef(q));
    }
    q = 0;
    for(i = 1; i  <= n ; i++)
    {
      if (i != m)
      {
        a[i] = random(-e, e);
      }
    }
    nbtrial++;
    if (nbtrial >= 10 * n)
    {
      e++;
      nbtrial = 0;
    }
  }
}
////////////////////////////////////////////////////////////////////
static proc RothsteinTrager(list g, poly p, int m, int expectedDegQ)
"USAGE:  RothsteinTrager(g,p,m,d);  g list, p poly, m,d int
RETURN: list L consisting of two entries of type poly
NOTE:   the return value is the Rothstein-Trager partial fraction
        decomposition of the quotient s/p, where s is a generic linear
        combination of the elements of g. The genericity via d
        (the expected degree of L[1]).
"
{
  def MPz = basering;
  int n = nvars(MPz) - 1;
  poly dp = diff(p, var(m));
  int r = size(g);
  list a;
  int i;
  for(i=1;i<=r;i++){a[i] = 0;}
  a[r] = 1;
  int nbtrial = 0;
  int e = 1;
  poly s;
  poly q;
  while (1)
  {
    s = 0;
    for(i = 1; i <= r; i++){s = s + a[i] * g[i];}
    q = RothsteinTragerResultant(s, p, m);
    q = smallestProperFactor(q);
    if (deg(q) == expectedDegQ)
    {
      // Go into the quotient by q(z)=0
      ring MP_z = (0,var(n+1)), (x(1..n)), dp;
      list lMP_z = ringlist(MP_z);
      lMP_z[1][4] = ideal(imap(MPz,q));
      list tmp = ringlist(MPz)[2];
      lMP_z[2] = list(tmp[1..n]);
      def MPq = ring(lMP_z);
      setring(MPq);
      poly f = gcd(imap(MPz, p), par(1) * imap(MPz, dp) - imap(MPz, s));
      f = f / leadcoef(f);
      setring(MPz);
      return(list(q, imap(MPq, f)));
    }
    for(i = 1; i  <= r ; i++)
    {
      a[i] = random(-e, e);
    }
    nbtrial++;
    if (nbtrial >= 10 * r)
    {
      e++;
      nbtrial = 0;
    }
  }
}
////////////////////////////////////////////////////////////////////
static proc absFactorizeIrreducible(poly p)
"USAGE:   absFactorizeIrreducible(p);   p poly
ASSUME:  p is an irreducible polynomial that does not depend on the last
         variable @z of the basering.
RETURN:  list L of two polynomials: q=L[1] is an irreducible polynomial of
         minimal degree in @z such that p has an absolute factor
         over K[@z]/<q>, and f represents such an absolute factor.
"
{
  int dblevel = printlevel - voice + 2;
  dbprint(dblevel,"Entering absfact.lib::absFactorizeIrreducible with ",p);
  def MPz = basering;
  int d = deg(p);
  int n = nvars(MPz) - 1;

  if (d < 1)
  {
    return(list(var(n+1), p));
  }
  int m = variableWithSmallestPositiveDegree(p);
  // var(m) is now considered as the main variable.

  poly q = extensionContainingSmoothPoint(p, m);
  int r = deg(q);
  if (r == 1)
  {
    return(list(var(n+1), p));
  }

  list tmp = ringlist(MPz)[2];
  // Go into the quotient by q(z)=0
  ring MP_z = (0,var(n+1)), (x(1..n)), dp;
  list lMP_z = ringlist(MP_z);
  lMP_z[1][4] = ideal(imap(MPz,q));
  lMP_z[2] = list(tmp[1..n]);
  def MPq = ring(lMP_z);
  setring(MPq);
  dbprint(dblevel-1,"Factoring in algebraic extension");
  // "Factoring p in the algebraic extension...";
  poly p_loc = imap(MPz, p);
  poly f = smallestProperSimpleFactor(p_loc);
  int degf = deg(f);

  if (degf == d)
  {
    setring(MPz);
    return(list(var(n+1), p));
  }

  if (degf * r == d)
  {
    setring(MPz);
    return(list(q, imap(MPq, f)));
  }
  dbprint(dblevel-1,"Absolutely irreducible factor found");
  dbprint(dblevel,"Minimizing field extension");
  // "Need to find a minimal extension";
  poly co_f = p_loc / f;
  poly e = diff(f, var(m)) * co_f;
  setring(MPz);
  poly e = imap(MPq, e);
  list g;
  int i;
  for(i = 1; i <= r; i++)
  {
    g[i] = subst(e, var(n+1), 0);
    e = diff(e, var(n+1));
  }
  return(RothsteinTrager(g, p, m, d div degf));
}

////////////////////////////////////////////////////////////////////
proc absFactorize(poly p, list #)
"USAGE:  absFactorize(p [,s]);   p poly, s string
ASSUME: coefficient field is the field of rational numbers or a
        transcendental extension thereof
RETURN: ring @code{R} which is obtained from the current basering
        by adding a new parameter (if a string @code{s} is given as a
        second input, the new parameter gets this string as name). The ring
        @code{R} comes with a list @code{absolute_factors} with the
        following entries:
@format
    absolute_factors[1]: ideal   (the absolute factors)
    absolute_factors[2]: intvec  (the multiplicities)
    absolute_factors[3]: ideal   (the minimal polynomials)
    absolute_factors[4]: int     (total number of nontriv. absolute factors)
@end format
        The entry @code{absolute_factors[1][1]} is a constant, the
        entry @code{absolute_factors[3][1]} is the parameter added to the
        current ring.@*
        Each of the remaining entries @code{absolute_factors[1][j]} stands for
        a class of conjugated absolute factors. The corresponding entry
        @code{absolute_factors[3][j]} is the minimal polynomial of the
        field extension over which the factor is minimally defined (its degree
        is the number of conjugates in the class). If the entry
        @code{absolute_factors[3][j]} coincides with @code{absolute_factors[3][1]},
        no field extension was necessary for the @code{j}th (class of)
        absolute factor(s).
NOTE:   All factors are presented denominator- and content-free. The constant
        factor (first entry) is chosen such that the product of all (!) the
        (denominator- and content-free) absolute factors of @code{p} equals
        @code{p / absolute_factors[1][1]}.
SEE ALSO: factorize, absPrimdecGTZ
EXAMPLE: example absFactorize; shows an example
"
{
  int dblevel = printlevel - voice + 2;
  dbprint(dblevel,"Entering absfact.lib::absFactorize with ",p);
  def MP = basering;
  int i;
  if (char(MP) != 0)
  {
    ERROR("// absfact.lib::absFactorize is only implemented for "+
          "characteristic 0");
  }
  if(minpoly!=0)
  {
    ERROR("// absfact.lib::absFactorize is not implemented for algebraic "
          +"extensions");
  }

  int n = nvars(MP);
  int pa=npars(MP);
  list lMP= ringlist(MP);
  list buflMP= lMP;
  intvec vv,vk;
  for(i=1;i<=n;i++){vv[i]=1;}
  vk=vv,1;

  //if the basering has parameters, add the parameters to the variables
  //takes care about coefficients and possible denominators
  if(pa>0)
  {
    poly qh=cleardenom(p);
    if (p==0)
    {
      number cok=0;
    }
    else
    {
      number cok=leadcoef(p)/leadcoef(qh);
    }
    p=qh;
    string sp;
    for(i=1;i<=npars(basering);i++)
    {
      sp=string(par(i));
      sp=sp[2..size(sp)-1];
      lMP[2][n+i]=sp;
      vv=vv,1;
    }
    lMP[1]=0;
    n=n+npars(MP);
  }

  // MPz is obtained by adding the new variable @z to MP
  // ordering is wp(1...1)
  // All the above subroutines work in MPz
  string newvar;
  if(size(#)>0)
  {
    if(typeof(#[1])=="string")
    {
      newvar=#[1];
    }
    else
    {
      newvar = "a";
    }
  }
  else
  {
    newvar = "a";
  }
  if (newvar=="a")
  {
    if(belongTo(newvar, lMP[2])||defined(a)){newvar = "b";}
    if(belongTo(newvar, lMP[2])||defined(b)){newvar = "c";}
    if(belongTo(newvar, lMP[2])||defined(c)){newvar = "@c";}
    while(belongTo(newvar, lMP[2]))
    {
       newvar = "@" + newvar;
    }
  }
  lMP[2][n+1] = newvar;

  // new ordering
  vv=vv,1;
  list orst;
  orst[1]=list("wp",vv);
  orst[2]=list("C",0);
  lMP[3]=orst;

  def MPz = ring(lMP);
  setring(MPz);
  poly p=imap(MP,p);

  // special treatment in the homogeneous case, dropping one variable
  // by substituting the first variable by 1
  int ho=homog(p);
  if(ho)
  {
    int dh=deg(p);
    p=subst(p,var(1),1);
    int di=deg(p);
  }
  list rat_facts = factorize(p);
  int s = size(rat_facts[1]);
  list tmpf; // absolute factors
  intvec tmpm; // respective multiplicities
  tmpf[1] = list(var(n+1), leadcoef(imap(MP,p)));
  tmpm[1] = 1;
  poly tmp;
  for(i = 2; i <= s; i++)
  {
    tmp = rat_facts[1][i];
    tmp = tmp / leadcoef(tmp);
    tmpf[i] = absFactorizeIrreducible(tmp);
    tmpm[i] = rat_facts[2][i];
  }
  // the homogeneous case, homogenizing the result
  // the new variable has to have degree 0
  // need to change the ring
  if(ho)
  {
    list ll=ringlist(MPz);
    vv[size(vv)]=0;
    ll[3][1][2]=vv;
    def MPhelp=ring(ll);
    setring(MPhelp);
    list tmpf=imap(MPz,tmpf);
    for(i=2;i<=size(tmpf);i++)
    {
      tmpf[i][2]=homog(tmpf[i][2],var(1));
    }
    if(dh>di)
    {
      tmpf[size(tmpf)+1]=list(var(n+1),var(1));
      tmpm[size(tmpm)+1]=dh-di;
    }
    setring(MPz);
    tmpf=imap(MPhelp,tmpf);
  }
  // in case of parameters we have to go back to the old ring
  // taking care about constant factors
  if(pa)
  {
    setring(MP);
    n=nvars(MP);
    list lM=ringlist(MP);
    orst[1]=list("wp",vk);
    orst[2]=list("C",0);
    lM[2][n+1] = newvar;
    lM[3]=orst;
    def MPout=ring(lM);
    setring(MPout);
    list tmpf=imap(MPz,tmpf);
    number cok=imap(MP,cok);
    tmpf[1][2]=cok*tmpf[1][2];
  }
  else
  {
    def MPout=MPz;
  }
  // if we want the output as string
  if(size(#)>0)
  {
    if(typeof(#[1])=="int")
    {
      if(#[1]==77)
      {  // undocumented feature for Gerhard's absPrimdecGTZ
        if (size(tmpf)<2){ list abs_fac = list(var(n+1),poly(1)); }
        else { list abs_fac=tmpf[2..size(tmpf)]; }
        abs_fac=abs_fac,newvar;
        return(string(abs_fac));
      }
    }
  }
  // preparing the output for SINGULAR standard
  // a list: factors(ideal),multiplicities(intvec),minpolys(ideal),
  // number of factors in the absolute factorization
  // the output(except the coefficient) should have no denominators
  // and no content
  ideal facts,minpols;
  intvec mults;
  int nfacts;
  number co=1;
  minpols[1]=tmpf[1][1];
  facts[1]=tmpf[1][2];     //the coefficient
  for(i=2;i<=size(tmpf);i++)
  {
    minpols[i]=cleardenom(tmpf[i][1]);
    facts[i]=cleardenom(tmpf[i][2]);
    co=co*(leadcoef(tmpf[i][2])/leadcoef(facts[i]))^(deg(minpols[i])*tmpm[i]);
  }
  facts[1]=facts[1]*co;
  for(i=1;i<=size(tmpm);i++)
  {
    mults[i]=tmpm[i];
  }
  for(i=2;i<=size(mults);i++)
  {
    nfacts=nfacts+mults[i]*deg(minpols[i]);
  }
  list absolute_factors=facts,mults,minpols,nfacts;

  // create ring with extra parameter `newvar` for output:
  setring(MP);
  list Lout=ringlist(MP);
  if(!pa)
  {
    list Lpar=list(char(MP),list(newvar),list(list("lp",intvec(1))),ideal(0));
  }
  else
  {
    list Lpar=Lout[1];
    Lpar[2][size(Lpar[2])+1]=newvar;
    vv=Lpar[3][1][2];
    vv=vv,1;
    Lpar[3][1][2]=vv;
  }
  Lout[1]=Lpar;
  def MPo=ring(Lout);
  setring(MPo);
  list absolute_factors=imap(MPout,absolute_factors);
  export absolute_factors;
  setring(MP);

  dbprint( printlevel-voice+3,"
// 'absFactorize' created a ring, in which a list absolute_factors (the
// absolute factors) is stored.
// To access the list of absolute factors, type (if the name S was assigned
// to the return value):
        setring(S); absolute_factors; ");
  return(MPo);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring R = (0), (x,y), lp;
  poly p = (-7*x^2 + 2*x*y^2 + 6*x + y^4 + 14*y^2 + 47)*(5x2+y2)^3*(x-y)^4;
  def S = absFactorize(p) ;
  setring(S);
  absolute_factors;
}

/*
  ring r=0,(x,t),dp;
  poly p=x^4+(t^3-2t^2-2t)*x^3-(t^5-2t^4-t^2-2t-1)*x^2
         -(t^6-4t^5+t^4+6t^3+2t^2)*x+(t^6-4t^5+2t^4+4t^3+t^2);
  def S = absFactorize(p,"s");
  setring(S);
  absolute_factors;

  ring r1=(0,a,b),(x,y),dp;
  poly p=(a3-a2b+27ab3-27b4)/(a+b5)*x2+(a2+27b3)*y;
  def S = absFactorize(p);
  setring(S);
  absolute_factors;

  ring r2=0,(x,y,z,w),dp;
  poly f=(x2+y2+z2)^2+w4;
  def S =absFactorize(f);
  setring(S);
  absolute_factors;

ring r=0,(x),dp;
poly p=0;
def S = absFactorize(p);
setring(S);
absolute_factors;

ring r=0,(x),dp;
poly p=7/11;
def S = absFactorize(p);
setring(S);
absolute_factors;

ring r=(0,a,b),(x,y),dp;
poly p=0;
def S = absFactorize(p);
setring(S);
absolute_factors;

ring r=(0,a,b),(x,y),dp;
poly p=(a+1)/b;
def S = absFactorize(p);
setring(S);
absolute_factors;

ring r=(0,a,b),(x,y),dp;
poly p=(a+1)/b*x;
def S = absFactorize(p,"s");
setring(S);
absolute_factors;

ring r=(0,a,b),(x,y),dp;
poly p=(a+1)/b*x + 1;
def S = absFactorize(p,"s");
setring(S);
absolute_factors;

ring r=(0,a,b),(x,y),dp;
poly p=(a+1)/b*x + y;
def S = absFactorize(p,"s");
setring(S);
absolute_factors;

ring r=0,(x,t),dp;
poly p=x^4+(t^3-2t^2-2t)*x^3-(t^5-2t^4-t^2-2t-1)*x^2
       -(t^6-4t^5+t^4+6t^3+2t^2)*x+(t^6-4t^5+2t^4+4t^3+t^2);
def S = absFactorize(p,"s");
setring(S);
absolute_factors;

ring r1=(0,a,b),(x,y),dp;
poly p=(a3-a2b+27ab3-27b4)/(a+b5)*x2+(a2+27b3)*y;
def S = absFactorize(p);
setring(S);
absolute_factors;

ring r2=0,(x,y,z,w),dp;
poly f=(x2+y2+z2)^2+w4;
def S =absFactorize(f);
setring(S);
absolute_factors;

ring r3=0,(x,y,z,w),dp;
poly f=(x2+y2+z2)^4+w8;
def S =absFactorize(f);
setring(S);
absolute_factors;

ring r4=0,(x,y),dp;
poly f=y6-(2x2-2x-14)*y4-(4x3+35x2-6x-47)*y2+14x4-12x3-94x2;
def S=absFactorize(f);
setring(S);
absolute_factors;

ring R1 = 0, x, dp;
def S1 = absFactorize(x4-2);
setring(S1);
absolute_factors;

ring R3 = 0, (x,y), dp;
poly f = x2y4+y6+2x3y2+2xy4-7x4+7x2y2+14y4+6x3+6xy2+47x2+47y2;
def S3 = absFactorize(f);
setring(S3);
absolute_factors;

ring R4 = 0, (x,y), dp;
poly f = y4+2*xy2-7*x2+14*y2+6*x+47;
def S4 = absFactorize(f);
setring(S4);
absolute_factors;

*/