axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / GBF.spad

)abbrev package GBF GroebnerFactorizationPackage
++ Author: H. Michael Moeller, Johannes Grabmeier
++ Date Created: 24 August 1989
++ Date Last Updated: 01 January 1992
++ References:
++ Normxx Notes 13: How to Compute a Groebner Basis
++ Description:
++ \spadtype{GroebnerFactorizationPackage} provides the function
++ groebnerFactor" which uses the factorization routines of Axiom to
++ factor each polynomial under consideration while doing the groebner basis
++ algorithm. Then it writes the ideal as an intersection of ideals
++ determined by the irreducible factors. Note that the whole ring may
++ occur as well as other redundancies. We also use the fact, that from the
++ second factor on we can assume that the preceding factors are
++ not equal to 0 and we divide all polynomials under considerations
++ by the elements of this list of "nonZeroRestrictions".
++ The result is a list of groebner bases, whose union of solutions
++ of the corresponding systems of equations is the solution of
++ the system of equation corresponding to the input list.
++ The term ordering is determined by the polynomial type used.
++ Suggested types include
++ \spadtype{DistributedMultivariatePolynomial},
++ \spadtype{HomogeneousDistributedMultivariatePolynomial},
++ \spadtype{GeneralDistributedMultivariatePolynomial}.

