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

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra syssolp.spad}
\author{Patrizia Gianni}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package SYSSOLP SystemSolvePackage}
<<package SYSSOLP SystemSolvePackage>>=
)abbrev package SYSSOLP SystemSolvePackage
++ Author: P. Gianni
++ Date Created: summer 1988
++ Date Last Updated: summer 1990
++ Basic Functions:
++ Related Constructors: Fraction, Polynomial, FloatingRealPackage,
++ FloatingComplexPackage, RadicalSolvePackage
++ Also See: LinearSystemMatrixPackage, GroebnerSolve
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ Symbolic solver for systems of rational functions with coefficients
++ in an integral domain R.
++ The systems are solved in the field of rational functions over R.
++ Solutions are exact of the form variable = value when the value is
++ a member of the coefficient domain R. Otherwise the solutions
++ are implicitly expressed as roots of univariate polynomial equations over R.
++ Care is taken to guarantee that the denominators of the input
++ equations do not vanish on the solution sets.
++ The arguments to solve can either be given as equations or
++ as rational functions interpreted as equal
++ to zero. The user can specify an explicit list of symbols to
++ be solved for, treating all other symbols appearing as parameters
++ or omit the list of symbols in which case the system tries to
++ solve with respect to all symbols appearing in the input.

NNI      ==> NonNegativeInteger
P        ==> Polynomial
EQ       ==> Equation
L        ==> List
V        ==> Vector
M        ==> Matrix
UP       ==> SparseUnivariatePolynomial
SE       ==> Symbol
IE       ==> IndexedExponents Symbol
SUP      ==> SparseUnivariatePolynomial

