aboutsummaryrefslogtreecommitdiff
path: root/src/algebra/galfact.spad.pamphlet
diff options
context:
space:
mode:
authordos-reis <gdr@axiomatics.org>2007-08-14 05:14:52 +0000
committerdos-reis <gdr@axiomatics.org>2007-08-14 05:14:52 +0000
commitab8cc85adde879fb963c94d15675783f2cf4b183 (patch)
treec202482327f474583b750b2c45dedfc4e4312b1d /src/algebra/galfact.spad.pamphlet
downloadopen-axiom-ab8cc85adde879fb963c94d15675783f2cf4b183.tar.gz
Initial population.
Diffstat (limited to 'src/algebra/galfact.spad.pamphlet')
-rw-r--r--src/algebra/galfact.spad.pamphlet862
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}