listgcd.spad

)abbrev package HEUGCD HeuGcd
++ Author: P.Gianni
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This package provides the functions for the heuristic integer gcd.
++ Geddes's algorithm, for univariate polynomials with integer coefficients
HeuGcd (BP) : C == T
 where
  BP       :   UnivariatePolynomialCategory Integer
  Z        ==> Integer
  ContPrim ==> Record(cont : Z, prim : BP)


  C == with
     gcd          : List BP  -> BP
       ++ gcd([f1, .., fk]) = gcd of the polynomials fi.
     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

  T == add

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

    --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)
    bigPrime := 4010021

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

    constNotZero(f : BP ) : Boolean == (degree f = 0) and not (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

    smallest_to_front(l : List Z) : List Z ==
        ress := []
        smallest := first(l)
        for p in rest l repeat
            if p < smallest then
                ress := cons(smallest, ress)
                smallest := p
            else
                ress := cons(p, ress)
        cons(smallest, ress)

    --gcd of a list of integers
    lintgcd(lval : List(Z)) : Z ==
      empty? lval => 0$Z
      member?(1, lval) => 1$Z
      lval := smallest_to_front(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| not zero? f]
      answr := 2+2*hgt
      minf := "mindegpol"/[f for f in listf| not 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])

    PolyV ==> U32VectorPolynomialOperations

    --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
          fa := to_mod_pa(makeSUP(f), bigPrime)$PolyV
          for g in rest listf  repeat
              lcg := leadingCoefficient(g)
              if positiveRemainder(lcg, bigPrime) = 0 then
                  break
              ga := to_mod_pa(makeSUP(g), bigPrime)$PolyV
              fa := gcd(fa, ga, bigPrime)$PolyV
              if degree(fa)$PolyV = 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)

--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.