Initial population.

1 files changed, 691 insertions, 0 deletions

diff --git a/src/algebra/gpgcd.spad.pamphlet b/src/algebra/gpgcd.spad.pamphlet
new file mode 100644
index 00000000..bf758915
--- /dev/null
+++ b/src/algebra/gpgcd.spad.pamphlet
@@ -0,0 +1,691 @@
+\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
+        (g = 1) => 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
+        (g = 1) => 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
+           (g = 1) => 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 = 1) => 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
+         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
+        (cp1 = 1) => 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}