\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra listgcd.spad}
\author{Patrizia Gianni}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package HEUGCD HeuGcd}
<<package HEUGCD HeuGcd>>=
)abbrev package HEUGCD HeuGcd
++ Author: P.Gianni
++ Date Created:
++ Date Last Updated: 13 September 94
++ 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")
    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 not (zero? f)

    negShiftz(n:Z,Modulus:PI):Z ==
      negative? n => n:= n+Modulus
      n > (Modulus quo 2) => n-Modulus
      n

    --compute the height of a polynomial
    height(f:BP):PI ==
      k:PI:=1
      while not zero? f 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 not zero? val 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(#1<#2,lval)
      val:=lval.first
      for val1 in lval.rest while not (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/[positive? n 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])

    --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 positive? mdeg 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
      for g in rest listf  repeat
        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)

@
\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 HEUGCD HeuGcd>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}