mfinfact.spad

)abbrev package MFINFACT MultFiniteFactorize
++ Author: P. Gianni
++ Date Created: Summer 1990
++ Basic Functions:
++ Related Constructors: PrimeField, FiniteField, Polynomial
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description: Package for factorization of multivariate polynomials
++ over finite fields.


MultFiniteFactorize(OV, E, F, PG) : C == T where
  F          :   FiniteFieldCategory
  OV         :   OrderedSet
  E          :   OrderedAbelianMonoidSup
  PG         :   PolynomialCategory(F, E, OV)
  SUP        ==> SparseUnivariatePolynomial
  R          ==> SUP F
  P          ==> SparseMultivariatePolynomial(R, OV)
  Z          ==> Integer
  NNI ==> NonNegativeInteger
  FFPOLY     ==> FiniteFieldPolynomialPackage(F)
  MParFact   ==> Record(irr : P, pow : NNI)
  MFinalFact ==> Record(contp : R, factors : List MParFact)
  SUParFact    ==> Record(irr : SUP P, pow : NNI)
  SUPFinalFact ==> Record(contp : R, factors : List SUParFact)

                 --  contp   =  content,
                 --  factors =  List of irreducible factors with exponent

  C == with

    factor      : PG     ->  Factored PG
      ++ factor(p) produces the complete factorization of the multivariate
      ++ polynomial p over a finite field.
    factor      : SUP PG     ->  Factored SUP PG
      ++ factor(p) produces the complete factorization of the multivariate
      ++ polynomial p over a finite field. p is represented as a univariate
      ++ polynomial with multivariate coefficients over a finite field.

  T == add

    import from LeadingCoefDetermination(OV, IndexedExponents OV, R, P)
    import from MultivariateLifting(IndexedExponents OV, OV, R, P)
    import from FactoringUtilities(IndexedExponents OV, OV, R, P)
    import from FactoringUtilities(E, OV, F, PG)
    import from GenExEuclid(R, SUP R)

    L         ==> List
    UPCF2     ==> UnivariatePolynomialCategoryFunctions2
    LeadFact  ==> Record(polfac : L P, correct : R, corrfact : L SUP R)
    ContPrim  ==> Record(cont : P, prim : P)
    ParFact   ==> Record(irr : SUP R, pow : NNI)
    FinalFact ==> Record(contp : R, factors : L ParFact)
    NewOrd    ==> Record(npol : SUP P, nvar : L OV, newdeg : L NNI)
    Valuf     ==> Record(inval : L L R, unvfact : L SUP R, lu : R, complead : L R)

                   ----  Local Functions  ----
    ran       :                   Z              -> R
    mFactor : (P, OV, Z) -> MFinalFact
    supFactor : (SUP P, OV, Z) -> SUPFinalFact
    mfconst : (SUP P, Z, L OV, L NNI) -> L SUP P
    mfpol : (SUP P, OV, Z, L OV, L NNI) -> L SUP P
    varChoose :           (P, L OV, L NNI)         -> NewOrd
    simplify : (P, OV, Z, L OV, L NNI) -> MFinalFact
    intChoose :        (SUP P, L OV, R, L P, L L R)  -> Valuf
    pretest   :         (P, NNI, L OV, L R)         -> FinalFact
    pushdcoef :                  PG              -> P
    pushdown  :                (PG, OV)           -> P
    pushupconst :               (R, OV)           -> PG
    pushup    :                 (P, OV)           -> PG
    norm      :               L SUP R            -> Integer
    constantCase :        (P, L MParFact)         -> MFinalFact
    pM          :             L SUP R            -> R
    intfact     :     (SUP P, L OV, L NNI, MFinalFact, L L R) -> L SUP P


    convertPUP(lfg : MFinalFact) : SUPFinalFact ==
      [lfg.contp, [[lff.irr ::SUP P, lff.pow]$SUParFact
                    for lff in lfg.factors]]$SUPFinalFact

    supFactor(um : SUP P, var : OV, dx : Z) : SUPFinalFact ==
        degree(um) = 0 => convertPUP(mFactor(ground um, var, dx))
        lvar:L OV := "setUnion"/[variables cf for cf in coefficients um]
        lcont : SUP P
        lf : L SUP P

        flead : SUPFinalFact := [0, empty()]$SUPFinalFact
        factorlist : L SUParFact := empty()

        mdeg := minimumDegree um     ---- is the Mindeg > 0? ----
        if mdeg>0 then
          f1 : SUP P := monomial(1, mdeg)
          um := (um exquo f1)::SUP P
          factorlist := cons([monomial(1, 1), mdeg], factorlist)
          if degree um = 0 then return
            lfg := convertPUP mFactor(ground um, var, dx)
            [lfg.contp, append(factorlist, lfg.factors)]


        om := map((p1 : P) : PG +-> pushup(p1, var), um
                 )$UPCF2(P, SUP P, PG, SUP PG)
        sqfacs := squareFree(om)
        lcont := map((p1 : PG) : P +-> pushdown(p1, var), unit(sqfacs)
                    )$UPCF2(PG, SUP PG, P, SUP P)

                                   ----   Factorize the content  ----
        if ground? lcont then
          flead := convertPUP constantCase(ground lcont, empty())
        else
          flead := supFactor(lcont, var, dx)

        factorlist := append(factorlist, flead.factors)

                                 ----  Make the polynomial square-free  ----
        sqqfact := [[map((p1 : PG) : P +-> pushdown(p1, var), ff.factor),
                      ff.exponent] for ff in factorList sqfacs]

                        ---  Factorize the primitive square-free terms ---
        for fact in sqqfact repeat
          ffactor : SUP P := fact.irr
          ffexp := fact.pow
          coefs := coefficients ffactor
          ldeg := ["max"/[degree(fc,xx) for fc in coefs] for xx in lvar]
          if ground?(leadingCoefficient ffactor) then
             lf := mfconst(ffactor, dx, lvar, ldeg)
          else lf := mfpol(ffactor, var, dx, lvar, ldeg)
          auxfl := [[lfp, ffexp]$SUParFact  for lfp in lf]
          factorlist := concat!(auxfl, factorlist)
        lcfacs := */[leadingCoefficient leadingCoefficient(f.irr)^(f.pow)
                             for f in factorlist]
        [(leadingCoefficient leadingCoefficient(um) exquo lcfacs)::R,
                       factorlist]$SUPFinalFact

    fUnion ==> Union("nil", "sqfr", "irred", "prime")
    supFF ==> Record(flag : fUnion, factor : SUP PG, exponent : NNI)

    sup_prime_factor(pol : SUP PG, pow : NNI) : supFF ==
        ["prime", pol, pow]$supFF

    factor_sup(m : R) : Factored SUP PG ==
        fm := sufactor(m)$UnivariateFiniteFieldFactorize(F, R)
        fr1 := [sup_prime_factor(map(c +-> c::PG, p.factor
                                 )$UPCF2(F, R, PG, SUP PG), p.exponent)
                        for p in fm.factors]
        makeFR(fm.cont::PG::SUP(PG), fr1)

    factor(um : SUP PG) : Factored SUP PG ==
        lv : List OV := variables um
        empty?(lv) =>
            factor_sup(map((p1 : PG) : F +-> ground(p1)@F, um
                          )$UPCF2(PG, SUP PG, F, R))
        ld := degree(um, lv)
        dx := "min"/ld
        var := lv.position(dx, ld)
        cm := map((p1 : PG) : P +-> pushdown(p1, var), um
                 )$UPCF2(PG, SUP PG, P, SUP P)
        flist := supFactor(cm, var, dx)
        fr1 := [sup_prime_factor(map((p1 : P) : PG +-> pushup(p1, var), u.irr
                                )$UPCF2(P, SUP P, PG, SUP PG),
                          u.pow) for u in flist.factors]
        makeFR(pushupconst(flist.contp, var)::SUP(PG), fr1)



    mFactor(m : P, var : OV, dx : Z) : MFinalFact ==
      ground?(m) => constantCase(m, empty())
      lvar : L OV := variables m
      lcont : P
      lf : L SUP P
      flead : MFinalFact := [1, empty()]$MFinalFact
      factorlist : L MParFact := empty()
                                  ---- is the Mindeg > 0? ----
      lmdeg := minimumDegree(m, lvar)
      or/[n > 0 for n in lmdeg] => simplify(m, var, dx, lvar, lmdeg)
                              ----  Make the polynomial square-free  ----
      om := pushup(m, var)
      sqfacs := squareFree(om)
      lcont := pushdown(unit(sqfacs), var)

                                  ----  Factorize the content  ----
      if ground? lcont then
        flead := constantCase(lcont, empty())
      else
        flead := mFactor(lcont, var, dx)
      factorlist := flead.factors
      sqqfact : List Record(factor : P, exponent : NonNegativeInteger)
      sqqfact := [[pushdown(ff.factor, var), ff.exponent]
                                              for ff in factorList sqfacs]
                       ---  Factorize the primitive square-free terms ---
      for fact in sqqfact repeat
        ffactor : P := fact.factor
        ffexp := fact.exponent
        ground? ffactor =>
          for lterm in constantCase(ffactor, empty()).factors repeat
            factorlist := cons([lterm.irr, lterm.pow * ffexp], factorlist)
        lvar := variables ffactor
        x : OV := lvar.1
        ldeg := degree(ffactor, lvar)
             ---  Is the polynomial linear in one of the variables ? ---
        member?(1, ldeg) =>
          x : OV := lvar.position(1, ldeg)
          lcont := gcd coefficients(univariate(ffactor, x))
          ffactor := (ffactor exquo lcont)::P
          factorlist := cons([ffactor, ffexp]$MParFact, factorlist)
          for lcterm in mFactor(lcont, var, dx).factors repeat
              factorlist := cons([lcterm.irr, lcterm.pow*ffexp], factorlist)

        varch := varChoose(ffactor, lvar, ldeg)
        um := varch.npol


        ldeg := ldeg.rest
        lvar := lvar.rest
        if varch.nvar.1 ~= x then
          lvar := varch.nvar
          x := lvar.1
          lvar := lvar.rest
          pc := gcd coefficients um
          if pc ~= 1 then
            um := (um exquo pc)::SUP P
            ffactor := multivariate(um, x)
            for lcterm in mFactor(pc, var, dx).factors repeat
              factorlist := cons([lcterm.irr, lcterm.pow*ffexp], factorlist)
          ldeg := degree(ffactor, lvar)

        -- should be unitNormal if unified, but for now it is easier
        lcum : F := leadingCoefficient leadingCoefficient
                leadingCoefficient um
        if lcum ~= 1  then
          um := ((inv lcum)::R::P) * um
          flead.contp := (lcum::R) *flead.contp

        if ground?(leadingCoefficient um)
        then lf := mfconst(um, dx, lvar, ldeg)
        else lf := mfpol(um, var, dx, lvar, ldeg)
        auxfl := [[multivariate(lfp, x), ffexp]$MParFact  for lfp in lf]
        factorlist := append(factorlist, auxfl)
      flead.factors := factorlist
      flead


    pM(lum : L SUP R) : R ==
      x := monomial(1, 1)$R
      for i in 1..size()$F repeat
         p := x + (index(i::PositiveInteger)$F) ::R
         testModulus(p, lum) => return p
      for e in 2.. repeat
          p :=  (createIrreduciblePoly(e::PositiveInteger))$FFPOLY
          testModulus(p, lum) => return p
          while not((q := nextIrreduciblePoly(p)$FFPOLY) case "failed") repeat
             p := q::SUP(F)
             if testModulus(p, lum)$GenExEuclid(R, SUP R) then return p

      ----  push x in the coefficient domain for a term ----
    pushdcoef(t : PG) : P ==
       map((f1 : F) : R +-> coerce(f1)$R, t)$MPolyCatFunctions2(OV, E,
                                           IndexedExponents OV, F, R, PG, P)


              ----  internal function, for testing bad cases  ----
    intfact(um : SUP P, lvar : L OV, ldeg : L NNI,
            tleadpol : MFinalFact, ltry : L L R) :  L SUP P ==
      polcase : Boolean := (not empty? tleadpol.factors )
      vfchoo : Valuf :=
        polcase =>
          leadpol : L P := [ff.irr for ff in tleadpol.factors]
          intChoose(um, lvar, tleadpol.contp, leadpol, ltry)
        intChoose(um, lvar, 1, empty(), empty())
      unifact : List SUP R := vfchoo.unvfact
      nfact : NNI := #unifact
      nfact = 1 => [um]
      ltry : L L R := vfchoo.inval
      lval : L R := first ltry
      dd := vfchoo.lu
      lpol : List P := empty()
      leadval : List R := empty()
      if polcase then
        leadval := vfchoo.complead
        distf := distFact(vfchoo.lu, unifact, tleadpol, leadval, lvar, lval)
        distf case "failed" =>
             return intfact(um, lvar, ldeg, tleadpol, ltry)
        dist := distf@LeadFact
          -- check the factorization of leading coefficient
        lpol := dist.polfac
        dd := dist.correct
        unifact := dist.corrfact
      if dd ~= 1 then
        unifact := [dd*unifact.i for i in 1..nfact]
        um := ((dd^(nfact-1)::NNI)::P)*um
      (ffin := lifting(um, lvar, unifact, lval, lpol, ldeg, pM(unifact)))
           case "failed" => intfact(um,lvar,ldeg,tleadpol,ltry)
      factfin : L SUP P := ffin
      if dd ~= 1 then
        factfin := [primitivePart ff  for ff in  factfin]
      factfin