SystemSolvePackage(R): Cat == Cap where
    R : IntegralDomain
    F   ==> Fraction Polynomial R
    PP2 ==> PolynomialFunctions2(R,F)
    PPR ==> Polynomial Polynomial R

    Cat == with
       solve:    (L F,    L SE) -> L L EQ F
         ++ solve(lp,lv) finds the solutions of the list lp of
         ++ rational functions with respect to the list of symbols lv.

       solve:    (L EQ F, L SE) -> L L EQ F
         ++ solve(le,lv) finds the solutions of the
         ++ list le of equations of rational functions
         ++ with respect to the list of symbols lv.

       solve:         L F       -> L L EQ F
         ++ solve(lp) finds the solutions of the list lp of rational
         ++ functions with respect to all symbols appearing in lp.

       solve:        L EQ F     -> L L EQ F
         ++ solve(le) finds the solutions of the list le of equations of
         ++ rational functions with respect to all symbols appearing in le.

       solve:    (F, SE) ->  L EQ F
         ++ solve(p,v) solves the equation p=0, where p is a rational function
         ++ with respect to the variable v.

       solve:    (EQ F,SE) -> L EQ F
         ++ solve(eq,v) finds the solutions of the equation
         ++ eq with respect to the variable v.

       solve:         F      -> L EQ F
         ++ solve(p) finds the solution of a rational function p = 0
         ++ with respect to the unique variable appearing in p.

       solve:         EQ F     -> L EQ F
         ++ solve(eq) finds the solutions of the equation eq
         ++ with respect to the unique variable appearing in eq.

       triangularSystems: (L F,    L SE) -> L L P R
         ++ triangularSystems(lf,lv) solves the system of equations
         ++ defined by lf with respect to the list of symbols lv;
         ++ the system of equations is obtaining
         ++ by equating to zero the list of rational functions lf.
         ++ The output is a list of solutions where
         ++ each solution is expressed as a "reduced" triangular system of
         ++ polynomials.

    Cap == add

       import MPolyCatRationalFunctionFactorizer(IE,SE,R,P F)

                     ---- Local Functions ----
       linSolve: (L F,    L SE) -> Union(L EQ F, "failed")
       makePolys :   L EQ F     ->  L F

       makeR2F(r : R) : F == r :: (P R) :: F

       makeP2F(p:P F):F ==
         lv:=variables p
         lv = [] => retract p
         for v in lv repeat p:=pushdown(p,v)
         retract p
                     ---- Local Functions ----
       makeEq(p:P F,lv:L SE): EQ F ==
         z:=last lv
         np:=numer makeP2F p
         lx:=variables np
         x : SE
         for x: free in lv repeat if member?(x,lx) then leave x
         up:=univariate(np,x)
         (degree up)=1 =>
           equation(x::P(R)::F,-coefficient(up,0)/leadingCoefficient up)
         equation(np::F,0$F)

       varInF(v: SE): F == v::P(R) :: F

       newInF(n: Integer):F==varInF new()$SE

       testDegree(f :P R , lv :L SE) : Boolean ==
         "or"/[positive? degree(f,vv) for vv in lv]
                    ---- Exported Functions ----

       -- solve a system of rational functions
       triangularSystems(lf: L F,lv:L SE) : L L P R ==
           empty? lv => empty()
           empty? lf => empty()
           #lf = 1 =>
              p:= numer(first lf)
              fp:=(factor p)$GeneralizedMultivariateFactorize(SE,IE,R,R,P R)
              [[ff.factor] for ff in factors fp | testDegree(ff.factor,lv)]
           dmp:=DistributedMultivariatePolynomial(lv,P R)
           OV:=OrderedVariableList(lv)
           DP:=DirectProduct(#lv, NonNegativeInteger)
           push:=PushVariables(R,DP,OV,dmp)
           lq : L dmp
           lvv:L OV:=[variable(vv)::OV for vv in lv]
           lq:=[pushup(df::dmp,lvv)$push for f in lf|not one?(df:=denom f)]
           lp:=[pushup(numer(f)::dmp,lvv)$push for f in lf]
           parRes:=groebSolve(lp,lvv)$GroebnerSolve(lv,P R,R)
           if lq~=[] then
             gb:=GroebnerInternalPackage(P R,DirectProduct(#lv,NNI),OV,dmp)
             parRes:=[pr for pr in parRes|
                       and/[not zero?(redPol(fq,pr pretend List(dmp))$gb)
                         for fq in lq]]
           [[retract pushdown(pf,lvv)$push for pf in pr] for pr in parRes]

      -- One polynomial. Implicit variable --
       solve(pol : F) ==
         zero? pol =>
            error "equation is always satisfied"
         lv:=removeDuplicates
             concat(variables numer pol, variables denom pol)
         empty? lv => error "inconsistent equation"
         #lv>1 => error "too many variables"
         solve(pol,first lv)

       -- general solver. Input in equation style. Implicit variables --
       solve(eq : EQ F) ==
         pol:= lhs eq - rhs eq
         zero? pol =>
            error "equation is always satisfied"
         lv:=removeDuplicates
             concat(variables numer pol, variables denom pol)
         empty? lv => error "inconsistent equation"
         #lv>1 => error "too many variables"
         solve(pol,first lv)

       -- general solver. Input in equation style  --
       solve(eq:EQ F,var:SE)  == solve(lhs eq - rhs eq,var)

       -- general solver. Input in polynomial style  --
       solve(pol:F,var:SE) ==
         if R has GcdDomain then
           p:=primitivePart(numer pol,var)
           fp:=(factor p)$GeneralizedMultivariateFactorize(SE,IE,R,R,P R)
           [makeEq(map(makeR2F,ff.factor)$PP2,[var]) for ff in factors fp]
         else empty()

       -- Convert a list of Equations in a list of Polynomials
       makePolys(l: L EQ F):L F == [lhs e - rhs e for e in l]

       -- linear systems solver. Input as list of polynomials  --
       linSolve(lp:L F,lv:L SE) ==
           rec:Record(particular:Union(V F,"failed"),basis:L V F)
           lr : L P R:=[numer f for f in lp]
           rec:=linSolve(lr,lv)$LinearSystemPolynomialPackage(R,IE,SE,P R)
           rec.particular case "failed" => "failed"
           rhs := rec.particular :: V F
           zeron:V F:=zero(#lv)
           for p in rec.basis | p ~= zeron repeat
               sym := newInF(1)
               for i in 1..#lv repeat
                   rhs.i := rhs.i + sym*p.i
           eqs: L EQ F := []
           for i in 1..#lv repeat
             eqs := append(eqs,[(lv.i)::(P R)::F = rhs.i])
           eqs

      -- general solver. Input in polynomial style. Implicit variables --
       solve(lr : L F) ==
         lv :="setUnion"/[setUnion(variables numer p, variables denom p)
               for p in lr]
         solve(lr,lv)

       -- general solver. Input in equation style. Implicit variables --
       solve(le : L EQ F) ==
         lr:=makePolys le
         lv :="setUnion"/[setUnion(variables numer p, variables denom p)
               for p in lr]
         solve(lr,lv)

       -- general solver. Input in equation style  --
       solve(le:L EQ F,lv:L SE)  == solve(makePolys le, lv)

       checkLinear(lr:L F,vl:L SE):Boolean ==
         ld:=[denom pol for pol in lr]
         for f in ld repeat
           if (or/[member?(x,vl) for x in variables f]) then return false
         and/[totalDegree(numer pol,vl) < 2 for pol in lr]

       -- general solver. Input in polynomial style  --
       solve(lr:L F,vl:L SE) ==
           empty? vl => empty()
           checkLinear(lr,vl) =>
                            -- linear system --
               soln := linSolve(lr, vl)
               soln case "failed" => []
               eqns: L EQ F := []
               for i in 1..#vl repeat
                   lhs := (vl.i::(P R))::F
                   rhs :=  rhs soln.i
                   eqns := append(eqns, [lhs = rhs])
               [eqns]

                         -- polynomial system --
           if R has GcdDomain then
             parRes:=triangularSystems(lr,vl)
             [[makeEq(map(makeR2F,f)$PP2,vl) for f in pr]
                                                        for pr in parRes]
           else [[]]

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