\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra groebsol.spad}
\author{Patrizia Gianni}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package GROEBSOL GroebnerSolve}
<<package GROEBSOL GroebnerSolve>>=
)abbrev package GROEBSOL GroebnerSolve
++ Author : P.Gianni, Summer '88, revised November '89
++ Solve systems of polynomial equations using Groebner bases
++ Total order Groebner bases are computed and then converted to lex ones
++ This package is mostly intended for internal use.
GroebnerSolve(lv,F,R) : C == T

  where
   R      :   GcdDomain
   F      :   GcdDomain
   lv     :   List Symbol

   NNI    ==>  NonNegativeInteger
   I      ==>  Integer
   S      ==>  Symbol

   OV     ==>  OrderedVariableList(lv)
   IES    ==>  IndexedExponents Symbol

   DP     ==>  DirectProduct(#lv,NonNegativeInteger)
   DPoly  ==>  DistributedMultivariatePolynomial(lv,F)

   HDP    ==>  HomogeneousDirectProduct(#lv,NonNegativeInteger)
   HDPoly ==>  HomogeneousDistributedMultivariatePolynomial(lv,F)

   SUP    ==>  SparseUnivariatePolynomial(DPoly)
   L      ==>  List
   P      ==>  Polynomial

   C == with
      groebSolve   : (L DPoly,L OV)  -> L L DPoly
        ++ groebSolve(lp,lv) reduces the polynomial system lp in variables lv
        ++ to triangular form. Algorithm based on groebner bases algorithm
        ++ with linear algebra for change of ordering.
        ++ Preprocessing for the general solver.
        ++ The polynomials in input are of type \spadtype{DMP}.

      testDim     : (L HDPoly,L OV)  -> Union(L HDPoly,"failed")
        ++ testDim(lp,lv) tests if the polynomial system lp
        ++ in variables lv is zero dimensional.

      genericPosition : (L DPoly, L OV) -> Record(dpolys:L DPoly, coords: L I)
        ++ genericPosition(lp,lv) puts a radical zero dimensional ideal
        ++ in general position, for system lp in variables lv.

   T == add
     import Boolean
     import PolToPol(lv,F)
     import GroebnerPackage(F,DP,OV,DPoly)
     import GroebnerInternalPackage(F,DP,OV,DPoly)
     import GroebnerPackage(F,HDP,OV,HDPoly)
     import LinGroebnerPackage(lv,F)

     nv:NNI:=#lv

          ---- test if f is power of a linear mod (rad lpol) ----
                     ----  f is monic  ----
     testPower(uf:SUP,x:OV,lpol:L DPoly) : Union(DPoly,"failed") ==
       df:=degree(uf)
       trailp:DPoly := coefficient(uf,(df-1)::NNI)
       (testquo := trailp exquo (df::F)) case "failed" => "failed"
       trailp := testquo::DPoly
       gg:=gcd(lc:=leadingCoefficient(uf),trailp)
       trailp := (trailp exquo gg)::DPoly
       lc := (lc exquo gg)::DPoly
       linp:SUP:=monomial(lc,1$NNI)$SUP + monomial(trailp,0$NNI)$SUP
       g:DPoly:=multivariate(uf-linp**df,x)
       redPol(g,lpol) ~= 0 => "failed"
       multivariate(linp,x)

            -- is the 0-dimensional ideal I in general position ?  --
                     ----  internal function  ----
     testGenPos(lpol:L DPoly,lvar:L OV):Union(L DPoly,"failed") ==
       rlpol:=reverse lpol
       f:=rlpol.first
       #lvar=1 => [f]
       rlvar:=rest reverse lvar
       newlpol:List(DPoly):=[f]
       for f: local in rlpol.rest repeat
         x:=first rlvar
         fi:= univariate(f,x)
         if (mainVariable leadingCoefficient fi case "failed") then
           if ((g:= testPower(fi,x,newlpol)) case "failed")
           then return "failed"
           newlpol :=concat(redPol(g::DPoly,newlpol),newlpol)
           rlvar:=rest rlvar
         else if not zero? redPol(f,newlpol) then return"failed"
       newlpol


        -- change coordinates and out the ideal in general position  ----
     genPos(lp:L DPoly,lvar:L OV): Record(polys:L HDPoly, lpolys:L DPoly,
                                           coord:L I, univp:HDPoly) ==
           rlvar:=reverse lvar
           lnp:=[dmpToHdmp(f) for f in lp]
           x := first rlvar;rlvar:=rest rlvar
           testfail:=true
           ranvals: L I
           gb: L HDPoly
           gbt: L DPoly
           gb1: Union(L DPoly,"failed")
           for count in 1.. while testfail repeat
             ranvals := [1+(random()$I rem (count*(# lvar))) for vv in rlvar]
             val:=+/[rv*(vv::HDPoly)
                        for vv in rlvar for rv in ranvals]
             val:=val+x::HDPoly
             gb := [elt(univariate(p,x),val) for p in lnp]
             gb:=groebner gb
             gbt:=totolex gb
             (gb1:=testGenPos(gbt,lvar)) case "failed"=>"try again"
             testfail:=false
           [gb,gbt,ranvals,dmpToHdmp(last (gb1::L DPoly))]

     genericPosition(lp:L DPoly,lvar:L OV) ==
        nans:=genPos(lp,lvar)
        [nans.lpolys, nans.coord]

        ---- select  the univariate factors
     select(lup:L L HDPoly) : L L HDPoly ==
       lup=[] => list []
       [:[cons(f,lsel) for lsel in select lup.rest] for f in lup.first]

        ---- in the non generic case, we compute the prime ideals ----
           ---- associated to leq, basis is the algebra basis  ----
     findCompon(leq:L HDPoly,lvar:L OV):L L DPoly ==
       teq:=totolex(leq)
       #teq = #lvar => [teq]
      -- not ((teq1:=testGenPos(teq,lvar)) case "failed") => [teq1::L DPoly]
       gp:=genPos(teq,lvar)
       lgp:= gp.polys
       g:HDPoly:=gp.univp
       fg:=(factor g)$GeneralizedMultivariateFactorize(OV,HDP,R,F,HDPoly)
       lfact:=[ff.factor for ff in factors(fg::Factored(HDPoly))]
       result: L L HDPoly := []
       #lfact=1 => [teq]
       for tfact in lfact repeat
         tlfact:=concat(tfact,lgp)
         result:=concat(tlfact,result)
       ranvals:L I:=gp.coord
       rlvar:=reverse lvar
       x:=first rlvar
       rlvar:=rest rlvar
       val:=+/[rv*(vv::HDPoly) for vv in rlvar for rv in ranvals]
       val:=(x::HDPoly)-val
       ans:=[totolex groebner [elt(univariate(p,x),val) for p in lp]
                           for lp in result]
       [ll for ll in ans | ll~=[1]]

     zeroDim?(lp: List HDPoly,lvar:L OV) : Boolean ==
       empty? lp => false
       n:NNI := #lvar
       #lp < n => false
       lvint1 := lvar
       for f in lp while not empty?(lvint1) repeat
          g:= f - reductum f
          x:=mainVariable(g)::OV
          if ground?(leadingCoefficient(univariate(g,x))) then
               lvint1 := remove(x, lvint1)
       empty? lvint1

     -- general solve, gives an error if the system not 0-dimensional
     groebSolve(leq: L DPoly,lvar:L OV) : L L DPoly ==
       lnp:=[dmpToHdmp(f) for f in leq]
       leq1:=groebner lnp
       #(leq1) = 1 and first(leq1) = 1 => list empty()
       not (zeroDim?(leq1,lvar)) =>
         error "system does not have a finite number of solutions"
       -- add computation of dimension, for a more useful error
       basis:=computeBasis(leq1)
       lup:L HDPoly:=[]
       llfact:L Factored(HDPoly):=[]
       for x in lvar repeat
         g:=minPol(leq1,basis,x)
         fg:=(factor g)$GeneralizedMultivariateFactorize(OV,HDP,R,F,HDPoly)
         llfact:=concat(fg::Factored(HDPoly),llfact)
         if degree(g,x) = #basis then leave "stop factoring"
       result: L L DPoly := []
       -- selecting a factor from the lists of the univariate factors
       lfact:=select [[ff.factor for ff in factors llf]
                       for llf in llfact]
       for tfact in lfact repeat
         tfact:=groebner concat(tfact,leq1)
         tfact=[1] => "next value"
         result:=concat(result,findCompon(tfact,lvar))
       result

     -- test if the system is zero dimensional
     testDim(leq : L HDPoly,lvar : L OV) : Union(L HDPoly,"failed") ==
       leq1:=groebner leq
       #(leq1) = 1 and first(leq1) = 1 => empty()
       not (zeroDim?(leq1,lvar)) => "failed"
       leq1

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