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

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra numsolve.spad}
\author{Patrizia Gianni}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package INFSP InnerNumericFloatSolvePackage}
<<package INFSP InnerNumericFloatSolvePackage>>=
)abbrev package INFSP InnerNumericFloatSolvePackage
++ Author: P. Gianni
++ Date Created: January 1990
++ Date Last Updated:
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This is an internal package
++ for computing approximate solutions to systems of polynomial equations.
++ The parameter K specifies the coefficient field of the input polynomials
++ and must be either \spad{Fraction(Integer)} or \spad{Complex(Fraction Integer)}.
++ The parameter F specifies where the solutions must lie and can
++ be one of the following: \spad{Float}, \spad{Fraction(Integer)}, \spad{Complex(Float)},
++ \spad{Complex(Fraction Integer)}. The last parameter specifies the type
++ of the precision operand and must be either \spad{Fraction(Integer)} or \spad{Float}.
InnerNumericFloatSolvePackage(K,F,Par): Cat == Cap where
    F    :   Field  -- this is the field where the answer will be
    K    :   GcdDomain  -- type of the  input
    Par  :   Join(Field, OrderedRing ) -- it will be NF or RN

    I        ==> Integer
    NNI      ==> NonNegativeInteger
    P        ==> Polynomial
    EQ       ==> Equation
    L        ==> List
    SUP      ==> SparseUnivariatePolynomial
    RN       ==> Fraction Integer
    NF       ==> Float
    CF       ==> Complex Float
    GI       ==> Complex Integer
    GRN      ==> Complex RN
    SE       ==> Symbol
    RFI      ==> Fraction P I

    Cat == with

       innerSolve1   :  (SUP K,Par)  -> L F
          ++ innerSolve1(up,eps) returns the list of the zeros
          ++ of the univariate polynomial up with precision eps.
       innerSolve1   :  (P K,Par)  -> L F
          ++ innerSolve1(p,eps) returns the list of the zeros
          ++ of the polynomial p with precision eps.
       innerSolve    : (L P K,L P K,L SE,Par) -> L L F
          ++ innerSolve(lnum,lden,lvar,eps) returns a list of
          ++ solutions of the system of polynomials lnum, with
          ++ the side condition that none of the members of lden
          ++ vanish identically on any solution. Each solution
          ++ is expressed as a list corresponding to the list of
          ++ variables in lvar and with precision specified by eps.
       makeEq        : (L F,L SE)     -> L EQ P F
          ++ makeEq(lsol,lvar) returns a list of equations formed
          ++ by corresponding members of lvar and lsol.

    Cap == add

                  ------  Local Functions  ------
       isGeneric? : (L P K,L SE) -> Boolean
       evaluate  : (P K,SE,SE,F) ->  F
       numeric   :     K          -> F
       oldCoord      : (L F,L I) -> L F
       findGenZeros  : (L P K,L SE,Par) -> L L F
       failPolSolve  : (L P K,L SE)  -> Union(L L P K,"failed")

       numeric(r:K):F ==
         K is I =>
           F is Float => r::I::Float
           F is RN    => r::I::RN
           F is CF    => r::I::CF
           F is GRN   => r::I::GRN
         K is GI =>
           gr:GI := r::GI
           F is GRN => complex(real(gr)::RN,imag(gr)::RN)$GRN
           F is CF  => convert(gr)
         error "case not handled"

       -- construct the equation
       makeEq(nres:L F,lv:L SE) : L EQ P F ==
           [equation(x::(P F),r::(P F)) for x in lv for r in nres]

       evaluate(pol:P K,xvar:SE,zvar:SE,z:F):F ==
         rpp:=map(numeric,pol)$PolynomialFunctions2(K,F)
         rpp := eval(rpp,zvar,z)
         upol:=univariate(rpp,xvar)
         retract(-coefficient(upol,0))/retract(leadingCoefficient upol)

       myConvert(eps:Par) : RN ==
         Par is RN => eps
         Par is NF => retract(eps)$NF

       innerSolve1(pol:P K,eps:Par) : L F == innerSolve1(univariate pol,eps)

       innerSolve1(upol:SUP K,eps:Par) : L F ==
         K is GI and (Par is RN or Par is NF) =>
             (complexZeros(upol,
                        eps)$ComplexRootPackage(SUP K,Par)) pretend L(F)
         K is I =>
           F is Float =>
             z:= realZeros(upol,myConvert eps)$RealZeroPackage(SUP I)
             [convert((1/2)*(x.left+x.right))@Float for x in z] pretend L(F)

           F is RN =>
             z:= realZeros(upol,myConvert eps)$RealZeroPackage(SUP I)
             [(1/2)*(x.left + x.right) for x in z] pretend L(F)
           error "improper arguments to INFSP"
         error "improper arguments to INFSP"


       -- find the zeros of components in "generic" position --
       findGenZeros(lp:L P K,rlvar:L SE,eps:Par) : L L F ==
         rlp:=reverse lp
         f:=rlp.first
         zvar:= rlvar.first
         rlp:=rlp.rest
         lz:=innerSolve1(f,eps)
         [reverse cons(z,[evaluate(pol,xvar,zvar,z) for pol in rlp
                         for xvar in rlvar.rest]) for z in lz]

       -- convert to the old coordinates --
       oldCoord(numres:L F,lval:L I) : L F ==
         rnumres:=reverse numres
         rnumres.first:= rnumres.first +
            (+/[n*nr for n in lval for nr in rnumres.rest])
         reverse rnumres

       -- real zeros of a system of 2 polynomials lp (incomplete)
       innerSolve2(lp:L P K,lv:L SE,eps: Par):L L F ==
          mainvar := first lv
          up1:=univariate(lp.1, mainvar)
          up2:=univariate(lp.2, mainvar)
          vec := subresultantVector(up1,up2)$SubResultantPackage(P K,SUP P K)
          p0 := primitivePart multivariate(vec.0, mainvar)
          p1 := primitivePart(multivariate(vec.1, mainvar),mainvar)
          zero? p1 or
            not one? gcd(p0, leadingCoefficient(univariate(p1,mainvar))) =>
              innerSolve(cons(0,lp),empty(),lv,eps)
          findGenZeros([p1, p0], reverse lv, eps)

       -- real zeros of the system of polynomial lp --
       innerSolve(lp:L P K,ld:L P K,lv:L SE,eps: Par) : L L F ==
          -- empty?(ld) and (#lv = 2) and (# lp = 2) => innerSolve2(lp, lv, eps)
           lnp:= [pToDmp(p)$PolToPol(lv,K) for p in lp]
           OV:=OrderedVariableList(lv)
           lvv:L OV:= [variable(vv)::OV for vv in lv]
           DP:=DirectProduct(#lv,NonNegativeInteger)
           dmp:=DistributedMultivariatePolynomial(lv,K)
           lq:L dmp:=[]
           if ld~=[] then
             lq:= [(pToDmp(q1)$PolToPol(lv,K)) pretend dmp for q1 in ld]
           partRes:=groebSolve(lnp,lvv)$GroebnerSolve(lv,K,K) pretend (L L dmp)
           partRes=list [] => []
           -- remove components where denominators vanish
           if lq~=[] then
             gb:=GroebnerInternalPackage(K,DirectProduct(#lv,NNI),OV,dmp)
             partRes:=[pr for pr in partRes|
                       and/[not zero?(redPol(fq,pr pretend List(dmp))$gb)
                         for fq in lq]]

           -- select the components in "generic" form
           rlv:=reverse lv
           rrlvv:= rest reverse lvv

           listGen:L L dmp:=[]
           for res in partRes repeat
             res1:=rest reverse res
             "and"/[("max"/degree(f,rrlvv))=1  for f in res1] =>
                      listGen:=concat(res pretend (L dmp),listGen)
           result:L L F := []
           if listGen~=[] then
             listG :L L P K:=
               [[dmpToP(pf)$PolToPol(lv,K) for pf in pr] for pr in listGen]
             result:=
               "append"/[findGenZeros(res,rlv,eps) for res in listG]
             for gres in listGen repeat
                partRes:=delete(partRes,position(gres,partRes))
           -- adjust the non-generic components
           for gres in partRes repeat
               genRecord := genericPosition(gres,lvv)$GroebnerSolve(lv,K,K)
               lgen := genRecord.dpolys
               lval := genRecord.coords
               lgen1:=[dmpToP(pf)$PolToPol(lv,K) for pf in lgen]
               lris:=findGenZeros(lgen1,rlv,eps)
               result:= append([oldCoord(r,lval) for r in lris],result)
           result

@
\section{package FLOATRP FloatingRealPackage}
<<package FLOATRP FloatingRealPackage>>=
)abbrev package FLOATRP FloatingRealPackage
++ Author: P. Gianni
++ Date Created: January 1990
++ Date Last Updated:
++ Basic Functions:
++ Related Constructors: SystemSolvePackage, RadicalSolvePackage,
++ FloatingComplexPackage
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++    This is a package for the approximation of real solutions for
++ systems of polynomial equations over the rational numbers.
++ The results are expressed as either rational numbers or floats
++ depending on the type of the precision parameter which can be
++ either a rational number or a floating point number.
FloatingRealPackage(Par): Cat == Cap where
    I        ==> Integer
    NNI      ==> NonNegativeInteger
    P        ==> Polynomial
    EQ       ==> Equation
    L        ==> List
    SUP      ==> SparseUnivariatePolynomial
    RN       ==> Fraction Integer
    NF       ==> Float
    CF       ==> Complex Float
    GI       ==> Complex Integer
    GRN      ==> Complex RN
    SE       ==> Symbol
    RFI      ==> Fraction P I
    INFSP ==> InnerNumericFloatSolvePackage

    Par : Join(OrderedRing, Field)  -- RN or NewFloat

    Cat == with

       solve:    (L RFI,Par) -> L L EQ P Par
         ++ solve(lp,eps) finds all of the real solutions of the
         ++ system lp of rational functions over the rational numbers
         ++ with respect to all the variables appearing in lp,
         ++ with precision eps.

       solve:    (L EQ RFI,Par) -> L L EQ P Par
         ++ solve(leq,eps) finds all of the real solutions of the
         ++ system leq of equationas of rational functions
         ++ with respect to all the variables appearing in lp,
         ++ with precision eps.

       solve:    (RFI,Par) ->  L EQ P Par
         ++ solve(p,eps) finds all of the real solutions of the
         ++ univariate rational function p with rational coefficients
         ++ with respect to the unique variable appearing in p,
         ++ with precision eps.

       solve:    (EQ RFI,Par) ->  L EQ P Par
         ++ solve(eq,eps) finds all of the real solutions of the
         ++ univariate equation eq of rational functions
         ++ with respect to the unique variables appearing in eq,
         ++ with precision eps.

       realRoots:    (L RFI,L SE,Par) -> L L Par
         ++ realRoots(lp,lv,eps) computes the list of the real
         ++ solutions of the list lp of rational functions with rational
         ++ coefficients with respect to the variables in lv,
         ++ with precision eps. Each solution is expressed as a list
         ++ of numbers in order corresponding to the variables in lv.

       realRoots : (RFI,Par) -> L Par
         ++ realRoots(rf, eps) finds the real zeros of a univariate
         ++ rational function with precision given by eps.

    Cap == add

       makeEq(nres:L Par,lv:L SE) : L EQ P Par ==
           [equation(x::(P Par),r::(P Par)) for x in lv for r in nres]

       -- find the real zeros of an univariate rational polynomial --
       realRoots(p:RFI,eps:Par) : L Par ==
         innerSolve1(numer p,eps)$INFSP(I,Par,Par)

       -- real zeros of the system of polynomial lp --
       realRoots(lp:L RFI,lv:L SE,eps: Par) : L L Par ==
         lnum:=[numer p for p in lp]
         lden:=[dp for p in lp |not one?(dp:=denom p)]
         innerSolve(lnum,lden,lv,eps)$INFSP(I,Par,Par)

       solve(lp:L RFI,eps : Par) : L L EQ  P Par ==
         lnum:=[numer p for p in lp]
         lden:=[dp for p in lp |not one?(dp:=denom p)]
         lv:="setUnion"/[variables np for np in lnum]
         if lden~=[] then
          lv:=setUnion(lv,"setUnion"/[variables dp for dp in lden])
         [makeEq(numres,lv) for numres
            in innerSolve(lnum,lden,lv,eps)$INFSP(I,Par,Par)]

       solve(le:L EQ RFI,eps : Par) : L L EQ  P Par ==
         lp:=[lhs ep - rhs ep for ep in le]
         lnum:=[numer p for p in lp]
         lden:=[dp for p in lp |not one?(dp:=denom p)]
         lv:="setUnion"/[variables np for np in lnum]
         if lden~=[] then
          lv:=setUnion(lv,"setUnion"/[variables dp for dp in lden])
         [makeEq(numres,lv) for numres
           in innerSolve(lnum,lden,lv,eps)$INFSP(I,Par,Par)]

       solve(p : RFI,eps : Par) :  L EQ  P Par ==
         (mvar := mainVariable numer p ) case "failed" =>
              error "no variable found"
         x:P Par:=mvar::SE::(P Par)
         [equation(x,val::(P Par)) for val in realRoots(p,eps)]

       solve(eq : EQ RFI,eps : Par) :  L EQ  P Par ==
         solve(lhs eq - rhs eq,eps)

@
\section{package FLOATCP FloatingComplexPackage}
<<package FLOATCP FloatingComplexPackage>>=
)abbrev package FLOATCP FloatingComplexPackage
++ Author: P. Gianni
++ Date Created: January 1990
++ Date Last Updated:
++ Basic Functions:
++ Related Constructors: SystemSolvePackage, RadicalSolvePackage,
++ FloatingRealPackage
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++    This is a package for the approximation of complex solutions for
++ systems of equations of rational functions with complex rational
++ coefficients. The results are expressed as either complex rational
++ numbers or complex floats depending on the type of the precision
++ parameter which can be either a rational number or a floating point number.
FloatingComplexPackage(Par): Cat == Cap where
    Par : Join(Field, OrderedRing)
    K   ==> GI
    FPK ==> Fraction P K
    C   ==> Complex
    I        ==> Integer
    NNI      ==> NonNegativeInteger
    P        ==> Polynomial
    EQ       ==> Equation
    L        ==> List
    SUP      ==> SparseUnivariatePolynomial
    RN       ==> Fraction Integer
    NF       ==> Float
    CF       ==> Complex Float
    GI       ==> Complex Integer
    GRN      ==> Complex RN
    SE       ==> Symbol
    RFI      ==> Fraction P I
    INFSP ==> InnerNumericFloatSolvePackage


    Cat == with

       complexSolve:    (L FPK,Par) -> L L EQ P C Par
         ++ complexSolve(lp,eps) finds all the complex solutions to
         ++ precision eps of the system lp of rational functions
         ++ over the complex rationals with respect to all the
         ++ variables appearing in lp.

       complexSolve:    (L EQ FPK,Par) -> L L EQ P C Par
         ++ complexSolve(leq,eps) finds all the complex solutions
         ++ to precision eps of the system leq of equations
         ++ of rational functions over complex rationals
         ++ with respect to all the variables appearing in lp.

       complexSolve:    (FPK,Par) -> L EQ P C Par
         ++ complexSolve(p,eps) find all the complex solutions of the
         ++ rational function p with complex rational coefficients
         ++ with respect to all the variables appearing in p,
         ++ with precision eps.

       complexSolve:    (EQ FPK,Par) ->  L EQ P C Par
         ++ complexSolve(eq,eps) finds all the complex solutions of the
         ++ equation eq of rational functions with rational rational coefficients
         ++ with respect to all the variables appearing in eq,
         ++ with precision eps.

       complexRoots : (FPK,Par) -> L C Par
         ++ complexRoots(rf, eps) finds all the complex solutions of a
         ++ univariate rational function with rational number coefficients.
         ++ The solutions are computed to precision eps.

       complexRoots : (L FPK,L SE,Par) -> L L C Par
         ++ complexRoots(lrf, lv, eps) finds all the complex solutions of a
         ++ list of rational functions with rational number coefficients
         ++ with respect the the variables appearing in lv.
         ++ Each solution is computed to precision eps and returned as
         ++ list corresponding to the order of variables in lv.

    Cap == add

       -- find the complex zeros of an univariate polynomial --
       complexRoots(q:FPK,eps:Par) : L C Par ==
         p:=numer q
         complexZeros(univariate p,eps)$ComplexRootPackage(SUP GI, Par)

       -- find the complex zeros of an univariate polynomial --
       complexRoots(lp:L FPK,lv:L SE,eps:Par) : L L C Par ==
         lnum:=[numer p for p in lp]
         lden:=[dp for p in lp |not one?(dp:=denom p)]
         innerSolve(lnum,lden,lv,eps)$INFSP(K,C Par,Par)

       complexSolve(lp:L FPK,eps : Par) : L L EQ  P C Par ==
         lnum:=[numer p for p in lp]
         lden:=[dp for p in lp |not one?(dp:=denom p)]
         lv:="setUnion"/[variables np for np in lnum]
         if lden~=[] then
          lv:=setUnion(lv,"setUnion"/[variables dp for dp in lden])
         [[equation(x::(P C Par),r::(P C Par)) for x in lv for r in nres]
           for nres in innerSolve(lnum,lden,lv,eps)$INFSP(K,C Par,Par)]

       complexSolve(le:L EQ FPK,eps : Par) : L L EQ  P C Par ==
         lp:=[lhs ep - rhs ep for ep in le]
         lnum:=[numer p for p in lp]
         lden:=[dp for p in lp |not one?(dp:=denom p)]
         lv:="setUnion"/[variables np for np in lnum]
         if lden~=[] then
          lv:=setUnion(lv,"setUnion"/[variables dp for dp in lden])
         [[equation(x::(P C Par),r::(P C Par)) for x in lv for r in nres]
           for nres in innerSolve(lnum,lden,lv,eps)$INFSP(K,C Par,Par)]

       complexSolve(p : FPK,eps : Par) : L EQ  P C Par ==
         (mvar := mainVariable numer p ) case "failed" =>
                 error "no variable found"
         x:P C Par:=mvar::SE::(P C Par)
         [equation(x,val::(P C Par)) for val in complexRoots(p,eps)]

       complexSolve(eq : EQ FPK,eps : Par) : L EQ  P C Par ==
         complexSolve(lhs eq - rhs eq,eps)

@
\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 INFSP InnerNumericFloatSolvePackage>>
<<package FLOATRP FloatingRealPackage>>
<<package FLOATCP FloatingComplexPackage>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}