syssolp.spad

)abbrev package SYSSOLP SystemSolvePackage
++ Author: P. Gianni
++ Date Created: summer 1988
++ 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
SY       ==> 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 SY) -> 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 SY) -> 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, SY) ->  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, SY) -> 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 SY) -> 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 from MPolyCatRationalFunctionFactorizer(IE, SY, R, P F)

                     ---- Local Functions ----
       linSolve : (L F,    L SY) -> 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 SY) : EQ F ==
         np := numer makeP2F p
         lx := variables np
         for x in lv repeat if member?(x, lx) then break
         up := univariate(np, x)
         (degree up)=1 =>
           equation(x::P(R)::F, -coefficient(up, 0)/leadingCoefficient up)
         equation(np::F, 0$F)

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

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

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

       -- solve a system of rational functions
       triangularSystems(lf : L F, lv : L SY) : L L P R ==
           empty? lv => empty()
           empty? lf => empty()
           #lf = 1 =>
              p := numer(first lf)
              fp := (factor p)$GeneralizedMultivariateFactorize(SY, IE, R, R, P R)
              [[ff.factor] for ff in factorList 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|(df := denom f) ~= 1]
           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), dmp)
             parRes := [pr for pr in parRes|
                       and/[(redPol(fq, pr pretend List(dmp))$gb) ~= 0
                         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 : SY)  == solve(lhs eq - rhs eq, var)

       -- general solver. Input in polynomial style  --
       solve(pol : F, var : SY) ==
         if R has GcdDomain then
           p := primitivePart(numer pol, var)
           fp := (factor p)$GeneralizedMultivariateFactorize(SY, IE, R, R, P R)
           [makeEq(map(makeR2F, ff.factor)$PP2, [var]) for ff in factorList 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 SY) ==
           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, SY, 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 SY)  == solve(makePolys le, lv)

       checkLinear(lr : L F, vl : L SY) : 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 SY) ==
           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 [[]]

--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.