path: root/src/algebra/gpgcd.spad.pamphlet

\documentclass{article}
\usepackage{axiom}
\begin{document}
\title{\$SPAD/src/algebra gpgcd.spad}
\author{The Axiom Team}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package GENPGCD GeneralPolynomialGcdPackage}
<<package GENPGCD GeneralPolynomialGcdPackage>>=
)abbrev package GENPGCD GeneralPolynomialGcdPackage
++ Description:
++ This package provides operations for GCD computations
++ on polynomials 
GeneralPolynomialGcdPackage(E,OV,R,P):C == T where
    R     :  PolynomialFactorizationExplicit
    P     :  PolynomialCategory(R,E,OV)
    OV    :  OrderedSet
    E     :  OrderedAbelianMonoidSup
 
    SUPP      ==> SparseUnivariatePolynomial P
--JHD    ContPrim  ==> Record(cont:P,prim:P)
 
    C == with
           gcdPolynomial     :   (SUPP,SUPP) -> SUPP
		++ gcdPolynomial(p,q) returns the GCD of p and q
           randomR        : ()                ->R
		++ randomR() should be local but conditional
--JHD      gcd               :   (P,P)    -> P
--JHD      gcd               :   List P   -> P
--JHD      gcdprim           :   (P,P)    -> P
--JHD      gcdprim           :   List P   -> P
 
--JHD      gcdcofact         :   List P   -> List P
--JHD      gcdcofactprim     :   List P   -> List P
 
--JHD      primitate         :   (P,OV)   -> P
--JHD      primitate         :    SUPP    -> SUPP
 
--JHD      content           :     P      -> P
--JHD      content           :   List P   -> List P
--JHD      contprim          :   List P   -> List ContPrim
 
--JHD      monomContent      :   (P,OV)   -> P
--JHD      monomContent      :    SUPP    -> SUPP
 
 
    T == add
 
      SUPR     ==> SparseUnivariatePolynomial R
--JHD      SUPLGcd  ==> Record(locgcd:SUPP,goodint:List R)
--JHD      LGcd     ==> Record(locgcd:P,goodint:List R)
--JHD      UTerm    ==> Record(lpol:List SUPR,lint:List R,mpol:P)
--JHD--JHD      pmod:R   :=  (prevPrime(2**26)$IntegerPrimesPackage(Integer))::R
 
--JHD      import MultivariateLifting(E,OV,R,P,pmod)
      import UnivariatePolynomialCategoryFunctions2(R,SUPR,P,SUPP)
      import UnivariatePolynomialCategoryFunctions2(P,SUPP,R,SUPR)
                 --------  Local  Functions  --------
 
--JHD      abs             :      P       -> P
      better          :    (P,P)     -> Boolean
--JHD      failtest        :   (P,P,P)    -> Boolean
--JHD      gcdMonom        :  (P,P,OV)    -> P
--JHD      gcdTermList     :    (P,P)     -> P
--JHD      gcdPrim         :  (P,P,OV)    -> P
--JHD      gcdSameMainvar  :  (P,P,OV)    -> P
--JHD      internal        :  (P,P,OV)    -> P
--JHD      good            :  (P,List OV) -> Record(upol:SUPR,inval:List R)
--JHD      gcdPrs          : (P,P,NNI,OV) -> Union(P,"failed")
--JHD
--JHD      chooseVal       :     (P,P,List OV)          -> UTerm
--JHD      localgcd        :     (P,P,List OV)          -> LGcd
--JHD      notCoprime      :  (P,P, List NNI,List OV)   -> P
--JHD      imposelc        : (List SUPR,List OV,List R,List P)  -> List SUPR
 
--JHD      lift? :(P,P,UTerm,List NNI,List OV)  -> Union("failed",P)
--      lift  :(P,SUPR,SUPR,P,List OV,List NNI,List R) -> P
      lift  :  (SUPR,SUPP,SUPR,List OV,List R) -> Union(SUPP,"failed")
         -- lifts first and third arguments as factors of the second
         -- fourth is number of variables.
--JHD      monomContent      :   (P,OV)   -> P
      monomContentSup   :    SUPP    -> SUPP
