path: root/src/algebra/twofact.spad.pamphlet

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra twofact.spad}
\author{Patrizia Gianni, James Davenport}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package NORMRETR NormRetractPackage}
<<package NORMRETR NormRetractPackage>>=
)abbrev package NORMRETR NormRetractPackage
++ Description:
++ This package \undocumented
NormRetractPackage(F, ExtF, SUEx, ExtP, n):C  == T where
  F          :   FiniteFieldCategory
  ExtF       :   FiniteAlgebraicExtensionField(F)
  SUEx       :   UnivariatePolynomialCategory ExtF
  ExtP       :   UnivariatePolynomialCategory SUEx
  n          :   PositiveInteger
  SUP       ==>  SparseUnivariatePolynomial
  R         ==>  SUP F
  P         ==>  SUP R

  C  ==> with
      normFactors : ExtP -> List ExtP
	++ normFactors(x) \undocumented
      retractIfCan : ExtP -> Union(P, "failed")
	++ retractIfCan(x) \undocumented
      Frobenius    : ExtP -> ExtP
	++ Frobenius(x) \undocumented

  T  ==> add

      normFactors(p:ExtP):List ExtP ==
          facs : List ExtP := [p]
          for i in 1..n-1 repeat 
             member?((p := Frobenius p), facs) => return facs
             facs := cons(p, facs)
          facs

      Frobenius(ff:ExtP):ExtP ==
         fft:ExtP:=0
         while ff~=0 repeat
           fft:=fft + monomial(map(Frobenius, leadingCoefficient ff),
                               degree ff)
           ff:=reductum ff
         fft

      retractIfCan(ff:ExtP):Union(P, "failed") ==          
         fft:P:=0
         while ff ~= 0 repeat
           lc : SUEx := leadingCoefficient ff
           plc: SUP F := 0
           while lc ~= 0 repeat
              lclc:ExtF := leadingCoefficient lc
              (retlc := retractIfCan lclc) case "failed" => return "failed"
              plc := plc + monomial(retlc::F, degree lc)
              lc := reductum lc
           fft:=fft+monomial(plc, degree ff)
           ff:=reductum ff
         fft

@
\section{package TWOFACT TwoFactorize}
<<package TWOFACT TwoFactorize>>=
)abbrev package TWOFACT TwoFactorize
++ Authors : P.Gianni, J.H.Davenport
++ Date Created : May 1990
++ Date Last Updated: March 1992
++ Description:
++ A basic package for the factorization of bivariate polynomials 
++ over a finite field.
++ The functions here represent the base step for the multivariate factorizer.
 
TwoFactorize(F) : C == T
 where
  F          :   FiniteFieldCategory
  SUP       ==>  SparseUnivariatePolynomial
  R         ==>  SUP F
  P         ==>  SUP R
  UPCF2     ==>  UnivariatePolynomialCategoryFunctions2
 
  C == with
    generalTwoFactor    : (P)  ->  Factored P 
      ++ generalTwoFactor(p) returns the factorisation of polynomial p,
      ++ a sparse univariate polynomial (sup) over a
      ++ sup over F.
    generalSqFr    : (P)  ->  Factored P 
      ++ generalSqFr(p) returns the square-free factorisation of polynomial p,
      ++ a sparse univariate polynomial (sup) over a
      ++ sup over F.
    twoFactor    : (P,Integer)  ->  Factored P 
      ++ twoFactor(p,n) returns the factorisation of polynomial p,
      ++ a sparse univariate polynomial (sup) over a
      ++ sup over F. 
      ++ Also, p is assumed primitive and square-free and n is the 
      ++ degree of the inner variable of p (maximum of the degrees
      ++ of the coefficients of p).
 
  T == add
    PI ==> PositiveInteger
    NNI ==> NonNegativeInteger
    import CommuteUnivariatePolynomialCategory(F,R,P)

                   ----  Local Functions  ----
    computeDegree  :  (P,Integer,Integer) -> PI
    exchangeVars   :           P          -> P
    exchangeVarTerm:        (R, NNI)      -> P
    pthRoot        :     (R, NNI, NNI)    -> R
 
    -- compute the degree of the extension to reduce the polynomial to a
    -- univariate one
    computeDegree(m : P,mx:Integer,q:Integer): PI ==
      my:=degree m
      n1:Integer:=length(10*mx*my)
      n2:Integer:=length(q)-1
      n:=(n1 quo n2)+1
      n::PI
