Initial population.

1 files changed, 862 insertions, 0 deletions

diff --git a/src/algebra/galfact.spad.pamphlet b/src/algebra/galfact.spad.pamphlet
new file mode 100644
index 00000000..8d9a6a02
--- /dev/null
+++ b/src/algebra/galfact.spad.pamphlet
@@ -0,0 +1,862 @@
+\documentclass{article}
+\usepackage{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
+        if (p.exponent = 1) 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 parts 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
+      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)) => 
+      not (abs(content(p)) = 1) => 
+       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
+      i: Z
+
+      -- Musser, see [3] --
+      diffp := differentiate p
+      for i 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
+                if (degf = 1) 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 parts 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 leadingCoefficient(f)<0 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 leadingCoefficient(lf.i)<0 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
+      if not (d = 1) 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 leadingCoefficient f < 0 then
+        c := -c
+        f := -f
+
+      -- is x**d factor of f
+      if (d := minimumDegree f) > 0 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
+      (d = 1) => [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 leadingCoefficient fac < 0 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 = 1) => 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 leadingCoefficient f < 0 then
+        c := -c
+        f := -f
+
+      -- is x**d factor of f
+      if (maxd := minimumDegree f) > 0 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
+      (d = 1) => [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 leadingCoefficient fac < 0 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 
+      if ((#(sqfflist)) = 1) 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 => 
+        (d = 1) => 
+          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"
+      ((d::N) = 1) 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}