axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / PGCD.spad

)abbrev package PGCD PolynomialGcdPackage
++ Author: Michael Lucks, P. Gianni, Frederic Lehobey
++ Date Last Updated: 17 June 1996
++ Description:
++ This package computes multivariate polynomial gcd's using
++ a hensel lifting strategy. The constraint on the coefficient
++ domain is imposed by the lifting strategy. It is assumed that
++ the coefficient domain has the property that almost all specializations
++ preserve the degree of the gcd.
 
PolynomialGcdPackage(E,OV,R,P) : SIG == CODE where
  E : OrderedAbelianMonoidSup
  OV : OrderedSet
  R : EuclideanDomain
  P : PolynomialCategory(R,E,OV)

  I    ==> Integer
  NNI  ==> NonNegativeInteger
  PI   ==> PositiveInteger
 
  SUPP ==> SparseUnivariatePolynomial P
 
  SIG ==> with

    gcd : (P,P) -> P
      ++ gcd(p,q) computes the gcd of the two polynomials p and q.
      ++
      ++X p1:=(x+1)*(x+6)
      ++X p2:=(x+1)*(x-6)
      ++X gcd(p1,p2)

    gcd : List P -> P
      ++ gcd(lp) computes the gcd of the list of polynomials lp.

    gcd : (SUPP,SUPP) -> SUPP
      ++ gcd(p,q) computes the gcd of the two polynomials p and q.

    gcd : List SUPP -> SUPP
      ++ gcd(lp) computes the gcd of the list of polynomials lp.

    gcdPrimitive : (P,P) -> P
      ++ gcdPrimitive(p,q) computes the gcd of the primitive polynomials
      ++ p and q.

    gcdPrimitive : (SUPP,SUPP) -> SUPP
      ++ gcdPrimitive(p,q) computes the gcd of the primitive polynomials
      ++ p and q.

    gcdPrimitive : List P -> P
      ++ gcdPrimitive lp computes the gcd of the list of primitive
      ++ polynomials lp.

  CODE ==> add
 
      SUP      ==> SparseUnivariatePolynomial R
 
      LGcd     ==> Record(locgcd:SUPP,goodint:List List R)
      UTerm    ==> Record(lpol:List SUP,lint:List List R,mpol:SUPP)
      pmod:R   :=  (prevPrime(2**26)$IntegerPrimesPackage(Integer))::R
 
      import MultivariateLifting(E,OV,R,P)
      import FactoringUtilities(E,OV,R,P)
 
                 --------  Local  Functions  --------
 
      myran           :    Integer   -> Union(R,"failed")
      better          :    (P,P)     -> Boolean
      failtest        :   (SUPP,SUPP,SUPP)    -> Boolean
      monomContent    :   (SUPP)   -> SUPP
      gcdMonom        :  (SUPP,SUPP)    -> SUPP
      gcdTermList     :    (P,P)     -> P
      good :  (SUPP,List OV,List List R) -> Record(upol:SUP,inval:List List R)
 
      chooseVal :  (SUPP,SUPP,List OV,List List R) -> Union(UTerm,"failed")
      localgcd        :  (SUPP,SUPP,List OV,List List R) -> LGcd
      notCoprime      :  (SUPP,SUPP, List NNI,List OV,List List R)   -> SUPP
      imposelc        : (List SUP,List OV,List R,List P) -> 
                          Union(List SUP, "failed")
 
      lift? :(SUPP,SUPP,UTerm,List NNI,List OV) -> _
               Union(s:SUPP,failed:"failed",notCoprime:"notCoprime")
      lift  :(SUPP,SUP,SUP,P,List OV,List NNI,List R) -> Union(SUPP,"failed")
 
                     ---- Local  functions ----
      -- test if something wrong happened in the gcd
      failtest(f:SUPP,p1:SUPP,p2:SUPP) : Boolean ==
        (p1 exquo f) case "failed" or (p2 exquo f) case "failed"
 
      -- Choose the integers
      chooseVal(p1:SUPP,p2:SUPP,lvr:List OV,_
                ltry:List List R):Union(UTerm,"failed") ==
        d1:=degree(p1)
        d2:=degree(p2)
        dd:NNI:=0$NNI
        nvr:NNI:=#lvr
        lval:List R :=[]
        range:I:=8
        repeat
          range:=2*range
          lval:=[ran(range) for i in 1..nvr]
          member?(lval,ltry) => "new point"
          ltry:=cons(lval,ltry)
          uf1:SUP:=completeEval(p1,lvr,lval)
          degree uf1 ^= d1 => "new point"
          uf2:SUP:= completeEval(p2,lvr,lval)
          degree uf2 ^= d2 => "new point"
          u:=gcd(uf1,uf2)
          du:=degree u
          --the univariate gcd is 1
          if du=0 then return [[1$SUP],ltry,0$SUPP]$UTerm
          ugcd:List SUP:=[u,(uf1 exquo u)::SUP,(uf2 exquo u)::SUP]
          uterm:=[ugcd,ltry,0$SUPP]$UTerm
          dd=0 => dd:=du
          --the degree is not changed
          du=dd =>
           --test if one of the polynomials is the gcd
            dd=d1 =>
              if ^((f:=p2 exquo p1) case "failed") then
                return [[u],ltry,p1]$UTerm
              if dd^=d2 then dd:=(dd-1)::NNI
 
            dd=d2 =>
              if ^((f:=p1 exquo p2) case "failed") then
                return [[u],ltry,p2]$UTerm
              dd:=(dd-1)::NNI
            return uterm
          --the new gcd has degree less
          du<dd => dd:=du
 
      good(f:SUPP,lvr:List OV, _
           ltry:List List R):Record(upol:SUP,inval:List List R) ==
        nvr:NNI:=#lvr
        range:I:=1
        while true repeat
          range:=2*range
          lval:=[ran(range) for i in 1..nvr]
          member?(lval,ltry) => "new point"
          ltry:=cons(lval,ltry)
          uf:=completeEval(f,lvr,lval)
          if degree gcd(uf,differentiate uf)=0 then return [uf,ltry]

      -- impose the right leading condition, check for failure.
      imposelc(lipol:List SUP, lvar:List OV, lval:List R,
               leadc:List P): Union(List SUP, "failed") ==
        result:List SUP :=[]
        for pol in lipol for leadpol in leadc repeat
           p1 := univariate eval(leadpol, lvar, lval) * pol
           p1u := p1 exquo leadingCoefficient pol
           p1u case "failed" => return "failed"
           result := cons(p1u::SUP, result)
        reverse result
 
      --Compute the gcd between not coprime polynomials
      notCoprime(g:SUPP, p2:SUPP, ldeg:List NNI,_
                 lvar1:List OV, ltry:List List R) : SUPP ==
        g1:=gcd(g,differentiate g)
        l1 := (g exquo g1)::SUPP
        lg:LGcd:=localgcd(l1,p2,lvar1,ltry)
        (l,ltry):=(lg.locgcd,lg.goodint)
        lval:=ltry.first
        p2l:=(p2 exquo l)::SUPP
        (gd1,gd2):=(l,l)
        ul:=completeEval(l,lvar1,lval)
        dl:=degree ul
        if degree gcd(ul,differentiate ul) ^=0 then
          newchoice:=good(l,lvar1,ltry)
          ul:=newchoice.upol
          ltry:=newchoice.inval
          lval:=ltry.first
        ug1:=completeEval(g1,lvar1,lval)
        ulist:=[ug1,completeEval(p2l,lvar1,lval)]
        lcpol:List P:=[leadingCoefficient g1, leadingCoefficient p2]
        while true repeat
          d:SUP:=gcd(cons(ul,ulist))
          if degree d =0 then return gd1
          lquo:=(ul exquo d)::SUP
          if degree lquo ^=0 then
            lgcd:=gcd(cons(leadingCoefficient l,lcpol))
            (gdl:=lift(l,d,lquo,lgcd,lvar1,ldeg,lval)) case "failed" =>
               return notCoprime(g,p2,ldeg,lvar1,ltry)
            l:=gd2:=gdl::SUPP
            ul:=completeEval(l,lvar1,lval)
            dl:=degree ul
          gd1:=gd1*gd2
          ulist:=[(uf exquo d)::SUP for uf in ulist]
 
      gcdPrimitive(p1:SUPP,p2:SUPP) : SUPP ==
        if (d1:=degree(p1)) > (d2:=degree(p2)) then
            (p1,p2):= (p2,p1)
            (d1,d2):= (d2,d1)
        degree p1 = 0 =>
            p1 = 0 => unitCanonical p2
            unitCanonical p1
        lvar:List OV:=
          sort((a:OV,b:OV):Boolean+->a>b,setUnion(variables p1,variables p2))
        empty? lvar =>
           raisePolynomial(gcd(lowerPolynomial p1,lowerPolynomial p2))
        (p2 exquo p1) case SUPP => unitCanonical p1
        ltry:List List R:=empty()
        totResult:=localgcd(p1,p2,lvar,ltry)
        result: SUPP:=totResult.locgcd
        -- special cases
        result=1 => 1$SUPP
        while failtest(result,p1,p2) repeat
          ltry:=totResult.goodint
          totResult:=localgcd(p1,p2,lvar,ltry)
          result:=totResult.locgcd
        result
 
       --local function for the gcd : it returns the evaluation point too
      localgcd(p1:SUPP,p2:SUPP,lvar:List(OV),ltry:List List R) : LGcd ==
        uterm:=chooseVal(p1,p2,lvar,ltry)::UTerm
        ltry:=uterm.lint
        listpol:= uterm.lpol
        ud:=listpol.first
        dd:= degree ud
        --the univariate gcd is 1
        dd=0 => [1$SUPP,ltry]$LGcd
        --one of the polynomials is the gcd
        dd=degree(p1) or dd=degree(p2) =>
                         [uterm.mpol,ltry]$LGcd
        ldeg:List NNI:=map(min,degree(p1,lvar),degree(p2,lvar))
        -- if there is a polynomial g s.t. g/gcd and gcd are coprime ...
        -- I can lift
        (h:=lift?(p1,p2,uterm,ldeg,lvar)) case notCoprime =>
          [notCoprime(p1,p2,ldeg,lvar,ltry),ltry]$LGcd
        h case failed => localgcd(p1,p2,lvar,ltry) -- skip bad values?
        [h.s,ltry]$LGcd
 
      -- content, internal functions return the poly if it is a monomial
      monomContent(p:SUPP):SUPP ==
        degree(p)=0 => 1
        md:= minimumDegree(p)
        monomial(gcd sort(better,coefficients p),md)
 
      -- Ordering for gcd purposes
      better(p1:P,p2:P):Boolean ==
        ground? p1 => true
        ground? p2 => false
        degree(p1,mainVariable(p1)::OV) < degree(p2,mainVariable(p2)::OV)
 
      best_to_front(l : List P) : List P ==
          ress := []
          best := first(l)
          for p in rest l repeat
              if better(p, best) then
                  ress := cons(best, ress)
                  best := p
              else
                  ress := cons(p, ress)
          cons(best, ress)

      -- Gcd between polynomial p1 and p2 with
      -- mainVariable p1 < x=mainVariable p2
      gcdTermList(p1:P,p2:P) : P ==
        termList := best_to_front(
           cons(p1,coefficients univariate(p2,(mainVariable p2)::OV)))
        q:P:=termList.first
        for term in termList.rest until q = 1$P repeat q:= gcd(q,term)
        q
 
      -- Gcd between polynomials with the same mainVariable
      gcd(p1:SUPP,p2:SUPP): SUPP ==
        if degree(p1) > degree(p2) then (p1,p2):= (p2,p1)
        degree p1 = 0 =>
           p1 = 0 => unitCanonical p2
           p1 = 1 => unitCanonical p1
           gcd(leadingCoefficient p1, content p2)::SUPP
        reductum(p1)=0 => gcdMonom(p1,monomContent p2)
        c1:= monomContent(p1)
        reductum(p2)=0 => gcdMonom(c1,p2)
        c2:= monomContent(p2)
        p1:= (p1 exquo c1)::SUPP
        p2:= (p2 exquo c2)::SUPP
        gcdPrimitive(p1,p2) * gcdMonom(c1,c2)
 
      -- gcd between 2 monomials
      gcdMonom(m1:SUPP,m2:SUPP):SUPP ==
        monomial(gcd(leadingCoefficient(m1),leadingCoefficient(m2)),
                 min(degree(m1),degree(m2)))

       --If there is a pol s.t. pol/gcd and gcd are coprime I can lift
      lift?(p1:SUPP,p2:SUPP,uterm:UTerm,ldeg:List NNI, _
            lvar:List OV) : _
              Union(s:SUPP,failed:"failed",notCoprime:"notCoprime") ==
        (listpol, lval) := (uterm.lpol, first(uterm.lint))
        d := first(listpol)
        listpol := rest(listpol)
        uf := listpol(1)
        f := p1
        --note uf and d not necessarily primitive
        if degree gcd(uf, d) ~= 0 then
          uf := listpol(2)
          f := p2
          if degree gcd(uf, d) ~= 0 then return ["notCoprime"]
        lgcd := gcd(leadingCoefficient p1, leadingCoefficient p2)
        l := lift(f, d, uf, lgcd, lvar, ldeg, lval)
        l case "failed" => ["failed"]
        [l :: SUPP]
 
       -- interface with the general "lifting" function
      lift(f:SUPP,d:SUP,uf:SUP,lgcd:P,lvar:List OV,
                        ldeg:List NNI,lval:List R):Union(SUPP,"failed") ==
        leadpol : Boolean := false
        lcf : P
        lcf := leadingCoefficient f
        df := degree f
        leadlist : List(P) := []

        if lgcd ^= 1 then
          leadpol := true
          f := lgcd*f
          ldeg := [n0+n1 for n0 in ldeg for n1 in degree(lgcd, lvar)]
          lcd : R := leadingCoefficient d
          lgcd1 :=
              degree(lgcd) = 0 => retract lgcd
              retract(eval(lgcd, lvar, lval))
          du := (lgcd1*d) exquo lcd
          du case "failed" => "failed"
          d := du::SUP
          uf := lcd*uf
        leadlist := [lgcd, lcf]
        lgu := imposelc([d, uf], lvar, lval, leadlist)
        lgu case "failed" => "failed"
        lg := lgu::List(SUP)
        (pl := lifting(f,lvar,lg,lval,leadlist,ldeg,pmod)) case "failed" =>
               "failed"
        plist := pl :: List SUPP
        (p0 : SUPP, p1 : SUPP) := (plist.first, plist.2)
        if completeEval(p0, lvar, lval) ^= lg.first then
           (p0, p1) := (p1, p0)
        not leadpol => p0
        p0 exquo content(p0)
 
      -- Gcd for two multivariate polynomials
      gcd(p1:P,p2:P) : P ==
        ground? p1 =>
          p1 := unitCanonical p1
          p1 = 1$P => p1
          p1 = 0$P => unitCanonical p2
          ground? p2 => gcd((retract p1)@R,(retract p2)@R)::P
          gcdTermList(p1,p2)
        ground? p2 =>
          p2 := unitCanonical p2
          p2 = 1$P => p2
          p2 = 0$P => unitCanonical p1
          gcdTermList(p2,p1)
        (p1:= unitCanonical(p1)) = (p2:= unitCanonical(p2)) => p1
        mv1:= mainVariable(p1)::OV
        mv2:= mainVariable(p2)::OV
        mv1 = mv2 => multivariate(gcd(univariate(p1,mv1),
                                      univariate(p2,mv1)),mv1)
        mv1 < mv2 => gcdTermList(p1,p2)
        gcdTermList(p2,p1)
 
      -- Gcd for a list of multivariate polynomials
      gcd(listp:List P) : P ==
        lf := best_to_front(listp)
        f:=lf.first
        for g in lf.rest repeat
          f:=gcd(f,g)
          if f=1$P then return f
        f

      gcd(listp:List SUPP) : SUPP ==
        lf:=sort((z1:SUPP,z2:SUPP):Boolean +-> degree(z1)<degree(z2),listp)
        f:=lf.first
        for g in lf.rest repeat
          f:=gcd(f,g)
          if f=1 then return f
        f

 
      -- Gcd for primitive polynomials
      gcdPrimitive(p1:P,p2:P):P ==
        (p1:= unitCanonical(p1)) = (p2:= unitCanonical(p2)) => p1
        ground? p1 =>
          ground? p2 => gcd((retract p1)@R,(retract p2)@R)::P
          p1 = 0$P => p2
          1$P
        ground? p2 =>
          p2 = 0$P => p1
          1$P
        mv1:= mainVariable(p1)::OV
        mv2:= mainVariable(p2)::OV
        mv1 = mv2 =>
          md:=min(minimumDegree(p1,mv1),minimumDegree(p2,mv2))
          mp:=1$P
          if md>1 then
            mp:=(mv1::P)**md
            p1:=(p1 exquo mp)::P
            p2:=(p2 exquo mp)::P
          up1 := univariate(p1,mv1)
          up2 := univariate(p2,mv2)
          mp*multivariate(gcdPrimitive(up1,up2),mv1)
        1$P
 
      -- Gcd for a list of primitive multivariate polynomials
      gcdPrimitive(listp:List P) : P ==
        lf:=sort(better,listp)
        f:=lf.first
        for g in lf.rest repeat
          f:=gcdPrimitive(f,g)
          if f=1$P then return f
        f