--      n=1 => 1$PositiveInteger
--      (nextPrime(max(n,min(mx,my)))$IntegerPrimesPackage(Integer))::PI
 
    exchangeVars(p : P) : P ==
       p = 0 => 0
       exchangeVarTerm(leadingCoefficient p, degree p) +
          exchangeVars(reductum p)

    exchangeVarTerm(c:R, e:NNI) : P ==
       c = 0 => 0
       monomial(monomial(leadingCoefficient c, e)$R, degree c)$P + 
          exchangeVarTerm(reductum c, e)

    pthRoot(poly:R,p:NonNegativeInteger,PthRootPow:NonNegativeInteger):R ==
       tmp:=divideExponents(map((#1::F)**PthRootPow,poly),p)
       tmp case "failed" => error "consistency error in TwoFactor"
       tmp
 
    fUnion ==> Union("nil", "sqfr", "irred", "prime")
    FF     ==> Record(flg:fUnion, fctr:P, xpnt:Integer)

    generalSqFr(m:P): Factored P ==
       m = 0 => 0
       degree m = 0 =>
         l:=squareFree(leadingCoefficient m)
         makeFR(unit(l)::P,[[u.flg,u.fctr::P,u.xpnt] for u in factorList l])
       cont := content m
       m := (m exquo cont)::P
       sqfrm := squareFree m
       pfaclist : List FF := empty()
       unitPart := unit sqfrm
       for u in factorList sqfrm repeat
          u.flg = "nil" =>
             uexp:NNI:=(u.xpnt):NNI
             nfacs:=squareFree(exchangeVars u.fctr)
             for v in factorList nfacs repeat
                pfaclist:=cons([v.flg, exchangeVars v.fctr, v.xpnt*uexp],
                              pfaclist)
             unitPart := unit(nfacs)**uexp * unitPart
          pfaclist := cons(u,pfaclist)
       cont ~= 1 =>
           sqp := squareFree cont
           contlist:=[[w.flg,(w.fctr)::P,w.xpnt] for w in factorList sqp]
           pfaclist:= append(contlist, pfaclist)
           makeFR(unit(sqp)*unitPart,pfaclist)
       makeFR(unitPart,pfaclist)

        
    generalTwoFactor(m:P): Factored P ==
       m = 0 => 0
       degree m = 0 =>
         l:=factor(leadingCoefficient m)$DistinctDegreeFactorize(F,R)
         makeFR(unit(l)::P,[[u.flg,u.fctr::P,u.xpnt] for u in factorList l])
       ll:List FF
       ll:=[]
       unitPart:P
       cont:=content m
       if positive? degree(cont) then 
          m1:=m exquo cont
          m1 case "failed" => error "content doesn't divide"
          m:=m1
          contfact:=factor(cont)$DistinctDegreeFactorize(F,R)
          unitPart:=(unit contfact)::P
          ll:=[[w.flg,(w.fctr)::P,w.xpnt] for w in factorList contfact]
       else
          unitPart:=cont::P
       sqfrm:=squareFree m
       for u in factors sqfrm repeat
           expo:=u.exponent
           if negative? expo then error "negative exponent in a factorisation"
           expon:NonNegativeInteger:=expo::NonNegativeInteger
           fac:=u.factor
           degree fac = 1 => ll:=[["irred",fac,expon],:ll]
           differentiate fac = 0 =>      
              -- the polynomial is  inseparable w.r.t. its main variable
              map(differentiate,fac) = 0 =>
                p:=characteristic$F
                PthRootPow:=(size()$F exquo p)::NonNegativeInteger
                m1:=divideExponents(map(pthRoot(#1,p,PthRootPow),fac),p)
                m1 case "failed" => error "consistency error in TwoFactor"
                res:=generalTwoFactor m1
                unitPart:=unitPart*unit(res)**((p*expon)::NNI)
                ll:=[:[[v.flg,v.fctr,expon *p*v.xpnt] for v in factorList res],:ll]
              m2:=generalTwoFactor swap fac
              unitPart:=unitPart*unit(m2)**(expon::NNI)
              ll:=[:[[v.flg,swap v.fctr,expon*v.xpnt] for v in factorList m2],:ll]
           ydeg:="max"/[degree w for w in coefficients fac]
           twoF:=twoFactor(fac,ydeg)
           unitPart:=unitPart*unit(twoF)**expon
           ll:=[:[[v.flg,v.fctr,expon*v.xpnt] for v in factorList twoF],
                :ll]
       makeFR(unitPart,ll)
 
    -- factorization of a primitive square-free bivariate polynomial --
    twoFactor(m:P,dx:Integer):Factored P ==
       -- choose the degree for the extension
       n:PI:=computeDegree(m,dx,size()$F)
       -- extend the field
       -- find the substitution for x
       look:Boolean:=true
       dm:=degree m
       tryCount:Integer:=min(5,size()$F)
       i:Integer:=0
       lcm := leadingCoefficient m
       umv : R
       vval : F
       while look and i < tryCount repeat
          vval := random()$F
          i:=i+1
          zero? elt(lcm, vval) => "next value"
          umv := map(elt(#1,vval), m)$UPCF2(R, P, F, R)
          degree(gcd(umv,differentiate umv))~=0 => "next val"
          n := 1
          look := false
       extField:=FiniteFieldExtension(F,n)
       SUEx:=SUP extField
       TP:=SparseUnivariatePolynomial SUEx
       mm:TP:=0
       m1:=m
       while m1~=0 repeat
         mm:=mm+monomial(map(coerce,leadingCoefficient m1)$UPCF2(F,R,
                extField,SUEx),degree m1)
         m1:=reductum m1
       lcmm := leadingCoefficient mm
       val : extField
       umex : SUEx
       if not look then
          val := vval :: extField
          umex := map(coerce, umv)$UPCF2(F, R, extField, SUEx)
       while look repeat
         val:=random()$extField
         i:=i+1
         elt(lcmm,val)=0 => "next value"
         umex := map(elt(#1,val), mm)$UPCF2(SUEx, TP, extField, SUEx)
         degree(gcd(umex,differentiate umex))~=0 => "next val"
         look:=false
       prime:SUEx:=monomial(1,1)-monomial(val,0)
       fumex:=factor(umex)$DistinctDegreeFactorize(extField,SUEx)
       lfact1:=factors fumex

       #lfact1=1 => primeFactor(m,1)
       lfact : List TP :=
          [map(coerce,lf.factor)$UPCF2(extField,SUEx,SUEx,TP)
           for lf in lfact1]
       lfact:=cons(map(coerce,unit fumex)$UPCF2(extField,SUEx,SUEx,TP),
                   lfact)
       import GeneralHenselPackage(SUEx,TP)
       dx1:PI:=(dx+1)::PI
       lfacth:=completeHensel(mm,lfact,prime,dx1)
       lfactk: List P :=[]
       Normp := NormRetractPackage(F, extField, SUEx, TP, n)
      
       while not empty? lfacth repeat
         ff := first lfacth
         lfacth := rest lfacth
         if (c:=leadingCoefficient leadingCoefficient ff) ~=1 then
           ff:=((inv c)::SUEx)* ff
         not ((ffu:= retractIfCan(ff)$Normp) case "failed") =>
                    lfactk := cons(ffu::P, lfactk)
         normfacs := normFactors(ff)$Normp
         lfacth := [g for g in lfacth | not member?(g, normfacs)]
         ffn := */normfacs
         lfactk:=cons(retractIfCan(ffn)$Normp :: P, lfactk)
       */[primeFactor(ff1,1) for ff1 in lfactk]

@
\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 NORMRETR NormRetractPackage>>
<<package TWOFACT TwoFactorize>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}