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

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra multsqfr.spad}
\author{Patrizia Gianni}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package MULTSQFR MultivariateSquareFree}
<<package MULTSQFR MultivariateSquareFree>>=
)abbrev package MULTSQFR MultivariateSquareFree
++Author : P.Gianni
++ This package provides the functions for the computation of the square
++ free decomposition of a multivariate polynomial.
++ It uses the package GenExEuclid for the resolution of
++ the equation \spad{Af + Bg = h} and its generalization to n polynomials
++ over an integral domain and the package \spad{MultivariateLifting}
++ for the "multivariate" lifting.
 
MultivariateSquareFree (E,OV,R,P) : C == T where
 Z   ==> Integer
 NNI ==> NonNegativeInteger
 R   : EuclideanDomain
 OV  : OrderedSet
 E   : OrderedAbelianMonoidSup
 P   : PolynomialCategory(R,E,OV)
 SUP ==> SparseUnivariatePolynomial P
 
 BP          ==> SparseUnivariatePolynomial(R)
 fUnion      ==> Union("nil","sqfr","irred","prime")
 ffSUP       ==> Record(flg:fUnion,fctr:SUP,xpnt:Integer)
 ffP         ==> Record(flg:fUnion,fctr:P,xpnt:Integer)
 FFE         ==> Record(factor:BP,exponent:Z)
 FFEP        ==> Record(factor:P,exponent:Z)
 FFES        ==> Record(factor:SUP,exponent:Z)
 Choice      ==> Record(upol:BP,Lval:List(R),Lfact:List FFE,ctpol:R)
 squareForm  ==> Record(unitPart:P,suPart:List FFES)
 Twopol      ==> Record(pol:SUP,polval:BP)
 UPCF2       ==> UnivariatePolynomialCategoryFunctions2

 
 C == with
   squareFree      :      P     -> Factored P
     ++ squareFree(p) computes the square free 
     ++ decomposition of a multivariate polynomial p.
   squareFree      :      SUP     -> Factored SUP
     ++ squareFree(p) computes the square free 
     ++ decomposition of a multivariate polynomial p presented as
     ++ a univariate polynomial with multivariate coefficients.
   squareFreePrim :      P     -> Factored P
     ++ squareFreePrim(p) compute the square free decomposition 
     ++ of a primitive multivariate polynomial p.


 
                    ----  local functions  ----
   compdegd   :             List FFE                   -> Z
	++ compdegd should be local
   univcase   :              (P,OV)                    -> Factored(P)
	++ univcase should be local
   consnewpol :            (SUP,BP,Z)                  -> Twopol
	++ consnewpol should be local
   nsqfree    :           (SUP,List(OV), List List R)  -> squareForm
	++ nsqfree should be local
   intChoose  :           (SUP,List(OV),List List R)   -> Choice
	++ intChoose should be local
   coefChoose :           (Z,Factored P)               -> P
	++ coefChoose should be local
   check      :     (List(FFE),List(FFE))              -> Boolean
	++ check should be local
   lift       : (SUP,BP,BP,P,List(OV),List(NNI),List(R)) -> Union(List(SUP),"failed")
	++ lift should be local
   myDegree   :       (SUP,List OV,NNI)                -> List NNI
	++ myDegree should be local
   normDeriv2 :            (BP,Z)                      ->  BP
	++ normDeriv2 should be local


 
 T == add
 
   pmod:R   :=  (prevPrime(2**26)$IntegerPrimesPackage(Integer))::R
 
 
   import GenExEuclid()
   import MultivariateLifting(E,OV,R,P)
   import PolynomialGcdPackage(E,OV,R,P)
   import FactoringUtilities(E,OV,R,P)
   import IntegerCombinatoricFunctions(Z)
 
 
    ----  Are the univariate square-free decompositions consistent?  ----
 
     ----  new square-free algorithm for primitive polynomial  ----
   nsqfree(oldf:SUP,lvar:List(OV),ltry:List List R) : squareForm ==
     f:=oldf
     univPol := intChoose(f,lvar,ltry)
