\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra galfact.spad}
\author{Frederic Lehobey}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package GALFACT GaloisGroupFactorizer}
<<package GALFACT GaloisGroupFactorizer>>=
)abbrev package GALFACT GaloisGroupFactorizer
++ Author: Frederic Lehobey
++ Date Created: 28 June 1994
++ Date Last Updated: 11 July 1997
++ Basic Operations: factor
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: factorization
++ Examples:
++ References:
++ [1] Bernard Beauzamy, Vilmar Trevisan and Paul S. Wang, Polynomial 
++ Factorization: Sharp Bounds, Efficient Algorithms,
++ J. Symbolic Computation (1993) 15, 393-413
++ [2] John Brillhart, Note on Irreducibility Testing,
++ Mathematics of Computation, vol. 35, num. 35, Oct. 1980, 1379-1381
++ [3] David R. Musser, On the Efficiency of a Polynomial Irreducibility Test,
++ Journal of the ACM, Vol. 25, No. 2, April 1978, pp. 271-282
++ Description: \spadtype{GaloisGroupFactorizer} provides functions
++ to factor resolvents.
-- improvements to do :
--   + reformulate the lifting problem in completeFactor -- See [1] (hard)
--   + implement algorithm RC -- See [1] (easy)
--   + use Dedekind's criterion to prove sometimes irreducibility (easy)
--     or even to improve early detection of true factors (hard)
--   + replace Sets by Bits
GaloisGroupFactorizer(UP): Exports == Implementation where
  Z ==> Integer
  UP: UnivariatePolynomialCategory Z
  N ==> NonNegativeInteger
  P ==> PositiveInteger
  CYC ==> CyclotomicPolynomialPackage()
  SUPZ ==> SparseUnivariatePolynomial Z

  ParFact ==> Record(irr: UP, pow: Z)
  FinalFact ==> Record(contp: Z, factors: List ParFact)
  DDRecord ==> Record(factor: UP, degree: Z) -- a Distinct-Degree factor
  DDList ==> List DDRecord
  MFact ==> Record(prime: Z,factors: List UP) -- Modular Factors
  LR ==> Record(left: UP, right: UP) -- Functional decomposition

  Exports ==> with
    makeFR: FinalFact -> Factored UP
      ++ makeFR(flist) turns the final factorization of henselFact into a
      ++ \spadtype{Factored} object.
    degreePartition: DDList -> Multiset N
      ++ degreePartition(ddfactorization) returns the degree partition of
      ++ the polynomial f modulo p where ddfactorization is the distinct
      ++ degree factorization of f computed by 
      ++ \spadfunFrom{ddFact}{ModularDistinctDegreeFactorizer}
      ++ for some prime p.
    musserTrials: () -> P
      ++ musserTrials() returns the number of primes that are tried in
      ++ \spadfun{modularFactor}.
    musserTrials: P -> P
      ++ musserTrials(n) sets to n the number of primes to be tried in
      ++ \spadfun{modularFactor} and returns the previous value.
    stopMusserTrials: () -> P
      ++ stopMusserTrials() returns the bound on the number of factors for
      ++ which \spadfun{modularFactor} stops to look for an other prime. You
      ++ will have to remember that the step of recombining the extraneous
      ++ factors may take up to \spad{2**stopMusserTrials()} trials. 
    stopMusserTrials: P -> P
      ++ stopMusserTrials(n) sets to n the bound on the number of factors for
      ++ which \spadfun{modularFactor} stops to look for an other prime. You
      ++ will have to remember that the step of recombining the extraneous
      ++ factors may take up to \spad{2**n} trials. Returns the previous
      ++ value.
    numberOfFactors: DDList -> N
      ++ numberOfFactors(ddfactorization) returns the number of factors of 
      ++ the polynomial f modulo p where ddfactorization is the distinct
      ++ degree factorization of f computed by 
      ++ \spadfunFrom{ddFact}{ModularDistinctDegreeFactorizer}
      ++ for some prime p.
    modularFactor: UP -> MFact
      ++ modularFactor(f) chooses a "good" prime and returns the factorization
      ++ of f modulo this prime in a form that may be used by
      ++ \spadfunFrom{completeHensel}{GeneralHenselPackage}. If prime is zero
      ++ it means that f has been proved to be irreducible over the integers
      ++ or that f is a unit (i.e. 1 or -1).
      ++ f shall be primitive (i.e. content(p)=1) and square free (i.e.
      ++ without repeated factors).
    useSingleFactorBound?: () -> Boolean
      ++ useSingleFactorBound?() returns \spad{true} if algorithm with single
      ++ factor bound is used for factorization, \spad{false} for algorithm
      ++ with overall bound.
    useSingleFactorBound: Boolean -> Boolean
      ++ useSingleFactorBound(b) chooses the algorithm to be used by the
      ++ factorizers: \spad{true} for algorithm with single
      ++ factor bound, \spad{false} for algorithm with overall bound.
      ++ Returns the previous value.
    useEisensteinCriterion?: () -> Boolean
      ++ useEisensteinCriterion?() returns \spad{true} if factorizers
      ++ check Eisenstein's criterion before factoring.
    useEisensteinCriterion: Boolean -> Boolean
      ++ useEisensteinCriterion(b) chooses whether factorizers check
      ++ Eisenstein's criterion before factoring: \spad{true} for
      ++ using it, \spad{false} else. Returns the previous value.
    eisensteinIrreducible?: UP -> Boolean
      ++ eisensteinIrreducible?(p) returns \spad{true} if p can be
      ++ shown to be irreducible by Eisenstein's criterion,
      ++ \spad{false} is inconclusive.
    tryFunctionalDecomposition?: () -> Boolean
      ++ tryFunctionalDecomposition?() returns \spad{true} if
      ++ factorizers try functional decomposition of polynomials before
      ++ factoring them.
    tryFunctionalDecomposition: Boolean -> Boolean
      ++ tryFunctionalDecomposition(b) chooses whether factorizers have
      ++ to look for functional decomposition of polynomials
      ++ (\spad{true}) or not (\spad{false}). Returns the previous value.
    factor: UP -> Factored UP
      ++ factor(p) returns the factorization of p over the integers.
    factor: (UP,N) -> Factored UP
      ++ factor(p,r) factorizes the polynomial p using the single factor bound
      ++ algorithm and knowing that p has at least r factors.
    factor: (UP,List N) -> Factored UP
      ++ factor(p,listOfDegrees) factorizes the polynomial p using the single
      ++ factor bound algorithm and knowing that p has for possible 
      ++ splitting of its degree listOfDegrees.
    factor: (UP,List N,N) -> Factored UP
      ++ factor(p,listOfDegrees,r) factorizes the polynomial p using the single
      ++ factor bound algorithm, knowing that p has for possible 
      ++ splitting of its degree listOfDegrees and that p has at least r
      ++ factors.
    factor: (UP,N,N) -> Factored UP
      ++ factor(p,d,r) factorizes the polynomial p using the single
      ++ factor bound algorithm, knowing that d divides the degree of all 
      ++ factors of p and that p has at least r factors.
    factorSquareFree: UP -> Factored UP
      ++ factorSquareFree(p) returns the factorization of p which is supposed
      ++ not having any repeated factor (this is not checked).
    factorSquareFree: (UP,N) -> Factored UP
      ++ factorSquareFree(p,r) factorizes the polynomial p using the single
      ++ factor bound algorithm and knowing that p has at least r factors.
      ++ f is supposed not having any repeated factor (this is not checked).
    factorSquareFree: (UP,List N) -> Factored UP
      ++ factorSquareFree(p,listOfDegrees) factorizes the polynomial p using
      ++ the single factor bound algorithm and knowing that p has for possible 
      ++ splitting of its degree listOfDegrees.
      ++ f is supposed not having any repeated factor (this is not checked).
    factorSquareFree: (UP,List N,N) -> Factored UP
      ++ factorSquareFree(p,listOfDegrees,r) factorizes the polynomial p using
      ++ the single factor bound algorithm, knowing that p has for possible 
      ++ splitting of its degree listOfDegrees and that p has at least r
      ++ factors.
      ++ f is supposed not having any repeated factor (this is not checked).
    factorSquareFree: (UP,N,N) -> Factored UP
      ++ factorSquareFree(p,d,r) factorizes the polynomial p using the single
      ++ factor bound algorithm, knowing that d divides the degree of all 
      ++ factors of p and that p has at least r factors.
      ++ f is supposed not having any repeated factor (this is not checked).
    factorOfDegree: (P,UP) -> Union(UP,"failed")
      ++ factorOfDegree(d,p) returns a factor of p of degree d.
    factorOfDegree: (P,UP,N) -> Union(UP,"failed")
      ++ factorOfDegree(d,p,r) returns a factor of p of degree
      ++ d knowing that p has at least r factors.
    factorOfDegree: (P,UP,List N) -> Union(UP,"failed")
      ++ factorOfDegree(d,p,listOfDegrees) returns a factor 
      ++ of p of degree d knowing that p has for possible splitting of its
      ++ degree listOfDegrees.
    factorOfDegree: (P,UP,List N,N) -> Union(UP,"failed")
      ++ factorOfDegree(d,p,listOfDegrees,r) returns a factor 
      ++ of p of degree d knowing that p has for possible splitting of its
      ++ degree listOfDegrees, and that p has at least r factors.
    factorOfDegree: (P,UP,List N,N,Boolean) -> Union(UP,"failed")
      ++ factorOfDegree(d,p,listOfDegrees,r,sqf) returns a
      ++ factor of p of degree d knowing that p has for possible splitting of
      ++ its degree listOfDegrees, and that p has at least r factors.
      ++ If \spad{sqf=true} the polynomial is assumed to be square free (i.e. 
      ++ without repeated factors).
    henselFact: (UP,Boolean) -> FinalFact
      ++ henselFact(p,sqf) returns the factorization of p, the result
      ++ is a Record such that \spad{contp=}content p,
      ++ \spad{factors=}List of irreducible factors of p with exponent.
      ++ If \spad{sqf=true} the polynomial is assumed to be square free (i.e. 
      ++ without repeated factors).
    btwFact: (UP,Boolean,Set N,N) -> FinalFact
      ++ btwFact(p,sqf,pd,r) returns the factorization of p, the result
      ++ is a Record such that \spad{contp=}content p,
      ++ \spad{factors=}List of irreducible factors of p with exponent.
      ++ If \spad{sqf=true} the polynomial is assumed to be square free (i.e. 
      ++ without repeated factors).
      ++ pd is the \spadtype{Set} of possible degrees. r is a lower bound for
      ++ the number of factors of p. Please do not use this function in your
      ++ code because its design may change.

  Implementation ==> add

    fUnion ==> Union("nil", "sqfr", "irred", "prime")
    FFE ==> Record(flg:fUnion, fctr:UP, xpnt:Z) -- Flag-Factor-Exponent
    DDFact ==> Record(prime:Z, ddfactors:DDList) -- Distinct Degree Factors
    HLR ==> Record(plist:List UP, modulo:Z) -- HenselLift Record

    mussertrials: P := 5
    stopmussertrials: P := 8
    usesinglefactorbound: Boolean := true
    tryfunctionaldecomposition: Boolean := true
    useeisensteincriterion: Boolean := true

    useEisensteinCriterion?():Boolean == useeisensteincriterion

    useEisensteinCriterion(b:Boolean):Boolean ==
      (useeisensteincriterion,b) := (b,useeisensteincriterion)
      b

    tryFunctionalDecomposition?():Boolean == tryfunctionaldecomposition

    tryFunctionalDecomposition(b:Boolean):Boolean ==
      (tryfunctionaldecomposition,b) := (b,tryfunctionaldecomposition)
      b

    useSingleFactorBound?():Boolean == usesinglefactorbound

    useSingleFactorBound(b:Boolean):Boolean ==
      (usesinglefactorbound,b) := (b,usesinglefactorbound)
      b

    stopMusserTrials():P == stopmussertrials

    stopMusserTrials(n:P):P ==
      (stopmussertrials,n) := (n,stopmussertrials)
      n

    musserTrials():P == mussertrials

    musserTrials(n:P):P ==
      (mussertrials,n) := (n,mussertrials)
      n

    import GaloisGroupFactorizationUtilities(Z,UP,Float)

    import GaloisGroupPolynomialUtilities(Z,UP)

    import IntegerPrimesPackage(Z)
    import IntegerFactorizationPackage(Z)

    import ModularDistinctDegreeFactorizer(UP)

    eisensteinIrreducible?(f:UP):Boolean ==
      rf := reductum f
      c: Z := content rf
      zero? c => false
      unit? c => false
      lc := leadingCoefficient f
      tc := lc
      while not zero? rf repeat
        tc := leadingCoefficient rf
        rf := reductum rf
      for p in factors(factor c)$Factored(Z) repeat
        if (one? p.exponent) and (not zero? (lc rem p.factor)) and
         (not zero? (tc rem ((p.factor)**2))) then return true
      false

    numberOfFactors(ddlist:DDList):N ==
      n: N := 0
      d: Z := 0
      for dd in ddlist repeat
        n := n +
          zero? (d := degree(dd.factor)::Z) => 1
          (d quo dd.degree)::N
      n

    -- local function, returns the a Set of shifted elements
    shiftSet(s:Set N,shift:N):Set N == set [ e+shift for e in members s ]

    -- local function, returns the "reductum" of an Integer (as chain of bits)
    reductum(n:Z):Z == n-shift(1,length(n)-1)

    -- local function, returns an integer with level lowest bits set to 1
    seed(level:Z):Z == shift(1,level)-1

    -- local function, returns the next number (as a chain of bit) for
    -- factor reconciliation of a given level (which is the number of
    -- extraneaous factors involved) or "End of level" if not any
    nextRecNum(levels:N,level:Z,n:Z):Union("End of level",Z) ==
      if (l := length n)<levels then return(n+shift(1,l-1))
      (n=shift(seed(level),levels-level)) => "End of level"
      b: Z := 1
      lr : Z
      while ((l-b) = (lr := length(n := reductum n)))@Boolean repeat b := b+1
      reductum(n)+shift(seed(b+1),lr)

    -- local function, return the set of N, 0..n
    fullSet(n:N):Set N == set [ i for i in 0..n ]

    modularFactor(p:UP):MFact ==
      not one? abs(content(p)) => 
       error "modularFactor: the polynomial is not primitive."
      zero? (n := degree p) => [0,[p]]

      -- declarations --
      cprime: Z := 2
      trials: List DDFact := empty()
      d: Set N := fullSet(n)
      dirred: Set N := set [0,n]
      s: Set N := empty()
      ddlist: DDList := empty()
      degfact: N := 0
      nf: N := stopmussertrials+1

      -- Musser, see [3] --
      diffp := differentiate p
      for i: Z in 1..mussertrials | nf>stopmussertrials repeat
        -- test 1: cprime divides leading coefficient
        -- test 2: "bad" primes: (in future: use Dedekind's Criterion)
        while (zero? ((leadingCoefficient p) rem cprime)) or
         (not zero? degree gcd(p,diffp,cprime)) repeat
          cprime := nextPrime(cprime)
        ddlist := ddFact(p,cprime)
        -- degree compatibility: See [3] --
        s := set [0]
        for f in ddlist repeat
          degfact := f.degree::N
          if not zero? degfact then 
            for j in 1..(degree(f.factor) quo degfact) repeat
              s := union(s, shiftSet(s,degfact))
        trials := cons([cprime,ddlist]$DDFact,trials)
        d := intersect(d, s)
        d = dirred => return [0,[p]] -- p is irreducible
        cprime := nextPrime(cprime)
        nf := numberOfFactors ddlist

      -- choose the one with the smallest number of factors
      choice := first trials
      nfc := numberOfFactors(choice.ddfactors)
      for t in rest trials repeat
        nf := numberOfFactors(t.ddfactors)
        if nf<nfc or ((nf=nfc) and (t.prime>choice.prime)) then
          nfc := nf
          choice := t
      cprime := choice.prime
      -- HenselLift$GHENSEL expects the degree 0 factor first 
      [cprime,separateFactors(choice.ddfactors,cprime)]

    degreePartition(ddlist:DDList):Multiset N ==
      dp: Multiset N := empty()
      d: N := 0
      dd: N := 0
      for f in ddlist repeat
        zero? (d := degree(f.factor)) => dp := insert!(0,dp)
        dd := f.degree::N
        dp := insert!(dd,dp,d quo dd)
      dp

    import GeneralHenselPackage(Z,UP)
    import UnivariatePolynomialDecompositionPackage(Z,UP)
    import BrillhartTests(UP) -- See [2]

    -- local function, finds the factors of f primitive, square-free, with
    -- positive leading coefficient and non zero trailing coefficient,
    -- using the overall bound technique. If pdecomp is true then look
    -- for a functional decomposition of f.
    henselfact(f:UP,pdecomp:Boolean):List UP ==
      if brillhartIrreducible? f or
       (useeisensteincriterion => eisensteinIrreducible? f ; false)
      then return [f]
      cf: Union(LR,"failed")
      if pdecomp and tryfunctionaldecomposition then
        cf := monicDecomposeIfCan f
      else
        cf := "failed"
      cf case "failed" =>
        m := modularFactor f
        zero? (cprime := m.prime) => m.factors
        b: P := (2*leadingCoefficient(f)*beauzamyBound(f)) :: P
        completeHensel(f,m.factors,cprime,b)
      lrf := cf::LR
      "append"/[ henselfact(g(lrf.right),false) for g in
       henselfact(lrf.left,true) ]

    -- local function, returns the complete factorization of its arguments,
    -- using the single-factor bound technique 
    completeFactor(f:UP,lf:List UP,cprime:Z,pk:P,r:N,d:Set N):List UP ==
      lc := leadingCoefficient f
      f0 := coefficient(f,0)
      ltrue: List UP := empty()
      found? := true
      degf: N := 0
      degg: N := 0
      g0: Z := 0
      g: UP := 0
      rg: N := 0
      nb: Z := 0
      lg: List UP := empty()
      b: P := 1
      dg: Set N := empty()
      llg: HLR := [empty(),0]
      levels: N := #lf
      level: Z := 1
      ic: Union(Z,"End of level") := 0
      i: Z := 0
      while level<levels repeat
        -- try all possible factors with degree in d
        ic := seed(level)
        while ((not found?) and (ic case Z)) repeat
          i := ic::Z
          degg := 0
          g0 := 1 -- LC algorithm
          for j in 1..levels repeat
            if bit?(i,j-1) then
              degg := degg+degree lf.j
              g0 := g0*coefficient(lf.j,0) -- LC algorithm
          g0 := symmetricRemainder(lc*g0,pk) -- LC algorithm
          if member?(degg,d) and (((lc*f0) exquo g0) case Z) then 
            --                       LC algorithm
            g := lc::UP -- build the possible factor -- LC algorithm
            for j in 1..levels repeat if bit?(i,j-1) then g := g*lf.j
            g := primitivePart reduction(g,pk)
            f1 := f exquo g
            if f1 case UP then -- g is a true factor
              found? := true
              -- remove the factors of g from lf
              nb := 1
              for j in 1..levels repeat
                if bit?(i,j-1) then 
                  swap!(lf,j,nb)
                  nb := nb+1
              lg := lf
              lf := rest(lf,level::N)
              setrest!(rest(lg,(level-1)::N),empty()$List(UP))
              f := f1::UP
              lc := leadingCoefficient f
              f0 := coefficient(f,0)
              -- is g irreducible?
              dg := select(#1<=degg,d)
              if not(dg=set [0,degg]) then -- implies degg >= 2
                rg := max(2,r+level-levels)::N
                b := (2*leadingCoefficient(g)*singleFactorBound(g,rg)) :: P
                if b>pk and (not brillhartIrreducible?(g)) and
                  (useeisensteincriterion => not eisensteinIrreducible?(g) ;
                  true)
                then
                  -- g may be reducible
                  llg := HenselLift(g,lg,cprime,b)
                  gpk: P := (llg.modulo)::P 
                  -- In case exact factorisation has been reached by
                  -- HenselLift before coefficient bound.
                  if gpk<b then
                    lg := llg.plist
                  else
                    lg := completeFactor(g,llg.plist,cprime,gpk,rg,dg)
                else lg := [ g ] -- g irreducible
              else lg := [ g ] -- g irreducible
              ltrue := append(ltrue,lg)
              r := max(2,(r-#lg))::N
              degf := degree f
              d := select(#1<=degf,d)
              if degf<=1 then -- lf exhausted
                if one? degf then
                  ltrue := cons(f,ltrue)
                return ltrue -- 1st exit, all factors found
              else -- can we go on with the same pk?
                b := (2*lc*singleFactorBound(f,r)) :: P
                if b>pk then -- unlucky: no we can't
                  llg := HenselLift(f,lf,cprime,b) -- I should reformulate
                   -- the lifting probleme, but hadn't time for that.
                   -- In any case, such case should be quite exceptional.
                  lf := llg.plist
                  pk := (llg.modulo)::P
                  -- In case exact factorisation has been reached by
                  -- HenselLift before coefficient bound.
                  if pk<b then return append(lf,ltrue) -- 2nd exit
                  level := 1
          ic := nextRecNum(levels,level,i)
        if found? then
          levels := #lf
          found? := false
        if not (ic case Z) then level := level+1
      cons(f,ltrue) -- 3rd exit, the last factor was irreducible but not "true"

    -- local function, returns the set of elements "divided" by an integer
    divideSet(s:Set N, n:N):Set N ==
      l: List N := [ 0 ]
      for e in members s repeat
        if (ee := (e exquo n)$N) case N then l := cons(ee::N,l)
      set(l)

    -- Beauzamy-Trevisan-Wang FACTOR, see [1] with some refinements
    -- and some differences. f is assumed to be primitive, square-free
    -- and with positive leading coefficient. If pdecomp is true then
    -- look for a functional decomposition of f. 
    btwFactor(f:UP,d:Set N,r:N,pdecomp:Boolean):List UP ==
      df := degree f
      not (max(d) = df) => error "btwFact: Bad arguments"
      reverse?: Boolean := false
      negativelc?: Boolean := false

      (d = set [0,df]) => [ f ]
      if abs(coefficient(f,0))<abs(leadingCoefficient(f)) then
        f := reverse f
        reverse? := true
      brillhartIrreducible? f or 
       (useeisensteincriterion => eisensteinIrreducible?(f) ; false) =>
        if reverse? then [ reverse f ] else [ f ]
      if negative? leadingCoefficient(f) then
        f := -f
        negativelc? := true
      cf: Union(LR,"failed")
      if pdecomp and tryfunctionaldecomposition then
        cf := monicDecomposeIfCan f
      else
        cf := "failed"
      if cf case "failed" then
        m := modularFactor f
        zero? (cprime := m.prime) => 
          if reverse? then
            if negativelc? then return [ -reverse f ]
            else return [ reverse f ]
          else if negativelc? then return [ -f ]
               else return [ f ]
        if noLinearFactor? f then d := remove(1,d)
        lc := leadingCoefficient f
        f0 := coefficient(f,0)
        b: P := (2*lc*singleFactorBound(f,r)) :: P -- LC algorithm
        lm := HenselLift(f,m.factors,cprime,b)
        lf := lm.plist
        pk: P := (lm.modulo)::P
        if ground? first lf then lf := rest lf
        -- in case exact factorisation has been reached by HenselLift
        -- before coefficient bound
        if not pk < b then lf := completeFactor(f,lf,cprime,pk,r,d)
      else
        lrf := cf::LR
        dh := degree lrf.right
        lg := btwFactor(lrf.left,divideSet(d,dh),2,true)
        lf: List UP := empty()
        for i in 1..#lg repeat
          g := lg.i
          dgh := (degree g)*dh
          df := subtractIfCan(df,dgh)::N
          lfg := btwFactor(g(lrf.right),
           select(#1<=dgh,d),max(2,r-df)::N,false)
          lf := append(lf,lfg)
          r := max(2,r-#lfg)::N
      if reverse? then lf := [ reverse(fact) for fact in lf ]
      for i in 1..#lf repeat
        if negative? leadingCoefficient(lf.i) then lf.i := -lf.i
        -- because we assume f with positive leading coefficient
      lf

    makeFR(flist:FinalFact):Factored UP ==
      ctp := factor flist.contp
      fflist: List FFE := empty()
      for ff in flist.factors repeat
        fflist := cons(["prime", ff.irr, ff.pow]$FFE, fflist)
      for fc in factorList ctp repeat
        fflist := cons([fc.flg, fc.fctr::UP, fc.xpnt]$FFE, fflist)
      makeFR(unit(ctp)::UP, fflist)

    import IntegerRoots(Z)

    -- local function, factorizes a quadratic polynomial
    quadratic(p:UP):List UP ==
      a := leadingCoefficient p
      b := coefficient(p,1)
      d := b**2-4*a*coefficient(p,0)
      r := perfectSqrt(d)
      r case "failed" => [p]
      b := b+(r::Z)
      a := 2*a
      d := gcd(a,b)
      if not one? d then
        a := a quo d
        b := b quo d
      f: UP := monomial(a,1)+monomial(b,0)
      cons(f,[(p exquo f)::UP])

    isPowerOf2(n:Z): Boolean ==
       n = 1 => true
       qr: Record(quotient: Z, remainder: Z) := divide(n,2)
       qr.remainder = 1 => false
       isPowerOf2 qr.quotient

    subMinusX(supPol: SUPZ): UP ==
       minusX: SUPZ := monomial(-1,1)$SUPZ
       unmakeSUP(elt(supPol,minusX)$SUPZ)

    henselFact(f:UP, sqf:Boolean):FinalFact ==
      factorlist: List(ParFact) := empty()

      -- make m primitive
      c: Z := content f
      f := (f exquo c)::UP

      -- make the leading coefficient positive
      if negative? leadingCoefficient f then
        c := -c
        f := -f

      -- is x**d factor of f
      if positive?(d := minimumDegree f) then
        f := monicDivide(f,monomial(1,d)).quotient
        factorlist := [[monomial(1,1),d]$ParFact]

      d := degree f

      -- is f constant?
      zero? d => [c,factorlist]$FinalFact

      -- is f linear?
      one? d => [c,cons([f,1]$ParFact,factorlist)]$FinalFact

      lcPol: UP := leadingCoefficient(f) :: UP

      -- is f cyclotomic (x**n - 1)?
      -lcPol = reductum(f) =>    -- if true, both will = 1
        for fac in map(unmakeSUP(#1)$UP,
         cyclotomicDecomposition(d)$CYC)$ListFunctions2(SUPZ,UP) repeat 
          factorlist := cons([fac,1]$ParFact,factorlist)
        [c,factorlist]$FinalFact

      -- is f odd cyclotomic (x**(2*n+1) + 1)?
      odd?(d) and (lcPol = reductum(f)) =>
        for sfac in cyclotomicDecomposition(d)$CYC repeat
           fac := subMinusX sfac
           if negative? leadingCoefficient fac then fac := -fac
           factorlist := cons([fac,1]$ParFact,factorlist)
        [c,factorlist]$FinalFact

      -- is the poly of the form x**n + 1 with n a power of 2?
      -- if so, then irreducible
      isPowerOf2(d) and (lcPol = reductum(f)) =>
        factorlist := cons([f,1]$ParFact,factorlist)
        [c,factorlist]$FinalFact

      -- other special cases to implement...

      -- f is square-free :
      sqf => [c, append([[pf,1]$ParFact for pf in henselfact(f,true)],
       factorlist)]$FinalFact

      -- f is not square-free :
      sqfflist := factors squareFree f
      for sqfr in sqfflist repeat
        mult := sqfr.exponent
        sqff := sqfr.factor
        d := degree sqff
        one? d => factorlist := cons([sqff,mult]$ParFact,factorlist)
        d=2 =>
          factorlist := append([[pf,mult]$ParFact for pf in quadratic(sqff)],
           factorlist)
        factorlist := append([[pf,mult]$ParFact for pf in
         henselfact(sqff,true)],factorlist) 
      [c,factorlist]$FinalFact

    btwFact(f:UP, sqf:Boolean, fd:Set N, r:N):FinalFact ==
      d := degree f
      not(max(fd)=d) => error "btwFact: Bad arguments"
      factorlist: List(ParFact) := empty()

      -- make m primitive
      c: Z := content f
      f := (f exquo c)::UP

      -- make the leading coefficient positive
      if negative? leadingCoefficient f then
        c := -c
        f := -f

      -- is x**d factor of f
      if positive?(maxd := minimumDegree f) then
        f := monicDivide(f,monomial(1,maxd)).quotient
        factorlist := [[monomial(1,1),maxd]$ParFact]
        r := max(2,r-maxd)::N
        d := subtractIfCan(d,maxd)::N
        fd := select(#1<=d,fd)

      -- is f constant?
      zero? d => [c,factorlist]$FinalFact

      -- is f linear?
      one? d => [c,cons([f,1]$ParFact,factorlist)]$FinalFact

      lcPol: UP := leadingCoefficient(f) :: UP

      -- is f cyclotomic (x**n - 1)?
      -lcPol = reductum(f) =>    -- if true, both will = 1
        for fac in map(unmakeSUP(#1)$UP,
         cyclotomicDecomposition(d)$CYC)$ListFunctions2(SUPZ,UP) repeat 
          factorlist := cons([fac,1]$ParFact,factorlist)
        [c,factorlist]$FinalFact

      -- is f odd cyclotomic (x**(2*n+1) + 1)?
      odd?(d) and (lcPol = reductum(f)) =>
        for sfac in cyclotomicDecomposition(d)$CYC repeat
           fac := subMinusX sfac
           if negative? leadingCoefficient fac then fac := -fac
           factorlist := cons([fac,1]$ParFact,factorlist)
        [c,factorlist]$FinalFact

      -- is the poly of the form x**n + 1 with n a power of 2?
      -- if so, then irreducible
      isPowerOf2(d) and (lcPol = reductum(f)) =>
        factorlist := cons([f,1]$ParFact,factorlist)
        [c,factorlist]$FinalFact

      -- other special cases to implement...

      -- f is square-free :
      sqf => [c, append([[pf,1]$ParFact for pf in btwFactor(f,fd,r,true)],
       factorlist)]$FinalFact

      -- f is not square-free :
      sqfflist := factors squareFree(f)
      if one?(#(sqfflist)) then -- indeed f was a power of a square-free 
        r := max(r quo ((first sqfflist).exponent),2)::N
      else
        r := 2
      for sqfr in sqfflist repeat
        mult := sqfr.exponent
        sqff := sqfr.factor
        d := degree sqff
        one? d => 
          factorlist := cons([sqff,mult]$ParFact,factorlist)
          maxd := (max(fd)-mult)::N
          fd := select(#1<=maxd,fd)
        d=2 =>
          factorlist := append([[pf,mult]$ParFact for pf in quadratic(sqff)],
           factorlist)
          maxd := (max(fd)-2*mult)::N
          fd := select(#1<=maxd,fd)
        factorlist := append([[pf,mult]$ParFact for pf in 
         btwFactor(sqff,select(#1<=d,fd),r,true)],factorlist)
        maxd := (max(fd)-d*mult)::N
        fd := select(#1<=maxd,fd)
      [c,factorlist]$FinalFact

    factor(f:UP):Factored UP ==
      makeFR
        usesinglefactorbound => btwFact(f,false,fullSet(degree f),2)
        henselFact(f,false)

    -- local function, returns true if the sum of the elements of the list
    -- is not the degree.
    errorsum?(d:N,ld:List N):Boolean == not (d = +/ld)

    -- local function, turns list of degrees into a Set
    makeSet(ld:List N):Set N ==
      s := set [0]
      for d in ld repeat s := union(s,shiftSet(s,d))
      s
      
    factor(f:UP,ld:List N,r:N):Factored UP ==
      errorsum?(degree f,ld) => error "factor: Bad arguments"
      makeFR btwFact(f,false,makeSet(ld),r)

    factor(f:UP,r:N):Factored UP == makeFR btwFact(f,false,fullSet(degree f),r)
    
    factor(f:UP,ld:List N):Factored UP == factor(f,ld,2)

    factor(f:UP,d:N,r:N):Factored UP ==
      n := (degree f) exquo d
      n case "failed" => error "factor: Bad arguments"
      factor(f,new(n::N,d)$List(N),r)

    factorSquareFree(f:UP):Factored UP ==
      makeFR
        usesinglefactorbound => btwFact(f,true,fullSet(degree f),2)
        henselFact(f,true)

    factorSquareFree(f:UP,ld:List(N),r:N):Factored UP ==
      errorsum?(degree f,ld) => error "factorSquareFree: Bad arguments"
      makeFR btwFact(f,true,makeSet(ld),r)

    factorSquareFree(f:UP,r:N):Factored UP ==
      makeFR btwFact(f,true,fullSet(degree f),r)
    
    factorSquareFree(f:UP,ld:List N):Factored UP == factorSquareFree(f,ld,2)

    factorSquareFree(f:UP,d:N,r:N):Factored UP ==
      n := (degree f) exquo d
      n case "failed" => error "factorSquareFree: Bad arguments"
      factorSquareFree(f,new(n::N,d)$List(N),r)

    factorOfDegree(d:P,p:UP,ld:List N,r:N,sqf:Boolean):Union(UP,"failed") ==
      dp := degree p
      errorsum?(dp,ld) => error "factorOfDegree: Bad arguments"
      (one? (d::N)) and noLinearFactor?(p) => "failed"
      lf := btwFact(p,sqf,makeSet(ld),r).factors
      for f in lf repeat
        degree(f.irr)=d => return f.irr
      "failed"

    factorOfDegree(d:P,p:UP,ld:List N,r:N):Union(UP,"failed") ==
      factorOfDegree(d,p,ld,r,false)

    factorOfDegree(d:P,p:UP,r:N):Union(UP,"failed") ==
      factorOfDegree(d,p,new(degree p,1)$List(N),r,false)

    factorOfDegree(d:P,p:UP,ld:List N):Union(UP,"failed") ==
      factorOfDegree(d,p,ld,2,false)

    factorOfDegree(d:P,p:UP):Union(UP,"failed") ==
      factorOfDegree(d,p,new(degree p,1)$List(N),2,false)

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