GroebnerFactorizationPackage(Dom, Expon, VarSet, Dpol) : SIG == CODE where
  Dom : Join(EuclideanDomain,CharacteristicZero)
  Expon : OrderedAbelianMonoidSup
  VarSet : OrderedSet
  Dpol : PolynomialCategory(Dom, Expon, VarSet)

  MF       ==>   MultivariateFactorize(VarSet,Expon,Dom,Dpol)
  sugarPol ==>   Record(totdeg: NonNegativeInteger, pol : Dpol)
  critPair ==>   Record(lcmfij: Expon,totdeg: NonNegativeInteger, poli: Dpol, polj: Dpol )
  L        ==>   List
  B        ==>   Boolean
  NNI      ==>   NonNegativeInteger
  OUT      ==>   OutputForm

  SIG ==> with

   factorGroebnerBasis : L Dpol -> L L Dpol
     ++ factorGroebnerBasis(basis) checks whether the basis contains
     ++ reducible polynomials and uses these to split the basis.

   factorGroebnerBasis : (L Dpol, Boolean) -> L L Dpol
     ++ factorGroebnerBasis(basis,info) checks whether the basis contains
     ++ reducible polynomials and uses these to split the basis.
     ++ If argument info is true, information is printed about
     ++ partial results.

   groebnerFactorize : (L Dpol, L Dpol) -> L L Dpol
     ++ groebnerFactorize(listOfPolys, nonZeroRestrictions) returns
     ++ a list of groebner basis. The union of their solutions
     ++ is the solution of the system of equations given by listOfPolys
     ++ under the restriction that the polynomials of nonZeroRestrictions
     ++ don't vanish.
     ++ At each stage the polynomial p under consideration (either from
     ++ the given basis or obtained from a reduction of the next S-polynomial)
     ++ is factorized. For each irreducible factors of p, a
     ++ new createGroebnerBasis is started
     ++ doing the usual updates with the factor
     ++ in place of p.

   groebnerFactorize : (L Dpol, L Dpol, Boolean) -> L L Dpol
     ++ groebnerFactorize(listOfPolys, nonZeroRestrictions, info) returns
     ++ a list of groebner basis. The union of their solutions
     ++ is the solution of the system of equations given by listOfPolys
     ++ under the restriction that the polynomials of nonZeroRestrictions
     ++ don't vanish.
     ++ At each stage the polynomial p under consideration (either from
     ++ the given basis or obtained from a reduction of the next S-polynomial)
     ++ is factorized. For each irreducible factors of p a
     ++ new createGroebnerBasis is started 
     ++ doing the usual updates with the factor in place of p.
     ++ If argument info is true, information is printed about
     ++ partial results.

   groebnerFactorize : L Dpol -> L L Dpol
     ++ groebnerFactorize(listOfPolys) returns 
     ++ a list of groebner bases. The union of their solutions
     ++ is the solution of the system of equations given by listOfPolys.
     ++ At each stage the polynomial p under consideration (either from
     ++ the given basis or obtained from a reduction of the next S-polynomial)
     ++ is factorized. For each irreducible factors of p, a
     ++ new createGroebnerBasis is started 
     ++ doing the usual updates with the factor 
     ++ in place of p.
     ++
     ++X mfzn : SQMATRIX(6,DMP([x,y,z],Fraction INT)) := 
     ++X   [ [0,1,1,1,1,1], [1,0,1,8/3,x,8/3], [1,1,0,1,8/3,y], 
     ++X   [1,8/3,1,0,1,8/3], [1,x,8/3,1,0,1], [1,8/3,y,8/3,1,0] ]
     ++X eq := determinant mfzn
     ++X groebnerFactorize 
     ++X   [eq,eval(eq, [x,y,z],[y,z,x]), eval(eq,[x,y,z],[z,x,y])] 

   groebnerFactorize : (L Dpol, Boolean) -> L L Dpol
     ++ groebnerFactorize(listOfPolys, info) returns
     ++ a list of groebner bases. The union of their solutions
     ++ is the solution of the system of equations given by listOfPolys.
     ++ At each stage the polynomial p under consideration (either from
     ++ the given basis or obtained from a reduction of the next S-polynomial)
     ++ is factorized. For each irreducible factors of p, a
     ++ new createGroebnerBasis is started 
     ++ doing the usual updates with the factor
     ++ in place of p.
     ++ If info is true, information is printed about partial results.

  CODE ==> add

   import GroebnerInternalPackage(Dom,Expon,VarSet,Dpol)
   -- next to help compiler to choose correct signatures:
   info: Boolean
   -- signatures of local functions

   newPairs : (L sugarPol, Dpol) -> L critPair
     ++ newPairs(lp, p) constructs list of critical pairs from the list of
     ++ lp of input polynomials and a given further one p.
     ++ It uses criteria M and T to reduce the list.
   updateCritPairs : (L critPair, L critPair, Dpol) -> L critPair
     ++ updateCritPairs(lcP1,lcP2,p) applies criterion B to lcP1 using
     ++ p. Then this list is merged with lcP2.
   updateBasis : (L sugarPol, Dpol, NNI) -> L sugarPol
     ++ updateBasis(li,p,deg) every polynomial in li is dropped if
     ++ its leading term is a multiple of the leading term of p.
     ++ The result is this list enlarged by p.
   createGroebnerBases : (L sugarPol, L Dpol, L Dpol, L Dpol, L critPair,_
                          L L Dpol, Boolean) -> L L Dpol
     ++ createGroebnerBases(basis, redPols, nonZeroRestrictions, inputPolys,
     ++   lcP,listOfBases): This function is used to be called from
     ++ groebnerFactorize.
     ++ basis: part of a Groebner basis, computed so far
     ++ redPols: Polynomials from the ideal to be used for reducing,
     ++   we don't throw away polynomials
     ++ nonZeroRestrictions: polynomials not zero in the common zeros
     ++   of the polynomials in the final (Groebner) basis
     ++ inputPolys: assumed to be in descending order
     ++ lcP: list of critical pairs built from polynomials of the
     ++   actual basis
     ++ listOfBases: Collects the (Groebner) bases constructed by this
     ++   recursive algorithm at different stages.
     ++   we print info messages if info is true
   createAllFactors: Dpol -> L Dpol
     ++ factor reduced critpair polynomial

   -- implementation of local functions


   createGroebnerBases(basis, redPols, nonZeroRestrictions, inputPolys,_
       lcP, listOfBases, info) ==
     doSplitting? : B := false
     terminateWithBasis : B := false
     allReducedFactors : L Dpol := []
     nP : Dpol  -- actual polynomial under consideration
     p :  Dpol  -- next polynomial from input list
     h :  Dpol  -- next polynomial from critical pairs
     stopDividing : Boolean
     --    STEP 1   do the next polynomials until a splitting is possible
     -- In the first step we take the first polynomial of "inputPolys"
     -- if empty, from list of critical pairs "lcP" and do the following:
     -- Divide it, if possible, by the polynomials from "nonZeroRestrictions".
     -- We factorize it and reduce each irreducible  factor with respect to
     -- "basis". If 0$Dpol occurs in the list we update the list and continue
     -- with next polynomial.
     -- If there are at least two (irreducible) factors
     -- in the list of factors we finish STEP 1 and set a boolean variable
     -- to continue with STEP 2, the splitting step.
     -- If there is just one of it, we do the following:
     -- If it is 1$Dpol we stop the whole calculation and put
     -- [1$Dpol] into the listOfBases
     -- Otherwise we update the "basis" and the other lists and continue
     -- with next polynomial.

     while (not doSplitting?) and (not terminateWithBasis) repeat
       terminateWithBasis := (null inputPolys and null lcP)
       not terminateWithBasis =>  -- still polynomials left
         -- determine next polynomial "nP"
         nP :=
           not null inputPolys =>
             p := first inputPolys
             inputPolys := rest inputPolys
             -- we know that p is not equal to 0 or 1, but, although,
             -- the inputPolys and the basis are ordered, we cannot assume
             -- that p is reduced w.r.t. basis, as the ordering is only quasi
             -- and we could have equal leading terms, and due to factorization
             -- polynomials of smaller leading terms, hence reduce p first:
             hMonic redPol(p,redPols)
           -- now we have inputPolys empty and hence lcP is not empty:
           -- create S-Polynomial from first critical pair:
           h := sPol first lcP
           lcP := rest lcP
           hMonic redPol(h,redPols)

         nP = 1$Dpol =>
           basis := [[0,1$Dpol]$sugarPol]
           terminateWithBasis := true

         -- if "nP" ^= 0, then  we continue, otherwise we determine next "nP"
         nP ^= 0$Dpol =>
           -- now we divide "nP", if possible, by the polynomials
           -- from "nonZeroRestrictions"
           for q in nonZeroRestrictions repeat
             stopDividing := false
             until stopDividing repeat
               nPq := nP exquo q
               stopDividing := (nPq case "failed")
               if not stopDividing then nP := autoCoerce nPq
               stopDividing := stopDividing or zero? degree nP

           zero? degree nP =>
             basis := [[0,1$Dpol]$sugarPol]
             terminateWithBasis := true  -- doSplitting? is still false

           -- a careful analysis has to be done, when and whether the
           -- following reduction and case nP=1 is necessary

           nP := hMonic redPol(nP,redPols)
           zero? degree nP =>
             basis := [[0,1$Dpol]$sugarPol]
             terminateWithBasis := true  -- doSplitting? is still false

           -- if "nP" ^= 0, then  we continue, otherwise we determine next "nP"
           nP ^= 0$Dpol =>
             -- now we factorize "nP", which is not constant
             irreducibleFactors : L Dpol := createAllFactors(nP)
             -- if there are more than 1 factors we reduce them and split
             (doSplitting? := not null rest irreducibleFactors) =>
               -- and reduce and normalize the factors
               for fnP in irreducibleFactors repeat
                 fnP := hMonic redPol(fnP,redPols)
                 -- no factor reduces to 0, as then "fP" would have been
                 -- reduced to zero,
                 -- but 1 may occur, which we will drop in a later version.
                 allReducedFactors := cons(fnP, allReducedFactors)
               -- end of "for fnP in irreducibleFactors repeat"

               -- we want that the smaller factors are dealt with first
               allReducedFactors := reverse allReducedFactors
             -- now the case of exactly 1 factor, but certainly not
             -- further reducible with respect to "redPols"
             nP := first irreducibleFactors
             -- put "nP" into "basis" and update "lcP" and "redPols":
             lcP : L critPair := updateCritPairs(lcP,newPairs(basis,nP),nP)
             basis := updateBasis(basis,nP,virtualDegree nP)
             redPols := concat(redPols,nP)
     -- end of "while not doSplitting? and not terminateWithBasis repeat"

     --    STEP 2  splitting step

     doSplitting? =>
       for fnP in allReducedFactors repeat
         if fnP ^= 1$Dpol
           then
             newInputPolys : L Dpol  := _
               sort((x,y) +-> degree x > degree y ,cons(fnP,inputPolys))
             listOfBases := createGroebnerBases(basis, redPols, _
               nonZeroRestrictions,newInputPolys,lcP,listOfBases,info)
             -- update "nonZeroRestrictions"
             nonZeroRestrictions := cons(fnP,nonZeroRestrictions)
           else
             if info then
               messagePrint("we terminated with [1]")$OUT
             listOfBases := cons([1$Dpol],listOfBases)

       -- we finished with all the branches on one level and hence
       -- finished this call of createGroebnerBasis. Therefore
       -- we terminate with the actual "listOfBasis" as
       -- everything is done in the recursions
       listOfBases
     -- end of "doSplitting? =>"

     --    STEP 3 termination step

     --  we found a groebner basis and put it into the list "listOfBases"
     --  (auto)reduce each basis element modulo the others
     newBasis := 
       minGbasis(sort((x,y)+->degree x > degree y,[p.pol for p in basis]))
     -- now check whether the normalized basis again has reducible
     -- polynomials, in this case continue splitting!
     if info then
       messagePrint("we found a groebner basis and check whether it ")$OUT
       messagePrint("contains reducible polynomials")$OUT
       print(newBasis::OUT)$OUT
       -- here we should create an output form which is reusable by the system
       -- print(convert(newBasis::OUT)$InputForm :: OUT)$OUT
     removeDuplicates append(factorGroebnerBasis(newBasis, info), listOfBases)

   createAllFactors(p: Dpol) ==
     loF : L Dpol := [el.fctr for el in factorList factor(p)$MF]
     sort((x,y) +-> degree x < degree y, loF)
   newPairs(lp : L sugarPol,p : Dpol) ==
     totdegreeOfp : NNI := virtualDegree p
     -- next list lcP contains all critPair constructed from
     -- p and and the polynomials q in lp
     lcP: L critPair := _
       --[[sup(degree q, degreeOfp), q, p]$critPair for q in lp]
       [makeCrit(q, p, totdegreeOfp) for q in lp]
     -- application of the criteria to reduce the list lcP
     critMTonD1 sort(critpOrder,lcP)
   updateCritPairs(oldListOfcritPairs, newListOfcritPairs, p)==
     updatD (newListOfcritPairs, critBonD(p,oldListOfcritPairs))
   updateBasis(lp, p, deg) == updatF(p,deg,lp)

   -- exported functions

   factorGroebnerBasis basis == factorGroebnerBasis(basis, false)

   factorGroebnerBasis (basis, info) ==
     foundAReducible : Boolean := false
     for p in basis while not foundAReducible repeat
       -- we use fact that polynomials have content 1
       foundAReducible := 1 < #[el.fctr for el in factorList factor(p)$MF]
     not foundAReducible =>
       if info then 
         messagePrint(_
          "factorGroebnerBasis: no reducible polynomials in this basis")$OUT
       [basis]
     -- improve! Use the fact that the irreducible ones already
     -- build part of the basis, use the done factorizations, etc.
     if info then  messagePrint("factorGroebnerBasis:_
        we found reducible polynomials and continue splitting")$OUT
     createGroebnerBases([],[],[],basis,[],[],info)

   groebnerFactorize(basis, nonZeroRestrictions) ==
     groebnerFactorize(basis, nonZeroRestrictions, false)

   groebnerFactorize(basis, nonZeroRestrictions, info) ==
     basis = [] => [basis]
     basis := remove((x:Dpol):Boolean +->(x = 0$Dpol),basis)
     basis = [] => [[0$Dpol]]
     -- normalize all input polynomial
     basis := [hMonic p for p in basis]
     member?(1$Dpol,basis) => [[1$Dpol]]
     basis :=  sort((x,y) +-> degree x > degree y, basis)
     createGroebnerBases([],[],nonZeroRestrictions,basis,[],[],info)

   groebnerFactorize(basis) == groebnerFactorize(basis, [], false)

   groebnerFactorize(basis,info) == groebnerFactorize(basis, [], info)