-- the following functions are used to "push" x in the coefficient ring -
               ----  push back the variable  ----
    pushup(f : P, x : OV) : PG ==
       ground? f => pushupconst((retract f)@R, x)
       rr : PG := 0
       while f ~= 0 repeat
         lf := leadingMonomial f
         cf := pushupconst(leadingCoefficient f, x)
         lvf := variables lf
         rr := rr+monomial(cf, lvf, degree(lf, lvf))$PG
         f := reductum f
       rr

        ----  push x in the coefficient domain for a polynomial ----
    pushdown(g : PG, x : OV) : P ==
       ground? g => ((retract g)@F)::R::P
       rf : P := 0$P
       ug := univariate(g, x)
       while ug ~= 0 repeat
         cf := monomial(1, degree ug)$R
         rf := rf+cf*pushdcoef(leadingCoefficient ug)
         ug := reductum ug
       rf

      ----  push x back from the coefficient domain ----
    pushupconst(r : R, x : OV) : PG ==
       ground? r => (retract r)@F ::PG
       rr : PG := 0
       while r ~= 0 repeat
         rr := rr+monomial((leadingCoefficient r)::PG, x, degree r)$PG
         r := reductum r
       rr

    -- This function has to be added to Eucliden domain
    ran(k1 : Z) : R ==
      --if R case Integer then random()$R rem (2*k1)-k1
      --else
      +/[monomial(random()$F, i)$R for i in 0..k1]

              ---  Choose the variable of least degree  ---
    varChoose(m : P, lvar : L OV, ldeg : L NNI) : NewOrd ==
      k := "min"/ldeg
      k = degree(m, first lvar) =>
                             [univariate(m, first lvar), lvar, ldeg]$NewOrd
      i := position(k, ldeg)
      x : OV := lvar.i
      ldeg := cons(k, delete(ldeg, i))
      lvar := cons(x, delete(lvar, i))
      [univariate(m, x), lvar, ldeg]$NewOrd


    norm(lum : L SUP R) : Integer == "max"/[degree lup for lup in lum]

          ---  Choose the values to reduce to the univariate case  ---
    intChoose(um : SUP P, lvar : L OV, clc : R, plist : L P, ltry : L L R) : Valuf ==
      -- declarations
      degum : NNI := degree um
      nvar1 := #lvar
      range : NNI := 0
      unifact : L SUP R
      testp : Boolean :=             -- polynomial leading coefficient
        plist = empty() => false
        true
      leadcomp, leadcomp1 : L R
      leadcomp := leadcomp1 := empty()
      nfatt : NNI := degum+1
      lffc : R := 1
      lffc1 := lffc
      leadtest := true --- the lc test with polCase has to be performed
      int : L R := empty()
      base_success : Z := 0
      base_failure : Z := 0

   --  New sets of values are chosen until we find twice the
   --  same number of "univariate" factors:the set smaller in modulo is
   --  is chosen.
      while true repeat
       lval := [ ran(range) for i in 1..nvar1]
       member?(lval, ltry) => range := 1+range
       ltry := cons(lval, ltry)
       leadcomp1 := [retract eval(pol, lvar, lval) for pol in plist]
       testp and or/[unit? epl for epl in leadcomp1] => range := range+1
       newm : SUP R := completeEval(um, lvar, lval)
       degum ~= degree newm or minimumDegree newm ~= 0 => range := range+1
       lffc1 := content newm
       newm := (newm exquo lffc1)::SUP R
       testp and leadtest and not polCase(lffc1*clc, #plist, leadcomp1)
                           => range := range+1
       Dnewm := differentiate newm
       D2newm := map(differentiate, newm)
       degree(gcd [newm, Dnewm, D2newm]) ~= 0 => range := range+1
      -- if R has Integer then luniv := henselFact(newm, false)$
      -- else
       lcnm : F := 1
        -- should be unitNormal if unified, but for now it is easier
       if (lcnm := leadingCoefficient leadingCoefficient newm) ~= 1 then
         newm := ((inv lcnm)::R)*newm
       luniv := tryTwoFactor(newm)$TwoFactorize(F)
       nfl := factorList luniv
       nf := #nfl
       nf = 0 or nf > nfatt => "next values"
       if nf = 1 and nfl(1).flag case "nil" then
           base_failure := base_failure + 1
           skip := base_failure < 3*base_success + 3
           luniv :=
               skip => luniv
               generalTwoFactor(newm)$TwoFactorize(F)
       else
           base_success := base_success + 1
           skip := false
       skip => iterate
       lunivf := factorList luniv
       nf := #lunivf

       nf = 0 or nf>nfatt => "next values"      ---  pretest failed ---

                        --- the univariate polynomial is irreducible ---
       if nf = 1 then (unifact := [newm]; break)

       lffc1 := lcnm * retract(unit(luniv))@R * lffc1

   --  the new integer give the same number of factors
       nfatt = nf =>
       -- if this is the first univariate factorization with polCase=true
       -- or if the last factorization has smaller norm and satisfies
       -- polCase
         if leadtest or
           ((norm unifact > norm [ff.factor for ff in lunivf]) and
             (not testp or polCase(lffc1*clc, #plist, leadcomp1))) then
                unifact := [uf.factor for uf in lunivf]
                int := lval
                lffc := lffc1
                if testp then leadcomp := leadcomp1
         -- found it
         break

   --  the first univariate factorization, initialize
       nfatt > degum =>
         unifact := [uf.factor for uf in lunivf]
         lffc := lffc1
         if testp then leadcomp := leadcomp1
         int := lval
         leadtest := false
         nfatt := nf

       nfatt>nf =>  -- for the previous values there were more factors
         if testp then leadtest := not polCase(lffc*clc, #plist, leadcomp)
         else leadtest := false
         -- if polCase=true we can consider the univariate decomposition
         if not leadtest then
           unifact := [uf.factor for uf in lunivf]
           lffc := lffc1
           if testp then leadcomp := leadcomp1
           int := lval
         nfatt := nf
      [cons(int, ltry), unifact, lffc, leadcomp]$Valuf


    constantCase(m : P, factorlist : List MParFact) : MFinalFact ==
      lunm := sufactor((retract m)@R)$UnivariateFiniteFieldFactorize(F, R)
      [(lunm.cont)::R, append(factorlist,
           [[(pp.factor)::P, pp.exponent] for pp in lunm.factors])]$MFinalFact

                ----  The polynomial has mindeg>0   ----

    simplify(m : P, var : OV, dm : Z, lvar : L OV, lmdeg : L NNI
            ) : MFinalFact ==
      factorlist : L MParFact := empty()
      pol1 : P := 1$P
      for x in lvar repeat
        i := lmdeg.(position(x, lvar))
        i = 0 => "next value"
        pol1 := pol1*monomial(1$P, x, i)
        factorlist := cons([x::P, i]$MParFact, factorlist)
      m := (m exquo pol1)::P
      ground? m => constantCase(m, factorlist)
      flead := mFactor(m, var, dm)
      flead.factors := append(factorlist, flead.factors)
      flead

                ----  m square-free, primitive, lc constant  ----
    mfconst(um : SUP P, dm : Z, lvar : L OV, ldeg : L NNI) : L SUP P ==
      nsign : Boolean
      empty? lvar =>
          um1 : SUP R := map(ground,
              um)$UPCF2(P, SUP P, R, SUP R)
          lum := generalTwoFactor(um1)$TwoFactorize(F)
          [map(coerce, lumf.factor)$UPCF2(R, SUP R, P, SUP P)
                for lumf in factorList lum]
      intfact(um, lvar, ldeg, [0, empty()]$MFinalFact, empty())

              --- m is square-free, primitive, lc is a polynomial  ---
    mfpol(um : SUP P, var : OV, dm : Z, lvar : L OV, ldeg : L NNI
         ) : L SUP P ==
      dist : LeadFact
      tleadpol := mFactor(leadingCoefficient um, var, dm)
      intfact(um, lvar, ldeg, tleadpol, empty())

    pFF ==> Record(flag : fUnion, factor : PG, exponent : NonNegativeInteger)

    factor(m : PG) : Factored PG ==
       lv := variables m
       lv = empty() => makeFR(m, empty() )
    -- reduce to multivariate over SUP
       ld := [degree(m, x) for x in lv]
       dx := "min"/ld
       var := lv(position(dx, ld))
       cm := pushdown(m, var)
       flist := mFactor(cm, var, dx)
       makeFR(pushupconst(flist.contp, var), [
           ["prime", pushup(u.irr, var), u.pow]$pFF
                                                 for u in flist.factors])

--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.