--
--JHD  gcdcofact         :   List P   -> List P
 
      gcdTrivial      :  (SUPP,SUPP)   -> SUPP
      gcdSameVariables:  (SUPP,SUPP,List OV)    -> SUPP
      recursivelyGCDCoefficients: (SUPP,List OV,SUPP,List OV) -> SUPP
      flatten         : (SUPP,List OV) -> SUPP
                        -- evaluates out all variables in the second
                        -- argument, leaving a polynomial of the same
                        -- degree
--    eval            : (SUPP,List OV,List R) -> SUPP
      variables       :  SUPP          -> List OV
                     ---- JHD's exported functions ---
      gcdPolynomial(p1:SUPP,p2:SUPP) ==
        zero? p1 => p2
        zero? p2 => p1
        0=degree p1  => gcdTrivial(p1,p2)
        0=degree p2  => gcdTrivial(p2,p1)
        if degree p1 < degree p2 then (p1,p2):=(p2,p1)
        p1 exquo p2 case SUPP  => (unitNormal p2).canonical
        c1:= monomContentSup(p1)
        c2:= monomContentSup(p2)
        p1:= (p1 exquo c1)::SUPP
        p2:= (p2 exquo c2)::SUPP
        (p1 exquo p2) case SUPP => (unitNormal p2).canonical * gcd(c1,c2)
        vp1:=variables p1
        vp2:=variables p2
        v1:=setDifference(vp1,vp2)
        v2:=setDifference(vp2,vp1)
        #v1 = 0 and #v2 = 0 => gcdSameVariables(p1,p2,vp1)*gcd(c1,c2)
                 -- all variables are in common
        v:=setDifference(vp1,v1)
        pp1:=flatten(p1,v1)
        pp2:=flatten(p2,v2)
        g:=gcdSameVariables(pp1,pp2,v)
        one? g => gcd(c1,c2)::SUPP
        (#v1 = 0 or not (p1 exquo g) case "failed") and
                     -- if #vi = 0 then pp1 = p1, so we know g divides
              (#v2 = 0 or not (p2 exquo g) case "failed")
            => g*gcd(c1,c2)  -- divdes them both, so is the gcd
        -- OK, so it's not the gcd: try again
        v:=variables g -- there can be at most these variables in answer
        v1:=setDifference(vp1,v)
        v2:=setDifference(vp2,v)
        if (#v1 = 0) then g:= gcdSameVariables(g,flatten(p2,v2),v)
        else if (#v2=0) then g:=gcdSameVariables(g,flatten(p1,v1),v)
        else g:=gcdSameVariables(g,flatten(p1,v1)-flatten(p2,v2),v)
        one? g => gcd(c1,c2)::SUPP
        (#v1 = 0 or not (p1 exquo g) case "failed") and
              (#v2 = 0 or not (p2 exquo g) case "failed")
            => g*gcd(c1,c2)::SUPP  -- divdes them both, so is the gcd
        v:=variables g -- there can be at most these variables in answer
        v1:=setDifference(vp1,v)
        if #v1 ~= 0 then
           g:=recursivelyGCDCoefficients(g,v,p1,v1)
           one? g => return gcd(c1,c2)::SUPP
           v:=variables g -- there can be at most these variables in answer
        v2:=setDifference(vp2,v)
        recursivelyGCDCoefficients(g,v,p2,v2)*gcd(c1,c2)
      if R has StepThrough then
         randomCount:R := init()
         randomR() ==
            (v:=nextItem(randomCount)) case R =>
                randomCount:=v
                v
            SAY("Taking next stepthrough range in GeneralPolynomialGcdPackage")$Lisp
            randomCount:=init()
            randomCount
      else
            randomR() == (random$Integer() rem 100)::R
                     ---- JHD's local functions ---
      gcdSameVariables(p1:SUPP,p2:SUPP,lv:List OV) ==
            -- two non-trivial primitive (or, at least, we don't care
            -- about content)
            -- polynomials with precisely the same degree
          #lv = 0 => map(#1::P,gcdPolynomial(map(ground,p1),
                                             map(ground,p2)))
          degree p2 = 1 =>
            p1 exquo p2 case SUPP => p2
            1
          gcdLC:=gcd(leadingCoefficient p1,leadingCoefficient p2)
          lr:=[randomR() for vv in lv]
          count:NonNegativeInteger:=0
          while count<10 repeat
	    while zero? eval(gcdLC,lv,lr) and count<10 repeat
		  lr:=[randomR() for vv in lv]
		  count:=count+1
	    count = 10 => error "too many evaluations in GCD code"
	    up1:SUPR:=map(ground eval(#1,lv,lr),p1)
	    up2:SUPR:=map(ground eval(#1,lv,lr),p2)
	    u:=gcdPolynomial(up1,up2)
	    degree u = 0 => return 1
            -- let's pick a second one, just to check
            lrr:=[randomR() for vv in lv]
	    while zero? eval(gcdLC,lv,lrr) and count<10 repeat
		  lrr:=[randomR() for vv in lv]
		  count:=count+1
	    count = 10 => error "too many evaluations in GCD code"
	    vp1:SUPR:=map(ground eval(#1,lv,lrr),p1)
	    vp2:SUPR:=map(ground eval(#1,lv,lrr),p2)
	    v:=gcdPolynomial(vp1,vp2)
	    degree v = 0 => return 1
            if degree v < degree u then
               u:=v
               up1:=vp1
               up2:=vp2
               lr:=lrr
	    up1:=(up1 exquo u)::SUPR
	    degree gcd(u,up1) = 0 =>
                ans:=lift(u,p1,up1,lv,lr)
                ans case SUPP => return ans
                "next"
	    up2:=(up2 exquo u)::SUPR
	    degree gcd(u,up2) = 0 =>
                ans:=lift(u,p2,up2,lv,lr)
                ans case SUPP => return ans
                "next"
	    -- so neither cofactor is relatively prime
	    count:=0
	    while count < 10 repeat
	       r:=randomR()
	       uu:=up1+r*up2
	       degree gcd(u,uu)=0 =>
                 ans:= lift(u,p1+r::P *p2,uu,lv,lr)
                 ans case SUPP => return ans
                 "next"
	    error "too many evaluations in GCD code"
          count >= 10 => error "too many evaluations in GCD code"
      lift(gR:SUPR,p:SUPP,cfR:SUPR,lv:List OV,lr:List R) ==
        -- lift the coprime factorisation gR*cfR = (univariate of p)
        -- where the variables lv have been evaluated at lr
        lcp:=leadingCoefficient p
        g:=monomial(lcp,degree gR)+map(#1::P,reductum gR)
        cf:=monomial(lcp,degree cfR)+map(#1::P,reductum cfR)
        p:=lcp*p       -- impose leaidng coefficient of p on each factor
        while lv ~= [] repeat
           v:=first lv
           r:=first lr
           lv:=rest lv
           lr:=rest lr
           thisp:=map(eval(#1,lv,lr),p)
           d:="max"/[degree(c,v) for c in coefficients p]
           prime:=v::P - r::P
           pn:=prime
           origFactors:=[g,cf]::List SUPP
           for n in 1..d repeat
              Ecart:=(thisp- g*cf) exquo pn
              Ecart case "failed" =>
                 error "failed lifting in hensel in Complex Polynomial GCD"
              zero? Ecart => leave
              step:=solveLinearPolynomialEquation(origFactors,
                                                  map(eval(#1,v,r),Ecart::SUPP))
              step case "failed" => return "failed"
              g:=g+pn*first step
              cf:=cf+pn*second step
              pn:=pn*prime
           thisp ~= g*cf => return "failed"
        g
      recursivelyGCDCoefficients(g:SUPP,v:List OV,p:SUPP,pv:List OV) ==
         mv:=first pv   -- take each coefficient w.r.t. mv
         pv:=rest pv    -- and recurse on pv as necessary
         d:="max"/[degree(u,mv) for u in coefficients p]
         for i in 0..d repeat
             p1:=map(coefficient(#1,mv,i),p)
             oldg:=g
             if pv = [] then g:=gcdSameVariables(g,p1,v)
             else g:=recursivelyGCDCoefficients(p,v,p1,pv)
             one? g => return 1
             g~=oldg =>
                oldv:=v
                v:=variables g
                pv:=setUnion(pv,setDifference(v,oldv))
         g
      flatten(p1:SUPP,lv:List OV) ==
         #lv = 0 => p1
         lr:=[ randomR() for vv in lv]
         dg:=degree p1
         ans : SUPP
         while dg ~= degree (ans:= map(eval(#1,lv,lr),p1)) repeat
           lr:=[ randomR() for vv in lv]
         ans
--      eval(p1:SUPP,lv:List OV,lr:List R) == map(eval(#1,lv,lr),p1)
      variables(p1:SUPP) ==
        removeDuplicates ("concat"/[variables u for u in coefficients p1])
      gcdTrivial(p1:SUPP,p2:SUPP) ==
        -- p1 is non-zero, but has degree zero
        -- p2 is non-zero
        cp1:=leadingCoefficient p1
        one? cp1 => 1
        degree p2 = 0 => gcd(cp1,leadingCoefficient p2)::SUPP
        un?:=unit? cp1
        while not zero? p2 and not un? repeat
           cp1:=gcd(leadingCoefficient p2,cp1)
           un?:=unit? cp1
           p2:=reductum p2
        un? => 1
        cp1::SUPP
 
                     ---- Local  functions ----
--JHD    -- test if something wrong happened in the gcd
--JHD      failtest(f:P,p1:P,p2:P) : Boolean ==
--JHD        (p1 exquo f) case "failed" or (p2 exquo f) case "failed"
--JHD
--JHD    -- Choose the integers
--JHD      chooseVal(p1:P,p2:P,lvar:List OV):UTerm ==
--JHD        x:OV:=lvar.first
--JHD        lvr:=lvar.rest
--JHD        d1:=degree(p1,x)
--JHD        d2:=degree(p2,x)
--JHD        dd:NNI:=0$NNI
--JHD        nvr:NNI:=#lvr
--JHD        lval:List R :=[]
--JHD        range:I:=8
--JHD        for i in 1..  repeat
--JHD          range:=2*range
--JHD          lval:=[(random()$I rem (2*range) - range)::R  for i in 1..nvr]
--JHD          uf1:SUPR:=univariate eval(p1,lvr,lval)
--JHD          degree uf1 ~= d1 => "new point"
--JHD          uf2:SUPR:=univariate eval(p2,lvr,lval)
--JHD          degree uf2 ~= d2 => "new point"
--JHD          u:=gcd(uf1,uf2)
--JHD          du:=degree u
--JHD         --the univariate gcd is 1
--JHD          if du=0 then return [[1$SUPR],lval,0$P]$UTerm
--JHD
--JHD          ugcd:List SUPR:=[u,(uf1 exquo u)::SUPR,(uf2 exquo u)::SUPR]
--JHD          uterm:=[ugcd,lval,0$P]$UTerm
--JHD          dd=0 => dd:=du
--JHD
--JHD        --the degree is not changed
--JHD          du=dd =>
--JHD
--JHD           --test if one of the polynomials is the gcd
--JHD            dd=d1 =>
--JHD              if ^((f:=p2 exquo p1) case "failed") then
--JHD                return [[u],lval,p1]$UTerm
--JHD              if dd~=d2 then dd:=(dd-1)::NNI
--JHD
--JHD            dd=d2 =>
--JHD              if ^((f:=p1 exquo p2) case "failed") then
--JHD                return [[u],lval,p2]$UTerm
--JHD              dd:=(dd-1)::NNI
--JHD            return uterm
--JHD
--JHD         --the new gcd has degree less
--JHD          du<dd => dd:=du
--JHD
--JHD      good(f:P,lvr:List OV):Record(upol:SUPR,inval:List R) ==
--JHD        nvr:NNI:=#lvr
--JHD        range:I:=1
--JHD        ltry:List List R:=[]
--JHD        while true repeat
--JHD          range:=2*range
--JHD          lval:=[(random()$I rem (2*range) -range)::R  for i in 1..nvr]
--JHD          member?(lval,ltry) => "new point"
--JHD          ltry:=cons(lval,ltry)
--JHD          uf:=univariate eval(f,lvr,lval)
--JHD          if degree gcd(uf,differentiate uf)=0 then return [uf,lval]
--JHD
--JHD      -- impose the right lc
--JHD      imposelc(lipol:List SUPR,
--JHD               lvar:List OV,lval:List R,leadc:List P):List SUPR ==
--JHD        result:List SUPR :=[]
--JHD        lvar:=lvar.rest
--JHD        for pol in lipol for leadpol in leadc repeat
--JHD           p1:= univariate eval(leadpol,lvar,lval) * pol
--JHD           result:= cons((p1 exquo leadingCoefficient pol)::SUPR,result)
--JHD        reverse result
--JHD
--JHD    --Compute the gcd between not coprime polynomials
--JHD      notCoprime(g:P,p2:P,ldeg:List NNI,lvar:List OV) : P ==
--JHD        x:OV:=lvar.first
--JHD        lvar1:List OV:=lvar.rest
--JHD        lg1:=gcdcofact([g,differentiate(g,x)])
--JHD        g1:=lg1.1
--JHD        lg:LGcd:=localgcd(g1,p2,lvar)
--JHD        (l,lval):=(lg.locgcd,lg.goodint)
--JHD        p2:=(p2 exquo l)::P
--JHD        (gd1,gd2):=(l,l)
--JHD        ul:=univariate(eval(l,lvar1,lval))
--JHD        dl:=degree ul
--JHD        if degree gcd(ul,differentiate ul) ~=0 then
--JHD          newchoice:=good(l,lvar.rest)
--JHD          ul:=newchoice.upol
--JHD          lval:=newchoice.inval
--JHD        ug1:=univariate(eval(g1,lvar1,lval))
--JHD        ulist:=[ug1,univariate eval(p2,lvar1,lval)]
--JHD        lcpol:=[leadingCoefficient univariate(g1,x),
--JHD                leadingCoefficient univariate(p2,x)]
--JHD        while true repeat
--JHD          d:SUPR:=gcd(cons(ul,ulist))
--JHD          if degree d =0 then return gd1
--JHD          lquo:=(ul exquo d)::SUPR
--JHD          if degree lquo ~=0 then
--JHD            lgcd:=gcd(cons(leadingCoefficient univariate(l,x),lcpol))
--JHD            gd2:=lift(l,d,lquo,lgcd,lvar,ldeg,lval)
--JHD            l:=gd2
--JHD            ul:=univariate(eval(l,lvar1,lval))
--JHD            dl:=degree ul
--JHD          gd1:=gd1*gd2
--JHD          ulist:=[(uf exquo d)::SUPR for uf in ulist]
--JHD
--JHD -- we suppose that the poly have the same mainvar, deg p1<deg p2 and the
--JHD -- polys primitive
--JHD      internal(p1:P,p2:P,x:OV) : P ==
--JHD        lvar:List OV:=sort(#1>#2,setUnion(variables p1,variables p2))
--JHD        d1:=degree(p1,x)
--JHD        d2:=degree(p2,x)
--JHD        result: P:=localgcd(p1,p2,lvar).locgcd
--JHD        -- special cases
--JHD        result=1 => 1$P
--JHD        (dr:=degree(result,x))=d1 or dr=d2  => result
--JHD        while failtest(result,p1,p2) repeat
--JHD          SAY$Lisp  "retrying gcd"
--JHD          result:=localgcd(p1,p2,lvar).locgcd
--JHD        result
--JHD
--JHD    --local function for the gcd : it returns the evaluation point too
--JHD      localgcd(p1:P,p2:P,lvar:List(OV)) : LGcd ==
--JHD        x:OV:=lvar.first
--JHD        uterm:=chooseVal(p1,p2,lvar)
--JHD        listpol:= uterm.lpol
--JHD        ud:=listpol.first
--JHD        dd:= degree ud
--JHD
--JHD        --the univariate gcd is 1
--JHD        dd=0 => [1$P,uterm.lint]$LGcd
--JHD
--JHD        --one of the polynomials is the gcd
--JHD        dd=degree(p1,x) or dd=degree(p2,x) =>
--JHD                                           [uterm.mpol,uterm.lint]$LGcd
--JHD        ldeg:List NNI:=map(min,degree(p1,lvar),degree(p2,lvar))
--JHD
--JHD       -- if there is a polynomial g s.t. g/gcd and gcd are coprime ...
--JHD        -- I can lift
--JHD        (h:=lift?(p1,p2,uterm,ldeg,lvar)) case "failed" =>
--JHD          [notCoprime(p1,p2,ldeg,lvar),uterm.lint]$LGcd
--JHD        [h::P,uterm.lint]$LGcd
--JHD
--JHD
--JHD  -- content, internal functions return the poly if it is a monomial
--JHD      monomContent(p:P,var:OV):P ==
--JHD        ground? p => 1$P
--JHD        md:= minimumDegree(p,var)
--JHD        ((var::P)**md)*(gcd sort(better,coefficients univariate(p,var)))
 
      monomContentSup(u:SUPP):SUPP ==
        degree(u) = 0$NonNegativeInteger => 1$SUPP
        md:= minimumDegree u
        gcd(sort(better,coefficients u)) * monomial(1$P,md)$SUPP
 
--JHD  -- change the polynomials to have positive lc
--JHD      abs(p:P): P == unitNormal(p).canonical
 
  -- 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)
 
   -- PRS algorithm
   -- gcdPrs(p1:P,p2:P,d:NNI,var:OV):Union(P,"failed") ==
   --   u1:= univariate(p1,var)
   --   u2:= univariate(p2,var)
   --   finished:Boolean:= false
   --   until finished repeat
   --     dd:NNI:=(degree u1 - degree u2)::NNI
   --     lc1:SUPP:=leadingCoefficient u2 * reductum u1
   --     lc2:SUPP:=leadingCoefficient u1 * reductum u2
   --     u3:SUPP:= primitate((lc1-lc2)*monomial(1$P,dd))$%
   --     (d3:=degree(u3)) <= d => finished:= true
   --     u1:= u2
   --     u2:= u3
   --     if d3 > degree(u1) then (u1,u2):= (u2,u1)
   --   g:= (u2 exquo u3)
   --   g case SUPP => abs multivariate(u3,var)
   --   "failed"
 
  -- Gcd between polynomial p1 and p2 with
  -- mainVariable p1 < x=mainVariable p2
--JHD      gcdTermList(p1:P,p2:P) : P ==
--JHD        termList:=sort(better,
--JHD           cons(p1,coefficients univariate(p2,(mainVariable p2)::OV)))
--JHD        q:P:=termList.first
--JHD        for term in termList.rest until q = 1$P repeat q:= gcd(q,term)
--JHD        q
--JHD
--JHD  -- Gcd between polynomials with the same mainVariable
--JHD      gcdSameMainvar(p1:P,p2:P,mvar:OV): P ==
--JHD        if degree(p1,mvar) < degree(p2,mvar) then (p1,p2):= (p2,p1)
--JHD        (p1 exquo p2) case P => abs p2
--JHD        c1:= monomContent(p1,mvar)$%
--JHD        c1 = p1 => gcdMonom(p1,p2,mvar)
--JHD        c2:= monomContent(p2,mvar)$%
--JHD        c2 = p2 => gcdMonom(p2,p1,mvar)
--JHD        p1:= (p1 exquo c1)::P
--JHD        p2:= (p2 exquo c2)::P
--JHD        if degree(p1,mvar) < degree(p2,mvar) then (p1,p2):= (p2,p1)
--JHD        (p1 exquo p2) case P => abs(p2) * gcd(c1,c2)
--JHD        abs(gcdPrim(p1,p2,mvar)) * gcd(c1,c2)
--JHD
--JHD   --  make the polynomial primitive with respect to var
--JHD      primitate(p:P,var:OV):P == (p exquo monomContent(p,var))::P
--JHD
--JHD      primitate(u:SUPP):SUPP == (u exquo monomContentSup u)::SUPP
--JHD
--JHD   -- gcd between primitive polynomials with the same mainVariable
--JHD      gcdPrim(p1:P,p2:P,mvar:OV):P ==
--JHD        vars:= removeDuplicates append(variables p1,variables p2)
--JHD        #vars=1 => multivariate(gcd(univariate p1,univariate p2),mvar)
--JHD        vars:=delete(vars,position(mvar,vars))
--JHD        --d:= degModGcd(p1,p2,mvar,vars)
--JHD        --d case "failed" => internal(p2,p1,mvar)
--JHD        --deg:= d:NNI
--JHD        --deg = 0$NNI => 1$P
--JHD        --deg = degree(p1,mvar) =>
--JHD        --  (p2 exquo p1) case P => abs(p1)  -- already know that
--JHD                                           -- ^(p1 exquo p2)
--JHD        --  internal(p2,p1,mvar)
--JHD        --cheapPrs?(p1,p2,deg,mvar) =>
--JHD        --  g:= gcdPrs(p1,p2,deg,mvar)
--JHD        --  g case P => g::P
--JHD        --  internal(p2,p1,mvar)
--JHD        internal(p2,p1,mvar)
--JHD
--JHD   -- gcd between a monomial and a polynomial
--JHD      gcdMonom(m:P,p:P,var:OV):P ==
--JHD        ((var::P) ** min(minimumDegree(m,var),minimumDegree(p,var))) *
--JHD          gcdTermList(leadingCoefficient(univariate(m,var)),p)
--JHD
--JHD    --If there is a pol s.t. pol/gcd and gcd are coprime I can lift
--JHD      lift?(p1:P,p2:P,uterm:UTerm,ldeg:List NNI,
--JHD                     lvar:List OV) : Union("failed",P) ==
--JHD        x:OV:=lvar.first
--JHD        leadpol:Boolean:=false
--JHD        (listpol,lval):=(uterm.lpol,uterm.lint)
--JHD        d:=listpol.first
--JHD        listpol:=listpol.rest
--JHD        nolift:Boolean:=true
--JHD        for uf in listpol repeat
--JHD              --note uf and d not necessarily primitive
--JHD          degree gcd(uf,d) =0 => nolift:=false
--JHD        nolift => "failed"
--JHD        f:P:=([p1,p2]$List(P)).(position(uf,listpol))
--JHD        lgcd:=gcd(leadingCoefficient univariate(p1,x),
--JHD                  leadingCoefficient univariate(p2,x))
--JHD        lift(f,d,uf,lgcd,lvar,ldeg,lval)
--JHD
--JHD   -- interface with the general "lifting" function
--JHD      lift(f:P,d:SUPR,uf:SUPR,lgcd:P,lvar:List OV,
--JHD                        ldeg:List NNI,lval:List R):P ==
--JHD        x:OV:=lvar.first
--JHD        leadpol:Boolean:=false
--JHD        lcf:P
--JHD        lcf:=leadingCoefficient univariate(f,x)
--JHD        df:=degree(f,x)
--JHD        leadlist:List(P):=[]
--JHD
--JHD        if lgcd~=1$P then
--JHD          leadpol:=true
--JHD          f:=lgcd*f
--JHD          ldeg:=[n0+n1 for n0 in ldeg for n1 in degree(lgcd,lvar)]
--JHD          lcd:R:=leadingCoefficient d
--JHD          if ground? lgcd then d:=((retract lgcd) *d exquo lcd)::SUPR
--JHD          else d:=(retract(eval(lgcd,lvar.rest,lval)) * d exquo lcd)::SUPR
--JHD          uf:=lcd*uf
--JHD        leadlist:=[lgcd,lcf]
--JHD        lg:=imposelc([d,uf],lvar,lval,leadlist)
--JHD        plist:=lifting(univariate(f,x),lvar,lg,lval,leadlist,ldeg)::List P
--JHD        (p0:P,p1:P):=(plist.first,plist.2)
--JHD        if univariate eval(p0,rest lvar,lval) ~= lg.first then
--JHD           (p0,p1):=(p1,p0)
--JHD        ^leadpol => p0
--JHD        cprim:=contprim([p0])
--JHD        cprim.first.prim
--JHD
--JHD  -- Gcd for two multivariate polynomials
--JHD      gcd(p1:P,p2:P) : P ==
--JHD        (p1:= abs(p1)) = (p2:= abs(p2)) => p1
--JHD        ground? p1 =>
--JHD          p1 = 1$P => p1
--JHD          p1 = 0$P => p2
--JHD          ground? p2 => gcd((retract p1)@R,(retract p2)@R)::P
--JHD          gcdTermList(p1,p2)
--JHD        ground? p2 =>
--JHD          p2 = 1$P => p2
--JHD          p2 = 0$P => p1
--JHD          gcdTermList(p2,p1)
--JHD        mv1:= mainVariable(p1)::OV
--JHD        mv2:= mainVariable(p2)::OV
--JHD        mv1 = mv2 => gcdSameMainvar(p1,p2,mv1)
--JHD        mv1 < mv2 => gcdTermList(p1,p2)
--JHD        gcdTermList(p2,p1)
--JHD
--JHD  -- Gcd for a list of multivariate polynomials
--JHD      gcd(listp:List P) : P ==
--JHD        lf:=sort(better,listp)
--JHD        f:=lf.first
--JHD        for g in lf.rest repeat
--JHD          f:=gcd(f,g)
--JHD          if f=1$P then return f
--JHD        f
--JHD   -- Gcd and cofactors for a list of polynomials
--JHD      gcdcofact(listp : List P) : List P ==
--JHD         h:=gcd listp
--JHD         cons(h,[(f exquo h) :: P  for f in listp])
--JHD
--JHD   -- Gcd for primitive polynomials
--JHD      gcdprim(p1:P,p2:P):P ==
--JHD        (p1:= abs(p1)) = (p2:= abs(p2)) => p1
--JHD        ground? p1 =>
--JHD          ground? p2 => gcd((retract p1)@R,(retract p2)@R)::P
--JHD          p1 = 0$P => p2
--JHD          1$P
--JHD        ground? p2 =>
--JHD          p2 = 0$P => p1
--JHD          1$P
--JHD        mv1:= mainVariable(p1)::OV
--JHD        mv2:= mainVariable(p2)::OV
--JHD        mv1 = mv2 =>
--JHD          md:=min(minimumDegree(p1,mv1),minimumDegree(p2,mv1))
--JHD          mp:=1$P
--JHD          if md>1 then
--JHD            mp:=(mv1::P)**md
--JHD            p1:=(p1 exquo mp)::P
--JHD            p2:=(p2 exquo mp)::P
--JHD          mp*gcdPrim(p1,p2,mv1)
--JHD        1$P
--JHD
--JHD  -- Gcd for a list of primitive multivariate polynomials
--JHD      gcdprim(listp:List P) : P ==
--JHD        lf:=sort(better,listp)
--JHD        f:=lf.first
--JHD        for g in lf.rest repeat
--JHD          f:=gcdprim(f,g)
--JHD          if f=1$P then return f
--JHD        f
--JHD   -- Gcd and cofactors for a list of primitive polynomials
--JHD      gcdcofactprim(listp : List P) : List P ==
--JHD         h:=gcdprim listp
--JHD         cons(h,[(f exquo h) :: P  for f in listp])
--JHD
--JHD   -- content of a polynomial (with respect to its main var)
--JHD      content(f:P):P ==
--JHD        ground? f => f
--JHD        x:OV:=(mainVariable f)::OV
--JHD        gcd sort(better,coefficients univariate(f,x))
--JHD
--JHD   -- contents of a list of polynomials
--JHD      content(listf:List P) : List P == [content f for f in listf]
--JHD
--JHD   -- contents and primitive parts of a list  of polynomials
--JHD      contprim(listf:List P) : List ContPrim ==
--JHD        prelim :List P := content listf
--JHD        [[q,(f exquo q)::P]$ContPrim for q in prelim for f in listf]
--JHD

@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>
 
<<package GENPGCD GeneralPolynomialGcdPackage>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}