--     debug msg
--     if not empty? ltry then output("ltry =", (ltry::OutputForm))$OutputPackage
     f0:=univPol.upol
     --the polynomial is square-free
     f0=1$BP => [1$P,[[f,1]$FFES]]$squareForm
     lfact:List(FFE):=univPol.Lfact
     lval:=univPol.Lval
     ctf:=univPol.ctpol
     leadpol:Boolean:=false
     sqdec:List FFES := empty()
     exp0:Z:=0
     unitsq:P:=1
     lcf:P:=leadingCoefficient f
     if not one? ctf then
       f0:=ctf*f0
       f:=(ctf::P)*f
       lcf:=ctf*lcf
     sqlead:List FFEP:= empty()
     sqlc:Factored P:=1
     if not one? lcf then
       leadpol:=true
       sqlc:=squareFree lcf
       unitsq:=unitsq*(unit sqlc)
       sqlead:= factors sqlc
     lfact:=sort(#1.exponent > #2.exponent,lfact)
     while lfact~=[] repeat
       pfact:=lfact.first
       (g0,exp0):=(pfact.factor,pfact.exponent)
       lfact:=lfact.rest
       lfact=[] and exp0 =1 =>
         f := (f exquo (ctf::P))::SUP
         gg := unitNormal leadingCoefficient f
         sqdec:=cons([gg.associate*f,exp0],sqdec)
         return  [gg.unit, sqdec]$squareForm
       if not one? ctf then g0:=ctf*g0
       npol:=consnewpol(f,f0,exp0)
       (d,d0):=(npol.pol,npol.polval)
       if leadpol then lcoef:=coefChoose(exp0,sqlc)
       else lcoef:=1$P
       ldeg:=myDegree(f,lvar,exp0::NNI)
       result:=lift(d,g0,(d0 exquo g0)::BP,lcoef,lvar,ldeg,lval)
       result case "failed" => return nsqfree(oldf,lvar,ltry)
       result0:SUP:= (result::List SUP).1
       r1:SUP:=result0**(exp0:NNI)
       if (h:=f exquo r1) case "failed" then return nsqfree(oldf,lvar,empty())
       sqdec:=cons([result0,exp0],sqdec)
       f:=h::SUP
       f0:=completeEval(h,lvar,lval)
       lcr:P:=leadingCoefficient result0
       if leadpol and not one? lcr then
         for lpfact in sqlead  while not one? lcr repeat
           ground? lcr =>
             unitsq:=(unitsq exquo lcr)::P
             lcr:=1$P
           (h1:=lcr exquo lpfact.factor) case "failed" => "next"
           lcr:=h1::P
           lpfact.exponent:=(lpfact.exponent)-exp0
     [((retract f) exquo ctf)::P,sqdec]$squareForm

 
   squareFree(f:SUP) : Factored SUP ==
     degree f =0 =>
       fu:=squareFree retract f
       makeFR((unit fu)::SUP,[["sqfr",ff.fctr::SUP,ff.xpnt]
               for ff in factorList fu])
     lvar:= "setUnion"/[variables cf for cf in coefficients f]
     empty? lvar =>  -- the polynomial is univariate
       upol:=map(ground,f)$UPCF2(P,SUP,R,BP)
       usqfr:=squareFree upol
       makeFR(map(coerce,unit usqfr)$UPCF2(R,BP,P,SUP),
              [["sqfr",map(coerce,ff.fctr)$UPCF2(R,BP,P,SUP),ff.xpnt]
                 for ff in factorList usqfr])

     lcf:=content f
     f:=(f exquo lcf) ::SUP
     lcSq:=squareFree lcf
     lfs:List ffSUP:=[["sqfr",ff.fctr ::SUP,ff.xpnt]
                        for ff in factorList lcSq]
     partSq:=nsqfree(f,lvar,empty())

     lfs:=append([["sqfr",fu.factor,fu.exponent]$ffSUP
                    for fu in partSq.suPart],lfs)
     makeFR((unit lcSq * partSq.unitPart) ::SUP,lfs)

   squareFree(f:P) : Factored P ==
     ground? f => makeFR(f,[])      ---   the polynomial is constant  ---
 
     lvar:List(OV):=variables(f)
     result1:List ffP:= empty()

     lmdeg :=minimumDegree(f,lvar)     ---       is the mindeg > 0 ? ---
     p:P:=1$P
     for im in 1..#lvar repeat
       (n:=lmdeg.im)=0 => "next im"
       y:=lvar.im
       p:=p*monomial(1$P,y,n)
       result1:=cons(["sqfr",y::P,n],result1)
     if not one? p then
       f := (f exquo p)::P
       if ground? f then return makeFR(f, result1)
       lvar:=variables(f)
 
 
     #lvar=1 =>                    ---  the polynomial is univariate ---
       result:=univcase(f,lvar.first)
       makeFR(unit result,append(result1,factorList result))
 
     ldeg:=degree(f,lvar)          ---  general case ---
     m:="min"/[j for j in ldeg|j~=0]
     i:Z:=1
     for j in ldeg while j>m repeat i:=i+1
     x:=lvar.i
     lvar:=delete(lvar,i)
     f0:=univariate (f,x)
     lcont:P:= content f0
     nsqfftot:=nsqfree((f0 exquo lcont)::SUP,lvar,empty())
     nsqff:List ffP:=[["sqfr",multivariate(fu.factor,x),fu.exponent]$ffP
                          for fu in nsqfftot.suPart]
     result1:=append(result1,nsqff)
     ground? lcont => makeFR(lcont*nsqfftot.unitPart,result1)
     sqlead:=squareFree(lcont)
     makeFR(unit sqlead*nsqfftot.unitPart,append(result1,factorList sqlead))
 
  -- Choose the integer for the evaluation.                        --
  -- If the polynomial is square-free the function returns upol=1. --
 
   intChoose(f:SUP,lvar:List(OV),ltry:List List R):Choice ==
     degf:= degree f
     tryCount:NNI:=0
     nvr:=#lvar
     range:Z:=10
     lfact1:List(FFE):=[]
     lval1:List R := []
     lfact:List(FFE)
     ctf1:R:=1
     f1:BP:=1$BP
     d1:Z
     while range < 10000000000 repeat
       range:=2*range
       lval:= [ran(range) for i in 1..nvr]
       member?(lval,ltry) => "new integer"
       ltry:=cons(lval,ltry)
       f0:=completeEval(f,lvar,lval)
       degree f0 ~=degf  => "new integer"
       ctf:=content f0
       lfact:List(FFE):=factors(squareFree((f0 exquo (ctf:R)::BP)::BP))
 
            ----  the univariate polynomial is square-free  ----
       if #lfact=1 and (lfact.1).exponent=1 then
         return [1$BP,lval,lfact,1$R]$Choice
 
       d0:=compdegd lfact
                 ----      inizialize lfact1      ----
       tryCount=0 =>
         f1:=f0
         lfact1:=lfact
         ctf1:=ctf
         lval1:=lval
         d1:=d0
         tryCount:=1
       d0=d1 =>
         return [f1,lval1,lfact1,ctf1]$Choice
       d0 < d1 =>
         tryCount:=1
         f1:=f0
         lfact1:=lfact
         ctf1:=ctf
         lval1:=lval
         d1:=d0
     error "intChoose$MULTQFR: fell off loop without value"
         
 
        ----  Choose the leading coefficient for the lifting  ----
   coefChoose(exp:Z,sqlead:Factored(P)) : P ==
     lcoef:P:=unit(sqlead)
     for term in factors(sqlead) repeat
       texp:=term.exponent
       texp<exp => "next term"
       texp=exp => lcoef:=lcoef*term.factor
       lcoef:=lcoef*(term.factor)**((texp quo exp)::NNI)
     lcoef
 
        ----  Construction of the polynomials for the lifting  ----
   consnewpol(g:SUP,g0:BP,deg:Z):Twopol ==
     deg=1 => [g,g0]$Twopol
     deg:=deg-1
     [normalDeriv(g,deg),normDeriv2(g0,deg)]$Twopol
 
         ----  lift the univariate square-free factor  ----
   lift(ud:SUP,g0:BP,g1:BP,lcoef:P,lvar:List(OV),
                        ldeg:List(NNI),lval:List(R)) : Union(List SUP,"failed") ==
     leadpol:Boolean:=false
     lcd:P:=leadingCoefficient ud
     leadlist:List(P):=empty()
 
     if not ground?(leadingCoefficient ud) then
       leadpol:=true
       ud:=lcoef*ud
       lcg0:R:=leadingCoefficient g0
       if ground? lcoef then g0:=retract(lcoef) quo lcg0 *g0
       else g0:=(retract(eval(lcoef,lvar,lval)) quo lcg0) * g0
       g1:=lcg0*g1
       leadlist:=[lcoef,lcd]
     plist:=lifting(ud,lvar,[g0,g1],lval,leadlist,ldeg,pmod)
     plist case "failed" => "failed" 
     (p0:SUP,p1:SUP):=((plist::List SUP).1,(plist::List SUP).2)
     if completeEval(p0,lvar,lval) ~= g0 then (p0,p1):=(p1,p0)
     [primitivePart p0,primitivePart p1]

                ----  the polynomial is univariate  ----
   univcase(f:P,x:OV) : Factored(P) ==
     uf := univariate f
     cf:=content uf
     uf :=(uf exquo cf)::BP
     result:Factored BP:=squareFree uf
     makeFR(multivariate(cf*unit result,x),
         [["sqfr",multivariate(term.factor,x),term.exponent]
           for term in factors result])

--   squareFreePrim(p:P) : Factored P ==
--     -- p is content free
--     ground? p => makeFR(p,[])      ---   the polynomial is constant  ---
-- 
--     lvar:List(OV):=variables p
--     #lvar=1 =>                    ---  the polynomial is univariate ---
--       univcase(p,lvar.first)
--     nsqfree(p,lvar,1)
 
   compdegd(lfact:List(FFE)) : Z ==
     ris:Z:=0
     for pfact in lfact repeat
       ris:=ris+(pfact.exponent -1)*degree pfact.factor
     ris
 
   normDeriv2(f:BP,m:Z) : BP ==
     (n1:Z:=degree f) < m => 0$BP
     n1=m => (leadingCoefficient f)::BP
     k:=binomial(n1,m)
     ris:BP:=0$BP
     n:Z:=n1
     while n>= m repeat
       while n1>n repeat
         k:=(k*(n1-m)) quo n1
         n1:=n1-1
       ris:=ris+monomial(k*leadingCoefficient f,(n-m)::NNI)
       f:=reductum f
       n:=degree f
     ris

   myDegree(f:SUP,lvar:List OV,exp:NNI) : List NNI==
     [n quo exp for n in degree(f,lvar)]

@
\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 MULTSQFR MultivariateSquareFree>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}