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

)abbrev package HEUGCD HeuGcd
++ Author: P.Gianni
++ Date Last Updated: 13 September 94
++ Description:
++ This package provides the functions for the heuristic integer gcd.
++ Geddes's algorithm,for univariate polynomials with integer coefficients

HeuGcd(BP) : SIG == CODE where
  BP : UnivariatePolynomialCategory Integer

  Z ==> Integer
  ContPrim ==> Record(cont:Z,prim:BP)

  SIG ==> with

    gcd : List BP -> BP
      ++ gcd([f1,..,fk]) = gcd of the polynomials fi.
      ++
      ++X gcd([671*671*x^2-1,671*671*x^2+2*671*x+1])
      ++X gcd([7*x^2+1,(7*x^2+1)^2])

    gcdprim : List BP -> BP
      ++ gcdprim([f1,..,fk]) = gcd of k PRIMITIVE univariate polynomials

    gcdcofact : List BP -> List BP
      ++ gcdcofact([f1,..fk]) = gcd and cofactors of k univariate polynomials.

    gcdcofactprim: List BP -> List BP
      ++ gcdcofactprim([f1,..fk]) = gcd and cofactors of k
      ++ primitive polynomials.

    content : List BP -> List Z
      ++ content([f1,..,fk]) = content of a list of univariate polynonials

    lintgcd : List Z -> Z
      ++ lintgcd([a1,..,ak]) = gcd of a list of integers

  CODE ==> add

    PI    ==> PositiveInteger
    NNI   ==> NonNegativeInteger
    Cases ==> Union("gcdprim","gcd","gcdcofactprim","gcdcofact")
    import ModularDistinctDegreeFactorizer BP

    --local functions
    localgcd     :        List BP       -> List BP
    constNotZero :           BP         -> Boolean
    height       :           BP         -> PI
    genpoly      :         (Z,PI)       -> BP
    negShiftz    :         (Z,PI)       -> Z
    internal     :     (Cases,List BP ) -> List BP
    constcase    : (List NNI ,List BP ) -> List BP
    lincase      : (List NNI ,List BP ) -> List BP
    myNextPrime  :        ( Z , NNI )   -> Z

    bigPrime:= prevPrime(2**26)$IntegerPrimesPackage(Integer)

    myNextPrime(val:Z,bound:NNI) : Z == nextPrime(val)$IntegerPrimesPackage(Z)

    constNotZero(f : BP ) : Boolean == (degree f = 0) and ^(zero? f)

    negShiftz(n:Z,Modulus:PI):Z ==
      n < 0 => n:= n+Modulus
      n > (Modulus quo 2) => n-Modulus
      n

    --compute the height of a polynomial
    height(f:BP):PI ==
      k:PI:=1
      while f^=0 repeat
           k:=max(k,abs(leadingCoefficient(f)@Z)::PI)
           f:=reductum f
      k

    --reconstruct the polynomial from the value-adic representation of
    --dval.
    genpoly(dval:Z,value:PI):BP ==
      d:=0$BP
      val:=dval
      for i in 0..  while (val^=0) repeat
        val1:=negShiftz(val rem value,value)
        d:= d+monomial(val1,i)
        val:=(val-val1) quo value
      d

    --gcd of a list of integers
    lintgcd(lval:List(Z)):Z ==
      empty? lval => 0$Z
      member?(1,lval) => 1$Z
      lval:=sort((z1,z2) +-> z1<z2,lval)
      val:=lval.first
      for val1 in lval.rest while ^(val=1) repeat val:=gcd(val,val1)
      val

    --content for a list of univariate polynomials
    content(listf:List BP ):List(Z) ==
      [lintgcd coefficients f for f in listf]

    --content of a list of polynomials with the relative primitive parts
    contprim(listf:List BP ):List(ContPrim) ==
       [[c:=lintgcd coefficients f,(f exquo c)::BP]$ContPrim  for f in listf]

    -- one polynomial is constant, remark that they are primitive
    -- but listf can contain the zero polynomial
    constcase(listdeg:List NNI ,listf:List BP ): List BP  ==
      lind:=select(constNotZero,listf)
      empty? lind =>
        member?(1,listdeg) => lincase(listdeg,listf)
        localgcd listf
      or/[n>0 for n in listdeg] => cons(1$BP,listf)
      lclistf:List(Z):= [leadingCoefficient f for f in listf]
      d:=lintgcd(lclistf)
      d=1 =>  cons(1$BP,listf)
      cons(d::BP,[(lcf quo d)::BP for lcf in lclistf])

    testDivide(listf: List BP, g:BP):Union(List BP, "failed") ==
      result:List BP := []
      for f in listf repeat
        if (f1:=f exquo g) case "failed" then return "failed"
        result := cons(f1::BP,result)
      reverse!(result)

    --one polynomial is linear, remark that they are primitive
    lincase(listdeg:List NNI ,listf:List BP ):List BP  ==
      n:= position(1,listdeg)
      g:=listf.n
      result:=[g]
      for f in listf repeat
        if (f1:=f exquo g) case "failed" then return cons(1$BP,listf)
        result := cons(f1::BP,result)
      reverse(result)

    IMG := InnerModularGcd(Z,BP,67108859,myNextPrime)

    mindegpol(f:BP, g:BP):BP ==
      degree(g) < degree (f) => g
      f

    --local function for the gcd among n PRIMITIVE univariate polynomials
    localgcd(listf:List BP ):List BP  ==
      hgt:="min"/[height(f) for f in listf|^zero? f]
      answr:=2+2*hgt
      minf := "mindegpol"/[f for f in listf|^zero? f]
      (result := testDivide(listf, minf)) case List(BP) =>
           cons(minf, result::List BP)
      if degree minf < 100 then for k in 1..10 repeat
        listval:=[f answr for f in listf]
        dval:=lintgcd(listval)
        dd:=genpoly(dval,answr)
        contd:=content(dd)
        d:=(dd exquo contd)::BP
        result:List BP :=[d]
        flag : Boolean := true
        for f in listf while flag repeat
          (f1:=f exquo d) case "failed" => flag:=false
          result := cons (f1::BP,result)
        if flag then return reverse(result)
        nvalue:= answr*832040 quo 317811
        if ((nvalue + answr) rem 2) = 0 then nvalue:=nvalue+1
        answr:=nvalue::PI
      gg:=modularGcdPrimitive(listf)$IMG
      cons(gg,[(f exquo gg) :: BP for f in listf])

    --internal function:it evaluates the gcd and avoids duplication of
    --code.
    internal(flag:Cases,listf:List BP ):List BP  ==
      --special cases
      listf=[] => [1$BP]
      (nlf:=#listf)=1 => [first listf,1$BP]
      minpol:=1$BP
      -- extract a monomial gcd
      mdeg:= "min"/[minimumDegree f for f in listf]
      if mdeg>0 then
        minpol1:= monomial(1,mdeg)
        listf:= [(f exquo minpol1)::BP for f in listf]
        minpol:=minpol*minpol1
      -- make the polynomials primitive
      Cgcd:List(Z):=[]
      contgcd : Z := 1
      if (flag case "gcd") or (flag case "gcdcofact") then
        contlistf:List(ContPrim):=contprim(listf)
        Cgcd:= [term.cont for term in contlistf]
        contgcd:=lintgcd(Cgcd)
        listf:List BP :=[term.prim for term in contlistf]
        minpol:=contgcd*minpol
      listdeg:=[degree f for f in listf ]
      f:= first listf
      if positiveRemainder(leadingCoefficient(f), bigPrime) ~= 0 then
        for g in rest listf repeat
          lcg := leadingCoefficient(g)
          if positiveRemainder(lcg, bigPrime) = 0 then
            leave
          f:=gcd(f,g,bigPrime)
          if degree f = 0 then return cons(minpol,listf)
      ans:List BP :=
         --one polynomial is constant
         member?(0,listdeg) => constcase(listdeg,listf)
         --one polynomial is linear
         member?(1,listdeg) => lincase(listdeg,listf)
         localgcd(listf)
      (result,ans):=(first ans*minpol,rest ans)
      if (flag case "gcdcofact") then
        ans:= [(p quo contgcd)*q for p in Cgcd for q in ans]
      cons(result,ans)

    --gcd among n PRIMITIVE univariate polynomials
    gcdprim (listf:List BP ):BP == first internal("gcdprim",listf)

    --gcd and cofactors for n PRIMITIVE univariate polynomials
    gcdcofactprim(listf:List BP ):List BP  == internal("gcdcofactprim",listf)

    --gcd for n generic univariate polynomials.
    gcd(listf:List BP ): BP  ==  first internal("gcd",listf)

    --gcd and cofactors for n generic univariate polynomials.
    gcdcofact (listf:List BP ):List BP == internal("gcdcofact",listf)